Annotation of /trunk/gadget/ecosystem.cc

Parent Directory | Revision Log

---

Revision 19 - (view) (download)

1 :	agomez	1	#include "ecosystem.h"
2 :			#include "runid.h"
3 :			#include "global.h"
4 :	ulcessvp	19	#include <ctime>
5 :	ulcessvp	16	#ifdef _OPENMP
6 :	ulcessvp	11	#include <omp.h>
7 :			#endif
8 :	agomez	1
9 :	ulcessvp	16	#ifdef _OPENMP
10 :	ulcessvp	11	extern Ecosystem** EcoSystems;
11 :			//extern StochasticData* data;
12 :			#endif
13 :	ulcessvp	17	extern volatile int interrupted_print;
14 :	ulcessvp	11
15 :	ulcessvp	19	time_t starttime, stoptime;
16 :
17 :	agomez	1	Ecosystem::Ecosystem(const MainInfo& main) : printinfo(main.getPI()) {
18 :
19 :			funceval = 0;
20 :	ulcessvp	17	//interrupted = 0;
21 :	agomez	1	likelihood = 0.0;
22 :			keeper = new Keeper;
23 :
24 :			// initialise counter used when printing output files
25 :			printcount = printinfo.getPrintIteration() - 1;
26 :
27 :			// read the model specification from the main file
28 :			char* filename = main.getMainGadgetFile();
29 :			ifstream infile;
30 :			infile.open(filename, ios::in);
31 :			CommentStream commin(infile);
32 :			handle.checkIfFailure(infile, filename);
33 :			handle.Open(filename);
34 :			this->readMain(commin, main);
35 :			handle.Close();
36 :			infile.close();
37 :			infile.clear();
38 :
39 :			// if this is an optimising run then read the optimisation parameters from file
40 :			if (main.runOptimise()) {
41 :			handle.logMessage(LOGMESSAGE, ""); //write blank line to log file
42 :			if (main.getOptInfoGiven()) {
43 :			filename = main.getOptInfoFile();
44 :			infile.open(filename, ios::in);
45 :			handle.checkIfFailure(infile, filename);
46 :			handle.Open(filename);
47 :	ulcessvp	11	unsigned seed;
48 :			main.getForcePrint();
49 :			this->readOptimisation(commin, main.getSeed());
50 :	agomez	1	handle.Close();
51 :			infile.close();
52 :			infile.clear();
53 :			} else {
54 :			handle.logMessage(LOGINFO, "Warning - no optimisation file specified, using default values");
55 :			optvec.resize(new OptInfoHooke());
56 :			}
57 :			}
58 :
59 :			if (main.runOptimise())
60 :			handle.logMessage(LOGINFO, "\nFinished reading model data files, starting to run optimisation");
61 :			else
62 :			handle.logMessage(LOGINFO, "\nFinished reading model data files, starting to run simulation");
63 :			handle.logMessage(LOGMESSAGE, ""); //write blank line to log file
64 :			}
65 :
66 :			Ecosystem::~Ecosystem() {
67 :			int i;
68 :			for (i = 0; i < optvec.Size(); i++)
69 :			delete optvec[i];
70 :			for (i = 0; i < printvec.Size(); i++)
71 :			delete printvec[i];
72 :			for (i = 0; i < likevec.Size(); i++)
73 :			delete likevec[i];
74 :			for (i = 0; i < tagvec.Size(); i++)
75 :			delete tagvec[i];
76 :			for (i = 0; i < basevec.Size(); i++)
77 :			delete basevec[i];
78 :
79 :			delete Area;
80 :			delete TimeInfo;
81 :			delete keeper;
82 :			}
83 :
84 :			void Ecosystem::writeStatus(const char* filename) const {
85 :			ofstream outfile;
86 :			outfile.open(filename, ios::out);
87 :			handle.checkIfFailure(outfile, filename);
88 :			handle.Open(filename);
89 :			RUNID.Print(outfile);
90 :			outfile << "The current simulation time is " << TimeInfo->getYear()
91 :			<< ", step " << TimeInfo->getStep() << endl;
92 :
93 :			int i;
94 :			for (i = 0; i < basevec.Size(); i++)
95 :			basevec[i]->Print(outfile);
96 :			for (i = 0; i < likevec.Size(); i++)
97 :			likevec[i]->Print(outfile);
98 :			for (i = 0; i < tagvec.Size(); i++)
99 :			tagvec[i]->Print(outfile);
100 :
101 :			handle.Close();
102 :			outfile.close();
103 :			outfile.clear();
104 :			}
105 :
106 :			void Ecosystem::Reset() {
107 :			int i;
108 :			TimeInfo->Reset();
109 :			for (i = 0; i < likevec.Size(); i++)
110 :			likevec[i]->Reset(keeper);
111 :			for (i = 0; i < basevec.Size(); i++)
112 :			basevec[i]->Reset(TimeInfo);
113 :			}
114 :
115 :			void Ecosystem::Optimise() {
116 :			int i;
117 :			for (i = 0; i < optvec.Size(); i++) {
118 :	ulcessvp	19	time(&starttime);
119 :	ulcessvp	16	#ifdef _OPENMP
120 :	ulcessvp	11	int j, numThr = omp_get_max_threads ( );
121 :			DoubleVector v = this->getValues();
122 :			for (j = 0; j < numThr; j++)
123 :			EcoSystems[j]->Update(v);
124 :			#endif
125 :	ulcessvp	15	#ifdef SPECULATIVE
126 :	ulcessvp	12	optvec[i]->OptimiseLikelihoodOMP();
127 :			#else
128 :			optvec[i]->OptimiseLikelihood();
129 :			#endif
130 :	ulcessvp	19	time(&stoptime);
131 :	ulcessvp	11	this->writeOptValues();
132 :			}
133 :			}
134 :
135 :	agomez	1	double Ecosystem::SimulateAndUpdate(const DoubleVector& x) {
136 :			int i, j;
137 :
138 :			if (funceval == 0) {
139 :			// JMB - only need to create these vectors once
140 :			initialval.resize(keeper->numVariables(), 0.0);
141 :			currentval.resize(keeper->numVariables(), 0.0);
142 :			optflag.resize(keeper->numVariables(), 0);
143 :			keeper->getOptFlags(optflag);
144 :			}
145 :
146 :			j = 0;
147 :			keeper->getCurrentValues(currentval);
148 :			keeper->getInitialValues(initialval);
149 :			for (i = 0; i < currentval.Size(); i++) {
150 :			if (optflag[i]) {
151 :			currentval[i] = x[j] * initialval[i];
152 :			j++;
153 :			}
154 :			}
155 :
156 :			keeper->Update(currentval);
157 :			this->Simulate(0); //dont print whilst optimising
158 :
159 :			if (printinfo.getPrint()) {
160 :			printcount++;
161 :			if (printcount == printinfo.getPrintIteration()) {
162 :			keeper->writeValues(likevec, printinfo.getPrecision());
163 :			printcount = 0;
164 :			}
165 :			}
166 :
167 :			funceval++;
168 :			return likelihood;
169 :			}
170 :
171 :			void Ecosystem::writeOptValues() {
172 :			int i;
173 :			DoubleVector tmpvec(likevec.Size(), 0.0);
174 :			for (i = 0; i < likevec.Size(); i++)
175 :			tmpvec[i] = likevec[i]->getUnweightedLikelihood();
176 :
177 :	ulcessvp	11	int iters = 0;
178 :	ulcessvp	16	#ifdef _OPENMP
179 :	ulcessvp	11	int numThr = omp_get_max_threads ( );
180 :			for (i = 0; i < numThr; i++)
181 :			iters += EcoSystems[i]->getFuncEval();
182 :			#endif
183 :			handle.logMessage(LOGINFO, "\nAfter a total of", funceval+iters, "function evaluations the best point found is");
184 :	agomez	1	keeper->writeBestValues();
185 :			handle.logMessage(LOGINFO, "\nThe scores from each likelihood component are");
186 :			handle.logMessage(LOGINFO, tmpvec);
187 :			if (!isZero(keeper->getBestLikelihoodScore())) // no better point has been found
188 :			handle.logMessage(LOGINFO, "\nThe overall likelihood score is", keeper->getBestLikelihoodScore());
189 :			else
190 :			handle.logMessage(LOGINFO, "\nThe overall likelihood score is", this->getLikelihood());
191 :	ulcessvp	19
192 :			handle.logMessage(LOGINFO, "\n Runtime for the algorithm was ", difftime(stoptime, starttime)," seconds\n");
193 :	agomez	1	}
194 :
195 :			void Ecosystem::writeInitialInformation(const char* const filename) {
196 :			keeper->openPrintFile(filename);
197 :			keeper->writeInitialInformation(likevec);
198 :			}
199 :
200 :			void Ecosystem::writeValues() {
201 :			keeper->writeValues(likevec, printinfo.getPrecision());
202 :			}
203 :
204 :			void Ecosystem::writeParams(const char* const filename, int prec) const {
205 :	ulcessvp	17	if ((funceval > 0) && (interrupted_print == 0)) {
206 :	agomez	1	//JMB - print the final values to any output files specified
207 :			//in case they have been missed by the -print value
208 :			if (printinfo.getPrint())
209 :			keeper->writeValues(likevec, printinfo.getPrecision());
210 :			}
211 :	ulcessvp	17	keeper->writeParams(optvec, filename, prec, interrupted_print);
212 :	agomez	1	}

---

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