cloudy trunk
Loading...
Searching...
No Matches
parse_optimize.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*ParseOptimize parse the optimize 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 "input.h"
11#include "prt.h"
12#include "energy.h"
13#include "predcont.h"
14#include "parser.h"
15#include "lines_service.h"
16
17static const realnum DEFERR = 0.05f;
18
19/*GetOptLineInt parse observed line intensities for optimization routines */
21
22/*GetOptColDen read observed column densities & errors for optimizer */
23STATIC void GetOptColDen(Parser &p );
24
25/*GetOptTemp read observed temperatures & errors for optimizer */
26STATIC void GetOptTemp(Parser &p );
27
28/* ParseOptimize - called from ParseCommands if OPTIMIZE command found */
30 /* command line, which was changed to all caps in main parsing routine */
31 Parser &p)
32{
33 DEBUG_ENTRY( "ParseOptimize()" );
34
35 /* this must come first so that contents of filename do not trigger wrong command */
36 if( p.nMatch("FILE") )
37 {
38 /* option to send final set of parameters to an input file
39 * get name within double quotes,
40 * and set to blanks in chCard and OrgCard */
41 /* chCard is all caps at this point. GetQuote will work with
42 * original version of command line, which preserves case of
43 * characters. Also removes string between quotes */
44 /* GetQuote will abort if it fails, so we can ignore the return value */
45 (void)p.GetQuote( chOptimFileName , true );
46 }
47
48 else if( p.nMatch("COLU") )
49 {
50 /* optimize column density */
51 /* read column densities to match */
52 GetOptColDen(p);
53
54 optimize.lgOptimize = true;
55 }
56
57 else if( p.nMatch("CONT") && p.nMatch("FLUX") )
58 {
59 /* optimize continuum flux at arbitrary wavelength */
60 double energy = p.FFmtRead();
61 if( p.lgEOL() )
62 p.NoNumb("energy");
63 const char* unit = p.StandardEnergyUnit();
64 long ind = t_PredCont::Inst().add( energy, unit );
65 Energy E( energy, unit );
66
67 double flux = p.FFmtRead();
68 if( p.lgEOL() )
69 p.NoNumb("flux");
70 if( flux <= 0. || p.nMatch( " LOG" ) )
71 flux = pow(10.,flux);
72 Flux F( E, flux, p.StandardFluxUnit() );
73
74 double relerr = p.FFmtRead();
75 if( p.lgEOL() )
76 relerr = DEFERR;
77 /* value is an upper limit only, use negative error to flag */
78 if( p.nMatch("<" ) )
79 relerr = -relerr;
80
81 optimize.ContIndex.push_back( ind );
82 optimize.ContEner.push_back( E );
83 optimize.ContNFnu.push_back( F );
84 optimize.ContNFnuErr.push_back( chi2_type(relerr) );
85 optimize.lgOptimize = true;
86 }
87
88 else if( p.nMatch("CONT") )
89 {
90 /* set flag saying that optimization should start from continue file */
91 optimize.lgOptCont = true;
92 optimize.lgOptimize = true;
93 }
94
95 else if( p.nMatch("DIAM") )
96 {
97 /* optimize angular diameter */
98 optimize.optDiam = chi2_type( p.FFmtRead() );
99 optimize.optDiamErr = chi2_type( p.FFmtRead() );
100 if( p.lgEOL() )
101 optimize.optDiamErr = chi2_type( DEFERR );
102 /* value is an upper limit only, use negative error to flag */
103 if( p.nMatch( "<" ) )
104 optimize.optDiamErr = -optimize.optDiamErr;
105
106 if( p.nMatch( " LOG" ) )
107 optimize.optDiam = chi2_type( pow(chi2_type(10.),optimize.optDiam) );
108
109 optimize.lgOptDiam = true;
110 /* angular diameter is default in arcsec, or in cm of keyword present */
111 optimize.lgDiamInCM = ( p.nMatch( " CM " ) != 0 );
112 optimize.lgOptimize = true;
113 }
114
115 else if( p.nMatch("INCR") )
116 {
117 /* scan off increments for the previously selected parameter */
118 if( optimize.nparm > 0 )
119 {
120 /* also called during optimization process, ignore then */
121 optimize.OptIncrm[optimize.nparm-1] =
122 (realnum)p.FFmtRead();
123 }
124 }
125
126 else if( p.nMatch("LUMI") || p.nMatch("INTE") )
127 {
128 /* scan off intensity or luminosity of normalization line */
129 optimize.optint = (realnum)p.FFmtRead();
130 optimize.optier = (realnum)p.FFmtRead();
131 if( p.lgEOL() )
132 optimize.optier = DEFERR;
133
134 // default intrinsic intensity, accept emergent
135 optimize.nOptLum = p.nMatch("EMER") ? 1 : 0;
136
137 /* set flag to say that intensity or luminosity of line set */
138 optimize.lgOptLum = true;
139 optimize.lgOptimize = true;
140 }
141
142 else if( p.nMatch("ITER") )
143 {
144 /* scan off number of iterations */
145 optimize.nIterOptim = (long)p.FFmtRead();
146 }
147
148 else if( p.nMatch("LINE") )
149 {
150 /* read lines to match */
151 GetOptLineInt(p);
152
153 optimize.lgOptimize = true;
154 }
155
156 else if( p.nMatch("PHYM") )
157 {
158 /* use PvH's PHYMIR to optimize parameters */
159 strcpy( optimize.chOptRtn, "PHYM" );
160# if defined(__unix) || defined(__APPLE__)
161 // turn parallel mode off if explicitly requested or if we are in single rank mode
162 // the latter means that each rank runs its own model in single-CPU mode
163 optimize.lgParallel = ! ( p.nMatch("SEQU") || cpu.i().lgMPISingleRankMode() );
164# else
165 optimize.lgParallel = false;
166# endif
167 if( optimize.lgParallel )
168 {
169 long dum = (long)p.FFmtRead();
170 /* default is the total number of available CPUs */
171 optimize.useCPU = p.lgEOL() ? cpu.i().nCPU() : dum;
172 }
173 else
174 {
175 optimize.useCPU = 1;
176 }
177 }
178
179 else if( p.nMatch("RANG") )
180 {
181 /* scan off range for the previously selected variable */
182 if( optimize.nparm > 0 )
183 {
184 bool lgFirstOneReal = false;
185
186 optimize.varang[optimize.nparm-1][0] =
187 (realnum)p.FFmtRead();
188 if( p.lgEOL() )
189 {
190 optimize.varang[optimize.nparm-1][0] = -FLT_MAX;
191 }
192 else
193 {
194 lgFirstOneReal = true;
195 }
196
197 optimize.varang[optimize.nparm-1][1] =
198 (realnum)p.FFmtRead();
199 if( p.lgEOL() )
200 {
201 optimize.varang[optimize.nparm-1][1] = FLT_MAX;
202 }
203 else if( lgFirstOneReal )
204 {
205 /* >>chng 06 aug 22, swap if second range parameter is less than the first,
206 * and always increment optimize.nRangeSet */
207 ++optimize.nRangeSet;
208 if( optimize.varang[optimize.nparm-1][1] < optimize.varang[optimize.nparm-1][0] )
209 {
210 realnum temp = optimize.varang[optimize.nparm-1][0];
211 optimize.varang[optimize.nparm-1][0] = optimize.varang[optimize.nparm-1][1];
212 optimize.varang[optimize.nparm-1][1] = temp;
213 }
214 }
215 }
216 }
217
218 else if( p.nMatch("SUBP") )
219 {
220 /* use subplex to optimize parameters */
221 strcpy( optimize.chOptRtn, "SUBP" );
222 }
223
224 /* match a temperature */
225 else if( p.nMatch("TEMP") )
226 {
227 /* read temperatures to match */
228 GetOptTemp(p);
229
230 optimize.lgOptimize = true;
231 }
232
233 else if( p.nMatch("TOLE") )
234 {
235 /* scan off tolerance of fit, sum of residuals must be smaller than this
236 * default is 0.10 */
237 optimize.OptGlobalErr = (realnum)p.FFmtRead();
238 }
239
240 else if( p.nMatch("TRAC") )
241 {
242 if( p.nMatch("STAR") )
243 {
244 /* trace start iteration number
245 * turn on trace printout starting on nth call to cloudy */
246 optimize.nTrOpt = (long)p.FFmtRead();
247 if( p.lgEOL() )
248 {
249 fprintf( ioQQQ, " optimize trace start command:\n" );
250 fprintf( ioQQQ, " The iteration number must appear.\n" );
252 }
253 optimize.lgTrOpt = true;
254 }
255 else if( p.nMatch("FLOW") )
256 {
257 /* trace flow
258 * follow logical flow within code */
259 optimize.lgOptimFlow = true;
260 }
261 else
262 {
263 fprintf( ioQQQ, " optimize trace flow command:\n" );
264 fprintf( ioQQQ, " One of the sub keys START or FLOW must appear.\n" );
266 }
267 }
268
269 else
270 {
271 p.PrintLine(ioQQQ);
272 fprintf( ioQQQ, " is unrecognized keyword, consult HAZY.\n" );
274 }
275 return;
276}
277
278/*GetOptColDen read observed column densities & errors for optimizer */
280{
281
282 DEBUG_ENTRY( "GetOptColDen()" );
283
284 /* read observed column densities & errors */
285
286 /* get first line */
287 p.getline();
288 if( p.m_lgEOF )
289 {
290 fprintf( ioQQQ, " Hit EOF while reading column density list; use END to end list.\n" );
292 }
293
294 while( !p.m_lgEOF )
295 {
296 /* order on line is element label (col 1-4), ionization stage, column density, err */
297 /* copy cap'd version of first 4 char of chCard to chColDen_label */
298 optimize.chColDen_label.push_back(p.getCommand(4));
299
300 /* now get the ion stage, this should be 1 for atom, up to element
301 * number plus one */
302 long ion = nint(p.FFmtRead());
303 if( p.lgEOL() )
304 {
305 p.PrintLine( ioQQQ );
306 fprintf( ioQQQ, " The ionization stage MUST appear on this line. Sorry.\n" );
308 }
309
310 /* the ion must be 1 or greater unless requesting a special,
311 * like a molecule or excited state population, in which
312 * case ion = 0
313 * can't check on upper limit yet since have not resolved element name */
314 if( ion < 0 )
315 {
316 p.PrintLine( ioQQQ );
317 fprintf( ioQQQ, " An ionization stage of %ld does not make sense. Sorry.\n", ion );
319 }
320 optimize.ion_ColDen.push_back(ion);
321
322 optimize.ColDen_Obs.push_back( realnum(pow(10.,p.FFmtRead())) );
323 if( p.lgEOL() )
324 {
325 p.PrintLine( ioQQQ );
326 fprintf( ioQQQ, " An observed column density MUST be entered. Sorry.\n" );
328 }
329
330 realnum err = realnum(p.FFmtRead());
331 if( err <= 0.0 )
332 {
333 /* this is the relative error allowed */
334 err = realnum(DEFERR);
335 }
336
337 /* check if number is a limit - if '<' appears on the line then it is */
338 if( p.nMatch( "<" ) )
339 {
340 /* value is an upper limit only, use negative error to flag */
341 err = -err;
342 }
343 optimize.ColDen_error.push_back(err);
344
345 p.getline();
346 if( p.m_lgEOF )
347 {
348 fprintf( ioQQQ, " Hit EOF while reading column density list; use END to end list.\n" );
350 }
351
352 if( p.strcmp( "END" ) == 0 )
353 {
354 p.m_lgEOF = true;
355 }
356 }
357
358 if( trace.lgTrace && optimize.lgTrOpt )
359 {
360 fprintf( ioQQQ, "%ld columns were entered, they were;\n",
361 (long int) optimize.ColDen_Obs.size() );
362 for( long i=0; i < long(optimize.ColDen_Obs.size()); i++ )
363 {
364 fprintf( ioQQQ, " %4.4s ion=%5ld%10.2e%10.2e\n",
365 optimize.chColDen_label[i].c_str(), optimize.ion_ColDen[i],
366 optimize.ColDen_Obs[i], optimize.ColDen_error[i] );
367 }
368 }
369 return;
370}
371
372/*GetOptLineInt parse observed line intensities for optimization routines */
374{
375 DEBUG_ENTRY( "GetOptLineInt()" );
376
377 // default intrinsic intensity, accept emergent
378 optimize.nEmergent = p.nMatch("EMER") ? 1 : 0;
379
380 /* read observed line fluxes & errors */
381 p.getline();
382 if( p.m_lgEOF )
383 {
384 fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" );
386 }
387
388 while( !p.m_lgEOF )
389 {
390 /* order on line is label (col 1-4), wavelength, flux, error */
391 optimize.chLineLabel.push_back( p.getCommand(4) );
392
393 /* next get the wavelength */
394 realnum wavl = realnum(p.getWaveOpt());
395 optimize.wavelength.push_back(wavl);
396
397 /* get the error associated with 4 significant figures */
398 optimize.errorwave.push_back( WavlenErrorGet(wavl) );
399
400 /* next get the observed intensity */
401 realnum xLineInt = realnum(p.FFmtRead());
402 if( p.lgEOL() )
403 {
404 fprintf( ioQQQ, " The wavelength and relative intensity MUST be entered on this line. Sorry.\n" );
405 fprintf( ioQQQ, " The command line is the following:\n " );
406 p.PrintLine(ioQQQ);
408 }
409
410 if( xLineInt <= 0. )
411 {
412 fprintf( ioQQQ, " An observed intensity of %.2e is not allowed. Sorry.\n",
413 xLineInt );
414 fprintf( ioQQQ, " The command line is the following:\n" );
415 p.PrintLine(ioQQQ);
417 }
418 optimize.xLineInt_Obs.push_back(xLineInt);
419
420 /* finally the optional error */
421 realnum err = realnum(p.FFmtRead());
422 /* most often will use the default */
423 if( err <= 0.0 )
424 err = realnum(DEFERR);
425
426 /* check if number is a limit - if '<' appears on the line then it is */
427 if( p.nMatch( "<" ) )
428 {
429 /* value is an upper limit only, use negative error to flag */
430 err = -err;
431 }
432 optimize.xLineInt_error.push_back(err);
433
434 /* get next line */
435 p.getline();
436 if( p.m_lgEOF )
437 {
438 fprintf( ioQQQ, " Hit EOF while reading line list for optimize command; use END to end list.\n" );
440 }
441
442 if( p.strcmp( "END" ) == 0 )
443 p.m_lgEOF = true;
444 }
445
446 if( trace.lgTrace && trace.lgTrOptm )
447 {
448 fprintf( ioQQQ, "%ld lines were entered, they were:\n",
449 (long int) optimize.xLineInt_Obs.size() );
450
451 for( long i=0; i < long(optimize.xLineInt_Obs.size()); i++ )
452 {
453 fprintf( ioQQQ, " %4.4s ", optimize.chLineLabel[i].c_str() );
454 prt_wl( ioQQQ, optimize.wavelength[i] );
455
456 fprintf( ioQQQ, " %10.2e%10.2e\n",
457 optimize.xLineInt_Obs[i],
458 optimize.xLineInt_error[i] );
459 }
460 }
461 return;
462}
463
464/*GetOptTemp parse observed line intensities for optimization routines */
466{
467 DEBUG_ENTRY( "GetOptTemp()" );
468
469 /* read observed line fluxes & errors - first set total number of observe temps */
470 p.getline();
471 if( p.m_lgEOF )
472 {
473 fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" );
475 }
476
477 /* make line caps so we can parse it */
478 while( !p.m_lgEOF )
479 {
480 /* order on line is label (col 1-4), ion, temperature, error */
481 optimize.chTempLab.push_back( p.getCommand(4) );
482
483 /* next get the ion stage */
484 optimize.ionTemp.push_back( nint(p.FFmtRead()) );
485
486 /* next get the observed temperature */
487 realnum temp_obs = realnum(p.FFmtRead());
488 if( p.lgEOL() )
489 {
490 fprintf( ioQQQ, " The ion stage and temperature MUST be entered on this line. Sorry.\n" );
491 fprintf( ioQQQ, " The command line is the following:\n " );
492 p.PrintLine(ioQQQ);
494 }
495 /* temperatures less than or equal to 10 are logs */
496 if( temp_obs <= 10. )
497 temp_obs = realnum(pow( 10., (double)temp_obs ) );
498 optimize.temp_obs.push_back(temp_obs);
499
500 /* finally the optional error */
501 realnum temp_error = realnum(p.FFmtRead());
502 /* most often will use the default */
503 if( temp_error <= 0.f )
504 temp_error = realnum(DEFERR);
505 /* check if number is a limit - if '<' appears on the line then it is */
506 if( p.nMatch( "<" ) )
507 temp_error = -temp_error;
508 optimize.temp_error.push_back(temp_error);
509
510 /* check for radius or volume for how to weight the mean temp
511 * this will be the default */
512 /* >>chng 05 dec 29, from chCard to chCap, unlike much of code in this file,
513 * chCard has been read in by this routine and contains the original form of
514 * the command line. It was converted to caps and stored in chCap above.
515 * As it was written only VOLUME was matched, would not match volume.
516 * Bug caught and corrected by Bohdan Melekh */
517 if( p.nMatch( "VOLUME" ) )
518 optimize.chTempWeight.push_back("volume");
519 else
520 optimize.chTempWeight.push_back("radius");
521
522 /* get next line */
523 p.getline();
524 if( p.m_lgEOF )
525 {
526 fprintf( ioQQQ, " Hit EOF while reading line list for optimize command; use END to end list.\n" );
528 }
529
530 if( p.strcmp( "END" ) == 0 )
531 p.m_lgEOF = true;
532 }
533
534 if( trace.lgTrace && trace.lgTrOptm )
535 {
536 fprintf( ioQQQ, "%ld temperatures were entered, they were;\n",
537 (long int) optimize.temp_obs.size() );
538
539 for( long i=0; i < long(optimize.temp_obs.size()); i++ )
540 {
541 fprintf( ioQQQ, " %4.4s ", optimize.chTempLab[i].c_str() );
542 fprintf( ioQQQ, " %li " , optimize.ionTemp[i] );
543
544 fprintf( ioQQQ, " %.2e %.2e\n",
545 optimize.temp_obs[i],
546 optimize.temp_error[i] );
547 }
548 }
549 return;
550}
FILE * ioQQQ
Definition cddefines.cpp:7
#define STATIC
Definition cddefines.h:97
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
long nint(double x)
Definition cddefines.h:719
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
Definition energy.h:8
Definition flux.h:10
bool getline(void)
Definition parser.cpp:164
double FFmtRead(void)
Definition parser.cpp:353
bool nMatch(const char *chKey) const
Definition parser.h:135
double getWaveOpt()
Definition parser.cpp:244
const char * StandardEnergyUnit(void) const
Definition parser.cpp:174
int strcmp(const char *s2)
Definition parser.h:177
bool lgEOL(void) const
Definition parser.h:98
string StandardFluxUnit(void) const
Definition parser.cpp:178
NORETURN void NoNumb(const char *chDesc) const
Definition parser.cpp:233
bool m_lgEOF
Definition parser.h:42
int GetQuote(char *chLabel, bool lgABORT)
Definition parser.h:209
string getCommand(long i)
Definition parser.h:215
int PrintLine(FILE *fp) const
Definition parser.h:204
static t_PredCont & Inst()
Definition cddefines.h:175
long add(double energy, const char *unit="Ryd")
Definition predcont.cpp:122
static t_cpu cpu
Definition cpu.h:355
realnum WavlenErrorGet(realnum wavelength)
char chOptimFileName[INPUT_LINE_LENGTH]
Definition optimize.cpp:6
t_optimize optimize
Definition optimize.cpp:5
double chi2_type
Definition optimize.h:49
STATIC void GetOptLineInt(Parser &p)
static const realnum DEFERR
void ParseOptimize(Parser &p)
STATIC void GetOptTemp(Parser &p)
STATIC void GetOptColDen(Parser &p)
void prt_wl(FILE *ioOUT, realnum wl)
Definition prt.cpp:13
t_trace trace
Definition trace.cpp:5