## Optimization algorithm again

Linux specific questions.
greenink
Posts: 200
Joined: Jan 28, 2016 15:45

### Optimization algorithm again

More scale free numeric optimization for Linux AMD64

Code: Select all

namespace EvoRnd

dim as long ffile

sub init()
if ffile<>0 then return
ffile=freefile()
Open "/dev/urandom" For Input As #ffile
end sub

sub finish()
close #ffile
ffile=0
end sub

function rand() as ulongint
dim as ulongint k
get #ffile,,k
return k
end function

end namespace

' returns a scale free mutation between (-1,1).  The precision defines the minimum magnitude and scaling range.
function mutateSingleSym naked (k as ulongint,precision as ulong) as single
asm
mov eax,edi
bswapq rdi
imul rax,rsi
mov esi,126
and edi,0x807fffff
shr rax,32
xorps xmm0,xmm0
sub esi,eax
jc endms
sal esi,23
or edi,esi
movd xmm0,edi
endms:
ret
end asm
end function

function mutateSym(precision as ulong) as single
return mutateSingleSym(EvoRnd.rand(),precision)
end function

function mutate(precision as ulong) as single
var r=mutateSingleSym(EvoRnd.rand(),precision)
return iif(r<0!,-r,r)
end function

'random uniform between 0 and 1 excludes 1
function uniform() as single
return EvoRnd.rand()*5.4210107E-20!
end function

'converts a random ulong k to a long between [min,max]. Includes min and max
function rndInt(k as ulong,min as longint,max as longint) as long
dim as ulongint r=k*(max-min+1)
return (r shr 32)+min
end function

sub permutate(x() as ulong,m as ulong)
for i as ulong=0 to m
var r=rndInt(EvoRnd.rand(),i,ubound(x))
swap x(i),x(r)
next
end sub

end namespace

namespace xfile

function openfile(filename as string) as long
var ff=freefile()
open filename for binary as #ff
seek #ff,1
return ff
end function

sub closefile(free_file as long)
close #free_file
end sub

end namespace

#define fnType function(x() as single) as single

type EvoSF
dimM as ulong
dimPrm(any) as ulong
precision1 as ulong
precision2 as ulong
parentIter as ulong
childIter as ulong
grandparentCost as single
grandparent(any) as single
parent(any) as single
child(any) as single
declare constructor(fnDim as ulong,precision1 as ulong,precision2 as ulong,parentIter as ulong,childIter as ulong)
declare destructor()
declare sub initMemory()
declare sub mutateVec(mutated() as single,in() as single)
declare sub optimise(fn as fnType)
declare function getResult(x() as single) as single
declare sub copyVec(x() as single,y() as single)
declare function load(free_file as long) as boolean
declare function save(free_file as long) as boolean

end type

constructor EvoSF(fnDim as ulong,precision1 as ulong,precision2 as ulong,parentIter as ulong,childIter as ulong)
this.precision1=precision1
this.precision2=precision2
this.parentIter=parentIter
this.childIter=childIter
this.dimM=fnDim-1
initMemory()
EvoRnd.init()
grandparentCost=1!/0!
for i as ulong=0 to dimM
next
end constructor

destructor EvoSF()
erase grandparent,parent,child,dimPrm
EvoRnd.finish()
end destructor

sub EvoSF.initMemory()
redim grandparent(dimM)
redim parent(dimM)
redim child(dimM)
redim dimPrm(dimM)
for i as ulong=0 to dimM
dimPrm(i)=i
next
end sub

