Blog, page 0

which fall from simple trig. There's one more equation, which constrains the side length to 61 cm:

I'm planning on using MATLAB's FMINCON, which means I need to formulate this as a minimization problem. This is accomplished by observing

$begin{split} min qquad & -L_1 - L_2 text{subject to} qquad & L_1^2 -T_1^2 - 51^2 = 0 & L_2^2 -T_2^2 - 61^2 = 0 & T_1 + T_2 - 61 = 0 end{split}$

function xsol = solveproblem()
f = @(x) -x(1)-x(2);
x0 = [0;0;0;0]; LB=[0;0;0;0];
settings = optimset('TolFun',1e-8,'Algorithm','interior-point');
xsol = fmincon(f,x0,[],[],[],[],LB,[],@nonlincon,settings);

function [c,ceq]=nonlincon(x)
l1=x(1);l2=x(2);t1=x(3);t2=x(4);
c=[];
ceq = [l1^2-t1^2-51^2;
        l2^2-t2^2-61^2;
        t1 + t2 - 61];
end
end

spamfunc for optimization in matlab (from 2012/01/30)

what spamfunc is

In developing optimization algorithms, one of the most tedious parts is trying different examples, each of which might have its own starting points or upper or lower bounds or other information. The tedium really starts when your algorithm requires first or second order information, which might be tricky to calculate correctly. These bugs can be pernicious, because it might be difficult to differentiate between a bug in your algorithm and a bug in your objective or constraint evaluation. Handily, Northwestern Professor Robert Fourer wrote a language called AMPL, which takes a programming-language specification of objective and constraints and calculates derivatives as needed. The official amplfunc/spamfunc reference is contained in Hooking Your Solver to AMPL, but I'm shooting for a more low-key introduction.

You'll need to compile spamfunc before using it, but I'm saving that for last–until after I've shown off it's power a little bit. You'll also need a .nl file for spamfunc; these are generated by AMPL's command line program from an AMPL .mod source file. I'll be using rosenb.nl as my example here, for the Rosenbrock function.

>> [x0,bl,bu,lam,cl,cu]=spamfunc('rosenb.nl');

This serves as the initialization. It's worth noting the mathematical form AMPL bends your problem into. In particular, AMPL's spamfunc returns your problem in the form

So we've initialized 6 variables for our environment: the initial x_0

we specified in the mod file, lower( x_L

) and upper bounds on those variables, the initial Lagrange multipler(lam, short for lambda

), and the lower( c_L

) and upper( c_U

) bounds on the constraints. To get zeroth order information(that is, f(x)

and

) for a particular point

, just call spamfunc again with the the

you'd like to evaluate at, as in

>> [f,c] = spamfunc(x,0);

To get first order information(the gradient of the objective and Jacobian of constraints), call spamfunc as

>> [g,A]=spamfunc(x,1);

so that the gradient nabla f(x)

is stored in g, and the $i^{th}$ row of A is the transposed-gradient (nabla c_i(x))^T

. That is, the entry $A_{ij}$ of the matrix

is the derivative $frac{partial c_i(x)}{partial x_j}$ .

Finally, to get the Hessian of the Lagrangian, call spamfunc with just the multipliers lambda

, as in

>> WL = spamfunc(lam);
>> WF = spamfunc(0*lam); % if you just want \nabla^2 f(x).

how to use spamfunc in your own code

Fine and good, but how do you use this in your code? I'll consider a very simple unconstrained optimization problem(and finally define the Rosenbrock function, if you've been waiting with bated breath). The Rosenbrock function–or, Rosenbrock's banana function, as I'll be calling it from now on, is defined as

If we wanted to use our spamfunc interface to solve this, we'd generate the .nl file as described below, then do a standard Newton method iteration. Since Newton steps are defined as p_k = -nabla^2 f(x_k) nabla f(x_k)

and then the iterate is updated as $x_{k+1}=x_k + p_k$ , the most basic possible form of the algorithm is

[x,bl,bu,lam,cl,cu]=spamfunc('rosenb.nl');

grad = inf;
while norm(grad) > 1e-6
	    [grad,A] = spamfunc(x,1);
	    W = spamfunc(lam);
	    pk = -W\grad;
	    x = x + pk
end

Note: This is a terrible algorithm; there're no bells and whistles like line search, Hessian modification, etc. But it demonstrates calling at least a couple of spamfunc options.

how to compile spamfunc

I recovered an outline for building spamfunc and amplfunc from a project I've worked on; it should get you a compiled spamfunc binary. You need to make sure to have mex (MATLAB's compiler wrapper) in your path and properly configured.

wget --ftp-user anonymous --ftp-password traviscj@traviscj.com ftp://netlib.org/ampl.tar
tar xf ampl.tar
cd ampl/solvers
sh configurehere
mv makefile makefile.dist
sed "s/CFLAGS = -O/CFLAGS = -fPIC -O/g" makefile.dist > Makefile
make
cd examples
# following line should reflect where mex(1) is.
export PATH=/meas/bin:$PATH
mv makefile makefile.dist
sed "s/- o/-o/g" makefile.dist | sed "s/mex -I/mex -ldl -I/g" > Makefile
make spamfunc.mex amplfunc.mex

It's worth noting the difference between amplfunc and spamfunc; amplfunc returns dense matrices while spamfunc returns them in sparse format. MATLAB can do linear algebra on sparse format matrices much more efficiently than dense format matrices(as long as the matrices themselves are sparse enough!)

how to generate nl files

rosenb.mod

var x {1..2};

minimize obj: 100*(x[2] - x[1]^2)^2 + (1-x[1])^2;

let x[1] := 1.2;
let x[2] := 1.2;

