1# Library of McStas runtime perl functions
2#
3#   This file is part of the McStas neutron ray-trace simulation package
4#   Copyright (C) 1997-2006, All rights reserved
5#   Risoe National Laborartory, Roskilde, Denmark
6#   Institut Laue Langevin, Grenoble, France
7#
8#   This program is free software; you can redistribute it and/or modify
9#   it under the terms of the GNU General Public License as published by
10#   the Free Software Foundation; version 2 of the License.
11#
12#   This program is distributed in the hope that it will be useful,
13#   but WITHOUT ANY WARRANTY; without even the implied warranty of
14#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15#   GNU General Public License for more details.
16#
17#   You should have received a copy of the GNU General Public License
18#   along with this program; if not, write to the Free Software
19#   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
20#
21use Config;
22
23use Math::Amoeba qw(MinimiseND);
24require "mccode_config.perl";
25
26# Overload with user's personal config
27if ($ENV{"HOME"} && -e $ENV{"HOME"}."/.mcstas/mccode_config.perl") {
28  require $ENV{"HOME"}."/.mcstas/mccode_config.perl";
29}
30
31our $optim_bestvalue=0;
32our $optim_variables;
33our @optim_datablock = ();
34our @optim_first;
35our @optim_youts=();
36
37# must run optimization steps as scans in oder to extract Detector lines
38# from output: (mcrun::do_scan) mcrun -N 1
39sub minimize_function {
40  my (@p) = @_; # array of instrument parameters at optimization step
41  my $j;
42  # function must return a value to be maximized
43  $y = 0;
44
45  my @scanned = ();
46  my @minval = ();
47  my @maxval = ();
48
49  my $out = "$optim_iterations ";
50  my @youts=();
51
52  # create a copy of the global $scan_info within limits
53  for($j = 0; $j < @{$scan_info->{VARS}}; $j++) {
54    if    ($p[$j] > $scan_info->{MAX}[$j]) { $p[$j] = $scan_info->{MAX}[$j]; }
55    elsif ($p[$j] < $scan_info->{MIN}[$j]) { $p[$j] = $scan_info->{MIN}[$j]; }
56    $minval[$j]  = $p[$j];
57    $maxval[$j]  = $p[$j];
58    $scanned[$j] = $scan_info->{VARS}[$j];
59    $out .= "$p[$j] "; # parameter values for this iteration step
60  }
61
62  # assemble $scan_info from optim_info and optimization step ($numpoints=1)
63  my $optim_info = { VARS => \@scanned, MIN => \@minval, MAX => \@maxval };
64
65  # execute single iteration step (do_scan). returns [params_val I err ... I err]
66  ($datablock, $variables, @youts) = do_scan($optim_info);
67  my @vars = split(" ", $variables);
68  my @vals = split(" ", $datablock);
69
70  # search optim_names in variables
71  # get list of monitors to maximize and loop
72  my $found_monitor=0;
73  for ($j = @{$scan_info->{VARS}}; $j < @vars; $j += 2) {
74    if ($optim_iterations == 0) { $optim_first[$j] = 0; }
75    my $value = $vals[$j];
76    if ($optim_first[$j] == 0 && $value != 0) {
77      $optim_first[$j] = abs($value);
78    }
79    if ($optim_first[$j] > 0) {
80      $value /= $optim_first[$j];
81    }
82    if ($optim_flag > 1) { # all monitors
83      $y = $y + $value; # add all values for criteria
84      $found_monitor++;
85    } else { # selected monitors
86      my $i;
87      for($i = 0; $i < @optim_names; $i++) {
88        my $this_name = $optim_names[$i];
89        # add each value of monitor to $y
90        if ($vars[$j] eq "$this_name" . "_I") {
91          $y = $y + $value; # add corresponding value to found name for criteria
92          $found_monitor++;
93        }
94      }
95    }
96    $out .= " $vals[$j] $vals[$j+1]";
97  } # end for j
98  die "optimization: ERROR: selected component is not a monitor\n" unless $found_monitor;
99  $out .= " $y 0";
100  push @optim_datablock, "$out\n";
101  if ($optim_iterations == 0) {
102    $optim_variables = "Point " . $variables;
103    @optim_youts = @youts;
104    $optim_variables .= " criteria null";
105    push @optim_youts, " (criteria,null)";
106  }
107  if ($optim_iterations == 0 || $y >= $optim_bestvalue) {
108    @optim_best      = @p;
109    $optim_bestvalue = $y;
110  }
111  $optim_iterations++;
112
113  # Handle writing of mcplot-compatible optim output file
114  $numpoints = $optim_iterations;
115  if (!$optfile) {
116    if ($MCSTAS::mcstas_config{'TEMP'} ne "no") {
117      require File::Temp;
118      ($OPT, $optfile) = File::Temp::tempfile("mcoptim"."_XXXX", SUFFIX => '.dat'); # throw file handle
119    } else {
120      $optfile = int(rand(9999));
121      $optfile = "mcoptim"."_$optfile.dat";
122    }
123  } else {
124    $OPT = new FileHandle;
125  }
126  open($OPT, ">$optfile");
127  autoflush $OPT 1;
128  push @opt_out, "$optim_iterations $y 0\n";
129  output_dat_header($OPT, "# ",$scan_info, '(criteria,null)', 'Point criteria null', $optfile);
130  print $OPT "@opt_out\n";
131  close($OPT);
132
133  print "Optimization $optim_iterations criteria=$y\n";
134
135  return -$y;
136}
137
138sub sighandler_optim {
139  my $signame = shift;
140  print "mcrun: Recieved signal $signame during optimization (iteration $optim_iterations).\n";
141  if ($signame eq "INT" || $signame eq "TERM") {
142    save_optimization(0);
143    exit(0);
144  } else {
145    print "# Optimization history:\n";
146    my $datablock = join(" ", @optim_datablock);  # mcoptimlib-> @optim_datablock, @optim_youts, $optim_variable
147    print "# variables: $optim_variables\n";
148    print "$datablock\n";
149  }
150}
151
152# function to save optimization history
153sub save_optimization {
154  my ($final_call) = @_;
155  my ($j);
156
157  # constrain within limits
158  for($j = 0; $j < @{$scan_info->{VARS}}; $j++) {
159    if    ($optim_best[$j] > $scan_info->{MAX}[$j]) { $optim_best[$j] = $scan_info->{MAX}[$j]; }
160    elsif ($optim_best[$j] < $scan_info->{MIN}[$j]) { $optim_best[$j] = $scan_info->{MIN}[$j]; }
161  }
162
163  # display result
164  print "Iteration $optim_iterations parameters:\n";
165
166  for($j = 0; $j < @{$scan_info->{VARS}}; $j++) {
167    $i = $scan_info->{VARS}[$j]; # Index of variable to be scanned
168    print "$params[$i]=$optim_best[$j] ";
169  }
170  print "\n";
171
172  # reactivate $data_dir and relaunch minimize_function to fill directory
173  # mcstas.sim is last iteration result
174  if ($data_dir_saved) {
175    $data_dir = $data_dir_saved;
176    print "Generating optimized configuration in --dir=$data_dir\n";
177    minimize_function(@optim_best); # this also sets scan-info by calling do_scan
178  }
179  # output optimization results as a scan
180  $numpoints = $optim_iterations;
181  my $prefix = "";
182  if($data_dir) { $prefix = "$data_dir/"; }
183  my $datfile = "mcstas.dat";
184  my $simfile = $datfile;
185  $simfile =~ s/\.dat$//;         # Strip any trailing ".dat" extension ...
186  $simfile .= $format_ext;        # ... and add ".sim|m|sci" extension.
187  my $datablock = join(" ", @optim_datablock);  # mcoptimlib-> @optim_datablock, @optim_youts, $optim_variable
188
189  # Only initialize / use $DAT datafile if format is PGPLOT
190  if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /PGPLOT|McStas/i) {
191    my $DAT = new FileHandle;
192    open($DAT, ">${prefix}$datfile");
193    output_dat_header($DAT, "# ", $scan_info, \@optim_youts, $optim_variables, $datfile);
194    print $DAT "$datablock\n";
195    close($DAT);
196  } else {
197    $datfile = $simfile;             # Any reference should be to the simfile
198  }
199
200  output_sim_file("${prefix}$simfile", $scan_info, \@optim_youts, $optim_variables, $datfile, $datablock);
201
202  print "Optimization file: '${prefix}$datfile'\nOptimized parameters: $optim_variables\n";
203  undef($data_dir);
204  $numpoints = 1;
205}
206
2071;
208