1 2# 3# GENERATED WITH PDL::PP! Don't modify! 4# 5package PDL::Minuit; 6 7@EXPORT_OK = qw( mn_init mn_def_pars mn_excm mn_pout mn_stat mn_err mn_contour mn_emat PDL::PP mninit PDL::PP mn_abre PDL::PP mn_cierra PDL::PP mnparm PDL::PP mnexcm PDL::PP mnpout PDL::PP mnstat PDL::PP mnemat PDL::PP mnerrs PDL::PP mncont ); 8%EXPORT_TAGS = (Func=>[@EXPORT_OK]); 9 10use PDL::Core; 11use PDL::Exporter; 12use DynaLoader; 13 14 15 16 17 @ISA = ( 'PDL::Exporter','DynaLoader' ); 18 push @PDL::Core::PP, __PACKAGE__; 19 bootstrap PDL::Minuit ; 20 21 22 23 24=head1 NAME 25 26PDL::Minuit -- a PDL interface to the Minuit library 27 28=head1 DESCRIPTION 29 30This package implements an interface to the Minuit minimization routines (part 31of the CERN Library) 32 33=head1 SYNOPSIS 34 35A basic fit with Minuit will call three functions in this package. First, a basic 36initialization is done with mn_init(). Then, the parameters are defined via 37the function mn_def_pars(), which allows setting upper and lower bounds. Then 38the function mn_excm() can be used to issue many Minuit commands, including simplex 39and migrad minimization algorithms (see Minuit manual for more details). 40 41See the test file minuit.t in the test (t/) directory for a basic example. 42 43 44 45 46 47 48 49=head1 FUNCTIONS 50 51 52 53=cut 54 55 56 57 58 59 60# Package variable 61my $mn_options; 62 63 64 65 66sub mn_init{ 67 my $fun_ref = shift; 68 69 $mn_options = { Log => undef, 70 Title => 'Minuit Fit', 71 N => undef, 72 Unit => undef, 73 Function => $fun_ref, 74 }; 75 76 if ( @_ ){ 77 my $args = $_[0]; 78 for my $key (qw/ Log Title Unit/){ 79 $mn_options->{$key} = $args->{$key} if exists $args->{$key}; 80 } 81 } 82 83 # Check if there was a valid F77 available and barf 84 # if there was not and the user is trying to pass Log 85 86 if (defined($mn_options->{Log})) { 87 $mn_options->{Unit} = 88 unless defined $mn_options->{Unit}; 88 } 89 else { $mn_options->{Unit} = 6; } 90 91 if (defined (my $logfile = $mn_options->{Log})){ 92 if (-e $logfile) { unlink $logfile; } 93 PDL::Minuit::mn_abre($mn_options->{Unit},$logfile,'new'); 94 print STDERR "# Opening file $logfile....\n"; 95 } 96 97 PDL::Minuit::mninit(5,$mn_options->{Unit},$mn_options->{Unit}); 98 PDL::Minuit::mnseti($mn_options->{Title}); 99 100 if (defined (my $logfile = $mn_options->{Log})){ 101 PDL::Minuit::mn_cierra($mn_options->{Unit}); 102 } 103 104} 105 106 107 108 109=head2 mninit 110 111=for sig 112 113 Signature: (int a();int b(); int c()) 114 115 116=for ref 117 118info not available 119 120 121=for bad 122 123mninit does not process bad values. 124It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. 125 126 127=cut 128 129 130 131 132 133 134*mninit = \&PDL::Minuit::mninit; 135 136 137 138 139 140=head2 mn_abre 141 142=for sig 143 144 Signature: (int l(); char* nombre; char* mode) 145 146 147=for ref 148 149info not available 150 151 152=for bad 153 154mn_abre does not process bad values. 155It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. 156 157 158=cut 159 160 161 162 163 164 165*mn_abre = \&PDL::Minuit::mn_abre; 166 167 168 169 170 171=head2 mn_cierra 172 173=for sig 174 175 Signature: (int l()) 176 177 178=for ref 179 180info not available 181 182 183=for bad 184 185mn_cierra does not process bad values. 186It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. 187 188 189=cut 190 191 192 193 194 195 196*mn_cierra = \&PDL::Minuit::mn_cierra; 197 198 199 200 201sub mn_def_pars{ 202 my $pars = shift; 203 my $steps = shift; 204 205 my $n = nelem($pars); 206 $mn_options->{N} = $n; 207 208 #print "Unit :".$mn_options->{Unit}."\n"; 209 210 my @names = (); 211 for (my $i=0; $i < $n; $i++) { $names[$i] = "Par_$i"; } 212 my $lo_bounds = zeroes($n); 213 my $up_bounds = zeroes($n); 214 215 if ( @_ ) { 216 my $opts = $_[0]; 217 $lo_bounds = $opts->{Lower_bounds} if defined $opts->{Lower_bounds}; 218 $up_bounds = $opts->{Upper_bounds} if defined $opts->{Upper_bounds}; 219 if (defined($opts->{Names})){ 220 $names_t = $opts->{Names}; 221 barf " Names has to be an array reference" unless ref($names_t) eq 'ARRAY'; 222 @names = @$names_t; 223 barf " Names has to have as many elements as there are parameters " unless ( @names == $n); 224 } 225 } 226 227 my $iflag = 0; 228 229 if (defined (my $logfile = $mn_options->{Log})){ 230 PDL::Minuit::mn_abre($mn_options->{Unit},$logfile,'old'); 231 } 232 233 foreach my $i ( 0..(nelem($pars)-1) ){ 234 my $ii = $i + 1; 235 $iflag = PDL::Minuit::mnparm($ii,$pars->slice("($i)"), 236 $steps->slice("($i)"), 237 $lo_bounds->slice("($i)"), 238 $up_bounds->slice("($i)"), 239 $names[$i]); 240 barf "Problem initializing parameter $i in Minuit " unless ($iflag == 0); 241 } 242 243 if (defined (my $logfile = $mn_options->{Log})){ 244 PDL::Minuit::mn_cierra($mn_options->{Unit}); 245 } 246} 247 248 249 250 251=head2 mnparm 252 253=for sig 254 255 Signature: (int a(); double b(); double c(); double d(); double e(); int [o] ia(); char* str) 256 257 258=for ref 259 260info not available 261 262 263=for bad 264 265mnparm does not process bad values. 266It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. 267 268 269=cut 270 271 272 273 274 275 276*mnparm = \&PDL::Minuit::mnparm; 277 278 279 280 281sub mn_excm{ 282 my $command = shift; 283 284 my $fun_ref = $mn_options->{Function}; 285 286 my ($arglis,$narg); 287 if ( @_ ) { $arglis = shift; $narg = nelem($arglis);} 288 else { $arglis = pdl(0); $narg = 0; } 289 290 if ( @_ ) { barf "Usage : mn_excm($command, [$arglis]) \n"; } 291 292 if (defined (my $logfile = $mn_options->{Log})){ 293 PDL::Minuit::mn_abre($mn_options->{Unit},$logfile,'old'); 294 } 295 296 my $iflag = pdl(0); 297 298 299 $iflag = PDL::Minuit::mnexcm($arglis, $narg, $command, $fun_ref,$mn_options->{N}); 300 warn "Problem executing command '$command' " unless ($iflag == 0); 301 302 if (defined (my $logfile = $mn_options->{Log})){ 303 PDL::Minuit::mn_cierra($mn_options->{Unit}); 304 } 305 306 return $iflag; 307} 308 309 310 311 312=head2 mnexcm 313 314=for sig 315 316 Signature: (double a(n); int ia(); int [o] ib(); char* str; SV* function; int numelem) 317 318 319=for ref 320 321info not available 322 323 324=for bad 325 326mnexcm does not process bad values. 327It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. 328 329 330=cut 331 332 333 334 335 336 337*mnexcm = \&PDL::Minuit::mnexcm; 338 339 340 341 342 sub mn_pout{ 343 barf "Usage: mn_pout(par_number)" unless ($#_ == 0); 344 my $par_num = shift; 345 my $n = $mn_options->{N}; 346 if (($par_num < 1) || ($par_num > $n)) { barf "Parameter numbers range from 1 to $n "; } 347 348 if (defined (my $logfile = $mn_options->{Log})){ 349 PDL::Minuit::mn_abre($mn_options->{Unit},$logfile,'old'); 350 } 351 352 my $val = pdl(0); 353 my $err = pdl(0); 354 my $bnd1 = pdl(0); 355 my $bnd2 = pdl(0); 356 my $ivarbl = pdl(0); 357 my $par_name = " "; 358 PDL::Minuit::mnpout($par_num,$val,$err,$bnd1,$bnd2,$ivarbl,\$par_name); 359 360 if (defined (my $logfile = $mn_options->{Log})){ 361 PDL::Minuit::mn_cierra($mn_options->{Unit}); 362 } 363 364 return ($val,$err,$bnd1,$bnd2,$ivarbl,$par_name); 365 } 366 367 368 369 370=head2 mnpout 371 372=for sig 373 374 Signature: (int ia(); double [o] a(); double [o] b(); double [o] c(); double [o] d();int [o] ib(); SV* str) 375 376 377=for ref 378 379info not available 380 381 382=for bad 383 384mnpout does not process bad values. 385It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. 386 387 388=cut 389 390 391 392 393 394 395*mnpout = \&PDL::Minuit::mnpout; 396 397 398 399 400 sub mn_stat{ 401 if (defined (my $logfile = $mn_options->{Log})){ 402 PDL::Minuit::mn_abre($mn_options->{Unit},$logfile,'old'); 403 } 404 405 406 my ($fmin,$fedm,$errdef,$npari,$nparx,$istat) = PDL::Minuit::mnstat(); 407 408 if (defined (my $logfile = $mn_options->{Log})){ 409 PDL::Minuit::mn_cierra($mn_options->{Unit}); 410 } 411 412 return ($fmin,$fedm,$errdef,$npari,$nparx,$istat); 413 } 414 415 416 417 418=head2 mnstat 419 420=for sig 421 422 Signature: (double [o] a(); double [o] b(); double [o] c(); int [o] ia(); int [o] ib(); int [o] ic()) 423 424 425=for ref 426 427info not available 428 429 430=for bad 431 432mnstat does not process bad values. 433It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. 434 435 436=cut 437 438 439 440 441 442 443*mnstat = \&PDL::Minuit::mnstat; 444 445 446 447 448 sub mn_emat{ 449 450 if (defined (my $logfile = $mn_options->{Log})){ 451 PDL::Minuit::mn_abre($mn_options->{Unit},$logfile,'old'); 452 } 453 454 my ($fmin,$fedm,$errdef,$npari,$nparx,$istat) = PDL::Minuit::mnstat(); 455 my $n = $npari->sum; 456 my $mat = zeroes($n,$n); 457 458 PDL::Minuit::mnemat($mat); 459 460 if (defined (my $logfile = $mn_options->{Log})){ 461 PDL::Minuit::mn_cierra($mn_options->{Unit}); 462 } 463 464 return $mat; 465 466 } 467 468 469 470 471=head2 mnemat 472 473=for sig 474 475 Signature: (double [o] mat(n,n)) 476 477 478=for ref 479 480info not available 481 482 483=for bad 484 485mnemat does not process bad values. 486It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. 487 488 489=cut 490 491 492 493 494 495 496*mnemat = \&PDL::Minuit::mnemat; 497 498 499 500 501 sub mn_err{ 502 503 barf "Usage: mn_err(par_number)" unless ($#_ == 0); 504 my $par_num = shift; 505 506 my $n = $mn_options->{N}; 507 if (($par_num < 1) || ($par_num > $n)) { barf "Parameter numbers range from 1 to $n "; } 508 509 if (defined (my $logfile = $mn_options->{Log})){ 510 PDL::Minuit::mn_abre($mn_options->{Unit},$logfile,'old'); 511 } 512 513 my ($eplus,$eminus,$eparab,$globcc) = PDL::Minuit::mnerrs($par_num); 514 515 if (defined (my $logfile = $mn_options->{Log})){ 516 PDL::Minuit::mn_cierra($mn_options->{Unit}); 517 } 518 519 return ($eplus,$eminus,$eparab,$globcc); 520 } 521 522 523 524 525=head2 mnerrs 526 527=for sig 528 529 Signature: (int ia(); double [o] a(); double [o] b(); double [o] c(); double [o] d()) 530 531 532=for ref 533 534info not available 535 536 537=for bad 538 539mnerrs does not process bad values. 540It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. 541 542 543=cut 544 545 546 547 548 549 550*mnerrs = \&PDL::Minuit::mnerrs; 551 552 553 554 555 sub mn_contour{ 556 barf "Usage: mn_contour(par_number_1,par_number_2,npt)" unless ($#_ == 2); 557 my $par_num_1 = shift; 558 my $par_num_2 = shift; 559 my $npt = shift; 560 561 my $fun_ref = $mn_options->{Function}; 562 563 my $n = $mn_options->{N}; 564 if (($par_num_1 < 1) || ($par_num_1 > $n)) { barf "Parameter numbers range from 1 to $n "; } 565 if (($par_num_2 < 1) || ($par_num_2 > $n)) { barf "Parameter numbers range from 1 to $n "; } 566 if ($npt < 5) { barf "Have to specify at least 5 points in routine contour "; } 567 568 my $xpt = zeroes($npt); 569 my $ypt = zeroes($npt); 570 my $nfound = pdl->new; 571 572 PDL::Minuit::mncont($par_num_1,$par_num_2,$npt,$xpt,$ypt,$nfound,$fun_ref,$n); 573 574 if (defined (my $logfile = $mn_options->{Log})){ 575 PDL::Minuit::mn_cierra($mn_options->{Unit}); 576 } 577 578 return ($xpt,$ypt,$nfound); 579 } 580 581 582 583 584=head2 mncont 585 586=for sig 587 588 Signature: (int ia(); int ib(); int ic(); double [o] a(n); double [o] b(n); int [o] id(); SV* function; int numelem) 589 590 591=for ref 592 593info not available 594 595 596=for bad 597 598mncont does not process bad values. 599It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. 600 601 602=cut 603 604 605 606 607 608 609*mncont = \&PDL::Minuit::mncont; 610 611 612 613 614=head2 mn_init() 615 616=for ref 617 618The function mn_init() does the basic initialization of the fit. The first argument 619has to be a reference to the function to be minimized. The function 620to be minimized has to receive five arguments 621($npar,$grad,$fval,$xval,$iflag). The first is the number 622of parameters currently variable. The second is the gradient 623of the function (which is not necessarily used, see 624the Minuit documentation). The third is the current value of the 625function. The fourth is a piddle with the values of the parameters. 626The fifth is an integer flag, which indicates what 627the function is supposed to calculate. The function has to 628return the values ($fval,$grad), the function value and 629the function gradient. 630 631There are three optional arguments to mn_init(). By default, the output of Minuit 632will come through STDOUT unless a filename $logfile is given 633in the Log option. Note that this will mercilessly erase $logfile 634if it already exists. Additionally, a title can be given to the fit 635by the Title option, the default is 'Minuit Fit'. If the output is 636written to a logfile, this is assigned Fortran unit number 88. If for 637whatever reason you want to have control over the unit number 638that Fortran associates to the logfile, you can pass the number 639through the Unit option. 640 641=for usage 642 643Usage: 644 645 mn_init($function_ref,{Log=>$logfile,Title=>$title,Unit=>$unit}) 646 647=for example 648 649Example: 650 651 mn_init(\&my_function); 652 653 #same as above but outputting to a file 'log.out'. 654 #title for fit is 'My fit' 655 mn_init(\&my_function, 656 {Log => 'log.out', Title => 'My fit'}); 657 658 659 sub my_function{ 660 # the five variables input to the function to be minimized 661 # xval is a piddle containing the current values of the parameters 662 my ($npar,$grad,$fval,$xval,$iflag) = @_; 663 664 665 # Here is code computing the value of the function 666 # and potentially also its gradient 667 # ...... 668 669 # return the two variables. If no gradient is being computed 670 # just return the $grad that came as input 671 return ($fval, $grad); 672 } 673 674=head2 mn_def_pars() 675 676=for ref 677 678The function mn_def_pars() defines the initial values of the parameters of the function to 679be minimized and the value of the initial steps around these values that the 680minimizer will use for the first variations of the parameters in the search for the minimum. 681There are several optional arguments. One allows assigning names to these parameters which 682otherwise get names (Par_0, Par_1,....,Par_n) by default. Another two arguments can give 683lower and upper bounds for the parameters via two piddles. If the lower and upper bound for a 684given parameter are both equal to 0 then the parameter is unbound. By default these lower and 685upper bound piddles are set to zeroes(n), where n is the number of parameters, i.e. the 686parameters are unbound by default. 687 688The function needs two input variables: a piddle giving the initial values of the 689parameters and another piddle giving the initial steps. An optional reference to a 690perl array with the variable names can be passed, as well as piddles 691with upper and lower bounds for the parameters (see example below). 692 693It returns an integer variable which is 0 upon success. 694 695=for usage 696 697Usage: 698 699 $iflag = mn_def_pars($pars, $steps,{Names => \@names, 700 Lower_bounds => $lbounds, 701 Upper_bounds => $ubounds}) 702 703 704=for example 705 706Example: 707 708 #initial parameter values 709 my $pars = pdl(2.5,3.0); 710 711 #steps 712 my $steps = pdl(0.3,0.5); 713 714 #parameter names 715 my @names = ('intercept','slope'); 716 717 #use mn_def_pars with default parameter names (Par_0,Par_1,...) 718 my $iflag = mn_def_pars($pars,$steps); 719 720 #use of mn_def_pars explicitly specify parameter names 721 $iflag = mn_def_pars($pars,$steps,{Names => \@names}); 722 723 # specify lower and upper bounds for the parameters. 724 # The example below leaves parameter 1 (intercept) unconstrained 725 # and constrains parameter 2 (slope) to be between 0 and 100 726 my $lbounds = pdl(0, 0); 727 my $ubounds = pdl(0, 100); 728 729 $iflag = mn_def_pars($pars,$steps,{Names => \@names, 730 Lower_bounds => $lbounds, 731 Upper_bounds => $ubounds}}); 732 733 #same as above because $lbounds is by default zeroes(n) 734 $iflag = mn_def_pars($pars,$steps,{Names => \@names, 735 Upper_bounds => $ubounds}}); 736 737 738 739=head2 mn_excm() 740 741The function mn_excm() executes a Minuit command passed as 742a string. The first argument is the command string and an optional 743second argument is a piddle with arguments to the command. 744The available commands are listed in Chapter 4 of the Minuit 745manual (see url below). 746 747It returns an integer variable which is 0 upon success. 748 749=for usage 750 751Usage: 752 753 $iflag = mn_excm($command_string, {$arglis}) 754 755=for example 756 757Example: 758 759 #start a simplex minimization 760 my $iflag = mn_excm('simplex'); 761 762 #same as above but specify the maximum allowed numbers of 763 #function calls in the minimization 764 my $arglist = pdl(1000); 765 $iflag = mn_excm('simplex',$arglist); 766 767 #start a migrad minimization 768 $iflag = mn_excm('migrad') 769 770 #set Minuit strategy in order to get the most reliable results 771 $arglist = pdl(2) 772 $iflag = mn_excm('set strategy',$arglist); 773 774 # each command can be specified by a minimal string that uniquely 775 # identifies it (see Chapter 4 of Minuit manual). The comannd above 776 # is equivalent to: 777 $iflag = mn_excm('set stra',$arglis); 778 779=head2 mn_pout() 780 781The function mn_pout() gets the current value of a parameter. It 782takes as input the parameter number and returns an array with the 783parameter value, the current estimate of its uncertainty (0 if 784parameter is constant), lower bound on the parameter, if any 785(otherwise 0), upper bound on the parameter, if any (otherwise 0), 786integer flag (which is equal to the parameter number if variable, 787zero if the parameter is constant and negative if parameter is 788not defined) and the parameter name. 789 790=for usage 791 792Usage: 793 794 ($val,$err,$bnd1,$bnd2,$ivarbl,$par_name) = mn_pout($par_number); 795 796=head2 mn_stat() 797 798The function mn_stat() gets the current status of the minimization. 799It returns an array with the best function value found so far, 800the estimated vertical distance remaining to minimum, the value 801of UP defining parameter uncertainties (default is 1), the number 802of currently variable parameters, the highest parameter defined 803and an integer flag indicating how good the covariance matrix is 804(0=not calculated at all; 1=diagonal approximation, not accurate; 8052=full matrix, but forced positive definite; 3=full accurate matrix) 806 807=for usage 808 809Usage: 810 811 ($fmin,$fedm,$errdef,$npari,$nparx,$istat) = mn_stat(); 812 813=head2 mn_emat() 814 815The function mn_emat returns the covariance matrix as a piddle. 816 817=for usage 818 819Usage: 820 821 $emat = mn_emat(); 822 823=head2 mn_err() 824 825The function mn_err() returns the current existing values for 826the error in the fitted parameters. It returns an array 827with the positive error, the negative error, the "parabolic" 828parameter error from the error matrix and the global correlation 829coefficient, which is a number between 0 and 1 which gives 830the correlation between the requested parameter and that linear 831combination of all other parameters which is most strongly 832correlated with it. Unless the command 'MINOS' has been issued via 833the function mn_excm(), the first three values will be equal. 834 835=for usage 836 837Usage: 838 839 ($eplus,$eminus,$eparab,$globcc) = mn_err($par_number); 840 841=head2 mn_contour() 842 843The function mn_contour() finds contours of the function being minimized 844with respect to two chosen parameters. The contour level is given by 845F_min + UP, where F_min is the minimum of the function and UP is the ERRordef 846specified by the user, or 1.0 by default (see Minuit manual). The contour 847calculated by this function is dynamic, in the sense that it represents the 848minimum of the function being minimized with respect to all the other NPAR-2 parameters 849(if any). 850 851The function takes as input the parameter numbers with respect to which the contour 852is to be determined (two) and the number of points $npt required on the contour (>4). 853It returns an array with piddles $xpt,$ypt containing the coordinates of the contour 854and a variable $nfound indicating the number of points actually found in the contour. 855If all goes well $nfound will be equal to $npt, but it can be negative if the input 856arguments are not valid, zero if less than four points have been found or <$npt if the 857program could not find $npt points. 858 859=for usage 860 861Usage: 862 863 ($xpt,$ypt,$nfound) = mn_contour($par_number_1,$par_number_2,$npt) 864 865=head1 SEE ALSO 866 867L<PDL> 868 869The Minuit documentation is online at 870 871 http://wwwasdoc.web.cern.ch/wwwasdoc/minuit/minmain.html 872 873 874=head1 AUTHOR 875 876This file copyright (C) 2007 Andres Jordan <ajordan@eso.org>. 877All rights reserved. There is no warranty. You are allowed to redistribute this 878software/documentation under certain conditions. For details, see the file 879COPYING in the PDL distribution. If this file is separated from the 880PDL distribution, the copyright notice should be included in the file. 881 882=cut 883 884 885 886; 887 888 889 890# Exit with OK status 891 8921; 893 894