tangetools/find-optimal/find-optimal

275 lines
7.1 KiB
Python
Executable file

#!/usr/bin/python3
'''=pod
=head1 NAME
find-optimal - Find optimal values
=head1 SYNOPSIS
B<find-optimal> I<command> [[options] I<value>] [options]
=head1 DESCRIPTION
B<find-optimal> will search for the optimal values for options.
It does this by calling I<command> with different values.
I<command> must as the last line output a value that says how bad it
was. B<find-optimal> will search for the smallest value.
=head1 EXAMPLE
Find the best value of pi - use 10 as a starting guess:
pi() { perl -e 'print(abs(shift()-3.14159265))' -- "$@"; }
export -f pi
find-optimal pi 10
Find the best values of pi and e - use 10 as starting guess for pi,
and 100 for e:
pie() { perl -e 'print(abs(shift()-3.14159265)+abs(shift()-2.718281828))' -- "$@"; }
export -f pie
find-optimal pie 10 100
=cut
#Find the fastest block size for dd
#
# mydd() { ( \time -f%e dd bs=$@ if=/dev/zero | head -c1G >/dev/null) 2>&1;}
# export -f mydd
# find-optimal mydd 10
# find-optimal -i --min 1 mydd 10
=pod
=head1 BUGS
B<find-optimal> used Nelder-Mead. Unfortunately, this only works well
for continous functions and not integers.
So currently B<find-optimal> will not work well for integers.
=head1 AUTHOR
Copyright (C) 2022 Ole Tange,
http://ole.tange.dk and Free Software Foundation, Inc.
=head1 LICENSE
Copyright (C) 2012 Free Software Foundation, Inc.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 3 of the License, or
at your option any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
=head1 DEPENDENCIES
B<vid> uses B<G>, and B<vlc>.
=head1 SEE ALSO
B<G>
=cut
'''
def add(x,y):
z = list(0 for i in x)
for i in range(0,len(x)):
z[i] = x[i] + y[i]
return (z)
def sub(x,y):
z = list(0 for i in x)
for i in range(0,len(x)):
z[i] = x[i] - y[i]
return (z)
def mul(x,y):
z = list(0 for i in y)
for i in range(0,len(y)):
z[i] = x * y[i]
return (z)
def nelder_mead(f, x_start,
step=0.001, no_improve_thr=10e-6,
no_improv_break=10, max_iter=0,
alpha=1., gamma=2., rho=-0.5, sigma=0.5):
'''
Pure Python implementation of the Nelder-Mead algorithm.
Reference: https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method
MIT Licence. Copyright (c) 2017 Matteo Maggioni
@param f (function): function to optimize, must return a scalar score
and operate over a numpy array of the same dimensions as x_start
@param x_start (numpy array): initial position
@param step (float): look-around radius in initial step
@no_improv_thr, no_improv_break (float, int): break after no_improv_break iterations with
an improvement lower than no_improv_thr
@max_iter (int): always break after this number of iterations.
Set it to 0 to loop indefinitely.
@alpha, gamma, rho, sigma (floats): parameters of the algorithm
(see Wikipedia page for reference)
return: tuple (best parameter array, best score)
'''
# init
dim = len(x_start)
prev_best = f(x_start)
no_improv = 0
res = [[x_start, prev_best]]
for i in range(dim):
x = mul(1,x_start)
x[i] = x[i] + step
score = f(x)
res.append([x, score])
# simplex iter
iters = 0
while 1:
# order
res.sort(key=lambda x: x[1])
best = res[0][1]
# break after max_iter
if max_iter and iters >= max_iter:
return res[0]
iters += 1
if best < prev_best - no_improve_thr:
no_improv = 0
prev_best = best
else:
no_improv += 1
if no_improv >= no_improv_break:
return res[0]
# centroid
x0 = [0.] * dim
for tup in res[:-1]:
for i, c in enumerate(tup[0]):
x0[i] += c / (len(res)-1)
# reflection
xr = add(x0, mul(alpha,(sub(x0, res[-1][0]))))
rscore = f(xr)
if res[0][1] <= rscore < res[-2][1]:
del res[-1]
res.append([xr, rscore])
continue
# expansion
if rscore < res[0][1]:
xe = add(x0, mul(gamma,(sub(x0, res[-1][0]))))
escore = f(xe)
if escore < rscore:
del res[-1]
res.append([xe, escore])
continue
else:
del res[-1]
res.append([xr, rscore])
continue
# contraction
xc = add(x0, mul(rho,(sub(x0, res[-1][0]))))
cscore = f(xc)
if cscore < res[-1][1]:
del res[-1]
res.append([xc, cscore])
continue
# reduction
x1 = res[0][0]
nres = []
for tup in res:
redx = add(x1, mul(sigma,(sub(tup[0],x1))))
score = f(redx)
nres.append([redx, score])
res = nres
def test():
import math
def f(x):
return math.sin(x[0]) * math.cos(x[1]) * (1. / (abs(x[2]) + 1))
print(nelder_mead(f, list([0., 20., 0.])))
if __name__ == "__main__":
import sys
import subprocess
import re
opt_verbose=True
def f(x):
# Run the shell command with the numeric values
# return the value of last line of output
run_arr=[]
i=0
for s in cmdtemplate:
if s == "0":
v = x[i];
try:
if opt_integer:
v = int(v)
if v < opt_min:
v = opt_min
except NameError:
True
run_arr.append(str(v))
i += 1
else:
run_arr.append(s)
run = ' '.join(run_arr)
try:
if opt_verbose:
print(run)
except NameError:
True
# Last line of stdout
try:
return float(subprocess.check_output(args=run,
shell=True,
executable='/bin/bash')
.splitlines()[-1])
except IndexError:
print("Command returns nothing.")
exit(1);
except subprocess.CalledProcessError:
exit(1);
cmdtemplate = []
values = []
# Convert:
# cmd -s 6 -f 7 -o foo
# =>
# [cmd -s 0 -f 0 -o foo],[6,7]
for s in sys.argv[1:]:
if re.match(r'^[-0-9.]+$',s):
values.append(s)
cmdtemplate.append("0")
else:
cmdtemplate.append(s)
print(nelder_mead(f, list([float(i) for i in values])))