This file tells AMPL what how many variables you need, how you'd like to name them, what objective function you're minimizing, and what initial values they should take. To compile it to a .nl file, use

$ ampl -P -ogrosenb rosenb.mod

VirtualBox and HRD/JT65 (from 2012/01/10)

One of the semi-frustrating aspects of the state of computing and ham radio is the unavailability of some software packages on OSX. In particular, I've been interested in experimenting with JT65 and also wanted to use my Yaesu FT-857d's CAT cable to program memories and also try out HRD. (I should add–I've tried the JT65 stuff and some Yaesu memory programmers, but wasn't satisfied with the Linux setup for those. I'll probably look into getting them working more later, but for now, this article.)

In any case, I use VirtualBox with a Windows install to get support for these packages. The difficulty, though, is getting the Signalink interfaced with Windows programs. I've been unable to get them to interface directly, but one workaround is to:

This works, but there's another difficulty now: Your windows sounds will be broadcast over the air, which is a bit of a ’'lid’’ mistake. So probably turn off oyur Windows sound effects as well.

Congratulations! You have an OSX box running Windows running ham radio software!

tcjblog.py (from 2011/09/09)

I've now (for the most part) finished a working version of some software I've wanted to tackle for a while. I call it tcjblog.py. It's a blog platform built using jemdoc.py to format simple text files into nice looking HTML. What're the benefits of this, and features of tcjblog.py in general?

There's still a fair bit of work for it, but probably it's time to focus a bit on the content and less on the software, at least until it seems more pressing. This has been a great couple-afternoon project for me–a nice break from the almost 5000 lines of numerical/algorithms stuff in C I've done in the last couple of weeks.

It's worth noting that the static approach generates a fair number of files for this approach: one for each entry, one for each category, and one for each month, plus a tenth of the total number again for the index pages. So I use a makefile to manage all of it.

Five Critical Textbooks for (Applied!) Math & Physics Students (from 2011/09/04)

In the course of working through my first year at grad school, I've come up with five favorites for the basics of an undergrad understanding of the essential topics for applied math and physics. Without ado and in the order I'd take them off my shelf:

I've spent probably the most time with Boas’ Mathematical Methods in the Physical Sciences–I've worked nearly 1000 problems out of the book to get ready for the preliminary exams my first year. It was completely worthwhile.

Getting GnuPlot to compile on OSX (from 2010/05/27)

I misunderstood an officemate’s question yesterday and set out to compile gnuplot. Stock GNUPlot will fail to compile on OSX due to an incompatibility in the readline library. So when you compile, let GNUPlot know to use the one included:

> tar xzvf gnuplot*
> cd gnuplot-4.2.39
> ./configure --prefix=/usr/local --with-readline=builtin
> make
> sudo make install

ZSH and scp/rsync globbing (from 2010/05/25)

tjohnson> scp tcj:*txt .
zsh: no matches found: tcj:*txt

I’d always fixed this by the up-CTRL-A-over-over-over-over-quote-CTRL-E-over-over-quote method, ie

tjohnson> scp "tcj:*txt" .

tjohnson> scp tcj:\*txt .

However, in the spirit of inspired overkill, there are a couple other options. You can define some new aliases with

alias scp='noglob scp'

autoload -U url-quote-magic
zle -N self-insert url-quote-magic
zstyle -e :urlglobber url-other-schema \
'[[ $words[1] == scp ]] && reply=("*") || reply=(http https ftp)'

Now, as I type scp tcj:* it replaces it with \*, negating the problem. Also, this fixes the problem where

wget http://traviscj.com/dynpage.php?var=1&importantvar=zyzyx

interprets the & as a ‘background this command’ and throws a couple of errors. Nice!

Pylab PDF (from 2010/05/25)

I’m in the process of kicking the MATLAB habit and replacing it with Python. Pylab is a reasonable approximation to the 2D plotting capabilities of MATLAB, but by default Pylab doesn’t run on headless(ie, non-GUI) boxes. So on my macbook, where I’m often running things interactively, I have

backend : macosx

backend : PDF

Spaces setup (from 2010/05/25)

My workflow is probably pretty specific, but I figured I’d document it for people looking for ideas on how to set it up. This one seems to work pretty well for me.

I have this arranged in a 3×2 grid, such that Textmate is top center. From there I can navigate to terminal, math programs, and Firefox directly(control-{left,right,down}). It takes two navigations to get to communication programs and entertainment, which is good for concentration while working.

Why We Know Less About the World than Ever (from 2010/01/11)

In a short break from studying for tomorrow's preliminary exam, I came across a TED talk about Why We Know Less About the World Than Ever. I've spent quite a lot of energy thinking about the main question Alisa Miller brings up, which is: how did we get news get so focused on things like Anna Nichole Smith's death while ignoring North Korean nuclear disarmament, the IPCC's global warming report, and other arguably more important news stories?

Miller is quick to blame the networks: It's just cheaper to cover Smith's death than it is to have foreign news offices sprinkled throughout the world. I'd like to believe her, but I'm not sure there'd be a market for the cheaper celebrity coverage if Americans still insisted on our news being about something else.

I was also a little discouraged to hear her quote the studies that Google News covered the same news stories as everyone else. Naturally, they have to; but the beauty of Google News is that people can hone in what they'd like covered. Put another way, I'd be more interested in hearing what the statistics on material reported indicated among GN users who had customized their coverage(remove celebrity news, increase world news, increase technology & science news, and soforth).

Anyways, I certainly agree with her that we should be served more diverse material for news, but her anti-everyone-else words seem more like an advertisement for her news service than anything, even if it's not mentioned in the TED talk. But I suppose exposing the ills of everything else is a better business move than telling people that they care about the wrong things.