cloudy trunk
Loading...
Searching...
No Matches
parse_grid.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/*ParseGrid parse the grid command lines */
4#include "cddefines.h"
5#include "optimize.h"
6#include "grid.h"
7#include "parser.h"
8
9/* ParseGrid - called from ParseCommands if GRID command found */
11 /* command line, which was changed to all caps in main parsing routine */
12 Parser &p)
13{
14 DEBUG_ENTRY( "ParseGrid()" );
15
16 /* RP fake optimizer to run a grid of calculations, also accepts
17 * keyword XSPEC */
18 strcpy( optimize.chOptRtn, "XSPE" );
19 grid.lgGrid = true;
20
21 if( p.nMatch("REPE") )
22 {
23 /* just keep repeating, don't actually change the values in the grid.
24 * useful for debugging unintentional crosstalk */
25 grid.lgStrictRepeat = true;
26 }
27
28 /* 06 aug 22, change to accept three parameters: lower and upper limit and number of points. */
29 /* scan off range for the previously selected variable */
30 if( optimize.nparm > 0 )
31 {
32 ASSERT( optimize.nparm <= LIMPAR );
33
34 grid.paramLimits[optimize.nparm-1][0] = (realnum)p.FFmtRead();
35 grid.paramLimits[optimize.nparm-1][1] = (realnum)p.FFmtRead();
36 grid.paramIncrements[optimize.nparm-1] = (realnum)p.FFmtRead();
37 grid.lgLinearIncrements[optimize.nparm-1] = p.nMatch("LINE") ? true : false ;
38 if( grid.paramIncrements[optimize.nparm-1] < realnum(0.) )
39 grid.lgNegativeIncrements = true;
40
41 /* the increase step should not be 0 */
42 if( grid.paramIncrements[optimize.nparm-1] == 0. )
43 {
44 fprintf( ioQQQ," The increment (third parameter) should not be zero.\n" );
45 fprintf( ioQQQ," Sorry.\n" );
47 }
48
49 if( p.lgEOL() )
50 {
51 fprintf( ioQQQ," This command has changed since the definition given in Porter et al. 2006, PASP, 118, 920.\n" );
52 fprintf( ioQQQ," The grid command now requires three parameters: lower limit, upper limit, and increment.\n" );
53 fprintf( ioQQQ," The keywords RANGE and STEPS are no longer necessary.\n" );
54 fprintf( ioQQQ," Sorry.\n" );
56 }
57 else
58 {
59 ++optimize.nRangeSet;
60 }
61
62 if( p.nMatch("CYCL") )
63 {
64 /* cycle through the grid multiple times, only used for testing */
65 grid.nCycle = nint(p.FFmtRead());
66 if( p.lgEOL() )
67 grid.nCycle = 2;
68 if( grid.nCycle < 2 )
69 {
70 fprintf( ioQQQ, " Invalid repetion number for cycle: %ld\n", grid.nCycle );
71 fprintf( ioQQQ, " Usage: grid <p1> <p2> <p3> cycle [ <n> ] with n >= 2.\n" );
73 }
74 }
75
76 realnum ratio = (grid.paramLimits[optimize.nparm-1][1] - grid.paramLimits[optimize.nparm-1][0])/
77 grid.paramIncrements[optimize.nparm-1];
78
79 /* Alert if the uplimit and lowlimit are wrong */
80 if( ratio < realnum(0.) )
81 {
82 fprintf( ioQQQ, "The increment (third parameter) has the wrong sign. It doesn't take you from the initial to the final grid value (first and second parameter, resp.).\n" );
83 fprintf( ioQQQ," Sorry.\n" );
85 }
86
87 // this takes care of the blowup in the error due to cancellation in limits[1]-limits[0]
88 // it assumes that limits[1]-limits[0] is accurate within 3*eps*(limits[0]+limits[1])/2
89 // which should be a very conservative estimate...
90 realnum feps = realnum(1.5)*
91 (grid.paramLimits[optimize.nparm-1][1] + grid.paramLimits[optimize.nparm-1][0])/
92 (grid.paramLimits[optimize.nparm-1][1] - grid.paramLimits[optimize.nparm-1][0]);
93 long eps = max(nint(abs(feps)),3);
94
95 // this will blow for pathologically narrow grid ranges
96 ASSERT( eps <= INT_MAX );
97
98 // take special care if step is integer fraction of max-min (which is nearly always the case)
99 if( fp_equal( ratio, realnum(nint(ratio)), int(eps) ) )
100 grid.numParamValues[optimize.nparm-1] = nint(ratio) + 1;
101 else
102 grid.numParamValues[optimize.nparm-1] = long(ratio) + 1;
103
104 if( grid.numParamValues[optimize.nparm-1] < 2 )
105 {
106 fprintf( ioQQQ, " There must be at least two grid points in each dimension.\n" );
107 fprintf( ioQQQ, " Sorry.\n" );
109 }
110
111 grid.numParamValues[optimize.nparm-1] = MAX2( 2, grid.numParamValues[optimize.nparm-1] );
112
113 // Create some buffer area in the allowed range of parameter values to prevent
114 // accidentally going over the limit due to roundoff error. The buffer is 1/10th
115 // of a step, so should still guard against doing too many steps due to bugs.
116 realnum safety = 0.001f*grid.paramIncrements[optimize.nparm-1];
117
118 if( grid.lgLinearIncrements[optimize.nparm-1] )
119 {
120 if( grid.paramLimits[optimize.nparm-1][0]-safety<=0. )
121 {
122 fprintf(ioQQQ,"The current implementation of the grid command works with log parameter values even when you specify LINEAR.\n");
123 fprintf(ioQQQ,"A non-positive value was entered. The grid command cannot deal with this.\n");
125 }
126 optimize.varang[optimize.nparm-1][0] = log10(grid.paramLimits[optimize.nparm-1][0]-safety);
127 optimize.varang[optimize.nparm-1][1] = log10(grid.paramLimits[optimize.nparm-1][1]+safety);
128 }
129 else
130 {
131 optimize.varang[optimize.nparm-1][0] = grid.paramLimits[optimize.nparm-1][0]-safety;
132 optimize.varang[optimize.nparm-1][1] = grid.paramLimits[optimize.nparm-1][1]+safety;
133 }
134 }
135
136 return;
137}
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
long nint(double x)
Definition cddefines.h:719
long max(int a, long b)
Definition cddefines.h:775
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
double FFmtRead(void)
Definition parser.cpp:353
bool nMatch(const char *chKey) const
Definition parser.h:135
bool lgEOL(void) const
Definition parser.h:98
t_grid grid
Definition grid.cpp:5
t_optimize optimize
Definition optimize.cpp:5
const long LIMPAR
Definition optimize.h:61
void ParseGrid(Parser &p)