heikki/fitting/fit.pl

   1 #!/usr/bin/perl -w
   2 # fit.c - experiments in curve fitting
   3 # for pazpar'2 ranking normalizing
   4
   5 # We have a number of data points ( position, score) from
   6 # different sources. The task is to normalize them so that
   7 # they all fall near the curve y=1/p, where p is the position
   8 # This is done by adjusting the ranks R so that Rn = aR+b
   9 # We need to find parameters a,b so as to minimize the chi-
  10 # squared difference from y=1/p
  11
  12
  13 my $plotnr = 1; # number the tmp files for plotting
  14 my $plotcmd = ""; # the plot commands for gnuplot
  15
  16 # Calculate the (squared) difference from the normalized rank to the 1/n function
  17 # Params
  18 #  p = position (x)
  19 #  r = rank, not normalized
  20 #  a,b normalizing params
  21 sub diff {
  22     my ( $p, $r, $a, $b ) = @_;
  23     my $rn = $r * $a + $b;
  24     my $f = 1.0 / $p;  # target value
  25     my $d = $rn - $f;
  26     return $d * $d;
  27 }
  28
  29 # Read and process one data file
  30 # Just one float number per line, nothing else
  31 sub onefile {
  32     my $fn = shift;
  33     my @d;
  34     open F, $fn or die "Could not open $fn: $!\n";
  35     my $n = 1; # number of data points
  36     my $first;
  37     my $last;
  38     my $title;
  39     while ( <F> ){
  40         chomp();
  41         $title = $_ unless defined($title);
  42         next unless /^[0-9]/; # skip comments etc
  43         my $v = 1.0 * $_ ;
  44         $first = $v unless defined($first);
  45         $last = $v;
  46         #print "Data $n is $v\n";
  47         $d[$n++] = $v;
  48     }
  49     $title =~ s/^[# ]+//; # clean the '#' and leading space
  50     print "$fn: '$title' $n points: $first - $last \n";
  51     # Initial guess Rn = a*R + b
  52     my $a = 1.0 / $first;
  53     my $b = - $last;
  54     # step sizes for a and b
  55     my $da = $a / 3;
  56     my $db = - $b / 3;
  57     my $iteration = 0;
  58     my $prev = 0.0;
  59     while (1) {
  60         $iteration++;
  61         # 5 sums: at (a,b) (a+,b), (a-,b), (a,b+), (a,b-)
  62         my $sab = 0.0; # at a,b
  63         my $sap = 0.0; # at a+da,b
  64         my $sam = 0.0; # at a-da,n
  65         my $sbp = 0.0; # at a, b+db
  66         my $sbm = 0.0; # at a, b-db
  67         for ( my $p = 1 ; $p < $n; $p++ ) {
  68             $sab += diff( $p, $d[$p], $a, $b );
  69             $sap += diff( $p, $d[$p], $a+$da, $b );
  70             $sam += diff( $p, $d[$p], $a-$da, $b );
  71             $sbp += diff( $p, $d[$p], $a, $b+$db );
  72             $sbm += diff( $p, $d[$p], $a, $b-$db );
  73         }
  74         my $dif = $sab - $prev;
  75         #print "iteration $iteration: a=$a +- $da   b=$b +- $db chisq=$sab dif=$dif\n";
  76         if ( (abs($da) < abs($a)/100.0 && abs($db) < abs($b)/100.0) ||
  77              ($iteration >= 100 ) ||
  78              (abs($dif) < 0.00001 ) ) {
  79             print "it-$iteration: a=$a +- $da   b=$b +- $db chisq=$sab dif=$dif\n";
  80             last;
  81         }
  82         $prev = $sab;
  83         # adjust a
  84         if ( $sap < $sab && $sap < $sam ) {
  85             $a += $da;
  86         } elsif ( $sam < $sab && $sam < $sap ) {
  87             $a -= $da;
  88         } else {
  89             $da = $da /2;
  90         }
  91         $da = $da * 0.99;
  92         # adjust b
  93         if ( $sbp < $sab && $sbp < $sbm ) {
  94             $b += $db;
  95         } elsif ( $sbm < $sab && $sbm < $sbp ) {
  96             $b -= $db;
  97         } else {
  98             $db = $db /2;
  99         }
 100         $db = $db * 0.99;
 101     }
 102
 103     # plot the file
 104     my $pf = "/tmp/plot.$plotnr.data";
 105     $plotnr++;
 106     open PF, ">$pf" or die "Could not open plot file $pf: $!\n";
 107     for ( my $p = 1 ; $p < $n; $p++ ) {
 108         my $rn = $d[$p] * $a + $b;
 109         print PF "$p $rn\n";
 110     }
 111     close PF;
 112     $plotcmd .= "," if ($plotcmd);
 113     $plotcmd .= "\"$pf\" using 1:2 with points title \"$title\"";
 114
 115
 116
 117 }
 118
 119 # main
 120
 121 if ( !defined($ARGV[0]) ) {
 122     die "Need at least one file to plot\n";
 123 }
 124 while ($ARGV[0]) {
 125   onefile( $ARGV[0] );
 126   shift(@ARGV);
 127 }
 128 my $cmd =
 129     "set term png\n" .
 130     "set out \"plot.png\" \n" .
 131     "plot $plotcmd \n";
 132
 133 print "$cmd \n";
 134
 135 open GP, "| gnuplot" or die "Could not open a pipe to gnuplot: $!\n";
 136 print GP $cmd;
 137 close GP;