Optimization algorithm again

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

Optimization algorithm again

Postby greenink » Oct 20, 2016 15:10

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
      
      
namespace RndAdapt   
   ' 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
      grandparent(i)=RndAdapt.uniform()
   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())
   var m=int(dimM*RndAdapt.mutate(precision2))
   RndAdapt.permutate(dimPrm(),m)
   for i as ulong=0 to m
      var idx=dimPrm(i)
      var r=in(idx)+RndAdapt.mutateSym(precision1)
      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

Postby MrSwiss » Oct 20, 2016 20:11

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

Postby greenink » Oct 21, 2016 1:31

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

Postby greenink » Oct 21, 2016 4:37

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

Postby greenink » Oct 21, 2016 5:01

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

Postby greenink » Oct 22, 2016 13:43

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

namespace RndAdapt   

   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
      grandparent(i)=RndAdapt.uniformSym()
   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

Postby greenink » Oct 24, 2016 5:55

https://drive.google.com/open?id=0BwsgMLjV0BnhWmZYQUg4X05sVzg
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

Postby greenink » Oct 24, 2016 13:59

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

Postby greenink » Oct 24, 2016 15:17

After a bit of evolution.
https://drive.google.com/open?id=0BwsgMLjV0BnhR2tvSlR2djYzMnc
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

Postby greenink » Oct 29, 2016 15:22

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

Postby greenink » Oct 31, 2016 16:19

Edit: Code link removed.
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

Postby greenink » Nov 02, 2016 15:31

I'll put this back since I can see it improve a little bit over an hour or so, on very small training sets:
https://drive.google.com/open?id=0BwsgMLjV0BnhQnhxVExidGZSSHc
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

Postby MichaelW » Nov 02, 2016 22:31

Why not try RDRAND, assuming your processor supports it.

Return to “Linux”

Who is online

Users browsing this forum: No registered users and 1 guest