Annotation of /trunk/gadget/ecosystem.cc

Parent Directory | Revision Log

---

Revision 17 - (view) (download)

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

---

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