Log In | Get Help   
Home My Page Projects Code Snippets Project Openings Mareframe
Summary Activity Forums Tracker Lists Tasks Docs Surveys News SCM Files
[mareframe] Annotation of /trunk/gadget/sionstep.cc
[mareframe] / trunk / gadget / sionstep.cc Repository:
ViewVC logotype

Annotation of /trunk/gadget/sionstep.cc

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (view) (download)

1 : agomez 1 #include "sionstep.h"
2 :     #include "areatime.h"
3 :     #include "errorhandler.h"
4 :     #include "readfunc.h"
5 :     #include "readword.h"
6 :     #include "gadget.h"
7 :     #include "global.h"
8 :    
9 :     SIOnStep::~SIOnStep() {
10 :     int i;
11 :     if (LgrpDiv != 0)
12 :     delete LgrpDiv;
13 :     if (LR != 0)
14 :     delete LR;
15 :     for (i = 0; i < areaindex.Size(); i++)
16 :     delete[] areaindex[i];
17 :     for (i = 0; i < colindex.Size(); i++)
18 :     delete[] colindex[i];
19 :     for (i = 0; i < obsIndex.Size(); i++) {
20 :     delete obsIndex[i];
21 :     delete modelIndex[i];
22 :     }
23 :     for (i = 0; i < weightIndex.Size(); i++)
24 :     delete weightIndex[i];
25 :     }
26 :    
27 :     SIOnStep::SIOnStep(CommentStream& infile, const char* datafilename, const CharPtrVector& aindex,
28 :     const TimeClass* const TimeInfo, const IntMatrix& areas, const CharPtrVector& charindex,
29 :     const char* givenname, int bio, SIType type) : HasName(givenname), Areas(areas), alptr(0) {
30 :    
31 :     biomass = bio;
32 :     sitype = type;
33 :     useweight = 0;
34 :     timeindex = 0;
35 :     slope = 0.0;
36 :     intercept = 0.0;
37 :    
38 :     int i;
39 :     for (i = 0; i < aindex.Size(); i++) {
40 :     areaindex.resize(new char[strlen(aindex[i]) + 1]);
41 :     strcpy(areaindex[i], aindex[i]);
42 :     }
43 :     for (i = 0; i < charindex.Size(); i++) {
44 :     colindex.resize(new char[strlen(charindex[i]) + 1]);
45 :     strcpy(colindex[i], charindex[i]);
46 :     }
47 :    
48 :     //read the fittype
49 :     readSIRegressionData(infile);
50 :    
51 :     //read the survey indices data from the datafile
52 :     ifstream datafile;
53 :     CommentStream subdata(datafile);
54 :     datafile.open(datafilename, ios::in);
55 :     handle.checkIfFailure(datafile, datafilename);
56 :     handle.Open(datafilename);
57 :     readSIData(subdata, TimeInfo);
58 :     handle.Close();
59 :     datafile.close();
60 :     datafile.clear();
61 :    
62 :     //resize to store the regression information
63 :     slopes.AddRows(areaindex.Size(), colindex.Size(), slope);
64 :     intercepts.AddRows(areaindex.Size(), colindex.Size(), intercept);
65 :     sse.AddRows(areaindex.Size(), colindex.Size(), 0.0);
66 :     if (useweight)
67 :     tmpWeight.resize(Years.Size(), 0.0);
68 :     tmpModel.resize(Years.Size(), 0.0);
69 :     tmpData.resize(Years.Size(), 0.0);
70 :     likelihoodValues.resize(areaindex.Size(), 0.0);
71 :     }
72 :    
73 :     void SIOnStep::readSIData(CommentStream& infile, const TimeClass* const TimeInfo) {
74 :    
75 :     double tmpnumber, tmpweight;
76 :     char tmparea[MaxStrLength], tmplabel[MaxStrLength];
77 :     strncpy(tmparea, "", MaxStrLength);
78 :     strncpy(tmplabel, "", MaxStrLength);
79 :     int keepdata, timeid, colid, areaid;
80 :     int i, year, step, count, reject;
81 :    
82 :     //Check the number of columns in the inputfile
83 :     infile >> ws;
84 :     if ((!useweight) && (countColumns(infile) != 5))
85 :     handle.logFileMessage(LOGFAIL, "wrong number of columns in inputfile - should be 5");
86 :     else if ((useweight) && (countColumns(infile) != 6))
87 :     handle.logFileMessage(LOGFAIL, "wrong number of columns in inputfile - should be 6");
88 :    
89 :     year = step = count = reject = 0;
90 :     while (!infile.eof()) {
91 :     keepdata = 1;
92 :    
93 :     if (useweight)
94 :     infile >> year >> step >> tmparea >> tmplabel >> tmpnumber >> tmpweight >> ws;
95 :     else
96 :     infile >> year >> step >> tmparea >> tmplabel >> tmpnumber >> ws;
97 :    
98 :     //crude check to see if something has gone wrong and avoid infinite loops
99 :     if (strlen(tmparea) == 0)
100 :     handle.logFileMessage(LOGFAIL, "failed to read data from file");
101 :    
102 :     //if tmparea is in areaindex keep data, else dont keep the data
103 :     areaid = -1;
104 :     for (i = 0; i < areaindex.Size(); i++)
105 :     if (strcasecmp(areaindex[i], tmparea) == 0)
106 :     areaid = i;
107 :    
108 :     if (areaid == -1)
109 :     keepdata = 0;
110 :    
111 :     //if tmplabel is in colindex find colid, else dont keep the data
112 :     colid = -1;
113 :     for (i = 0; i < colindex.Size(); i++)
114 :     if (strcasecmp(colindex[i], tmplabel) == 0)
115 :     colid = i;
116 :    
117 :     if (colid == -1)
118 :     keepdata = 0;
119 :    
120 :     //check if the year and step are in the simulation
121 :     timeid = -1;
122 :     if ((TimeInfo->isWithinPeriod(year, step)) && (keepdata == 1)) {
123 :     for (i = 0; i < Years.Size(); i++)
124 :     if ((Years[i] == year) && (Steps[i] == step))
125 :     timeid = i;
126 :    
127 :     if (timeid == -1) {
128 :     Years.resize(1, year);
129 :     Steps.resize(1, step);
130 :     timeid = (Years.Size() - 1);
131 :     obsIndex.resize(new DoubleMatrix(areaindex.Size(), colindex.Size(), 0.0));
132 :     modelIndex.resize(new DoubleMatrix(areaindex.Size(), colindex.Size(), 0.0));
133 :     if (useweight)
134 :     weightIndex.resize(new DoubleMatrix(areaindex.Size(), colindex.Size(), 0.0));
135 :     }
136 :    
137 :     } else
138 :     keepdata = 0;
139 :    
140 :     if (keepdata == 1) {
141 :     //survey indices data is required, so store it
142 :     count++;
143 :     (*obsIndex[timeid])[areaid][colid] = tmpnumber;
144 :     if (useweight)
145 :     (*weightIndex[timeid])[areaid][colid] = tmpweight;
146 :     } else
147 :     reject++; //count number of rejected data points read from file
148 :     }
149 :    
150 :     AAT.addActions(Years, Steps, TimeInfo);
151 :     if (count == 0)
152 :     handle.logMessage(LOGWARN, "Warning in surveyindex - found no data in the data file for", this->getName());
153 :     if (Years.Size() < 2)
154 :     handle.logMessage(LOGWARN, "Warning in surveyindex - insufficient data to fit a regression line for", this->getName());
155 :    
156 :     if ((handle.getLogLevel() >= LOGWARN) && (Steps.Size() > 0) && (this->getSIType() != SIEFFORT)) {
157 :     //JMB to be comparable, this should only take place on the same step in each year
158 :     //However the effort data will be on most timesteps so this check isnt required
159 :     step = Steps[0];
160 :     timeid = 0;
161 :     for (i = 1; i < Steps.Size(); i++)
162 :     if (Steps[i] != step)
163 :     timeid++;
164 :    
165 :     if (timeid != 0)
166 :     handle.logMessage(LOGWARN, "Warning in surveyindex - differing timesteps for", this->getName());
167 :     }
168 :    
169 :     if (reject != 0)
170 :     handle.logMessage(LOGMESSAGE, "Discarded invalid surveyindex data - number of invalid entries", reject);
171 :     handle.logMessage(LOGMESSAGE, "Read surveyindex data file - number of entries", count);
172 :     }
173 :    
174 :     void SIOnStep::readSIRegressionData(CommentStream& infile) {
175 :    
176 :     char text[MaxStrLength];
177 :     strncpy(text, "", MaxStrLength);
178 :    
179 :     infile >> ws >> text;
180 :     if (strcasecmp(text, "linearfit") == 0) {
181 :     fittype = LINEARFIT;
182 :     LR = new LinearRegression(FREE);
183 :     } else if (strcasecmp(text, "loglinearfit") == 0) {
184 :     fittype = LOGLINEARFIT;
185 :     LR = new LogLinearRegression(FREE);
186 :     } else if (strcasecmp(text, "weightlinearfit") == 0) {
187 :     useweight = 1;
188 :     fittype = WEIGHTLINEARFIT;
189 :     LR = new WeightRegression(FREE);
190 :     } else if (strcasecmp(text, "logweightlinearfit") == 0) {
191 :     useweight = 1;
192 :     fittype = LOGWEIGHTLINEARFIT;
193 :     LR = new LogWeightRegression(FREE);
194 :     } else if (strcasecmp(text, "fixedslopelinearfit") == 0) {
195 :     fittype = FIXEDSLOPELINEARFIT;
196 :     LR = new LinearRegression(FIXEDSLOPE);
197 :     } else if (strcasecmp(text, "fixedslopeloglinearfit") == 0) {
198 :     fittype = FIXEDSLOPELOGLINEARFIT;
199 :     LR = new LogLinearRegression(FIXEDSLOPE);
200 :     } else if (strcasecmp(text, "fixedslopeweightlinearfit") == 0) {
201 :     useweight = 1;
202 :     fittype = FIXEDSLOPEWEIGHTLINEARFIT;
203 :     LR = new WeightRegression(FIXEDSLOPE);
204 :     } else if (strcasecmp(text, "fixedslopelogweightlinearfit") == 0) {
205 :     useweight = 1;
206 :     fittype = FIXEDSLOPELOGWEIGHTLINEARFIT;
207 :     LR = new LogWeightRegression(FIXEDSLOPE);
208 :     } else if (strcasecmp(text, "fixedinterceptlinearfit") == 0) {
209 :     fittype = FIXEDINTERCEPTLINEARFIT;
210 :     LR = new LinearRegression(FIXEDINTERCEPT);
211 :     } else if (strcasecmp(text, "fixedinterceptloglinearfit") == 0) {
212 :     fittype = FIXEDINTERCEPTLOGLINEARFIT;
213 :     LR = new LogLinearRegression(FIXEDINTERCEPT);
214 :     } else if (strcasecmp(text, "fixedinterceptweightlinearfit") == 0) {
215 :     useweight = 1;
216 :     fittype = FIXEDINTERCEPTWEIGHTLINEARFIT;
217 :     LR = new WeightRegression(FIXEDINTERCEPT);
218 :     } else if (strcasecmp(text, "fixedinterceptlogweightlinearfit") == 0) {
219 :     useweight = 1;
220 :     fittype = FIXEDINTERCEPTLOGWEIGHTLINEARFIT;
221 :     LR = new LogWeightRegression(FIXEDINTERCEPT);
222 :     } else if (strcasecmp(text, "fixedlinearfit") == 0) {
223 :     fittype = FIXEDLINEARFIT;
224 :     LR = new LinearRegression(FIXED);
225 :     } else if (strcasecmp(text, "fixedloglinearfit") == 0) {
226 :     fittype = FIXEDLOGLINEARFIT;
227 :     LR = new LogLinearRegression(FIXED);
228 :     } else if (strcasecmp(text, "fixedweightlinearfit") == 0) {
229 :     useweight = 1;
230 :     fittype = FIXEDWEIGHTLINEARFIT;
231 :     LR = new WeightRegression(FIXED);
232 :     } else if (strcasecmp(text, "fixedlogweightlinearfit") == 0) {
233 :     useweight = 1;
234 :     fittype = FIXEDLOGWEIGHTLINEARFIT;
235 :     LR = new LogWeightRegression(FIXED);
236 :     } else
237 :     handle.logFileMessage(LOGFAIL, "\nError in surveyindex - unrecognised fittype", text);
238 :    
239 :     switch (fittype) {
240 :     case LINEARFIT:
241 :     case LOGLINEARFIT:
242 :     case WEIGHTLINEARFIT:
243 :     case LOGWEIGHTLINEARFIT:
244 :     break;
245 :     case FIXEDSLOPELINEARFIT:
246 :     case FIXEDSLOPELOGLINEARFIT:
247 :     case FIXEDSLOPEWEIGHTLINEARFIT:
248 :     case FIXEDSLOPELOGWEIGHTLINEARFIT:
249 :     readWordAndVariable(infile, "slope", slope);
250 :     LR->setSlope(slope);
251 :     break;
252 :     case FIXEDINTERCEPTLINEARFIT:
253 :     case FIXEDINTERCEPTLOGLINEARFIT:
254 :     case FIXEDINTERCEPTWEIGHTLINEARFIT:
255 :     case FIXEDINTERCEPTLOGWEIGHTLINEARFIT:
256 :     readWordAndVariable(infile, "intercept", intercept);
257 :     LR->setIntercept(intercept);
258 :     break;
259 :     case FIXEDLINEARFIT:
260 :     case FIXEDLOGLINEARFIT:
261 :     case FIXEDWEIGHTLINEARFIT:
262 :     case FIXEDLOGWEIGHTLINEARFIT:
263 :     readWordAndVariable(infile, "slope", slope);
264 :     readWordAndVariable(infile, "intercept", intercept);
265 :     LR->setSlope(slope);
266 :     LR->setIntercept(intercept);
267 :     break;
268 :     default:
269 :     handle.logFileMessage(LOGFAIL, "\nError in surveyindex - unrecognised fittype", text);
270 :     break;
271 :     }
272 :    
273 :     //JMB - check that the slope of the regression line is positive
274 :     if (slope < 0.0)
275 :     handle.logFileMessage(LOGFAIL, "\nError in surveyindex - slope of the regression line must be positive", slope);
276 :     }
277 :    
278 :     void SIOnStep::printLikelihood(ofstream& outfile, const TimeClass* const TimeInfo) {
279 :    
280 :     int a, i;
281 :     if (AAT.atCurrentTime(TimeInfo)) {
282 :     timeindex = -1;
283 :     for (i = 0; i < Years.Size(); i++)
284 :     if ((Years[i] == TimeInfo->getYear()) && (Steps[i] == TimeInfo->getStep()))
285 :     timeindex = i;
286 :     if (timeindex == -1)
287 :     handle.logMessage(LOGFAIL, "Error in surveyindex - invalid timestep");
288 :    
289 :     for (a = 0; a < areaindex.Size(); a++) {
290 :     for (i = 0; i < colindex.Size(); i++) {
291 :     outfile << setw(lowwidth) << Years[timeindex] << sep << setw(lowwidth)
292 :     << Steps[timeindex] << sep << setw(printwidth) << areaindex[a] << sep
293 :     << setw(printwidth) << colindex[i] << sep << setw(largewidth);
294 :    
295 :     //JMB crude filter to remove the 'silly' values from the output
296 :     if ((*modelIndex[timeindex])[a][i] < rathersmall)
297 :     outfile << 0;
298 :     else
299 :     outfile << setprecision(largeprecision) << (*modelIndex[timeindex])[a][i];
300 :    
301 :     if (useweight)
302 :     outfile << sep << setw(printwidth) << (*weightIndex[timeindex])[a][i];
303 :     outfile << endl;
304 :     }
305 :     }
306 :     }
307 :    
308 :     //JMB - this is nasty hack to output the regression information
309 :     if (TimeInfo->getTime() == TimeInfo->numTotalSteps()) {
310 :     for (a = 0; a < areaindex.Size(); a++) {
311 :     outfile << "; Regression information for area " << areaindex[a] << endl;
312 :     for (i = 0; i < colindex.Size(); i++)
313 :     outfile << "; " << colindex[i] << " intercept " << intercepts[a][i]
314 :     << " slope " << slopes[a][i] << " sse " << sse[a][i] << endl;
315 :     }
316 :     }
317 :     }
318 :    
319 :     double SIOnStep::calcSSE() {
320 :    
321 :     if (handle.getLogLevel() >= LOGMESSAGE)
322 :     handle.logMessage(LOGMESSAGE, "Calculating likelihood score for surveyindex component", this->getName());
323 :    
324 :     int a, i, j;
325 :     double score = 0.0;
326 :     for (a = 0; a < areaindex.Size(); a++) {
327 :     likelihoodValues[a] = 0.0;
328 :     for (i = 0; i < colindex.Size(); i++) {
329 :     for (j = 0; j < tmpData.Size(); j++) {
330 :     tmpData[j] = (*obsIndex[j])[a][i];
331 :     tmpModel[j] = (*modelIndex[j])[a][i];
332 :     }
333 :    
334 :     //if the weights are required then pass them to the regression line
335 :     if (useweight) {
336 :     for (j = 0; j < tmpWeight.Size(); j++)
337 :     tmpWeight[j] = (*weightIndex[j])[a][i];
338 :     LR->setWeights(tmpWeight);
339 :     }
340 :    
341 :     //fit the data to the (log) linear regression line
342 :     LR->storeVectors(tmpModel, tmpData);
343 :     LR->calcFit();
344 :    
345 :     //and then store the results
346 :     slopes[a][i] = LR->getSlope();
347 :     intercepts[a][i] = LR->getIntercept();
348 :     sse[a][i] = LR->getSSE();
349 :     likelihoodValues[a] += LR->getSSE();
350 :     }
351 :     score += likelihoodValues[a];
352 :     }
353 :    
354 :     if (handle.getLogLevel() >= LOGMESSAGE)
355 :     handle.logMessage(LOGMESSAGE, "The likelihood score from the regression line for this component is", score);
356 :     return score;
357 :     }
358 :    
359 :     void SIOnStep::printSummary(ofstream& outfile, const double weight) {
360 :     int a;
361 :     for (a = 0; a < areaindex.Size(); a++)
362 :     outfile << "all all " << setw(printwidth) << areaindex[a] << sep
363 :     << setw(largewidth) << this->getName() << sep << setprecision(smallprecision)
364 :     << setw(smallwidth) << weight << sep << setprecision(largeprecision)
365 :     << setw(largewidth) << likelihoodValues[a] << endl;
366 :     outfile.flush();
367 :     }
368 :    
369 :     void SIOnStep::Reset() {
370 :     int i;
371 :     for (i = 0; i < modelIndex.Size(); i++)
372 :     (*modelIndex[i]).setToZero();
373 :     }

root@forge.cesga.es
ViewVC Help
Powered by ViewVC 1.0.0  

Powered By FusionForge