cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
optimize_do.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*lgOptimize_do main driver for optimization runs*/
4 #include "cddefines.h"
5 #define NPLXMX (LIMPAR*(LIMPAR+6)+1)
6 #include "input.h"
7 #include "called.h"
8 #include "prt.h"
9 #include "optimize.h"
10 
11 /* called by cdDrive, this returns false if things went ok, true for disaster */
12 bool lgOptimize_do(void)
13 {
14  long int i,
15  iflag,
16  ii,
17  iworke[NPLXMX],
18  j,
19  mode,
20  need,
21  nfe,
22  np;
23  chi2_type fret;
24  realnum fx,
25  ptem[LIMPAR],
26  delta[LIMPAR],
27  toler,
28  worke[NPLXMX];
29 
30  DEBUG_ENTRY( "lgOptimize_do()" );
31 
32  /* main driver for optimization runs
33  * Drives cloudy to optimize variables;*/
34 
35  /* code originally written by R.F. Carswell, IOA Cambridge */
36 
37  toler = (realnum)log10(1. + optimize.OptGlobalErr);
38 
39  if( strcmp(optimize.chOptRtn,"PHYM") == 0 )
40  {
41  /* Phymir method */
42  for( j=0; j < optimize.nvary; j++ )
43  {
44  ptem[j] = optimize.vparm[0][j];
45  delta[j] = optimize.vincr[j];
46  }
47  /* >>chng 06 jan 02, fix uninitialized var problem detected by valgrind/purify, PvH */
48  for( j=optimize.nvary; j < LIMPAR; j++ )
49  {
50  ptem[j] = -FLT_MAX;
51  delta[j] = -FLT_MAX;
52  }
53  optimize_phymir(ptem,delta,optimize.nvary,&fret,toler);
54  for( j=0; j < optimize.nvary; j++ )
55  optimize.vparm[0][j] = ptem[j];
56  }
57 
58  else if( strcmp(optimize.chOptRtn,"SUBP") == 0 )
59  {
60  fprintf( ioQQQ, " Begin optimization with SUBPLEX\n" );
61  need = 2*optimize.nvary + optimize.nvary*(optimize.nvary + 4) + 1;
62  if( need > NPLXMX )
63  {
64  fprintf( ioQQQ, " Increase size of NPLXMX in parameter statements to handle this many variables.\n" );
65  fprintf( ioQQQ, " I need at least %5ld\n", need );
67  }
68  for( j=0; j < optimize.nvary; j++ )
69  ptem[j] = optimize.vparm[0][j];
70 
71  /* The subroutine SUBPLX input into cloudy 8/4/94.
72  * The program itself is very well commented.
73  * The mode must set to 0 for the default values.
74  * The switch iflag tells if the program terminated normally. */
75  mode = 0;
76 
77  /* >>chng 97 dec 08, remove first arg, optimize_func, since not used in routines */
79  /* the number of parameters to vary */
81  /* the relative error, single number */
82  toler,
83  /* maximum number of function evaluations before giving up */
85  /* mode of operation, we simply set to zero */
86  mode,
87  /* the initial changes in the guessed best coefficients, typically 0.2 to 1 */
89  /* a vector of nvary initial parameters that are the starting guesses for the parameters */
90  ptem,
91  /* a realnum, this is simply ignored */
92  &fx,
93  /* another parameter that is simply ignored, a long int */
94  &nfe,
95  /* a realnum that is NPLXMX long, used for working space by the routine */
96  /* an array that is 20*26 + 1 elements long, used for working space */
97  worke,
98  /* a long int that is NPLXMX long, used for working space by the routine */
99  /* an array that is 20*26 + 1 elements long, used for working space */
100  iworke,
101  /* a long int - says what happened, if -1 then exceeded nIterOptim iteration */
102  &iflag);
103 
104  if( iflag == -1 )
105  {
106  fprintf( ioQQQ, " SUBPLEX exceeding maximum iterations.\n"
107  " This can be reset with the OPTIMZE ITERATIONS command.\n" );
108  }
109 
110  for( j=0; j < optimize.nvary; j++ )
111  optimize.vparm[0][j] = ptem[j];
112 
113  if( optimize.lgOptimFlow )
114  {
115  fprintf( ioQQQ, " trace return optimize_subplex:\n" );
116  for( j=0; j < optimize.nvary; j++ )
117  {
118  fprintf( ioQQQ, " Values:" );
119  for( ii=1; ii <= optimize.nvarxt[j]; ii++ )
120  fprintf( ioQQQ, " %.2e", optimize.vparm[ii-1][j] );
121  fprintf( ioQQQ, "\n" );
122  }
123  }
124  }
125  else
126  TotalInsanity();
127 
128  // turn output back on for final model
129  called.lgTalk = cpu.i().lgMPI_talk();
131  prt.lgFaintOn = true;
132 
133  if( called.lgTalk )
134  {
135  fprintf( ioQQQ, " **************************************************\n" );
136  fprintf( ioQQQ, " **************************************************\n" );
137  fprintf( ioQQQ, " **************************************************\n" );
138  fprintf( ioQQQ, "\n Cloudy was called %4ld times.\n\n", optimize.nOptimiz );
139 
140  for( i=0; i < optimize.nvary; i++ )
141  {
142  np = optimize.nvfpnt[i];
143 
144  /* now generate the actual command with parameter */
145  if( optimize.nvarxt[i] == 1 )
146  {
147  sprintf( input.chCardSav[np], optimize.chVarFmt[i], optimize.vparm[0][i] );
148  }
149 
150  else if( optimize.nvarxt[i] == 2 )
151  {
152  sprintf( input.chCardSav[np], optimize.chVarFmt[i], optimize.vparm[0][i],
153  optimize.vparm[1][i]);
154  }
155 
156  else if( optimize.nvarxt[i] == 3 )
157  {
158  sprintf( input.chCardSav[np], optimize.chVarFmt[i],
159  optimize.vparm[0][i], optimize.vparm[1][i], optimize.vparm[2][i] );
160  }
161 
162  else if( optimize.nvarxt[i] == 4 )
163  {
164  sprintf( input.chCardSav[np], optimize.chVarFmt[i],
165  optimize.vparm[0][i], optimize.vparm[1][i], optimize.vparm[2][i],
166  optimize.vparm[3][i] );
167  }
168 
169  else if( optimize.nvarxt[i] == 5 )
170  {
171  sprintf( input.chCardSav[np], optimize.chVarFmt[i],
172  optimize.vparm[0][i], optimize.vparm[1][i], optimize.vparm[2][i],
173  optimize.vparm[3][i], optimize.vparm[4][i]);
174  }
175 
176  else
177  {
178  fprintf(ioQQQ,"The number of variable options on this line makes no sense to me.\n");
180  }
181 
182  strncpy( input.chCardStrip[np], input.chCardSav[np], INPUT_LINE_LENGTH );
183 
184  fprintf( ioQQQ, " Optimal command: %s\n", input.chCardSav[np]);
185  fprintf( ioQQQ, " Smallest value:%10.2e Largest value:%10.2e Allowed range %10.2e to %10.2e\n",
186  optimize.varmin[i], optimize.varmax[i], optimize.varang[i][0],
187  optimize.varang[i][1] );
188  }
189  }
190 
191  if( cpu.i().lgMaster() )
192  {
193  /* save optimal parameters */
194  FILE *ioOptim = open_data( chOptimFileName, "w" );
195  for( i=0; i <= input.nSave; i++ )
196  if( input.InclLevel[i] == 0 )
197  fprintf( ioOptim, "%s\n", input.chCardSav[i]);
198  fclose( ioOptim );
199 
200  /* recalculate values in cloudy for the best fit, and print out
201  * all the information */
202  fprintf( ioQQQ, "\f" );
203 
204  for( i=0; i < optimize.nvary; i++ )
205  ptem[i] = optimize.vparm[0][i];
206 
207  (void)optimize_func(ptem);
208  }
209 
210  return lgAbort;
211 }
long int nSave
Definition: input.h:102
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:765
realnum varmax[LIMPAR]
Definition: optimize.h:187
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
t_input input
Definition: input.cpp:12
long int nvfpnt[LIMPAR]
Definition: optimize.h:198
t_cpu_i & i()
Definition: cpu.h:419
chi2_type optimize_func(const realnum param[], int grid_index=-1)
long int nOptimiz
Definition: optimize.h:250
realnum varang[LIMPAR][2]
Definition: optimize.h:201
void optimize_phymir(realnum[], const realnum[], long, chi2_type *, realnum)
char chVarFmt[LIMPAR][FILENAME_PATH_LENGTH_2]
Definition: optimize.h:267
FILE * ioQQQ
Definition: cddefines.cpp:7
realnum vparm[LIMEXT][LIMPAR]
Definition: optimize.h:192
bool lgTalk
Definition: called.h:12
char chCardStrip[NKRD][INPUT_LINE_LENGTH]
Definition: input.h:77
void optimize_subplex(long int n, double tol, long int maxnfe, long int mode, realnum scale[], realnum x[], realnum *fx, long int *nfe, realnum work[], long int iwork[], long int *iflag)
long int nIterOptim
Definition: optimize.h:209
#define NPLXMX
Definition: optimize_do.cpp:5
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
const int INPUT_LINE_LENGTH
Definition: cddefines.h:301
bool lgMaster() const
Definition: cpu.h:396
bool lgFaintOn
Definition: prt.h:245
#define cdEXIT(FAIL)
Definition: cddefines.h:482
bool lgMPI_talk() const
Definition: cpu.h:397
const long LIMPAR
Definition: optimize.h:61
t_optimize optimize
Definition: optimize.cpp:6
realnum vincr[LIMPAR]
Definition: optimize.h:195
t_prt prt
Definition: prt.cpp:14
char chOptRtn[5]
Definition: optimize.h:268
double chi2_type
Definition: optimize.h:49
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
realnum varmin[LIMPAR]
Definition: optimize.h:188
bool lgOptimize_do(void)
Definition: optimize_do.cpp:12
char chCardSav[NKRD][INPUT_LINE_LENGTH]
Definition: input.h:74
realnum OptGlobalErr
Definition: optimize.h:247
static t_cpu cpu
Definition: cpu.h:427
bool lgOptimFlow
Definition: optimize.h:252
char chOptimFileName[INPUT_LINE_LENGTH]
Definition: optimize.cpp:7
long int nvarxt[LIMPAR]
Definition: optimize.h:198
t_called called
Definition: called.cpp:4
long int nvary
Definition: optimize.h:204
bool lgTalkIsOK
Definition: called.h:23
bool lgAbort
Definition: cddefines.cpp:10
int InclLevel[NKRD]
Definition: input.h:91