cloudy  trunk
 All Data Structures 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-2008 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 "punch.h"
10 #include "optimize.h"
11 
12 /* called by cdDrive, this returns false if things went ok, true for disaster */
13 bool lgOptimize_do(void)
14 {
15  long int i,
16  iflag,
17  ii,
18  iworke[NPLXMX],
19  j,
20  mode,
21  need,
22  nfe,
23  np;
24  realnum fret,
25  fx,
26  param[LIMPAR],
27  ptem[LIMPAR],
28  delta[LIMPAR],
29  toler,
30  worke[NPLXMX];
31 /* double _e0[210]; */
32 /* realnum (*const p)[21] = (realnum(*)[21])_e0; */
33 /* realnum (*const xi)[20] = (realnum(*)[20])_e0; */
34 
35  DEBUG_ENTRY( "lgOptimize_do()" );
36 
37  /* main driver for optimization runs
38  * Drives cloudy to optimize variables;*/
39 
40  /* code originally written by R.F. Carswell, IOA Cambridge */
41 
42  toler = (realnum)log10(1. + optimize.OptGlobalErr);
43 
44  /*>>chng 07 feb 19, rm powell's method
45  * this is phymir */
46  if( strcmp(optimize.chOptRtn,"PHYM") == 0 )
47  {
48  /* Phymir method */
49  for( j=0; j < optimize.nvary; j++ )
50  {
51  ptem[j] = optimize.vparm[0][j];
52  delta[j] = optimize.vincr[j];
53  }
54  /* >>chng 06 jan 02, fix uninitialized var problem detected by valgrind/purify, PvH */
55  for( j=optimize.nvary; j < LIMPAR; j++ )
56  {
57  ptem[j] = -FLT_MAX;
58  delta[j] = -FLT_MAX;
59  }
60  optimize_phymir(ptem,delta,optimize.nvary,&fret,toler);
61  for( j=0; j < optimize.nvary; j++ )
62  {
63  optimize.vparm[0][j] = ptem[j];
64  }
65 
66  }
67 
68  else if( strcmp(optimize.chOptRtn,"SUBP") == 0 )
69  {
70  fprintf( ioQQQ, " Begin optimization with SUBPLEX\n" );
71  need = 2*optimize.nvary + optimize.nvary*(optimize.nvary + 4) + 1;
72  if( need > NPLXMX )
73  {
74  fprintf( ioQQQ, " Increase size of NPLXMX in parameter statements to handle this many variables.\n" );
75  fprintf( ioQQQ, " I need at least %5ld\n", need );
76  cdEXIT(EXIT_FAILURE);
77  }
78  for( j=0; j < optimize.nvary; j++ )
79  {
80  ptem[j] = optimize.vparm[0][j];
81  }
82 
83  /* The subroutine SUBPLX input into cloudy 8/4/94.
84  * The program itself is very well commented.
85  * The mode must set to 0 for the default values.
86  * The switch iflag tells if the program terminated normally. */
87  mode = 0;
88 
89  /* >>chng 97 dec 08, remove first arg, optimize_func, since not used in routines */
91  /* the number of parameters to vary */
93  /* the relative error, single number */
94  toler,
95  /* maximum number of function evaluations before giving up */
97  /* mode of operation, we simply set to zero */
98  mode,
99  /* the initial changes in the guessed best coefficients, typically 0.2 to 1 */
100  optimize.vincr,
101  /* a vector of nvary initial parameters that are the starting guesses for the parameters */
102  ptem,
103  /* a realnum, this is simply ignored */
104  &fx,
105  /* another parameter that is simply ignored, a long int */
106  &nfe,
107  /* a realnum that is NPLXMX long, used for working space by the routine */
108  /* an array that is 20*26 + 1 elements long, used for working space */
109  worke,
110  /* a long int that is NPLXMX long, used for working space by the routine */
111  /* an array that is 20*26 + 1 elements long, used for working space */
112  iworke,
113  /* a long int - says what happened, if -1 then exceeded nIterOptim iteration */
114  &iflag);
115 
116  if( iflag == -1 )
117  {
118  fprintf( ioQQQ, " SUBPLEX exceeding maximum iterations.\n This can be reset with the OPTIMZE ITERATIONS command.\n" );
119  }
120 
121  for( j=0; j < optimize.nvary; j++ )
122  {
123  optimize.vparm[0][j] = ptem[j];
124  }
125 
126  if( optimize.lgOptimFlow )
127  {
128  fprintf( ioQQQ, " trace return optimize_subplex:\n" );
129  for( j=0; j < optimize.nvary; j++ )
130  {
131  fprintf( ioQQQ, " Values:" );
132  for( ii=1; ii <= optimize.nvarxt[j]; ii++ )
133  {
134  fprintf( ioQQQ, "%10.2e", optimize.vparm[ii-1][j] );
135  }
136  fprintf( ioQQQ, "\n" );
137  }
138  }
139  }
140  else
141  TotalInsanity();
142  /*>>chng 07 feb 08m rm optimize_amoeba due to Robin & Peter license
143  * concerns - R837*/
144 
145  fprintf( ioQQQ, " **************************************************\n" );
146  fprintf( ioQQQ, " **************************************************\n" );
147  fprintf( ioQQQ, " **************************************************\n" );
148  fprintf( ioQQQ, "\n Cloudy was called %4ld times.\n\n", optimize.nOptimiz );
149 
150  for( i=0; i < optimize.nvary; i++ )
151  {
152  optimize.vparm[0][i] = (realnum)MIN2(optimize.vparm[0][i],optimize.varang[i][1]);
153  optimize.vparm[0][i] = (realnum)MAX2(optimize.vparm[0][i],optimize.varang[i][0]);
154  param[i] = optimize.vparm[0][i];
155  np = optimize.nvfpnt[i];
156 
157  /* now generate the actual command with parameter,
158  * there will be from 1 to 3 numbers on the line */
159  if( optimize.nvarxt[i] == 1 )
160  {
161  /* case with 1 parameter */
162  sprintf( input.chCardSav[np] , optimize.chVarFmt[i], optimize.vparm[0][i] );
163  }
164 
165  else if( optimize.nvarxt[i] == 2 )
166  {
167  /* case with 2 parameter */
168  sprintf( input.chCardSav[np] , optimize.chVarFmt[i], optimize.vparm[0][i], optimize.vparm[1][i]);
169  }
170 
171  else if( optimize.nvarxt[i] == 3 )
172  {
173  /* case with 3 parameter */
174  sprintf( input.chCardSav[np] , optimize.chVarFmt[i],
175  optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i] );
176  }
177 
178  else if( optimize.nvarxt[i] == 4 )
179  {
180  /* case with 4 parameter */
181  sprintf( input.chCardSav[np] , optimize.chVarFmt[i],
182  optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i], optimize.vparm[3][i] );
183  }
184 
185  else if( optimize.nvarxt[i] == 5 )
186  {
187  /* case with 5 parameter */
188  sprintf( input.chCardSav[np] , optimize.chVarFmt[i],
189  optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i] ,
190  optimize.vparm[3][i] , optimize.vparm[4][i]);
191  }
192 
193  else
194  {
195  fprintf(ioQQQ,"The number of variable options on this line makes no sense to me.\n");
196  cdEXIT(EXIT_FAILURE);
197  }
198 
199 
200  fprintf( ioQQQ, " Optimal command: %s\n", input.chCardSav[np]);
201  fprintf( ioQQQ, " Smallest value:%10.2e Largest value:%10.2e Allowed range %10.2e to %10.2e\n",
202  optimize.varmin[i], optimize.varmax[i], optimize.varang[i][0],
203  optimize.varang[i][1] );
204  }
205 
206  called.lgTalk = true;
207  called.lgTalkIsOK = true;
208  prt.lgFaintOn = true;
209 
210  /* though a page eject */
211  fprintf( ioQQQ, "\f" );
212 
213  /* punch optimal parameters on unit ioOptim */
214  if( optimize.ioOptim == NULL )
215  {
216  /* open default file name, optimal.in */
218  }
219 
220  for( i=0; i <= input.nSave; i++ )
221  {
222  fprintf( optimize.ioOptim, "%s\n", input.chCardSav[i]);
223  }
224  fclose( optimize.ioOptim );
225 
226  /* recalculate values in cloudy for the best fit, and print out
227  * all the information */
228  fret = (realnum)optimize_func(param);
229 
230 
231  if( lgAbort )
232  {
233  /* busted set means there were serious problems somewhere */
234  return true;
235  }
236  else
237  {
238  /* return 0 if everything is ok */
239  return false;
240  }
241 }

Generated for cloudy by doxygen 1.8.1.2