Math-Evol

 view release on metacpan or  search on metacpan

lua/Evol.lua  view on Meta::CPAN

--------------------------------------------------------------------
--     This Lua5 module is Copyright (c) 2010, Peter J Billam     --
--                       www.pjb.com.au                           --
--                                                                --
--  This module is free software; you can redistribute it and/or  --
--         modify it under the same terms as Lua5 itself.         --
--------------------------------------------------------------------
local M = {} -- public interface
M.Version = '1.13'
M.VersionDate = '27aug2010'

-- Example usage:
--local MM = require 'Evol'
--MM.bar()

----------------- infrastructure for evol ----------------
local function arr2txt(a) -- neat printing of arrays for debug use
	local txt = {}
	for k,v in ipairs(a) do txt[#txt+1] = string.format('%g',v); end
	return table.concat(txt,' ')
end
local function warn(str)
	io.stderr:write(str,'\n')
end
local function die(str)
	io.stderr:write(str,'\n')
	os.exit(1)
end
math.randomseed(os.time())
local gaussn_a = math.random()  -- reject 1st call to rand in case it's zero
local gaussn_b
local gaussn_flag = false
local function gaussn(standdev)
	-- returns normal distribution around 0.0 by the Box-Muller rules
	if not gaussn_flag then
		gaussn_a = math.sqrt(-2.0 * math.log(0.999*math.random()+0.001))
		gaussn_b = 6.28318531 * math.random()
		gaussn_flag = true
		return (standdev * gaussn_a * math.sin(gaussn_b))
	else
		gaussn_flag = false
		return (standdev * gaussn_a * math.cos(gaussn_b))
	end
end
--------------------------------------------------------------

-- various epsilons used in convergence testing ...
M.ea = 0.00000000000001;   -- absolute stepsize
M.eb = 0.000000001;        -- relative stepsize
M.ec = 0.0000000000000001; -- absolute error
M.ed = 0.00000000001;      -- relative error

function M.evol (xb ,sm, func,constrain, tm)
	if (type(xb) ~= 'table') then
		die "Evol.evol 1st arg must be a table\n";
	elseif (type(sm) ~= 'table') then
		die "Evol.evol 2nd arg must be a table\n";
	elseif (type(func) ~= 'function') then
		die "Evol.evol 3rd arg must be a function\n";
	elseif constrain and (type(constrain) ~= 'function') then
		die "Evol.evol 4th arg must be a function\n";
	end
	if not tm then tm = 10.0 end
	tm = math.abs(tm)

	local debug = false
	local n = #xb;   -- number of variables
	local clock_o = os.clock()
	local fc; local rel_limit    -- used to test convergence

	-- step-size adjustment stuff, courtesy Rechenberg ...
	local l = {}; local ten = 10
	local i; for i=1,ten do l[i] = n*i/5.0 end -- 6

lua/Evol.lua  view on Meta::CPAN

		elseif (defined $max[$i]) then
			push @sub_constr,"\tif (\$_[$i]>$max[$i]) { \$_[$i]=$max[$i]; }\n";
		end
		if (defined $min[$i] || defined $max[$i]) then
			some_constraints = some_constraints + 1
		end
		i = i + 1
	end
	push @sub_constr, "\treturn \@_;\n}\n";
	if debug then warn join ('', @sub_constr)."\n" end

	sub choose_best {
		local $xbref; local $linenum; @texts = {};
		while ($xbref = shift @_) do
			local @newtext = @text; local $i = $[;
			foreach $linenum (@linenums) do
				$newtext[$linenum] = $firstbits[$i] . string.format ('%g', $$xbref[$i])
				. $middlebits[$i];
				i = i + 1
			end
			push @texts, join ("\n", @newtext);
		end
		return &$func(@texts);
	end

	local ($xbref, $smref);
	if ($some_constraints) then
		eval join '', @sub_constr; if ($@) then die "text_evol: $@\n"; end
		($xbref, $smref) =
		 &select_evol(\@xb, \@sm, \&choose_best, \&constrain, $nchoices);
	else
		($xbref, $smref) = &select_evol(\@xb,\@sm,\&choose_best,0,$nchoices);
	end

	local @new_text = @text; $i = $[;
	foreach $linenum (@linenums) do
		$new_text[$linenum] = $firstbits[$i] . string.format ('%g', $$xbref[$i])
		. $middlebits[$i] . ' evol step '. string.format ('%g', $$smref[$i]);
		if (defined $min[$i]) then $new_text[$linenum] .= " min $min[$i]"; end
		if (defined $max[$i]) then $new_text[$linenum] .= " max $max[$i]"; end
		i = i + 1
	end
	if debug then warn   join ("\n", @new_text)."\n" end
	return join ("\n", @new_text)."\n";
]]
end

-- warn('Evol: M = '..tostring(M))
return M

--[[

=pod

=head1 NAME

Evol - Evolution search optimisation

=head1 SYNOPSIS

 local EV = require 'Evol'
 xb, sm, fb, lf = evol(xb, sm, function, constrain, tm)
 -- or
 xb, sm = select_evol(xb, sm, choose_best, constrain)

 -- not yet implemented:
 -- new_text = text_evol(text, choose_best_text, nchoices );

=head1 DESCRIPTION

This module implements the evolution search strategy.  Derivatives of
the objective function are not required.  Constraints can be incorporated.
The caller must supply initial values for the variables and for the
initial step sizes.

This evolution strategy is a random strategy, and as such is
particularly robust and will cope well with large numbers of variables,
or rugged objective funtions.

Evol.pm works either automatically (evol) with an objective function to
be minimised, or interactively (select_evol) with a (suitably patient)
human who at each step will choose the better of two possibilities.
Another subroutine (text_evol) allows the evolution of numeric parameters
in a text file, the parameters to be varied being identified in the text
by means of special comments.  A script I<ps_evol> which uses that is
included for human-judgement-based fine-tuning of drawings in PostScript.

Version 1.13

=head1 FUNCTIONS

=over 3

=item I<evol>( xb, sm, minimise, constrain, tm);

Where the arguments are:
 I<xb> the initial values of variables,
 I<sm> the initial values of step sizes,
 I<minimise> the function to be minimised,
 I<constrain> a function constraining the values,
 I<tm> the max allowable cpu time in seconds

The step sizes and the caller-supplied functions
I<function> and I<constrain> are discussed below.
The default value of I<tm> is 10 seconds.

I<evol> returns a list of four things:
 I<xb> the best values of the variables,
 I<sm> the final values of step sizes,
 I<fb> the best value of objective function,
 I<lf> a success parameter

I<lf> is set false if the search ran out of cpu time before converging.
For more control over the convergence criteria, see the
CONVERGENCE CRITERIA section below.

=item I<select_evol>( xb, sm, choose_best, constrain, nchoices);

Where the arguments are:
 I<xb> the initial values of variables,
 I<sm> the initial values of step sizes,



( run in 0.620 second using v1.01-cache-2.11-cpan-39bf76dae61 )