The following optimizer code is almost a drop in replacement for the optimizers in the FBMath library. It is rather slow because it uses both global and local optimization. It should be able to solve some difficult problems. Oh yeh, I am pasting from linux, I think there was some issue about that before???? End of line characters?
http://code.google.com/p/lemontree/downloads/list
Code: Select all
'******************************************************************
' Algorithm name: Hanoi Implosion EDA
' Purpose: To locate the minimum of a function in a given domain
' Outline: Uses a global EDA to
' feed vectors to a local EDA. The local EDA creates a population
' of micro-mutated vectors from the given vector, which then
' explode into the basin of the nearest local minimum. The vectors
' then undergo an implosion to the actual local minimum itself.
' ******************************************************************
' Devised and written by: Sean O'Connor
' Email: readordream@yahoo.ie
' Date: 22 August 2008
' Programming Language: FreeBasic - Version 0.20.0 Beta
' ******************************************************************
' Walsh-Hadamard Transform
#macro WHT(wVec)
w0=ubound(wVec)+1
w1=1
while(w1<w0)
w2=w1+w1
w3=0
while(w3<w0)
for w4=w3 to w3+w1-1
wx0=wVec(w4)
wx1=wVec(w4+w1)
wVec(w4)=wx0+wx1
wVec(w4+w1)=wx0-wx1
next
w3+=w2
wend
w1=w2
wend
#endmacro
'Builds a heap to store the lowest cost atoms
#macro buildheap(wVec,wMark)
w0=ubound(wVec)+1
for w1=2 to w0
w2=w1
w3=w1
do
w3=w3 shr 1
if wVec(w2-1)wMark<=wVec(w3-1)wMark then exit do
swap wVec(w2-1),wVec(w3-1)
w2=w3
loop while w3<>1
next
#endmacro
'Inserts atoms into the heap if their cost is low enough
#macro insertHeap(wHeap,wVec,wMark)
w0=ubound(wHeap)+1
for w1=0 to ubound(wVec)
if wVec(w1)wMark>=wHeap(0)wMark then continue for
swap wVec(w1),wHeap(0)
w2=1
do
w3=w2+w2
if w3>w0 then exit do
if w3=w0 then
if wHeap(w2-1)wMark>=wHeap(w3-1)wMark then exit do
swap wHeap(w2-1),wHeap(w3-1)
exit do
end if
if wHeap(w3)wMark>wHeap(w3-1)wMark then w3+=1
if wHeap(w2-1)wMark>=wHeap(w3-1)wMark then exit do
swap wHeap(w2-1),wHeap(w3-1)
w2=w3
loop
next
#endmacro
#macro permute(wprm)
w0=ubound(wprm)+1
for w1=0 to ubound(wprm)
swap wprm(w1),wprm(int(rnd*w0))
next
#endmacro
#macro permute1(wprm)
w0=ubound(wprm)
for w1=1 to w0
swap wprm(w1),wprm(1+int(rnd*w0))
next
#endmacro
#macro clip(wArray)
for w0=0 to vecTop
wx0=xmin(w0+Lbx)
wx1=xmax(w0+Lbx)
for w1=0 to ubound(wArray)
wx2=wArray(w1).vec[w0]
if wx2<wx0 then wx2=wx0
if wx2>wx1 then wx2=wx1
wArray(w1).vec[w0]=wx2
next
next
#endmacro
#macro evaluate(wArray)
for w0=0 to ubound(wArray)
for w1=0 to vecTop
work(Lbx+w1)=wArray(w0).vec[w1]
next
wArray(w0).cost=func(work())
next
#endmacro
#macro rand(wArray,wstart,wtop)
for w0=0 to vecTop
wx0=xmin(w0+Lbx)
wx1=xmax(w0+Lbx)-wx0
for w1=wStart to wTop
wArray(w1).vec[w0]=wx0+wx1*rnd
next
next
#endmacro
'Applies an O'Connor transform on each dimension in turn converting it to
'a Gaussian (Normal) distribution. In this particular OCT the mean is invariant.
'converge stays zero if the population has collapsed to a common point.
#macro octmix(wSource,wDest,wMix,wScale,wScaleSD)
converge=0
for w5=0 to vecTop 'for each dimension
for w6=0 to ubound(wSource)
wMix(w6)=wSource(w6).vec[w5] 'gather dimension
next
permute(wMix)
wht(wMix)
wMix(0)*=wScale
for w6=1 to ubound(wSource)
wx0=wMix(w6)
wx0*=wScaleSD 'effect on standard deviation
wMix(w6)=wx0
converge+=wx0*wx0
next
permute1(wMix) 'all-sum term of previous WHT not moved
wht(wMix)
for w6=0 to ubound(wSource)
wDest(w6).vec[w5]=wMix(w6) 'scatter dimension
next
next
#endmacro
#macro implode(wArray)
for w7=0 to globalTop
for w8=0 to localTop
for w9=0 to vecTop
localAtoms(w8).vec[w9]=wArray(w7).vec[w9]+(.5-rnd)*delta
next
next
clip(localAtoms)
evaluate(localAtoms)
buildHeap(localAtoms,.cost)
for w8=0 to localGen-1
octmix(localAtoms,localTrial,localMix,localScale,localScale*(AlphaStart-w8*AlphaStep))
clip(localTrial)
evaluate(localTrial)
insertHeap(localAtoms,localTrial,.cost)
if converge=0 then exit for
next
wx0=localAtoms(0).cost
w8=0
for w9=1 to localtop
if localAtoms(w9).cost<wx0 then
wx0=localAtoms(w9).cost
w8=w9
end if
next
wArray(w7).cost=wx0
for w9=0 to vecTop
wArray(w7).vec[w9]=localAtoms(w8).vec[w9]
next
if wx0<fMin then
fMin=wx0
for w9=0 to vecTop
x(Lbx+w9)=localAtoms(w8).vec[w9]
next
end if
next
#endmacro
type atom
vec as double ptr
cost as double
end type
' L2Global: Log base 2 of the population size of the global EDA.
' L2Local: Log base 2 of the population size of the local EDA. Larger values
' increase the size of the explosion into the local minimum basin.
' GenGlobal: Number of generations the global EDA is to run for.
' GenLocal: Number of generations each local EDA is to run for.
' StartAlpha: Controls the size of the explosion into the local minimum basin
' The normal value is 1, larger values increase the size of the explosion.Using
' too high a value will cause the local EDA population to splatter into several
' local minima basins, making subsequent implosion more difficult. A value of
' 1.1 would be very high.
' EndAlpha: Controls the rate of implosion to the actual local minimum. 1 is the
' normal value. You can use a value of less than 1 to increase the rate of
' implosion.
' Perturb: The amount each vector presented to the local EDA is micro-mutated
' by to create the local population.
declare SUB InitEDAParams( L2Global AS UINTEGER = 5, L2Local AS UINTEGER = 5, _
GenGlobal AS UINTEGER = 300, GenLocal AS UINTEGER = 100, _
StartAlpha As DOUBLE = 1, EndAlpha As DOUBLE =1, _
Perturb AS DOUBLE = 1E-6)
' Func: The function to minimize
' X(): The result vector. It can also be used to pass a previous good result
' to the algorithm.
' Xmin(): Minimum domain values for each dimension.
' Xmax(): Maximum domain values for each dimension.
' Fmin: The minimum funcion value found.
declare SUB EDAalg( Func AS FUNCTION(X() AS DOUBLE) AS DOUBLE, _
X() AS DOUBLE, _
Xmin() AS DOUBLE, _
Xmax() AS DOUBLE, _
BYREF Fmin AS DOUBLE)
declare SUB EDA_CreateLogFile(LogFileName AS STRING = "Hanoi_Implosion_EDA.txt" )
declare SUB setEDACallback(cb AS SUB(x() as double,fMin as double,gen as uinteger))
' ------------------------------------------------------------------
' Global variable
' ------------------------------------------------------------------
COMMON SHARED ErrCode AS INTEGER
' Error code from the last function evaluation
' Strings used for printing results in log file
CONST Hdr = " Iter F"
CONST Mask = "##### #####.############"
CONST Hdr2 = " X(n) Value"
' Global variables
DIM SHARED AS INTEGER WriteLogFile = 0 ' Log file ?
DIM SHARED AS INTEGER FF ' File handle
DIM SHARED AS INTEGER GlobalL2 = 5
DIM SHARED AS INTEGER LocalL2 = 5
DIM SHARED AS INTEGER GlobalGen = 300
DIM SHARED AS INTEGER LocalGen = 100
DIM SHARED AS DOUBLE AlphaStart = 1
DIM SHARED AS DOUBLE AlphaEnd = 1
DIM SHARED AS DOUBLE Delta = 1E-6
DIM SHARED AS SUB(x() as double,Fmin as double,gen as uinteger) callBackSub
SUB InitEDAParams( L2Global AS UINTEGER = 5, L2Local AS UINTEGER = 5, _
GenGlobal AS UINTEGER = 300, GenLocal AS UINTEGER = 100, _
StartAlpha As DOUBLE = 1, EndAlpha As DOUBLE =1, _
Perturb AS DOUBLE = 1E-6)
' ------------------------------------------------------------------
' Initialize EDA Algorithm parameters
' ------------------------------------------------------------------
GlobalL2=L2Global
LocalL2=L2Local
GlobalGen=GenGlobal
LocalGen=GenLocal
AlphaStart=StartAlpha
AlphaEnd=EndAlpha
Delta=Perturb
END SUB
SUB EDAalg( Func AS FUNCTION(X() AS DOUBLE) AS DOUBLE, _
X() AS DOUBLE, _
Xmin() AS DOUBLE, _
Xmax() AS DOUBLE, _
BYREF Fmin AS DOUBLE)
' ------------------------------------------------------------------
' Minimization of a function of several variables
' by the MekongEDA algorithm
' ------------------------------------------------------------------
' Input parameters : Func = objective function to be minimized
' X = initial minimum coordinates
' Xmin = minimum value of X
' Xmax = maximum value of X
' ------------------------------------------------------------------
' Output parameters: X = refined minimum coordinates
' Fmin = function value at minimum
' ------------------------------------------------------------------
dim as integer Lbx,gen,i
dim as integer globalSize,globalTop,vecLen,vecTop
dim as integer localSize,localTop
dim as integer w0,w1,w2,w3,w4,w5,w6,w7,w8,w9
dim as double wx0,wx1,wx2,wx3,converge
dim as double oldFmin,GlobalScale,LocalScale,AlphaStep,cost
If lbound(Xmin)<>lbound(Xmax) or lbound(x)<>lbound(Xmin) or ubound(Xmin)<> _
ubound(Xmax) or ubound(x)<>ubound(Xmin) THEN
ErrCode = -3
EXIT SUB
END IF
ErrCode = 0
randomize timer , 3 '3=Mersenne Twister
globalSize=1 shl GlobalL2
globalTop=globalSize-1
GlobalScale=1.0/globalSize
localSize=1 shl LocalL2
localTop=localSize-1
localScale=1.0/localSize
AlphaStep=(AlphaStart-AlphaEnd)/(localGen-1)
Lbx=lbound(x)
vecTop=ubound(x)-Lbx
vecLen=vecTop+1
dim globalMix(globalTop) as double
dim localMix(localTop) as double
dim work(Lbx to ubound(x)) as double
dim globalAtoms(globalTop) as atom
dim globalTrial(globalTop) as atom
dim localAtoms(localTop) as atom
dim localTrial(localTop) as atom
' Allocate space for the atom vectors in the localAtoms and trial arrays
for i=0 to localTop
localAtoms(i).vec=allocate(vecLen*sizeof(double))
localTrial(i).vec=allocate(vecLen*sizeof(double))
next
' Allocate space for the atom vectors in the globalAtoms and trial arrays
for i=0 to globalTop
globalAtoms(i).vec=allocate(vecLen*sizeof(double))
globalTrial(i).vec=allocate(vecLen*sizeof(double))
next
' initalize lowest cost vector and make it part of the global population as well
for i=Lbx to ubound(x)
if x(i)<xmin(i) or x(i)>xmax(i) then x(i)=.5*(xmin(i)+xmax(i))
globalAtoms(0).vec[i-Lbx]=x(i)
next
rand(globalAtoms,1,globalTop)
evaluate(globalAtoms)
Fmin=globalAtoms(0).cost
oldFmin=fMin
buildHeap(globalAtoms,.cost)
'For each generation
for gen=1 to globalGen
octmix(globalAtoms,globalTrial,globalMix,globalScale,globalScale)
clip(globalTrial)
implode(globalTrial)
insertHeap(globalAtoms,globalTrial,.cost)
''' finish of current generation
if callBackSub<>0 and (gen=1 or gen=globalGen or Fmin<oldFmin) then
callbackSub(x(),Fmin,gen)
end if
if WriteLogFile then
PRINT #FF, USING Mask; gen; Fmin
end if
oldFmin=Fmin
next 'next generation
for i=0 to localTop
deallocate(localAtoms(i).vec)
deallocate(localTrial(i).vec)
next
for i=0 to globalTop
deallocate(globalAtoms(i).vec)
deallocate(globalTrial(i).vec)
next
callbackSub=0 'Don't use again by mistake
IF WriteLogFile THEN
print #FF,
print #FF,"Optimum Vector"
print #FF,Hdr2
for i=Lbx to ubound(x)
print #FF,using mask;i;x(i)
next
CLOSE #FF
WriteLogFile = 0
END IF
END SUB
SUB EDA_CreateLogFile(LogFileName AS STRING = "Hanoi_Implosion_EDA.txt" )
' ------------------------------------------------------------------
' Initialize log file
' ------------------------------------------------------------------
FF = FREEFILE
OPEN LogFileName FOR OUTPUT AS #FF
PRINT #FF, "Hanoi Implosion EDA"
PRINT #FF,
PRINT #FF, Hdr
WriteLogFile = -1
END SUB
sub setEDACallback(cb AS SUB(x() as double,fMin as double,gen as uinteger))
callBackSub=cb
end sub
Demo:
Code: Select all
#include "Hanoi_Implosion_EDA.bas"
#include "gl\gl.bi"
#include "gl\glu.bi"
#include "fbgfx.bi"
' Hanoi Implosion EDA Function cost minimiser
' If you want to maximize adjust your function to return minus it's original value
Const PI As Double = 3.1415926535897932
dim as integer i
dim as single smx, smy, viewdist = 30, viewangle = 0.75, viewheight = 25, lookheight = 2.5, aa
dim as single a1,a2,at
dim as string k, skipframe
print
print
print "OpenGL demo to solve a problem in 60 dimensions."
print "That's way to much for most Genetic Algorithms."
print "The problem involves repulsive forces between the centers of"
print "squares and attractive forces between the corners of the squares."
dim shared as double ox(20*4),oy(20*4),cx(20*4),cy(20*4)
declare sub optThreadSub(pram as any ptr)
declare function optCost(v() as double) as double
declare sub glUpdate(x() as double,fMin as double,gen as uinteger)
dim thr as any ptr
thr=ThreadCreate(@optThreadSub,0,100000)
screen 18, 32, , 2
glViewport 0, 0, 640, 480
glMatrixMode GL_PROJECTION
glLoadIdentity
gluPerspective 45.0, 64.0/48.0, 0.1, 10000.0
glMatrixMode GL_MODELVIEW
glLoadIdentity
glShadeModel GL_SMOOTH
glClearColor 0.0, 0.0, 0.0, 1.0
glClearDepth 1.0
glEnable GL_DEPTH_TEST
glDepthFunc GL_LEQUAL
glHint GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST
dim as integer a,b,c,d,e,f,hc,dx,dy
do
glClear GL_COLOR_BUFFER_BIT or GL_DEPTH_BUFFER_BIT
glLoadIdentity
gluLookAt viewdist * cos(viewangle), viewdist * sin(viewangle),viewheight,0,0,lookheight,0,0,1
glPushMatrix
for i = -16 to 16
glColor3f .5,0,0
glBegin( GL_LINES )
glVertex3f( i, 16, 0 )
glVertex3f( i, -16, 0 )
glVertex3f( 16, i, 0 )
glVertex3f( -16, i, 0 )
glEnd()
next
glPopMatrix
glPushMatrix
for i=0 to 19
glColor3f 1-i*.05,i*.05,1
glBegin( GL_QUADS)
glVertex3f(cx(i*4),cy(i*4), .1+i*.01)
glVertex3f(cx(i*4+1),cy(i*4+1), .1+i*.01)
glVertex3f(cx(i*4+2),cy(i*4+2), .1+i*.01)
glVertex3f(cx(i*4+3),cy(i*4+3), .1+i*.01)
glEnd
next
glPopMatrix
flip
ScreenSync
k = inkey
select case k
case "Q", "q"
lookheight = lookheight + 10
case "A", "a"
lookheight = lookheight - 10
case "W", "w"
viewheight = viewheight + 10
case "S", "s"
viewheight = viewheight - 10
case "X", "x"
viewangle = viewangle + 0.1
case "Z", "z"
viewangle = viewangle - 0.1
case "+"
viewdist = viewdist - 10
case "-"
viewdist = viewdist + 10
end select
loop until k = chr(27) or k=chr(255,107)
screen 0
end
sub optThreadSub(pram as any ptr)
dim as integer i
dim as double res
InitEDAParams(5,6,&h7fffffff,200,1.01,1,.1)
setedaCallback(@glUpdate)
dim optMin(20*3-1) as double
dim optMax(20*3-1) as double
dim optX(20*3-1) as double
for i=0 to ubound(optX) step 3
optMin(i)=-15
optMax(i)=15
optMin(i+1)=-15
optMax(i+1)=15
optMin(i+2)=-PI
optMax(i+2)=PI
next
EDAalg(@optCost,optX(),optmin(),optmax(),res)
end sub
function optCost(v() as double) as double
dim as integer i,j,k
dim as double rep,att,d1,d2
dim as double x,y,z,t1,t2,cA,sA,cB,sB,r,s
rep=0 'sum of repulsive forces
att=0 'sum of attractive forces
' the centers of the squares repel each other with app. inverse square law
for i=0 to 19
x=v(i*3)
y=v(i*3+1)
for j=0 to 19
if i<>j then
d1=v(j*3)-x
d2=v(j*3+1)-y
rep+=1/(.1+d1*d1+d2*d2)
end if
next
next
' set up all the 2D coordinates of the squares
for i=0 to 19
r=v(i*3)
s=v(i*3+1)
x=cos(v(i*3+2))
y=sin(v(i*3+2))
ox(i*4)=x+r
oy(i*4)=y+s
x=cos(v(i*3+2)+.5*PI)
y=sin(v(i*3+2)+.5*PI)
ox(i*4+1)=x+r
oy(i*4+1)=y+s
x=cos(v(i*3+2)+PI)
y=sin(v(i*3+2)+PI)
ox(i*4+2)=x+r
oy(i*4+2)=y+s
x=cos(v(i*3+2)+1.5*PI)
y=sin(v(i*3+2)+1.5*PI)
ox(i*4+3)=x+r
oy(i*4+3)=y+s
next
' Create an attractive force between the corners of all the squares
for i=0 to 19 'for each square
for j=0 to 3 'for each corner
x=ox(i*4+j)
y=oy(i*4+j)
for k=0 to 19
if(i<>k) then
d1=ox(k*4)-x
d2=oy(k*4)-y
att+=1/(.1+d1*d1+d2*d2)
d1=ox(k*4+1)-x
d2=oy(k*4+1)-y
att+=1/(.1+d1*d1+d2*d2)
d1=ox(k*4+2)-x
d2=oy(k*4+2)-y
att+=1/(.1+d1*d1+d2*d2)
d1=ox(k*4+3)-x
d2=oy(k*4+3)-y
att+=1/(.1+d1*d1+d2*d2)
end if
next
next
next
return 17*rep-att ' minimizing therefore the attactive force is negative
end function
sub glUpdate(x() as double,fMin as double,gen as uinteger)
dim as integer i
dim as double cost
cost=optCost(x())
for i=0 to 79
cx(i)=ox(i)
cy(i)=oy(i)
next
end sub