File Coverage

blib/lib/Math/Evol.pm

Criterion	Covered	Total	%
statement	132	234	56.4
branch	55	122	45.0
condition	8	21	38.1
subroutine	5	9	55.5
pod	4	7	57.1
total	204	393	51.9

# Math::Evol.pm
#########################################################################
#        This Perl module is Copyright (c) 2002, Peter J Billam         #
#               c/o P J B Computing, www.pjb.com.au                     #
#                                                                       #
#     This module is free software; you can redistribute it and/or      #
#            modify it under the same terms as Perl itself.             #
#########################################################################

package Math::Evol;
no strict;
$VERSION = '1.13';
# gives a -w warning, but I'm afraid $VERSION .= ''; would confuse CPAN
require Exporter;
@ISA = qw(Exporter);
@EXPORT = qw(evol select_evol text_evol);
@EXPORT_OK = qw(arr2txt gaussn);

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

sub new {
	my $arg1 = shift;
	my $class = ref($arg1) || $arg1; # can be used as class or instance method
	my $self  = {};   # ref to an empty hash
	bless $self, $class;
	$self->_initialise();
	return $self;
}


sub evol { my ( $xbref,$smref, $func_ref,$constrain_ref, $tm) = @_;
	if (ref $xbref ne 'ARRAY') {
		die "Math::Evol::evol 1st arg must be an array ref\n";
	} elsif (ref $smref ne 'ARRAY') {
		die "Math::Evol::evol 2nd arg must be an array ref\n";
	} elsif (ref $func_ref ne 'CODE') {
		die "Math::Evol::evol 3rd arg must be a function ref\n";
	} elsif ($constrain_ref && (ref $constrain_ref) ne 'CODE') {
		die "Math::Evol::evol 4th arg must be a function ref\n";
	}
	my @xb = @$xbref;
	my @sm = @$smref;
	if (! $tm) { $tm = 10.0; }
	$tm = abs $tm;

	my $debug = 0;
	my $n = scalar @xb;      # number of variables
	my $nm1 = $n + $[ - 1;   # last index
	my ($usr_o, $sys_o) = times;
	my ($fc, $rel_limit);    # used to test convergence

	# step-size adjustment stuff, courtesy Rechenberg ...
	my @l; foreach (1 .. 10) { push @l, $n*$_/5.0; } # 6
	my $le = $n+$n; my $lm = 0; my $lc = 0;

	if ($constrain_ref) { @xb = &$constrain_ref(@xb); }
	$fb = &$func_ref(@xb); $fc = $fb;
	while (1) {
		foreach $i ($[ .. $nm1) { $x[$i] = $xb[$i]+&gaussn($sm[$i]); }
		if ($constrain_ref) {
			@x = &$constrain_ref(@x);
			warn 'new @x is '.&arr2txt(@x) if $debug;
		}
		$ff = &$func_ref (@x);   # find new objective function
		if ($ff <= $fb) { $le++; $fb=$ff; @xb=@x; }
		$lm++;  if ($lm < $n) { next; }

   		# adjust the step sizes ...
		my $k = 0.85 ** ($n+$n <=> $le-$l[$[]);
		foreach $i ($[ .. $nm1) {
			$sm[$i] *=  $k;
			$rel_limit = abs ($eb * $xb[$i]);
			if ($sm[$i] < $rel_limit) { $sm[$i] = $rel_limit; }
			if ($sm[$i] < $ea) { $sm[$i] = $ea; }
		}
		# could i_l%10 the index as in lua, rather than shift and push as here
		shift @l; push (@l, $le); $lm = 0;
		if ($debug) {
			warn "le=$le l[0]=$l[$[] n=$n k=$k\n";
			warn "new step sizes sm = ".&arr2txt(@sm);
			warn "new l = ".&arr2txt(@l);
		}

		if (defined $ea) { $ea = abs $ea; }
		if (defined $eb) { $eb = abs $eb; }
		if (defined $ec) { $ec = abs $ec; }
		if (defined $ed) { $ed = abs $ed; }
		warn sprintf("fb=%g fc=%g ec=%g ed=%g\n",$fb,$fc,$ec,$ed) if $debug;
		$lc++;
		if ($lc < 25) {
			my ($usr_n, $sys_n) = times;
			if (($usr_n+$sys_n - $usr_o-$sys_o) > $tm) {   # out of time ?
				warn "exceeded $tm seconds\n" if $debug;
				return (\@xb, \@sm, $fb, 0);   # return best current value
			} else { next;
			}
		}
		if ((defined $ec) && (($fc-$fb) <= $ec)) {
			warn "converged absolutely\n" if $debug;
			return (\@xb, \@sm, $fb, 1);   # 29
		}
		if ((defined $ed) && (($fc-$fb)/$ed <= abs($fc))) {
			warn "converged relativelty\n" if $debug;
			return (\@xb, \@sm, $fb, 1);
		}
		$lc = 0; $fc = $fb; next;
	}
}

sub select_evol { my ($xbref,$smref,$func_ref,$constrain_ref,$nchoices) = @_;
	if (ref $xbref ne 'ARRAY') {
		die "Math::Evol::select_evol 1st arg must be an array ref\n";
	} elsif (ref $smref ne 'ARRAY') {
		die "Math::Evol::select_evol 2nd arg must be an array ref\n";
	} elsif (ref $func_ref ne 'CODE') {
		die "Math::Evol::select_evol 3rd arg must be a function ref\n";
	} elsif ($constrain_ref && (ref $constrain_ref) ne 'CODE') {
		die "Math::Evol::select_evol 4th arg must be a function ref\n";
	}
	my @xb = @$xbref;
	my @sm = @$smref;
	if (! $nchoices) { $nchoices = 1; }

	my $debug = 0;
	my $n = scalar @xb;      # number of variables
	# warn and return if $n==0,  else $test_every=0 and we get a div by 0
	if ($n==0) { die("Math::Evol::select_evol 1st arg was an empty array\n"); }
	my $nm1 = $n + $[ - 1;   # last index
	my ($choice, $continue);

	# step-size adjustment stuff, courtesy Rechenberg
	# modified from Schwefel's sliding flat window average to an EMA
	my $desired_success_rate = 1.0 - 0.8**$nchoices;
	my $success_rate = $desired_success_rate;
	# test every 5*$n binary choices equivalent - 10*$n is just too slow.
	my $lm = 0;
	my $test_every = 5 * $n * (log 2) / (log ($nchoices+1));
	my $ema = exp (-1.0 / $test_every);
	my $one_m_ema = 1.0 - $ema;
	warn sprintf
	("n=$n nchoices=$nchoices success_rate=%g test_every=%g ema=%g\n",
	 $success_rate, $test_every, $ema) if $debug;

	if ($constrain_ref) { @xb = &$constrain_ref(@xb); }
	my @func_args;  # array of refs to arrays
	while () {
		my @x;  #XXX?
		@func_args = ( \@xb );
		foreach (1 .. $nchoices) {
			foreach $i ($[ .. $nm1) {
				$x[$i] = $xb[$i]+&gaussn($sm[$i]);
			}
			if ($constrain_ref) { @x = &$constrain_ref(@x); }
			push @func_args, [ @x ];
		}
		($choice, $continue) = &$func_ref(@func_args);
		if ($choice > $[) {
			@xb = @{$func_args[$choice]};
			$success_rate = $one_m_ema + $ema*$success_rate;
		} else {
			$success_rate = $ema*$success_rate;
		}
		warn "success_rate=$success_rate\n" if $debug;
		if (!$continue) { return (\@xb, \@sm); }
		$lm++; if ($lm < $test_every) { next; }

   	# adjust the step sizes ...
		my $k = 0.85 ** ($desired_success_rate <=> $success_rate);
		warn "success_rate=$success_rate k=$k\n" if $debug;
		foreach $i ($[ .. $nm1) {
			$sm[$i] *=  $k;
			$rel_limit = abs ($eb * $xb[$i]);
			if ($sm[$i] < $rel_limit) { $sm[$i] = $rel_limit; }
			if ($sm[$i] < $ea) { $sm[$i] = $ea; }
		}
		warn "new step sizes sm = ".&arr2txt(@sm) if $debug;
		$lm = 0;
	}
}

sub text_evol { my ($text, $nchoices); local ($func_ref);
	($text, $func_ref, $nchoices) = @_;
	return unless $text;
	if (ref $func_ref ne 'CODE') {
		die "Math::Evol::text_evol 2nd arg must be a function ref\n";
	}
	if (! $nchoices) { $nchoices = 1; }

	my $debug = 0;
	local @text = split ("\n", $text);
	my ($linenum,$m,@xb,@sm,@min,@max) = ($[,0);
	local (@linenums, @firstbits, @middlebits, $lastbit);
	my $n = $[ - 1; foreach (@text) {
		if (/^(.*?)(-?[\d.]+)(\D+)evol\s+step\s+([\d.]+)(.*)$/) {
			$n++; $linenums[$n] = $linenum; $firstbits[$n]=$1;
			$xb[$n]=$2; $middlebits[$n]=$3; $sm[$n]=$4; $lastbit = $5;
			if ($lastbit =~ /min\s+([-\d.]+)/) { $min[$n] = $1; }
			if ($lastbit =~ /max\s+([-\d.]+)/) { $max[$n] = $1; }
		}
		$linenum++;
	}
	warn "xb = ".&arr2txt(@xb) if $debug;
	warn "sm = ".&arr2txt(@sm) if $debug;

	# construct the constraint routine
	my $some_constraints = 0;
	my @sub_constr = ("sub constrain {\n");
	my $i = $[; while ($i <= $n) {
		if (defined $min[$i] && defined $max[$i]) {
		  push @sub_constr,"\tif (\$_[$i]>$max[$i]) { \$_[$i]=$max[$i];\n";
		  push @sub_constr,"\t} elsif (\$_[$i]<$min[$i]) {\$_[$i]=$min[$i];\n";
		  push @sub_constr,"\t}\n";
		} elsif (defined $min[$i]) {
			push @sub_constr,"\tif (\$_[$i]<$min[$i]) { \$_[$i]=$min[$i]; }\n";
		} elsif (defined $max[$i]) {
			push @sub_constr,"\tif (\$_[$i]>$max[$i]) { \$_[$i]=$max[$i]; }\n";
		}
		if (defined $min[$i] || defined $max[$i]) { $some_constraints++; }
		$i++;
	}
	push @sub_constr, "\treturn \@_;\n}\n";
	warn join ('', @sub_constr)."\n" if $debug;

	sub choose_best {
		my $xbref; my $linenum; @texts = ();
		while ($xbref = shift @_) {
			my @newtext = @text; my $i = $[;
			foreach $linenum (@linenums) {
				$newtext[$linenum] = $firstbits[$i].sprintf('%g',$$xbref[$i])
				. $middlebits[$i];
				$i++;
			}
			push @texts, join ("\n", @newtext);
		}
		return &$func_ref(@texts);
	}

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

	my @new_text = @text; $i = $[;
	foreach $linenum (@linenums) {
		$new_text[$linenum] = $firstbits[$i] . sprintf ('%g', $$xbref[$i])
		. $middlebits[$i] . ' evol step '. sprintf ('%g', $$smref[$i]);
		if (defined $min[$i]) { $new_text[$linenum] .= " min $min[$i]"; }
		if (defined $max[$i]) { $new_text[$linenum] .= " max $max[$i]"; }
		$i++;
	}
	warn   join ("\n", @new_text)."\n" if $debug;
	return join ("\n", @new_text)."\n";
}

# --------------- infrastructure for evol ----------------

sub arr2txt { # neat printing of arrays for debug use
	my @txt; foreach (@_) { push @txt, sprintf('%g',$_); }
	return join (' ',@txt)."\n";
}
my $gaussn_a = rand;  # reject 1st call to rand in case it's zero
my $gaussn_b;
my $gaussn_flag;
sub gaussn {   my $standdev = $_[$[];
	# returns normal distribution around 0.0 by the Box-Muller rules
	if (! $gaussn_flag) {
		# $gaussn_a = sqrt(-2.0 * log(rand)); BUG #44777: Log of zero error
		$gaussn_a = sqrt(-2.0 * log(rand(0.999)+0.001));  # 1.12
		$gaussn_b = 6.28318531 * rand;
		$gaussn_flag = 1;
		return ($standdev * $gaussn_a * sin($gaussn_b));
	} else {
		$gaussn_flag = 0;
		return ($standdev * $gaussn_a * cos($gaussn_b));
	}
}
1;

__END__