earth is my favourite planet
life in the pedestrian lane: science, faith, ideas, politics, techArchive for March, 2009
Taming a Beast called UNIX
The year was 1988. Voices hushed in trepidation, the apprentices gathered in a dingy corner of the Engineering building. In the undergraduate computer lab we assembled to offer sacrifices to the mainframe, hunched over the fading amber or green glow of hesitant dumb terminals, all connected to the mighty UNIX (VAX/VMS) mainframe. Poor first year acolytes we were, learning various complex incantations to extract data from the Beast. We learned to navigate the perils of UNIX: using ed, vi, man, fortran, and other strange spells to tame the Beast. For some of those innocent souls, the Beast caused them to cringe in fear and retreat to a “lesser” academic career, or flee gibbering and drooling in fear. Alas, I was an overconfident youth whose world was shattered by the Beast. So it was that I gave up my dream of conquest and glory.
But fate is fickle, and the world ever turns. So it was that seven years later, after a journeyman’s life, I wandered accidentally back into the nether realms of UNIX. However this time, I was different — stronger in body and mind from years of hard labour — and was wary of the nature of the Beast. Also, Unix itself was a different creature, somewhat mellowed with age, more tractable, less hostile: even unto supporting lower-case symbols and X-windows. So it was that an uneasy truce was forged between Unix and I. We learned to work together. Thus it was that I saw the Beast in a different light: a creature of mankind designed to serve, but only those who prove themselves worthy by mastering its secrets.
One of my greatest masterpieces was written in g77, the GNU port of Fortran. This program applies the Revised Simplex Method in two phases to find the optimal solution to a set of linear equations. I can’t remember all the mathematical theory and strategies involved, but the code still looks cool. For posterity I now give it to the world.
program main
implicit none ! Set up variables..
integer o,p
parameter(o=20,p=100)
integer i,j,m,n,r,s
integer apos(o),bpos(o)
real A(o,p),Bi(o,o),b(o),xb(o),cost(p),ca(o+p)
! Read data..
open(unit=8,file='lp.dat',status='old')
open(unit=9,file='lp.out',status='new')
read(8,*) m,n
if (m.gt.o.or.n.gt.p) then
write(9,*)'Exceeded variable/constraint limits'
endif
read(8,*) (cost(j),j=1,n)
write(9,*) 'm ',m
write(9,*) 'n ',n
write(9,*) 'c'' ',(cost(j),j=1,n)
write(9,*)
write(9,*) 'A.x = b'
do i=1,m
read(8,*) (A(i,j),j=1,n),b(i)
write(9,*) (A(i,j),j=1,n),' =',b(i)
enddo
c ------------------extra bit..-----------------
do i=1,m
read(8,*) (Bi(i,j),j=1,m)
write(9,*) (Bi(i,j),j=1,m)
enddo
read(8,*) (bpos(i),i=1,m)
write(9,*) (bpos(i),i=1,m)
read(8,*) (apos(i),i=1,n)
write(9,*) (apos(i),i=1,n)
do i=1,m
do j=1,m
xb(i)=xb(i)+Bi(i,j)*b(j)
enddo
enddo
goto 5
c -----------------------miss out for now------------------------
write(9,*)
write(9,*) 'Phase 1'
do i=1,n ! Phase 1
ca(i)=0 ! Prepare artificial cost,basis,variables
apos(i)=0 ! all original variables are nonbasic
enddo
do i=1,m
do j=1,m
Bi(i,j)=0 ! initially Bi = identity matrix, I
enddo
Bi(i,i)=1
ca(n+i)=1 ! artificial cost vector now constructed
bpos(i)=n+i ! artificial variables fill the basis
xb(i)=b(i) ! initially xb = b (because xb=Bi.b, Bi=I)
enddo
call rsm(m,n,A,Bi,b,xb,ca,bpos,apos,1)
c ---------------------------------------------------------------
5 write(9,*)
write(9,*) 'Phase 2'
! Phase 2
! Now do RSM with the new basis
call rsm(m,n,A,Bi,b,xb,cost,bpos,apos,2)
end
C--------------------------------------------------
subroutine rsm(m,n,A,Bi,b,xb,c,bpos,apos,phase)
parameter(o=20,p=100)
integer apos(n),bpos(m),phase ! variables passed to rsm
real A(o,p),Bi(o,o),b(m),c(n+m) ! "
integer i,j,imax,r,s,note ! introducing new variables
real bestrc,bias(m),pi(m),xb(m),z ! "
character retval
imax=2*(m+n)
tiny=1E-10
i=0
call currentsoln(Bi,b,bpos,c,m,xb,z)
1 if (i.lt.imax) then ! Iterate Revised Simplex Method..
i=i+1
call computedual(Bi,c,bpos,m,n,pi) ! find pi (initially pi=e)
call reducedcost(A,c,pi,apos,m,n,s,bestrc)
if(s.eq.0) then
if(z.gt.tiny.and.phase.eq.1) then
write(9,*) ' Infeasibility detected.'
note=-2
else
write(9,*) ' Optimality reached.'
note=1
endif
goto 3
endif
call mivitwwbias(A,Bi,m,s,bias)
call lvratiotest(xb,bias,m,n,r,phase,bpos)
if (r.eq.0) then
write(9,*) ' Unboundedness detected.'
note=-1
goto 3
endif
write(9,*) 'basisupdate:'
write(9,*) ' bpos',(bpos(j),j=1,m),' apos',(apos(j),j=1,n)
call basisupdate(Bi,bias,c,pi,xb,bpos,apos,m,n,r,s)
write(9,*) ' ->',(bpos(j),j=1,m),' ->',(apos(j),j=1,n)
call currentsoln(Bi,b,bpos,c,m,xb,z)
else
goto 3
endif
2 goto 1
3 write(9,*) ' RSM iterations executed: ',i
end
C-----------------------------------------------------
subroutine currentsoln(Bi,b,bpos,c,m,xb,z)
parameter(o=20)
integer bpos(*) ! old
real Bi(o,o),b(*),c(*),xb(*),z ! old
integer i,j ! new
z=0
do i=1,m
z=z+c(bpos(i))*xb(i)
enddo
write(9,*) 'currentsoln: Bi, xb'
do i=1,m
write(9,*) ' ',(Bi(i,j),j=1,m),' x',bpos(i),' =',xb(i)
enddo
write(9,*)
write(9,*) ' --------------------------- Objective value z =',z
write(9,*)
end
C-----------------------------------------------------
subroutine computedual(Bi,c,bpos,m,n,pi)
parameter(o=20)
integer bpos(*) ! old
real Bi(o,o),c(*),pi(*) ! old
integer i,j,k ! new
do i=1,m
pi(i)=0
do j=1,m
k=bpos(j)
pi(i)=pi(i)+c(k)*Bi(j,i)
enddo
enddo
write(9,*) 'computedual: pi ',(pi(i),i=1,m)
end
C----------------------------------------------------
subroutine reducedcost(A,c,pi,apos,m,n,s,bestrc)
parameter(o=20,p=100)
integer apos(*),s ! old
real A(o,p),c(*),pi(*),bestrc ! old
integer i,j ! new
real rc,pian,tiny ! new
rc=0
bestrc=0
tiny=-1E-10
do i=1,n ! determine entering var x(s)
if(apos(i).eq.0) then !< look at nonbasic columns
pian=0 ! (never price artificials)
do j=1,m
pian=pian+pi(j)*A(j,i)
enddo
rc=c(i)-pian !< red. cost for nonbasic x(i)
if(rc.lt.bestrc) then
bestrc=rc
s=i
endif
endif
enddo
if(bestrc.ge.tiny) then !< This means no ev was found,
s=0 ! optimality has been reached,
endif ! so flag s to 0
write(9,*) 'reducedcost: x',s,' enters with rc',bestrc
end
C----------------------------------------------------
subroutine mivitwwbias(A,Bi,m,s,bias)
parameter(o=20,p=100)
real A(o,p),Bi(o,o),bias(*)
integer i,j,s
do i=1,m ! a quick calculation of bias,
bias(i)=0 ! aka. the most important vector
do j=1,m ! in the whole world.
bias(i)=bias(i)+Bi(i,j)*A(j,s)
enddo
enddo
write(9,*) 'mivitwwbias: bias',(bias(j),j=1,m)
end
C----------------------------------------------------
subroutine lvratiotest(xb,bias,m,n,r,phase,bpos)
integer bpos(*),r,phase ! old
real xb(*),bias(*) ! old
integer i ! new
real ratio,rmin,tiny ! new
r=0 ! r=0 returned if no lv can be found
tiny=1E-10
teeny=-1E-10
rmin=10E+10
do i=1,m
if(bpos(i).gt.n.and.phase.eq.2) then
! Artificial present in phase2 basis:
! Extended lv routine to force out
! artificial variables..
if(bias(i).lt.teeny.or.bias(i).gt.tiny) then
rmin=0
r=i
endif
else ! Usual lv routine..
if(bias(i).gt.tiny) then ! ratio = rate of change of xb(i), as
ratio=xb(i)/bias(i) ! xs increases from 0 (enters)
if(ratio.lt.rmin) then
rmin=ratio ! the leaving variable is xb(r), the
r=i ! basic variable which decreases most
endif ! rapidly as xs enters.
endif
endif
enddo
write(9,*) 'lvratiotest: xb',r,' leaves.'
end
C----------------------------------------------------
subroutine basisupdate(Bi,bias,c,pi,xb,bpos,apos,m,n,r,s)
parameter(o=20)
integer apos(*),bpos(*),m,n,r,s ! old
real Bi(o,o),bias(*),c(*),pi(*),xb(*) ! old
integer i,j ! new
do i=1,m ! do gauss-jordan pivot on rth element
if(i.eq.r) then ! of bias, and thus also on Bi and xb
do j=1,m
Bi(r,j)=Bi(r,j)/bias(r)
enddo
xb(r)=xb(r)/bias(r)
else
do j=1,m
Bi(i,j)=Bi(i,j)-bias(i)/bias(r)*Bi(r,j)
enddo
xb(i)=xb(i)-bias(i)/bias(r)*xb(r)
endif
enddo
if(bpos(r).le.n) then
apos(bpos(r))=0 ! xb(r) now nonbasic; 0 in apos
endif
bpos(r)=s ! xb(r) replaced by xs in bpos
apos(s)=r ! xs located in rth basic position
end
simplex.f
Note: the above code has not been compiled or executed since sometime around 1998.
Here are some problems solved using this glorious program.
Read the rest of this entry »
Fermat’s Last Theorem
Fermat’s Last Theorem states that:
There are no positive integersand
with
if
is an integer greater than 2.
For 
( 3, 4, 5) ( 5, 12, 13) ( 7, 24, 25) ( 8, 15, 17)
( 9, 40, 41) (11, 60, 61) (12, 35, 37) (13, 84, 85)
(16, 63, 65) (20, 21, 29) (28, 45, 53) (33, 56, 65)
(36, 77, 85) (39, 80, 89) (48, 55, 73) (65, 72, 97) … the Pythagorean triples.
In the margin of his copy of the Arithmetica of Diophantus the French jurist Fermat wrote circa 1637 that for greater
Cubum autem in duos cubos, aut quadrato-quadratum in duos quadrato-quadratos, et generaliter nullam in infinitum ultra quadratum potestatum in duos ejusdem nominis fas est dividere; cujus rei demonstrationem mirabilem sane detexi. Hanc marginis exiguitas non caparet.
Some of the greatest mathematical minds in history have applied their genius to this problem, none succeeded. Every other result which Fermat had announced in like manner had long ago been dealt with; only this one, the last, remained. Finally, 358 years after Fermat wrote his theorem, Andrew Wiles saw an application of elliptic curves and modular forms (20th century mathematical constructs, not known to earlier mathematicians) that could contribute to the proof.
Theorem:
Letwith
.
Ifthen
Proof:
The proof follows a program formulated around 1985 by Frey and Serre[F,S]. By classical results of Fermat, Euler, Dirichlet, Legendre, and Lamé, we may assume that, an odd prime ≥ 11. Suppose
and
. Without loss of generality we may assume that
and
. Frey [F] observed that the elliptic curve
has the following ‘remarkable’ properties:
(1)is semistable with conductor
; and
(2)is unramified outside
and is flat at
.
By the modularity theorem of Wiles and Taylor-Wiles [W,T-W], there is an eigenformsuch that
A theorem of Mazur implies that
is irreducible, so Ribet’s theorem [R] produces a Hecke eigenform
such that
for some
But
has genus zero, so
. This is a contradiction and Fermat’s Last Theorem follows. Q.E.D.
References: Read the rest of this entry »
RE: Darwin & The Deity
-
Should it be of concern to Christians that Darwin was never more than a nominal believer? – well, are you concerned whether your surgeon, mechanic, or hair stylist goes to church? Of course not. Your only concern is that she wields a scalpel, wrench, or scissors with know-how and dexterity. So too with a scientist: one’s only concern should be that he is an honest and skilled practitioner of his craft. And Darwin wasn’t just an able and meticulous biologist, he was a bloody genius. If his theory of evolution by natural selection is the best theory in town that explains the evidence (palaeontological, morphological, genetic) – and it is – deal with it. Of course refute it on empirical grounds if you can, but don’t rubbish it because you don’t like its theological or moral implications, or because you have a political agenda. Fight science with science – not with the pseudo-science of creationism or the bad science of ID (not to mention the bad theology of both).
Unix bash snippet: global replace text in multiple files
Unlike this overkill solution in python, bash is the best tool for this sort of task. The following script searches all files the current directory and subdirectories, and replaces all instances of text “abc” with text “xyz”.
grep and sed are your friends…
for file in `grep -r -l abc *`
do sed 's/abc/xyz/g' $file > $file.fixed
echo $file.fixed
done
USA is going crazy.
Rick Santelli vented a forceful and furious rant against Obama’s multi-trillion stimulus package — many people suspect it will only cause further distortions in a deeply troubled economy, and the only realistic solution is to let people suffer the consequences of vast overspending. Anything else is just an attempt to defy the law of gravity. Personally I don’t know — I hope America can reinvent itself, but there’s a real sense of panic emerging from the States right now.
Even if it was some kind of staged PR exercise, what Santelli actually said is entirely correct. I’ve been saving doggedly all these years, watching everyone jump on board the property gravy train, and there’s NO WAY that I would willingly subsidize them. As Peter Schiff said, the President should be defending the dollar, not embarking on mad inflationary schemes. Obama’s plan is fundamentally unjust.
I Bought an Expensive House. My Bad, Not Yours
The only people affected by plummeting real estate prices are the ones who bought a house that cost more than they could afford, hoping for a spike in value so they could sell at a profit or take out a new loan based on an increased value. Their home wasn’t just a place to live; it was an investment they thought they could liquefy at will. If we’re saving these poor souls from the 26.7% drop in their investment, we should give twice as much aid to everyone who has lost approximately 50% in the stock market since its peak.
Transcript over the jump.. Read the rest of this entry »
RE: THIS IS LIFE, Get Creative.
-
1. This is Life
Quite possibly the most obvious but most eye-opening advice I can ever offer people is the statement: this is life. We may have more than one shot, but nobody knows that for certain so we have to make the most of what we have. This is it. Right now. This is life. More likely than not you will never get this opportunity again.
Let go of all the grudges you hold against others, make sure all the people that are important in your life know how you feel, and don’t waste another minute complaining about something petty. Could you honestly say that if you knew you were dying tomorrow you’ve made the most of the opportunity you have?
If not (which I’m sure applies to most people), then look at areas that instantly popped into your mind when asking yourself that question. What is holding you back?
-
Programming is a mapper’s game. General Tips on Mapping:
Packers have a whole proceduralised culture that provides behavioural tramlines for just about everything. It’s so complete you don’t even notice it until you solve a problem perfectly effectively one day, by a method that’s not on the list. It might be trivial..
Mappers hardly ever get the upper hand on these cultural issues, but when it does happen it can be hilarious. A packer gave a dinner party and it so happened that over half of the guests were mapper types, IT workers and others. The host pulled a pile of warm plates from the oven, and started handing them to the guy on his left.”Just pass them around!”, he cried cheerfully. Everything went well until he realised he needed to shout “Stop!” [everyone was just passing plates around randomly in an endless loop!]Mappers don’t have a general cultural context to learn from, so we are almost entirely self taught.
.JPG)
