← Back to branch summary

~ubuntu-branches/ubuntu/lucid/pdl/lucid

« back to all changes in this revision

Viewing changes to Example/Simplex/tsimp2.pl

  • Committer: Bazaar Package Importer
  • Author(s): Ben Gertzfield
  • Date: 2002-04-08 18:47:16 UTC
  • Revision ID: james.westby@ubuntu.com-20020408184716-0hf64dc96kin3htp
Tags: upstream-2.3.2
Import upstream version 2.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Example/Simplex/tsimp2.pl
 
1
use PDL;
 
2
use PDL::Opt::Simplex;
 
3
#
 
4
# Simplex Demo - Alison Offer (aro@aaocbn.aao.gov.au)
 
5
#
 
6
#
 
7
# this is some sort of convergence criterion
 
8
$minsize = 1.e-6;
 
9
# max number of iterations ?
 
10
$maxiter = 100;
 
11
#
 
12
print " \n
 
13
        1: least squares gaussian fit to data + noise \n
 
14
           32 *exp (-((x-10)/6)^2) + noise
 
15
        2: minimise polynomial \n 
 
16
           (x-3)^2 + 2.*(x-3)*(y-2.5) + 3.*(y-2.5)^2 \n
 
17
        Please make your choice (1 or 2):";
 
18
chop($choice = <>);
 
19
if ($choice == 1) {
 
20
  print "Please enter noise factor (small number, 0-6):";
 
21
  chop($factor = <>);
 
22
#
 
23
# data : gaussian + noise
 
24
#
 
25
  foreach $j (0..19) {
 
26
       $data[$j] = 32*exp(-(($j-10)/6)**2) + 
 
27
                     $factor  * (rand() - 0.5);
 
28
  }
 
29
#
 
30
# initial guess - $initsize controls the size of the initial simplex (I think)
 
31
#
 
32
  $init = pdl [ 33, 9, 12 ];
 
33
  $initsize = 2;
 
34
#
 
35
#
 
36
  ($optimum,$ssize) = simplex($init,$initsize,$minsize,$maxiter,  
 
37
# this sub returns the function to be minimised.          
 
38
                            sub {my ($xv) =@_;
 
39
                                my $a = $xv->slice("(0)"); 
 
40
                                my $b = $xv->slice("(1)");                  
 
41
                                my $c = $xv->slice("(2)");   
 
42
                                my $sum = $a * 0.0;
 
43
                                foreach $j (0..19) {
 
44
                                    $sum += ($data[$j] -
 
45
                                     $a*exp(-(($j-$b)/$c)**2))**2;
 
46
                                }
 
47
                                return $sum;
 
48
                            });
 
49
                            
 
50
} else {
 
51
  $init = pdl [ 2 , 2 ];
 
52
  $initsize = 2;
 
53
  ($optimum,$ssize) = simplex($init,$initsize,$minsize,$maxiter,           
 
54
                            sub {my ($xv) =@_;
 
55
                                my $x = $xv->slice("(0)"); 
 
56
                                my $y = $xv->slice("(1)"); 
 
57
                                return ($x-3)**2 + 2.*($x-3)*($y-2.5) 
 
58
                                          + 3.*($y-2.5)**2; 
 
59
                            });
 
60
                             
 
61
 
 
62
}
 
63
print "OPTIMUM = $optimum \n";
 
64
print "SSIZE   = $ssize\n";
 
65