Annotation of /trunk/gadget/ecosystem.cc

Parent Directory | Revision Log

---

Revision 11 - (view) (download)

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

---

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