function EvoSF.save(free_file as long) as boolean
var e=put(#free_file,,precision1)
e or=put(#free_file,,precision2)
e or=put(#free_file,,parentIter)
e or=put(#free_file,,childIter)
e or=put(#free_file,,dimM)
e or=put(#free_file,,grandparentcost)
e or=put(#free_file,,grandparent())
return e=0
end function

function EvoSF.load(free_file as long) as boolean
var e=get(#free_file,,precision1)
e or=get(#free_file,,precision2)
e or=get(#free_file,,parentIter)
e or=get(#free_file,,childIter)
e or=get(#free_file,,dimM)
e or=get(#free_file,,grandparentcost)
initMemory()
e or=get(#free_file,,grandparent())
return e=0
end function

sub EvoSF.mutateVec(mutated() as single,in() as single)
copyVec(mutated(),in())
for i as ulong=0 to m
var idx=dimPrm(i)
if r>1! or r<0! then r=in(idx)
mutated(idx)=r
next
end sub

sub EvoSF.copyVec(x() as single,y() as single)
for i as ulong=0 to dimM
x(i)=y(i)
next
end sub

function EvoSF.getResult(x() as single) as single
copyVec(x(),grandparent())
return grandparentcost
end function

sub EvoSF.optimise(fn as fnType)
for i as ulong=0 to parentIter-1
mutateVec(parent(),grandparent())
var parentcost=fn(parent())
for j as ulong=0 to childIter-1
mutateVec(child(),parent())
var childcost=fn(child())
if childcost<parentcost then
parentcost=childcost
copyVec(parent(),child())
end if
next
if parentcost<grandparentcost then
grandparentcost=parentcost
copyVec(grandparent(),parent())
end if
next
end sub

/'f(-0.54719,-1.54719)=-1.9133
function test(u() as single) as single
var x=u(0)*5.5!-1.5!
var y=u(1)*7!-3!
return sin(x+y)+(x-y)*(x-y)-1.5*x+2.5*y+1
end function '/
' -391.16599
function test(u() as single) as single
dim as single res
for i as ulong=0 to ubound(u)
var x=u(i)*10!-5!
res+=0.5!*(x*x*x*x-16*x*x+5*x)
next
return res
end function

dim as EvoSF e=EvoSF(10,40,40,100,1000)
e.optimise(@test)
print e.grandparentcost
getkey

It's kinda ugly. I'll have to tidy it up later when I get a new keyboard. I wanted to use truly random numbers for the algorithm, The best I could do was to use /dev/urandom. I'll try using the sound card later to sample noise from the real world as a better source of entropy.
Last edited by greenink on Oct 22, 2016 5:49, edited 1 time in total.
MrSwiss
Posts: 3367
Joined: Jun 02, 2013 9:27
Location: Switzerland

### Re: Optimization algorithm again

I doubt it will work on AMD 64, since the ASM Routine:
mutateSingleSym naked (k as ulongint, precision as ulong) as single
is coded in 32 bit ASM ...
greenink
Posts: 200
Joined: Jan 28, 2016 15:45

### Re: Optimization algorithm again

mov eax,edi -move the lower 32 bits of parameter k into eax, upper 32 bits of rax get zeroed by this move
bswapq rdi - move the highest 32 bits of rdi (parameter k) to the lower 32 bits of rdi, it doesn't matter that they get permuted because k is supposed to be random
imul rax,rsi - the lower 32 bits of rax are random the upper 32 bits zero, multiply by rsi (precision)
mov esi,126 - this is the biased exponent for numbers between 0.5 and 0.9999999 in 32 bit floating point format
and edi,0x807fffff - edi contains a random number (from k) keep the random sign and 23 bit significand as understood as a 32 bit float
shr rax,32 - upper 32 bits of rax contain a random number from 0 to precision-1 (move to lower 32 bits)
xorps xmm0,xmm0 -zero floating point register xmm0 (return result)
sub esi,eax -subtract the random number from the exponent (126)
jc endms - an exponent below zero means a mutation so small it can't be represented in 32 bit float format so just return zero
sal esi,23 -move the exponent into its correct (bitwise) position for the 32 bit floating point format.
or edi,esi -or it in place with the random sign and random significand fractional bits
movd xmm0,edi -put it into the floating point register that is used for floating point returns in Linux AMD64
endms:
ret
Last edited by greenink on Oct 21, 2016 5:03, edited 1 time in total.
greenink
Posts: 200
Joined: Jan 28, 2016 15:45

### Re: Optimization algorithm again

Code: Select all

' http://xoroshiro.di.unimi.it/
namespace xor128
const as ulongint PHI= &h9E3779B97F4A7C15
dim as ulongint a,b

function rand() as ulongint
var result=a+b
b xor=a
a=((a shl 55) or (a shr (64-55))) xor b xor (b shl 14)
b=(b shl 36) or (b shr (64-36))
return result
end function

sub init() constructor
var t=timer()
a=*(@t)
b=a+PHI
end sub

sub seed(s as ulongint)
a=s+PHI
b=s*PHI
end sub
end namespace

union FI
i as ulong
f as single
end union

'scale free random mutation between -2 and 2 exclusive. Precision 127.
function mutateSymP127() as single
dim as FI x
x.i=xor128.rand()
x.i and=&hbfffffff
return x.f
end function

'Precision 63
function mutateSymP63() as single
dim as FI x
x.i=xor128.rand()
x.i and=&hbfffffff
x.i or= &h20000000
return x.f
end function

'Precision 31
function mutateSymP31() as single
dim as FI x
x.i=xor128.rand()
x.i and=&hbfffffff
x.i or= &h30000000
return x.f
end function

'scale free random mutation between 0 and 2. Precision 127.
function mutateP127() as single
dim as FI x
x.i=xor128.rand()
x.i and=&h3fffffff
return x.f
end function

'Precision 63
function mutateP63() as single
dim as FI x
x.i=xor128.rand()
x.i and=&h3fffffff
x.i or= &h20000000
return x.f
end function

'Precision 31
function mutateP31() as single
dim as FI x
x.i=xor128.rand()
x.i and=&h3fffffff
x.i or= &h30000000
return x.f
end function

'https://en.wikipedia.org/wiki/Test_functions_for_optimization
'Lévi function N.13
function levi(x as single,y as single) as single
var pi=4*atn(1)
return sin(3*pi*x)*sin(3*pi*x)+(x-1)*(x-1)*(1+sin(3*pi*y)*sin(3*pi*y))+(y-1)*(y-1)*(1+sin(2*pi*y)*sin(2*pi*y))
end function

screenres 512,512,32
'mutations equally likely to ocurr at any scale down to the resolution limit of
'the floating point representation.
for i as ulong=0 to 100000
var x=-log(0.5*abs(mutateSymP127()))
var y=-log(0.5*abs(mutateSymP127()))
pset (x*5,y*5)
next
getkey
cls
dim as single x,y,cost=levi(x,y) 'cost at (0,0)
for i as ulong=0 to 10000
var cx=x+mutateSymP127()
var cy=y+mutateSymP127()
if cx>1! or cx<-1! then cx=x
if cy>1! or cy<-1! then cy=y
var mcost=levi(cx,cy)
if mcost<cost then
x=cx
y=cy
cost=mcost
end if
next
print "cost",cost,"at",x,y
getkey

greenink
Posts: 200
Joined: Jan 28, 2016 15:45

### Re: Optimization algorithm again

Levi function over the -10 to 10 domain.

Code: Select all

getkey
cls
dim as single x,y,cost=levi(10*x,10*y) 'cost at (0,0)
for i as ulong=0 to 10000
var cx=x+mutateSymP127()
var cy=y+mutateSymP127()
if cx>1! or cx<-1! then cx=x
if cy>1! or cy<-1! then cy=y
var mcost=levi(10*cx,10*cy)
if mcost<cost then
x=cx
y=cy
cost=mcost
end if
next
print "cost",cost,"at",10*x,10*y
getkey
greenink
Posts: 200
Joined: Jan 28, 2016 15:45

### Re: Optimization algorithm again

I need to find some interesting visual problems to test this code. I also have to check some things like load/save.

Code: Select all

' http://xoroshiro.di.unimi.it/
namespace Xor128
const as ulongint PHI= &h9E3779B97F4A7C15
dim as ulongint a,b

function rand() as ulongint
var result=a+b
b xor=a
a=((a shl 55) or (a shr (64-55))) xor b xor (b shl 14)
b=(b shl 36) or (b shr (64-36))
return result
end function

sub init() constructor
var t=timer()
a=*(@t)
b=a+PHI
end sub

sub seed(s as ulongint)
a=s+PHI
b=s*PHI
end sub

end namespace

union FI
i as ulong
f as single
end union

function mutateSymP31() as single
dim as FI x
x.i=xor128.rand()
x.i and=&hbfffffff
x.i or= &h30000000
return x.f
end function

'returns a scale free mutation between 0 and 2. Precision 127
function mutateP127() as single
dim as FI x
x.i=xor128.rand()
x.i and=&h3fffffff
return x.f
end function

'random uniform between -1 and 1 excludes -1 and 1
function uniformSym() as single
dim as longint k=Xor128.rand()
return k*1.08420214e-19!
end function

'returns a long between [min,max]. Includes min and max
function rndInt(min as longint,max as longint) as long
dim as ulong k=Xor128.rand()
dim as ulongint r=k*(max-min+1)
return (r shr 32)+min
end function

'randomly permute the elements 0 to m of an array with the rest of the array
sub permutate(x() as ulong,m as ulong)
for i as ulong=0 to m
var r=rndInt(i,ubound(x))
swap x(i),x(r)
next
end sub

end namespace

namespace xfile

function openfile(filename as string) as long
var ff=freefile()
open filename for binary as #ff
seek #ff,1
return ff
end function

sub closefile(free_file as long)
close #free_file
end sub

end namespace

#define fnType function(x() as single) as single

type Opt5
dimM as ulong
parentIter as ulong
childIter as ulong
grandparentCost as single
dimPrm(any) as ulong
grandparent(any) as single
parent(any) as single
child(any) as single
declare constructor(fnDim as ulong,parentIter as ulong,childIter as ulong)
declare destructor()
declare sub initMemory()
declare function mutateVec(mutated() as single,in() as single) as ulong
declare sub copyVec(x() as single,y() as single)

declare sub optimise(fn as fnType)
declare function getResult(x() as single) as single
declare sub recostGrandparent(fn as fnType)
declare function load(free_file as long) as boolean
declare function save(free_file as long) as boolean
end type

constructor Opt5(fnDim as ulong,parentIter as ulong,childIter as ulong)
this.dimM=fnDim-1
this.parentIter=parentIter
this.childIter=childIter
initMemory()
grandparentCost=1!/0!
for i as ulong=0 to dimM
next
end constructor

destructor Opt5()
erase grandparent,parent,child,dimPrm
end destructor

sub Opt5.initMemory()
redim grandparent(dimM)
redim parent(dimM)
redim child(dimM)
redim dimPrm(dimM)
for i as ulong=0 to dimM
dimPrm(i)=i
next
end sub

function Opt5.save(free_file as long) as boolean
var e=put(#free_file,,parentIter)
e or=put(#free_file,,childIter)
e or=put(#free_file,,dimM)
e or=put(#free_file,,grandparentcost)
e or=put(#free_file,,grandparent())
return e=0
end function

function Opt5.load(free_file as long) as boolean
var e=get(#free_file,,parentIter)
e or=get(#free_file,,childIter)
e or=get(#free_file,,dimM)
e or=get(#free_file,,grandparentcost)
initMemory()
e or=get(#free_file,,grandparent())
return e=0
end function

function Opt5.mutateVec(mutated() as single,in() as single) as ulong
dim as ulong m,mcount
copyVec(mutated(),in())
m=int((dimM+1)*0.5!*RndAdapt.mutateP127())    'number of elements to mutate
RndAdapt.permutate(dimPrm(),m)            'randomly select which elements to mutate
for i as ulong=0 to m
var idx=dimPrm(i)
var inVal=mutated(idx)             'because in() is already copied to mutated()
var r=inVal+RndAdapt.mutateSymP31()    'this could result in no change if inVal is large and the mutation small
if r<=-1! or r>=1! then r=inVal    'the parameters supply to the function to optimise are between -1 and 1 exclusive
if r<>inVal then mcount+=1          'count genuine mutations
mutated(idx)=r
next
return mcount
end function

' to deal with optima that drift around
sub Opt5.recostGrandparent(fn as fnType)
grandparentcost=fn(grandparent())
end sub

sub Opt5.copyVec(x() as single,y() as single)
for i as ulong=0 to dimM
x(i)=y(i)
next
end sub

function Opt5.getResult(x() as single) as single
copyVec(x(),grandparent())
return grandparentcost
end function

sub Opt5.optimise(fn as fnType)
for i as ulong=0 to parentIter-1
var parentcost=iif(mutateVec(parent(),grandparent())=0,grandparentcost,fn(parent()))
for j as ulong=0 to childIter-1
if mutateVec(child(),parent())=0 then continue for 'child is the parent mutated in some scale free number of elements by a scale free mutation
var childcost=fn(child())         ' if the mutations are so small as to have no effect then skip this part
if childcost<parentcost then        'if the child is better than the parent then it takes it's place
parentcost=childcost
copyVec(parent(),child())
end if
next
if parentcost<grandparentcost then     'if the parent is better than the grandparent it takes it's place
grandparentcost=parentcost
copyVec(grandparent(),parent())
end if
next
end sub

' -391.16599
function test(u() as single) as single
dim as single res
for i as ulong=0 to ubound(u)
var x=u(i)*5!
res+=0.5!*(x*x*x*x-16*x*x+5*x)
next
return res
end function

dim as Opt5 e=Opt5(10,10,700)
e.optimise(@test)
print e.grandparentcost
getkey
greenink
Posts: 200
Joined: Jan 28, 2016 15:45

### Re: Optimization algorithm again

If you increase the precision to 1000 or so it works a lot faster. I've just started experimenting with it.
Last edited by greenink on Oct 24, 2016 14:01, edited 1 time in total.
greenink
Posts: 200
Joined: Jan 28, 2016 15:45

### Re: Optimization algorithm again

Code: Select all

const testpic="a.bmp"
const points=1000
const smooth=50
const brightness=100
const precision=20000
const parentIter=1
const childIter=25
greenink
Posts: 200
Joined: Jan 28, 2016 15:45

### Re: Optimization algorithm again

After a bit of evolution.
I guess you could create a real image compression algorithm out of it by combining points under different degrees of blurring.
Also, I think it is not generally known that it is possible to combine Gaussian blurred images to create a sharp image. Sort of Gaussian synthesis similar to Fourier synthesis.
greenink
Posts: 200
Joined: Jan 28, 2016 15:45

### Re: Optimization algorithm again

Deleted.
Last edited by greenink on Nov 09, 2016 7:27, edited 1 time in total.
greenink
Posts: 200
Joined: Jan 28, 2016 15:45

### Re: Optimization algorithm again

Last edited by greenink on Nov 04, 2016 14:19, edited 1 time in total.
greenink
Posts: 200
Joined: Jan 28, 2016 15:45

### Re: Optimization algorithm again

I'll put this back since I can see it improve a little bit over an hour or so, on very small training sets:
I think there is a problem in some of the current deep neural nets where they have to be trained on (mini) batches of examples. Very inconvenient. It will be interesting to see if the drop-in idea of using effectively random crossover breaks that curse in back-propagation neural nets. I'm sure you would pay a price though in extra training time.
MichaelW
Posts: 3500
Joined: May 16, 2006 22:34
Location: USA

### Re: Optimization algorithm again

Why not try RDRAND, assuming your processor supports it.

### Who is online

Users browsing this forum: No registered users and 1 guest