DownhillSimplex.pm

package DownhillSimplex;

use strict;
use warnings;
use 5.008;

sub simplex {
    my ($init_vertices, $params, $coes, $eval_func, $disp_func) = @_;

    my $alpha = defined($coes->{alpha}) ? $coes->{alpha} : 1;
    my $gamma = defined($coes->{gamma}) ? $coes->{gamma} : 2;
    my $rho   = defined($coes->{rho})   ? $coes->{rho}   : 0.5;
    my $sigma = defined($coes->{sigma}) ? $coes->{sigma} : 0.5;
    my $maxiter = defined($coes->{maxiter}) ? $coes->{maxiter} : 100;
    my $epsilon = defined($coes->{epsilon}) ? $coes->{epsilon} : 1e-8;

    my $n = @$init_vertices - 1;
    my @vertices = ();
    foreach my $vertex (@$init_vertices) {
        my $score = $eval_func->($vertex, $params);
        push @vertices, { vertex => $vertex, score => $score };
    }

    my $iter = 0;
    while ($iter < $maxiter) {
        $iter += 1;
        @vertices = sort { $a->{score} <=> $b->{score} } @vertices;
        $disp_func->($iter, \@vertices, $params);

        my $grav = get_gravity(\@vertices);

        last if calc_distance($grav, $vertices[0]->{vertex}) < $epsilon;

        my $refl = get_reflection(\@vertices, $grav, $alpha);
        my $refl_score = $eval_func->($refl, $params);

        my $best_score = $vertices[0]->{score};
        my $second_worst_score = $vertices[$n - 1]->{score};
        my $worst_score = $vertices[$n]->{score};

        if ($best_score <= $refl_score && $refl_score < $second_worst_score) {
            pop @vertices;
            push @vertices, { vertex => $refl,  score => $refl_score};
            next;
        } elsif ($refl_score < $best_score) {
            my $exp = get_reflection(\@vertices, $grav, $gamma);
            my $exp_score = $eval_func->($exp, $params);
            pop @vertices;
            if ($exp_score < $refl_score) {
                push @vertices, { vertex => $exp,  score => $exp_score };
            } else {
                push @vertices, { vertex => $refl, score => $refl_score };
            }
            next;
        } elsif ($second_worst_score <= $refl_score) {
            my $cont = get_contraction(\@vertices, $grav, $rho);
            my $cont_score = $eval_func->($cont, $params);
            if ($cont_score < $worst_score) {
                pop @vertices;
                push @vertices, { vertex => $cont, score => $cont_score };
                next;
            }
        }

        my $best = $vertices[0]->{vertex};
        foreach my $i (1..$n) {
            my $vertex = $vertices[$i]->{vertex};
            my $l = @$vertex;
            foreach my $j (0..$l-1) {
                $vertex->[$j] = (1 - $sigma) * $best->[$j]
                              + $sigma * $vertex->[$j]; 
            }
        }
    }

    @vertices = sort { $a->{score} <=> $b->{score} } @vertices;
    my $best = $vertices[0]->{vertex};
    my $best_score = $vertices[0]->{score};

    ($best, $best_score);
}


sub get_gravity {
    my $vertices = shift;
    my $grav = [];
    my $n = @$vertices - 1;
    foreach my $i (0..$n-1) {
        my $vertex = $vertices->[$i]{vertex};
        my $l = @$vertex;
        die "" unless $n == $l;
        foreach my $j (0..$l-1) {
            $grav->[$j] += $vertex->[$j] / $n;
        }
    }
    $grav;
}

sub get_reflection {
    my ($vertices, $grav, $alpha) = @_;
    my $vertex = $vertices->[int(@$vertices)-1]{vertex};
    my $refl = [];
    foreach my $i (0..int(@$vertex)-1) {
        $refl->[$i] = (1 + $alpha) * $grav->[$i]
                    - $alpha * $vertex->[$i];
    }
    $refl;
}

sub get_contraction {
    my ($vertices, $grav, $rho) = @_;
    my $vertex = $vertices->[int(@$vertices)-1]{vertex};
    my $cont = [];
    foreach my $i (0..int(@$vertex)-1) {
        $cont->[$i] = (1 - $rho) * $vertex->[$i]
                    + $rho * $grav->[$i];
    }
    $cont;
}

sub calc_distance {
    my ($vertex1, $vertex2) = @_;
    my $distance = 0;
    foreach my $i (0..int(@$vertex1)-1) {
        $distance += ($vertex1->[$i] - $vertex2->[$i]) ** 2;
    }
    $distance = sqrt($distance);

    $distance;
}

1;