cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
parse_optimize.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 /*ParseOptimize parse the optimize and grid command lines */
4 /*GetOptColDen read observed column densities & errors for optimizer */
5 /*GetOptLineInt parse observed line intensities for optimization routines */
6 /*GetOptTemp read observed temperatures & errors for optimizer */
7 #include "cddefines.h"
8 #include "trace.h"
9 #include "optimize.h"
10 #include "grid.h"
11 #include "input.h"
12 #include "prt.h"
13 #include "parse.h"
14 #include "lines_service.h"
15 
16 static const realnum DEFERR = 0.05f;
17 
18 /*GetOptLineInt parse observed line intensities for optimization routines */
19 STATIC void GetOptLineInt(char *chCard );
20 
21 /*GetOptColDen read observed column densities & errors for optimizer */
22 STATIC void GetOptColDen(char *chCard );
23 
24 /*GetOptTemp read observed temperatures & errors for optimizer */
25 STATIC void GetOptTemp(char *chCard );
26 
27 /* ParseOptimize - called from ParseCommands if GRID or OPTIMIZE command found */
29  /* command line, which was changed to all caps in main parsing routine */
30  char *chCard)
31 {
32  bool lgEOL;
33  long int i;
34 
35  DEBUG_ENTRY( "ParseOptimize()" );
36 
37  /* this must come first so that contents of filename do not trigger wrong command */
38  if( nMatch("FILE",chCard) )
39  {
40  /* option to send final set of parameters to an input file
41  * get name within double quotes,
42  * and set to blanks in chCard and OrgCard */
43  /* file handle set NULL in optimize_do when first set up,
44  * this routine called many times during optimization process,
45  * only want to open file one time */
46  if( optimize.ioOptim == NULL )
47  {
48  /* chCard is all caps at this point. GetQuote will work with
49  * original version of command line, which preserves case of
50  * characters. Also removes string between quotes */
51  if( GetQuote( chOptimFileName , chCard , true ) )
52  /* this can't happen since GetQuote exits if quote not
53  * found, since true sets abort mode */
54  TotalInsanity();
55 
56  /* open the file */
58  }
59  }
60 
61  else if( nMatch("COLU",chCard) )
62  {
63  /* optimize column density */
64  optimize.lgOptCol = true;
65  optimize.lgOptimize = true;
66 
67  /* read column densities to match */
68  GetOptColDen(chCard);
69  }
70 
71  else if( nMatch("CONT",chCard) )
72  {
73  /* set flag saying that optimization should start from continue file */
74  optimize.lgOptCont = true;
75  optimize.lgOptimize = true;
76  }
77 
78  /* >>chng 06 sep 4, exclude lines with "GRID", handled below. */
79  else if( nMatch("INCR",chCard) && !nMatch("GRID",chCard) )
80  {
81  /* scan off increments for the previously selected parameter */
82  if( optimize.nparm > 0 )
83  {
84  /* also called during optimization process, ignore then */
85  i = 5;
87  (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
88  }
89  }
90 
91  else if( nMatch("LUMI",chCard) || nMatch("INTE",chCard) )
92  {
93  /* scan off intensity or luminosity of normalization line */
94  i = 5;
95  optimize.optint = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
96  optimize.optier = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
97  if( lgEOL )
99 
100  /* set flag to say that intensity or luminosity of line set */
101  optimize.lgOptLum = true;
102  optimize.lgOptimize = true;
103  }
104 
105  else if( nMatch("ITER",chCard) )
106  {
107  /* scan off number of iterations */
108  i = 5;
109  optimize.nIterOptim = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
110  }
111 
112  else if( nMatch("LINE",chCard) )
113  {
114  /* read lines to match */
115  GetOptLineInt(chCard);
116 
117  /* set flag saying that something has been set */
118  optimize.lgOptLin = true;
119  optimize.lgOptimize = true;
120  }
121 
122  else if( nMatch("PHYM",chCard) )
123  {
124  /* use PvH's PHYMIR to optimize parameters */
125  strcpy( optimize.chOptRtn, "PHYM" );
126 # ifdef __unix
127  optimize.lgParallel = ! nMatch("SEQU",chCard);
128 # else
129  optimize.lgParallel = false;
130 # endif
131  if( optimize.lgParallel )
132  {
133  i = 5;
134  long dum = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
135  /* default is the total number of available CPUs */
136  optimize.useCPU = lgEOL ? cpu.nCPU() : dum;
137  }
138  else
139  {
140  optimize.useCPU = 1;
141  }
142  }
143 
144  /* >>chng 06 aug 22, grid now done differently, requires two range parameters
145  * and the number of points, parsed below as "GRID". */
146  else if( nMatch("RANG",chCard) && !nMatch("GRID",chCard) )
147  {
148  /* scan off range for the previously selected variable */
149  if( optimize.nparm > 0 )
150  {
151  bool lgFirstOneReal = false;
152 
153  i = 5;
154  optimize.varang[optimize.nparm-1][0] =
155  (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
156  if( lgEOL )
157  {
158  optimize.varang[optimize.nparm-1][0] = -1e38f;
159  }
160  else
161  {
162  lgFirstOneReal = true;
163  }
164 
165  optimize.varang[optimize.nparm-1][1] =
166  (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
167  if( lgEOL )
168  {
169  optimize.varang[optimize.nparm-1][1] = 1e38f;
170  }
171  else if( lgFirstOneReal )
172  {
173  /* >>chng 06 aug 22, swap if second range parameter is less than the first,
174  * and always increment optimize.nRangeSet */
177  {
178  realnum temp = optimize.varang[optimize.nparm-1][0];
180  optimize.varang[optimize.nparm-1][1] = temp;
181  }
182  }
183  }
184  }
185 
186  else if( nMatch("SUBP",chCard) )
187  {
188  /* use subplex to optimize parameters */
189  strcpy( optimize.chOptRtn, "SUBP" );
190  }
191 
192  /* match a temperature */
193  else if( nMatch("TEMP",chCard) )
194  {
195  /* read temperatures to match */
196  GetOptTemp(chCard);
197 
198  /* set flag saying that optimize temps has been set */
199  optimize.lgOptTemp = true;
200  optimize.lgOptimize = true;
201  }
202 
203  else if( nMatch("TOLE",chCard) )
204  {
205  /* scan off tolerance of fit, sum of residuals must be smaller than this
206  * default is 0.10 */
207  i = 5;
209  }
210 
211  else if( nMatch("TRAC",chCard) )
212  {
213  if( nMatch("STAR",chCard) )
214  {
215  /* trace start iteration number
216  * turn on trace printout starting on nth call to cloudy */
217  i = 5;
218  optimize.nTrOpt = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
219  if( lgEOL )
220  {
221  fprintf( ioQQQ, " optimize trace start command:\n" );
222  fprintf( ioQQQ, " The iteration number must appear.\n" );
223  cdEXIT(EXIT_FAILURE);
224  }
225  optimize.lgTrOpt = true;
226  }
227  else if( nMatch("FLOW",chCard) )
228  {
229  /* trace flow
230  * follow logical flow within code */
231  optimize.lgOptimFlow = true;
232  }
233  else
234  {
235  fprintf( ioQQQ, " optimize trace flow command:\n" );
236  fprintf( ioQQQ, " One of the sub keys START or FLOW must appear.\n" );
237  cdEXIT(EXIT_FAILURE);
238  }
239  }
240 
241  else
242  {
243  fprintf( ioQQQ, " ==%s== is unrecognized keyword, consult HAZY.\n", chCard );
244  cdEXIT(EXIT_FAILURE);
245  }
246  return;
247 }
248 
249 /*GetOptColDen read observed column densities & errors for optimizer */
250 STATIC void GetOptColDen(char *chCard )
251 {
252  char chCap[INPUT_LINE_LENGTH];
253  bool lgEOF,
254  lgEOL;
255  long int i;
256 
257  DEBUG_ENTRY( "GetOptColDen()" );
258 
259  /* read observed column densities & errors */
260  optimize.ncobs = 0;
261 
262  /* get first line */
263  input_readarray(chCard,&lgEOF);
264  if( lgEOF )
265  {
266  fprintf( ioQQQ, " Hit EOF while reading column density list; use END to end list.\n" );
267  cdEXIT(EXIT_FAILURE);
268  }
269 
270  /* now copy this line to chCap, then convert to caps */
271  strcpy( chCap, chCard );
272  caps(chCap);
273 
274  while( !lgEOF )
275  {
276  if( optimize.ncobs > NCOLLM )
277  {
278  fprintf( ioQQQ, " Too many column densities have been entered; the limit is%4ld. Increase variable NCOLLM.\n",
279  NCOLLM );
280  cdEXIT(EXIT_FAILURE);
281  }
282 
283  /* order on line is element label (col 1-4), ionization stage, column density, err */
284  /* copy cap'd version of first 4 char of chCard to chColDen_label */
285  cap4((char*)optimize.chColDen_label[optimize.ncobs] , chCard );
286 
287  /* now get the ion stage, this should be 1 for atom, up to element
288  * number plus one */
289  i = 5;
290  optimize.ion_ColDen[optimize.ncobs] = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
291  if( lgEOL )
292  {
293  fprintf( ioQQQ, " %s\n", chCard );
294  fprintf( ioQQQ, " The ionization stage MUST appear on this line. Sorry.\n" );
295  cdEXIT(EXIT_FAILURE);
296  }
297 
298  /* the ion must be 1 or greater unless requesting a special,
299  * like a molecule or excited state population, in which
300  * case ion = 0
301  * can't check on upper limit yet since have not resolved element name */
302  if( optimize.ion_ColDen[optimize.ncobs] < 0 )
303  {
304  fprintf( ioQQQ, " %s\n", chCard );
305  fprintf( ioQQQ, " An ionization stage of%4ld does not make sense. Sorry.\n",
307  cdEXIT(EXIT_FAILURE);
308  }
309 
310  optimize.ColDen_Obs[optimize.ncobs] = (realnum)pow(10.,FFmtRead(chCard,&i,
311  INPUT_LINE_LENGTH,&lgEOL));
312  if( lgEOL )
313  {
314  fprintf( ioQQQ, " %80.80s\n", chCard );
315  fprintf( ioQQQ, " An observed column density MUST be entered. Sorry.\n" );
316  cdEXIT(EXIT_FAILURE);
317  }
318 
320  if( optimize.chColDen_error[optimize.ncobs] <= 0.0 )
321  {
322  /* this is the relative error allowed */
324  }
325 
326  /* check if number is a limit - if '<' appears on the line then it is */
327  if( strchr( chCard , '<' ) != NULL )
328  {
329  /* value is an upper limit only, use negative error to flag */
331  }
332 
333  input_readarray(chCard,&lgEOF);
334  strcpy( chCap, chCard );
335  caps(chCap);
336  if( lgEOF )
337  {
338  fprintf( ioQQQ, " Hit EOF while reading column density list; use END to end list.\n" );
339  cdEXIT(EXIT_FAILURE);
340  }
341 
342  if( strncmp( chCap , "END" , 3) == 0 )
343  {
344  lgEOF = true;
345  }
346 
347  /* now increment the number of columns we have entered */
348  optimize.ncobs += 1;
349  }
350 
351  if( trace.lgTrace && optimize.lgTrOpt )
352  {
353  fprintf( ioQQQ, "%4ld columns were entered, they were;\n",
354  optimize.ncobs );
355  for( i=0; i < optimize.ncobs; i++ )
356  {
357  fprintf( ioQQQ, " %4.4s ion=%5ld%10.2e%10.2e\n",
360  }
361  }
362  return;
363 }
364 
365 /*GetOptLineInt parse observed line intensities for optimization routines */
366 STATIC void GetOptLineInt(char *chCard )
367 {
368  char chCap[INPUT_LINE_LENGTH];
369 
370  bool lgEOF,
371  lgEOL;
372  long int i;
373 
374  DEBUG_ENTRY( "GetOptLineInt()" );
375 
376  /* read observed line fluxes & errors */
377  optimize.nlobs = 0;
378 
379  input_readarray(chCard,&lgEOF);
380  if( lgEOF )
381  {
382  fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" );
383  cdEXIT(EXIT_FAILURE);
384  }
385 
386  strcpy( chCap, chCard );
387  caps(chCap);
388 
389  while( !lgEOF )
390  {
391  if( optimize.nlobs >= NOBSLM )
392  {
393  fprintf( ioQQQ,
394  " Too many lines have been entered; the limit is %ld. Increase variable NOBSLM.\n",
395  NOBSLM );
396  cdEXIT(EXIT_FAILURE);
397  }
398 
399  /* order on line is label (col 1-4), wavelength, flux, error */
400  strncpy( optimize.chLineLabel[optimize.nlobs], chCard , 4 );
401  /* null terminate the label*/
403 
404  i = 5;
405  /* next get the wavelength */
407 
408  /* now convert wavelength to angstroms */
409  /* line was entered, look for possible micron or cm label */
410  if( input.chCARDCAPS[i-1] == 'M' )
411  {
412  /* microns */
414  }
415  else if( input.chCARDCAPS[i-1] == 'C' )
416  {
417  /* microns */
419  }
420 
421  /* get the error associated with 4 significant figures */
424 
425  /* next get the observed intensity */
427  if( lgEOL )
428  {
429  fprintf( ioQQQ, " %s\n", chCard );
430  fprintf( ioQQQ, " The wavelength and relative intensity MUST be entered on this line. Sorry.\n" );
431  fprintf( ioQQQ, " The command line is the following:\n %s\n", chCard );
432  cdEXIT(EXIT_FAILURE);
433  }
434 
435  if( optimize.xLineInt_Obs[optimize.nlobs] <= 0. )
436  {
437  fprintf( ioQQQ, " An observed intensity of %.2e is not allowed. Sorry.\n",
439  fprintf( ioQQQ, " The command line is the following:\n %s\n", chCard );
440  cdEXIT(EXIT_FAILURE);
441  }
442 
443  /* finally the optional error */
445  /* most often will use the default */
446  if( optimize.xLineInt_error[optimize.nlobs] <= 0.0 )
447  {
448  /* this is the relative error allowed */
450  }
451 
452  /* check if number is a limit - if '<' appears on the line then it is */
453  if( strchr( chCard , '<' ) != NULL )
454  {
455  /* value is an upper limit only, use negative error to flag */
457  }
458 
459  /* get next line */
460  input_readarray(chCard,&lgEOF);
461  if( lgEOF )
462  {
463  fprintf( ioQQQ, " Hit EOF while reading line list for optimize command; use END to end list.\n" );
464  cdEXIT(EXIT_FAILURE);
465  }
466 
467  strcpy( chCap, chCard );
468  caps(chCap);
469  if( strncmp( chCap ,"END" , 3 ) == 0 )
470  lgEOF = true;
471 
472  /* finally increment the number of observed lines */
473  ++optimize.nlobs;
474  }
475 
476  if( trace.lgTrace && trace.lgTrOptm )
477  {
478  fprintf( ioQQQ, "%4ld lines were entered, they were;\n",
479  optimize.nlobs );
480 
481  for( i=0; i < optimize.nlobs; i++ )
482  {
483  fprintf( ioQQQ, " %4.4s ", optimize.chLineLabel[i] );
485 
486  fprintf( ioQQQ, " %10.2e%10.2e\n",
489  }
490  }
491  return;
492 }
493 
494 /*GetOptTemp parse observed line intensities for optimization routines */
495 STATIC void GetOptTemp(char *chCard )
496 {
497  char chCap[INPUT_LINE_LENGTH];
498 
499  bool lgEOF,
500  lgEOL;
501  long int i;
502 
503  DEBUG_ENTRY( "GetOptTemp()" );
504 
505  /* read observed line fluxes & errors - first set total number of observe temps */
506  optimize.nTempObs = 0;
507 
508  input_readarray(chCard,&lgEOF);
509  if( lgEOF )
510  {
511  fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" );
512  cdEXIT(EXIT_FAILURE);
513  }
514 
515  /* make line caps so we can parse it */
516  strcpy( chCap, chCard );
517  caps(chCap);
518 
519  while( !lgEOF )
520  {
521  if( optimize.nTempObs >= NOBSLM )
522  {
523  fprintf( ioQQQ,
524  " Too many temperatures have been entered; the limit is %ld. Increase variable NOBSLM.\n",
525  NOBSLM );
526  cdEXIT(EXIT_FAILURE);
527  }
528 
529  /* order on line is label (col 1-4), ion, temperature, error */
530  strncpy( optimize.chTempLab[optimize.nTempObs], chCard , 4 );
531  /* null terminate the label*/
533 
534  i = 5;
535  /* next get the ion stage */
536  optimize.ionTemp[optimize.nTempObs] = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
537 
538  /* next get the observed temperature */
540  if( lgEOL )
541  {
542  fprintf( ioQQQ, " %s\n", chCard );
543  fprintf( ioQQQ, " The ion stage and temperature MUST be entered on this line. Sorry.\n" );
544  fprintf( ioQQQ, " The command line is the following:\n %s\n", chCard );
545  cdEXIT(EXIT_FAILURE);
546  }
547 
548  /* temperatures less than or equal to 10 are logs */
549  if( optimize.temp_obs[optimize.nTempObs] <= 10. )
550  {
552  }
553 
554  /* finally the optional error */
556  /* most often will use the default */
557  if( optimize.temp_error[optimize.nTempObs] <= 0.0 )
558  {
559  /* this is the relative error allowed */
561  }
562 
563  /* check if number is a limit - if '<' appears on the line then it is */
564  if( strchr( chCard , '<' ) != NULL )
565  {
566  /* value is an upper limit only, use negative error to flag */
568  }
569 
570  /* check for radius or volume for how to weight the mean temp
571  * this will be the default */
572  strcpy( optimize.chTempWeight[optimize.nTempObs] , "radius" );
573  /* >>chng 05 dec 29, from chCard to chCap, unlike much of code in this file,
574  * chCard has been read in by this routine and contains the original form of
575  * the command line. It was converted to caps and stored in chCap above.
576  * As it was written only VOLUME was matched, would not match volume.
577  * Bug caught and corrected by Bohdan Melekh */
578  /*if( nMatch( "VOLUME" , chCard ) )*/
579  if( nMatch( "VOLUME" , chCap ) )
580  {
581  strcpy( optimize.chTempWeight[optimize.nTempObs] , "volume" );
582  }
583 
584  /* get next line */
585  input_readarray(chCard,&lgEOF);
586  if( lgEOF )
587  {
588  fprintf( ioQQQ, " Hit EOF while reading line list for optimize command; use END to end list.\n" );
589  cdEXIT(EXIT_FAILURE);
590  }
591 
592  strcpy( chCap, chCard );
593  caps(chCap);
594  if( strncmp( chCap ,"END" , 3 ) == 0 )
595  lgEOF = true;
596 
597  /* finally increment the number of observed lines */
598  ++optimize.nTempObs;
599  }
600 
601  if( trace.lgTrace && trace.lgTrOptm )
602  {
603  fprintf( ioQQQ, "%4ld temperatures were entered, they were;\n",
604  optimize.nTempObs );
605 
606  for( i=0; i < optimize.nTempObs; i++ )
607  {
608  fprintf( ioQQQ, " %4.4s ", optimize.chTempLab[i] );
609  fprintf( ioQQQ, " %li " , optimize.ionTemp[i] );
610 
611  fprintf( ioQQQ, " %.2e %.2e\n",
612  optimize.temp_obs[i],
613  optimize.temp_error[i] );
614  }
615  }
616  return;
617 }

Generated for cloudy by doxygen 1.8.4