Optimization algorithm again

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

Optimization algorithm again

Post by greenink »

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: 3910
Joined: Jun 02, 2013 9:27
Location: Switzerland

Re: Optimization algorithm again

Post by MrSwiss »

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

Post by greenink »

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

Post by greenink »

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

Post by greenink »

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

Post by greenink »

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

Post by greenink »

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

Post by greenink »

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

Post by greenink »

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

Post by greenink »

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

Post by greenink »

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

Post by greenink »

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=0BwsgM ... ExidGZSSHc
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

Post by MichaelW »

Why not try RDRAND, assuming your processor supports it.
Post Reply