1 /*
2 * cResourceCount.cc
3 * Avida
4 *
5 * Called "resource_count.cc" prior to 12/5/05.
6 * Copyright 1999-2011 Michigan State University. All rights reserved.
7 * Copyright 1993-2001 California Institute of Technology.
8 *
9 *
10 * This file is part of Avida.
11 *
12 * Avida is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
14 *
15 * Avida is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.
17 *
18 * You should have received a copy of the GNU Lesser General Public License along with Avida.
19 * If not, see <http://www.gnu.org/licenses/>.
20 *
21 */
22
23 #include "cResourceCount.h"
24 #include "cResource.h"
25 #include "cDynamicCount.h"
26 #include "cGradientCount.h"
27 #include "cWorld.h"
28 #include "cStats.h"
29
30 #include "nGeometry.h"
31
32 #include <cmath>
33
34 using namespace std;
35
36 const double cResourceCount::UPDATE_STEP(1.0 / 10000.0);
37 const double cResourceCount::EPSILON (1.0e-15);
38 const int cResourceCount::PRECALC_DISTANCE(100);
39
40
FlowMatter(cSpatialCountElem & elem1,cSpatialCountElem & elem2,double inxdiffuse,double inydiffuse,double inxgravity,double inygravity,int xdist,int ydist,double dist)41 void FlowMatter(cSpatialCountElem &elem1, cSpatialCountElem &elem2,
42 double inxdiffuse,
43 double inydiffuse, double inxgravity, double inygravity,
44 int xdist, int ydist, double dist) {
45
46 /* Routine to calculate the amount of flow from one Element to another.
47 Amount of flow is a function of:
48
49 1) Amount of material in each cell (will try to equalize)
50 2) Distance between each cell
51 3) x and y "gravity"
52
53 This method only effect the delta amount of each element. The State
54 method will need to be called at the end of each time step to complete
55 the movement of material.
56 */
57
58 double diff, flowamt, xgravity, xdiffuse, ygravity, ydiffuse;
59
60 if (((elem1.amount == 0.0) && (elem2.amount == 0.0)) && (dist < 0.0)) return;
61 diff = (elem1.amount - elem2.amount);
62 if (xdist != 0) {
63
64 /* if there is material to be effected by x gravity */
65
66 if (((xdist>0) && (inxgravity>0.0)) || ((xdist<0) && (inxgravity<0.0))) {
67 xgravity = elem1.amount * fabs(inxgravity)/3.0;
68 } else {
69 xgravity = -elem2.amount * fabs(inxgravity)/3.0;
70 }
71
72 /* Diffusion uses the diffusion constant x half the difference (as the
73 elements attempt to equalize) / the number of possible neighbors (8) */
74
75 xdiffuse = inxdiffuse * diff / 16.0;
76 } else {
77 xdiffuse = 0.0;
78 xgravity = 0.0;
79 }
80 if (ydist != 0) {
81
82 /* if there is material to be effected by y gravity */
83
84 if (((ydist>0) && (inygravity>0.0)) || ((ydist<0) && (inygravity<0.0))) {
85 ygravity = elem1.amount * fabs(inygravity)/3.0;
86 } else {
87 ygravity = -elem2.amount * fabs(inygravity)/3.0;
88 }
89 ydiffuse = inydiffuse * diff / 16.0;
90 } else {
91 ydiffuse = 0.0;
92 ygravity = 0.0;
93 }
94
95 flowamt = ((xdiffuse + ydiffuse + xgravity + ygravity)/
96 (fabs(xdist*1.0) + fabs(ydist*1.0)))/dist;
97 elem1.delta -= flowamt;
98 elem2.delta += flowamt;
99 }
100
cResourceCount(int num_resources)101 cResourceCount::cResourceCount(int num_resources)
102 : update_time(0.0)
103 , spatial_update_time(0.0)
104 , m_last_updated(0)
105 , m_spatial_update(0)
106 {
107 if(num_resources > 0) {
108 SetSize(num_resources);
109 }
110
111 return;
112 }
113
cResourceCount(const cResourceCount & rc)114 cResourceCount::cResourceCount(const cResourceCount &rc) {
115 *this = rc;
116
117 return;
118 }
119
operator =(const cResourceCount & rc)120 const cResourceCount &cResourceCount::operator=(const cResourceCount &rc) {
121 SetSize(rc.GetSize());
122 resource_name = rc.resource_name;
123 resource_initial = rc.resource_initial;
124 resource_count = rc.resource_count;
125 decay_rate = rc.decay_rate;
126 inflow_rate = rc.inflow_rate;
127 decay_precalc = rc.decay_precalc;
128 inflow_precalc = rc.inflow_precalc;
129 geometry = rc.geometry;
130
131 for (int i = 0; i < rc.spatial_resource_count.GetSize(); i++) {
132 *(spatial_resource_count[i]) = *(rc.spatial_resource_count[i]);
133 }
134
135 curr_grid_res_cnt = rc.curr_grid_res_cnt;
136 curr_spatial_res_cnt = rc.curr_spatial_res_cnt;
137 update_time = rc.update_time;
138 spatial_update_time = rc.spatial_update_time;
139 cell_lists = rc.cell_lists;
140
141 return *this;
142 }
143
SetSize(int num_resources)144 void cResourceCount::SetSize(int num_resources)
145 {
146 resource_name.ResizeClear(num_resources);
147 resource_initial.ResizeClear(num_resources);
148 resource_count.ResizeClear(num_resources);
149 decay_rate.ResizeClear(num_resources);
150 inflow_rate.ResizeClear(num_resources);
151 if(num_resources > 0) {
152 decay_precalc.ResizeClear(num_resources, PRECALC_DISTANCE+1);
153 inflow_precalc.ResizeClear(num_resources, PRECALC_DISTANCE+1);
154 }
155 geometry.ResizeClear(num_resources);
156
157 for(int i = 0; i < spatial_resource_count.GetSize(); i++){
158 delete spatial_resource_count[i];
159 }
160
161 spatial_resource_count.ResizeClear(num_resources);
162
163 for(int i = 0; i < num_resources; i++){
164 spatial_resource_count[i] = new cSpatialResCount();
165 }
166
167 curr_grid_res_cnt.ResizeClear(num_resources);
168 curr_spatial_res_cnt.ResizeClear(num_resources);
169 cell_lists.ResizeClear(num_resources);
170 resource_name.SetAll("");
171 resource_initial.SetAll(0.0);
172 resource_count.SetAll(0.0);
173 decay_rate.SetAll(0.0);
174 inflow_rate.SetAll(0.0);
175 decay_precalc.SetAll(1.0); // This is 1-inflow, so there should be no inflow by default, JEB
176 inflow_precalc.SetAll(0.0);
177 geometry.SetAll(nGeometry::GLOBAL);
178 curr_grid_res_cnt.SetAll(0.0);
179 //DO spacial resources need to be set to zero?
180 }
181
~cResourceCount()182 cResourceCount::~cResourceCount()
183 {
184 for(int i = 0; i < spatial_resource_count.GetSize(); i++){
185 delete spatial_resource_count[i];
186 }
187 }
188
SetCellResources(int cell_id,const tArray<double> & res)189 void cResourceCount::SetCellResources(int cell_id, const tArray<double> & res)
190 {
191 assert(resource_count.GetSize() == res.GetSize());
192
193 for (int i = 0; i < resource_count.GetSize(); i++) {
194 if (geometry[i] == nGeometry::GLOBAL || geometry[i]==nGeometry::PARTIAL) {
195 // Set global quantity of resource
196 } else {
197 spatial_resource_count[i]->SetCellAmount(cell_id, res[i]);
198
199 /* Ideally the state of the cell's resource should not be set till
200 the end of the update so that all processes (inflow, outflow,
201 diffision, gravity and organism demand) have the same weight. However
202 waiting can cause problems with negative resources so we allow
203 the organism demand to work immediately on the state of the resource */
204 }
205 }
206 }
207
Setup(cWorld * world,const int & res_index,const cString & name,const double & initial,const double & inflow,const double & decay,const int & in_geometry,const double & in_xdiffuse,const double & in_xgravity,const double & in_ydiffuse,const double & in_ygravity,const int & in_inflowX1,const int & in_inflowX2,const int & in_inflowY1,const int & in_inflowY2,const int & in_outflowX1,const int & in_outflowX2,const int & in_outflowY1,const int & in_outflowY2,tArray<cCellResource> * in_cell_list_ptr,tArray<int> * in_cell_id_list_ptr,const int & verbosity_level,const bool & isdynamic,const int & in_peaks,const double & in_min_height,const double & in_min_radius,const double & in_radius_range,const double & in_ah,const double & in_ar,const double & in_acx,const double & in_acy,const double & in_hstepscale,const double & in_rstepscale,const double & in_cstepscalex,const double & in_cstepscaley,const double & in_hstep,const double & in_rstep,const double & in_cstepx,const double & in_cstepy,const int & in_update_dynamic,const int & in_peakx,const int & in_peaky,const int & in_height,const int & in_spread,const double & in_plateau,const int & in_decay,const int & in_max_x,const int & in_min_x,const int & in_max_y,const int & in_min_y,const double & in_move_a_scaler,const int & in_updatestep,const int & in_halo,const int & in_halo_inner_radius,const int & in_halo_width,const int & in_halo_anchor_x,const int & in_halo_anchor_y,const int & in_move_speed,const double & in_plateau_inflow,const double & in_plateau_outflow,const double & in_cone_inflow,const double & in_cone_outflow,const double & in_gradient_inflow,const int & in_is_plateau_common,const double & in_floor,const int & in_habitat,const int & in_min_size,const int & in_max_size,const int & in_config,const int & in_count,const double & in_resistance,const double & in_init_plat,const double & in_threshold,const int & in_refuge,const bool & isgradient)208 void cResourceCount::Setup(cWorld* world, const int& res_index, const cString& name, const double& initial, const double& inflow, const double& decay,
209 const int& in_geometry, const double& in_xdiffuse, const double& in_xgravity,
210 const double& in_ydiffuse, const double& in_ygravity,
211 const int& in_inflowX1, const int& in_inflowX2, const int& in_inflowY1, const int& in_inflowY2,
212 const int& in_outflowX1, const int& in_outflowX2, const int& in_outflowY1,
213 const int& in_outflowY2, tArray<cCellResource> *in_cell_list_ptr,
214 tArray<int> *in_cell_id_list_ptr, const int& verbosity_level,
215 const bool& isdynamic, const int& in_peaks,
216 const double& in_min_height, const double& in_min_radius, const double& in_radius_range,
217 const double& in_ah, const double& in_ar,
218 const double& in_acx, const double& in_acy,
219 const double& in_hstepscale, const double& in_rstepscale,
220 const double& in_cstepscalex, const double& in_cstepscaley,
221 const double& in_hstep, const double& in_rstep,
222 const double& in_cstepx, const double& in_cstepy,
223 const int& in_update_dynamic, const int& in_peakx, const int& in_peaky,
224 const int& in_height, const int& in_spread, const double& in_plateau, const int& in_decay,
225 const int& in_max_x, const int& in_min_x, const int& in_max_y, const int& in_min_y, const double& in_move_a_scaler,
226 const int& in_updatestep, const int& in_halo, const int& in_halo_inner_radius, const int& in_halo_width,
227 const int& in_halo_anchor_x, const int& in_halo_anchor_y, const int& in_move_speed,
228 const double& in_plateau_inflow, const double& in_plateau_outflow, const double& in_cone_inflow, const double& in_cone_outflow,
229 const double& in_gradient_inflow, const int& in_is_plateau_common, const double& in_floor, const int& in_habitat,
230 const int& in_min_size, const int& in_max_size, const int& in_config, const int& in_count, const double& in_resistance,
231 const double& in_init_plat, const double& in_threshold, const int& in_refuge, const bool& isgradient)
232 {
233 assert(res_index >= 0 && res_index < resource_count.GetSize());
234 assert(initial >= 0.0);
235 assert(decay >= 0.0);
236 assert(inflow >= 0.0);
237 assert(spatial_resource_count[res_index]->GetSize() > 0);
238 int tempx = spatial_resource_count[res_index]->GetX();
239 int tempy = spatial_resource_count[res_index]->GetY();
240
241 cString geo_name;
242 if (in_geometry == nGeometry::GLOBAL) {
243 geo_name = "GLOBAL";
244 } else if (in_geometry == nGeometry::GRID) {
245 geo_name = "GRID";
246 } else if (in_geometry == nGeometry::TORUS) {
247 geo_name = "TORUS";
248 }
249 else if (in_geometry == nGeometry::PARTIAL) {
250 geo_name = "PARTIAL";
251 }
252 else {
253 cerr << "[cResourceCount::Setup] Unknown resource geometry " << in_geometry << ". Exiting.";
254 exit(2);
255 }
256
257
258 /* If the verbose flag is set print out information about resources */
259 verbosity = verbosity_level;
260 if (verbosity > VERBOSE_NORMAL) {
261 cout << "Setting up resource " << name
262 << "(" << geo_name
263 << ") with initial quatity=" << initial
264 << ", decay=" << decay
265 << ", inflow=" << inflow
266 << endl;
267 if ((in_geometry == nGeometry::GRID) || (in_geometry == nGeometry::TORUS)) {
268 cout << " Inflow rectangle (" << in_inflowX1
269 << "," << in_inflowY1
270 << ") to (" << in_inflowX2
271 << "," << in_inflowY2
272 << ")" << endl;
273 cout << " Outflow rectangle (" << in_outflowX1
274 << "," << in_outflowY1
275 << ") to (" << in_outflowX2
276 << "," << in_outflowY2
277 << ")" << endl;
278 cout << " xdiffuse=" << in_xdiffuse
279 << ", xgravity=" << in_xgravity
280 << ", ydiffuse=" << in_ydiffuse
281 << ", ygravity=" << in_ygravity
282 << endl;
283 }
284 }
285
286
287 /* recource_count gets only the values for global resources */
288
289 resource_name[res_index] = name;
290 resource_initial[res_index] = initial;
291 if (in_geometry == nGeometry::GLOBAL) {
292 resource_count[res_index] = initial;
293 spatial_resource_count[res_index]->RateAll(0);
294 }
295 else if (in_geometry == nGeometry::PARTIAL) {
296 resource_count[res_index]=initial;
297
298 spatial_resource_count[res_index]->RateAll(0);
299 // want to set list of cell ids here
300 cell_lists[res_index].Resize(in_cell_id_list_ptr->GetSize());
301 for (int i = 0; i < in_cell_id_list_ptr->GetSize(); i++)
302 cell_lists[res_index][i] = (*in_cell_id_list_ptr)[i];
303
304 }
305 else {
306 resource_count[res_index] = 0;
307 if(isdynamic){ //JW
308 delete spatial_resource_count[res_index];
309 spatial_resource_count[res_index] = new cDynamicCount(in_peaks, in_min_height, in_radius_range, in_min_radius, in_ah, in_ar,
310 in_acx, in_acy, in_hstepscale, in_rstepscale, in_cstepscalex, in_cstepscaley, in_hstep, in_rstep,
311 in_cstepx, in_cstepy, tempx, tempy, in_geometry, in_updatestep);
312 spatial_resource_count[res_index]->RateAll(0);
313 }
314
315 else if(isgradient){
316 delete spatial_resource_count[res_index];
317 spatial_resource_count[res_index] = new cGradientCount(world, in_peakx, in_peaky, in_height, in_spread, in_plateau, in_decay,
318 in_max_x, in_max_y, in_min_x, in_min_y, in_move_a_scaler, in_updatestep,
319 tempx, tempy, in_geometry, in_halo, in_halo_inner_radius,
320 in_halo_width, in_halo_anchor_x, in_halo_anchor_y, in_move_speed,
321 in_plateau_inflow, in_plateau_outflow, in_cone_inflow, in_cone_outflow,
322 in_gradient_inflow, in_is_plateau_common, in_floor, in_habitat,
323 in_min_size, in_max_size, in_config, in_count, in_init_plat);
324 spatial_resource_count[res_index]->RateAll(0);
325 }
326
327 else{
328 spatial_resource_count[res_index]->SetInitial(initial / spatial_resource_count[res_index]->GetSize());
329 spatial_resource_count[res_index]->RateAll(spatial_resource_count[res_index]->GetInitial());
330 }
331 }
332 spatial_resource_count[res_index]->StateAll();
333 decay_rate[res_index] = decay;
334 inflow_rate[res_index] = inflow;
335 geometry[res_index] = in_geometry;
336 spatial_resource_count[res_index]->SetGeometry(in_geometry);
337 spatial_resource_count[res_index]->SetPointers();
338 spatial_resource_count[res_index]->SetCellList(in_cell_list_ptr);
339
340 double step_decay = pow(decay, UPDATE_STEP);
341 double step_inflow = inflow * UPDATE_STEP;
342
343 decay_precalc(res_index, 0) = 1.0;
344 inflow_precalc(res_index, 0) = 0.0;
345 for (int i = 1; i <= PRECALC_DISTANCE; i++) {
346 decay_precalc(res_index, i) = decay_precalc(res_index, i-1) * step_decay;
347 inflow_precalc(res_index, i) = inflow_precalc(res_index, i-1) * step_decay + step_inflow;
348 }
349 spatial_resource_count[res_index]->SetXdiffuse(in_xdiffuse);
350 spatial_resource_count[res_index]->SetXgravity(in_xgravity);
351 spatial_resource_count[res_index]->SetYdiffuse(in_ydiffuse);
352 spatial_resource_count[res_index]->SetYgravity(in_ygravity);
353 spatial_resource_count[res_index]->SetInflowX1(in_inflowX1);
354 spatial_resource_count[res_index]->SetInflowX2(in_inflowX2);
355 spatial_resource_count[res_index]->SetInflowY1(in_inflowY1);
356 spatial_resource_count[res_index]->SetInflowY2(in_inflowY2);
357 spatial_resource_count[res_index]->SetOutflowX1(in_outflowX1);
358 spatial_resource_count[res_index]->SetOutflowX2(in_outflowX2);
359 spatial_resource_count[res_index]->SetOutflowY1(in_outflowY1);
360 spatial_resource_count[res_index]->SetOutflowY2(in_outflowY2);
361 }
362
SetGradientCount(cAvidaContext & ctx,cWorld * world,const int & res_id,const int & peakx,const int & peaky,const int & height,const int & spread,const double & plateau,const int & decay,const int & max_x,const int & min_x,const int & max_y,const int & min_y,const double & move_a_scaler,const int & updatestep,const int & halo,const int & halo_inner_radius,const int & halo_width,const int & halo_anchor_x,const int & halo_anchor_y,const int & move_speed,const double & plateau_inflow,const double & plateau_outflow,const double & cone_inflow,const double & cone_outflow,const double & gradient_inflow,const int & is_plateau_common,const double & floor,const int & habitat,const int & min_size,const int & max_size,const int & config,const int & count,const double & resistance,const double & plat_val,const double & threshold,const int & refuge)363 void cResourceCount::SetGradientCount(cAvidaContext& ctx, cWorld* world, const int& res_id, const int& peakx, const int& peaky,
364 const int& height, const int& spread, const double& plateau, const int& decay,
365 const int& max_x, const int& min_x, const int& max_y, const int& min_y, const double& move_a_scaler,
366 const int& updatestep, const int& halo, const int& halo_inner_radius, const int& halo_width,
367 const int& halo_anchor_x, const int& halo_anchor_y, const int& move_speed,
368 const double& plateau_inflow, const double& plateau_outflow, const double& cone_inflow, const double& cone_outflow,
369 const double& gradient_inflow, const int& is_plateau_common, const double& floor, const int& habitat,
370 const int& min_size, const int& max_size, const int& config, const int& count, const double& resistance,
371 const double& plat_val, const double& threshold, const int& refuge)
372 {
373 assert(res_id >= 0 && res_id < resource_count.GetSize());
374 assert(spatial_resource_count[res_id]->GetSize() > 0);
375 int worldx = spatial_resource_count[res_id]->GetX();
376 int worldy = spatial_resource_count[res_id]->GetY();
377
378 spatial_resource_count[res_id]->SetGradPeakX(peakx);
379 spatial_resource_count[res_id]->SetGradPeakY(peaky);
380 spatial_resource_count[res_id]->SetGradHeight(height);
381 spatial_resource_count[res_id]->SetGradSpread(spread);
382 spatial_resource_count[res_id]->SetGradPlateau(plateau);
383 spatial_resource_count[res_id]->SetGradInitialPlatVal(plat_val);
384 spatial_resource_count[res_id]->SetGradDecay(decay);
385 spatial_resource_count[res_id]->SetGradMaxX(max_x);
386 spatial_resource_count[res_id]->SetGradMaxY(max_y);
387 spatial_resource_count[res_id]->SetGradMinX(min_x);
388 spatial_resource_count[res_id]->SetGradMinY(min_y);
389 spatial_resource_count[res_id]->SetGradMoveScaler(move_a_scaler);
390 spatial_resource_count[res_id]->SetGradUpdateStep(updatestep);
391
392 spatial_resource_count[res_id]->SetGradIsHalo(halo);
393 spatial_resource_count[res_id]->SetGradHaloInnerRad(halo_inner_radius);
394 spatial_resource_count[res_id]->SetGradHaloWidth(halo_width);
395 spatial_resource_count[res_id]->SetGradHaloX(halo_anchor_x);
396 spatial_resource_count[res_id]->SetGradHaloY(halo_anchor_y);
397 spatial_resource_count[res_id]->SetGradMoveSpeed(move_speed);
398 spatial_resource_count[res_id]->SetGradPlatInflow(plateau_inflow);
399 spatial_resource_count[res_id]->SetGradPlatOutflow(plateau_outflow);
400 spatial_resource_count[res_id]->SetGradConeInflow(cone_inflow);
401 spatial_resource_count[res_id]->SetGradConeOutflow(cone_outflow);
402 spatial_resource_count[res_id]->SetGradientInflow(gradient_inflow);
403 spatial_resource_count[res_id]->SetGradPlatIsCommon(is_plateau_common);
404 spatial_resource_count[res_id]->SetGradFloor(floor);
405 spatial_resource_count[res_id]->SetGradHabitat(habitat);
406 spatial_resource_count[res_id]->SetGradMinSize(min_size);
407 spatial_resource_count[res_id]->SetGradMaxSize(max_size);
408 spatial_resource_count[res_id]->SetGradConfig(config);
409 spatial_resource_count[res_id]->SetGradCount(count);
410 spatial_resource_count[res_id]->SetGradResistance(resistance);
411 spatial_resource_count[res_id]->SetGradThreshold(threshold);
412 spatial_resource_count[res_id]->SetGradRefuge(refuge);
413
414 spatial_resource_count[res_id]->ResetGradRes(ctx, worldx, worldy);
415 }
416
SetGradientPlatInflow(const int & res_id,const double & inflow)417 void cResourceCount::SetGradientPlatInflow(const int& res_id, const double& inflow)
418 {
419 assert(res_id >= 0 && res_id < resource_count.GetSize());
420 assert(spatial_resource_count[res_id]->GetSize() > 0);
421 spatial_resource_count[res_id]->SetGradPlatInflow(inflow);
422 }
423
SetGradientPlatOutflow(const int & res_id,const double & outflow)424 void cResourceCount::SetGradientPlatOutflow(const int& res_id, const double& outflow)
425 {
426 assert(res_id >= 0 && res_id < resource_count.GetSize());
427 assert(spatial_resource_count[res_id]->GetSize() > 0);
428 spatial_resource_count[res_id]->SetGradPlatOutflow(outflow);
429 }
430
SetGradientConeInflow(const int & res_id,const double & inflow)431 void cResourceCount::SetGradientConeInflow(const int& res_id, const double& inflow)
432 {
433 assert(res_id >= 0 && res_id < resource_count.GetSize());
434 assert(spatial_resource_count[res_id]->GetSize() > 0);
435 spatial_resource_count[res_id]->SetGradConeInflow(inflow);
436 }
437
SetGradientConeOutflow(const int & res_id,const double & outflow)438 void cResourceCount::SetGradientConeOutflow(const int& res_id, const double& outflow)
439 {
440 assert(res_id >= 0 && res_id < resource_count.GetSize());
441 assert(spatial_resource_count[res_id]->GetSize() > 0);
442 spatial_resource_count[res_id]->SetGradConeOutflow(outflow);
443 }
444
SetGradientInflow(const int & res_id,const double & inflow)445 void cResourceCount::SetGradientInflow(const int& res_id, const double& inflow)
446 {
447 assert(res_id >= 0 && res_id < resource_count.GetSize());
448 assert(spatial_resource_count[res_id]->GetSize() > 0);
449 spatial_resource_count[res_id]->SetGradientInflow(inflow);
450 }
451
SetGradPlatVarInflow(const int & res_id,const double & mean,const double & variance,const int & type)452 void cResourceCount::SetGradPlatVarInflow(const int& res_id, const double& mean, const double& variance, const int& type)
453 {
454 assert(res_id >= 0 && res_id < resource_count.GetSize());
455 assert(spatial_resource_count[res_id]->GetSize() > 0);
456 spatial_resource_count[res_id]->SetGradPlatVarInflow(mean, variance, type);
457 }
458
SetPredatoryResource(const int & res_id,const double & odds,const int & juvsper)459 void cResourceCount::SetPredatoryResource(const int& res_id, const double& odds, const int& juvsper)
460 {
461 assert(res_id >= 0 && res_id < resource_count.GetSize());
462 assert(spatial_resource_count[res_id]->GetSize() > 0);
463 spatial_resource_count[res_id]->SetPredatoryResource(odds, juvsper);
464 }
465
SetProbabilisticResource(cAvidaContext & ctx,const int & res_id,const double & initial,const double & inflow,const double & outflow,const double & lambda,const double & theta,const int & x,const int & y,const int & count)466 void cResourceCount::SetProbabilisticResource(cAvidaContext& ctx, const int& res_id, const double& initial, const double& inflow,
467 const double& outflow, const double& lambda, const double& theta, const int& x, const int& y, const int& count)
468 {
469 assert(res_id >= 0 && res_id < resource_count.GetSize());
470 assert(spatial_resource_count[res_id]->GetSize() > 0);
471 spatial_resource_count[res_id]->SetProbabilisticResource(ctx, initial, inflow, outflow, lambda, theta, x, y, count);
472 }
473
474 /*
475 * This is unnecessary now that a resource has an index
476 * TODO:
477 * - Change name to GetResourceCountIndex
478 * - Fix anything that breaks by just using the index of the resource (not id)
479 * - Get rid of this function
480 */
GetResourceCountID(const cString & res_name)481 int cResourceCount::GetResourceCountID(const cString& res_name)
482 {
483 for (int i = 0; i < resource_name.GetSize(); i++) {
484 if (resource_name[i] == res_name) return i;
485 }
486 cerr << "Error: Unknown resource '" << res_name << "'." << endl;
487 return -1;
488 }
489
GetInflow(const cString & name)490 double cResourceCount::GetInflow(const cString& name)
491 {
492 int id = GetResourceCountID(name);
493 if (id == -1) return -1;
494
495 return inflow_rate[id];
496 }
497
SetInflow(const cString & name,const double _inflow)498 void cResourceCount::SetInflow(const cString& name, const double _inflow)
499 {
500 int id = GetResourceCountID(name);
501 if (id == -1) return;
502
503 inflow_rate[id] = _inflow;
504 double step_inflow = _inflow * UPDATE_STEP;
505 double step_decay = pow(decay_rate[id], UPDATE_STEP);
506
507 inflow_precalc(id, 0) = 0.0;
508 for (int i = 1; i <= PRECALC_DISTANCE; i++) {
509 inflow_precalc(id, i) = inflow_precalc(id, i-1) * step_decay + step_inflow;
510 }
511 }
512
GetDecay(const cString & name)513 double cResourceCount::GetDecay(const cString& name)
514 {
515 int id = GetResourceCountID(name);
516 if (id == -1) return -1;
517
518 return decay_rate[id];
519 }
520
SetDecay(const cString & name,const double _decay)521 void cResourceCount::SetDecay(const cString& name, const double _decay)
522 {
523 int id = GetResourceCountID(name);
524 if (id == -1) return;
525
526 decay_rate[id] = _decay;
527 double step_decay = pow(_decay, UPDATE_STEP);
528 decay_precalc(id, 0) = 1.0;
529 for (int i = 1; i <= PRECALC_DISTANCE; i++) {
530 decay_precalc(id, i) = decay_precalc(id, i-1) * step_decay;
531 }
532 }
533
Update(double in_time)534 void cResourceCount::Update(double in_time)
535 {
536 update_time += in_time;
537 spatial_update_time += in_time;
538 }
539
540
GetResources(cAvidaContext & ctx) const541 const tArray<double> & cResourceCount::GetResources(cAvidaContext& ctx) const
542 {
543 DoUpdates(ctx);
544 return resource_count;
545 }
546
GetCellResources(int cell_id,cAvidaContext & ctx) const547 const tArray<double> & cResourceCount::GetCellResources(int cell_id, cAvidaContext& ctx) const
548
549 // Get amount of the resource for a given cell in the grid. If it is a
550 // global resource pass out the entire content of that resource.
551
552 {
553 int num_resources = resource_count.GetSize();
554
555 DoUpdates(ctx);
556
557 for (int i = 0; i < num_resources; i++) {
558 if (geometry[i] == nGeometry::GLOBAL || geometry[i]==nGeometry::PARTIAL) {
559 curr_grid_res_cnt[i] = resource_count[i];
560 } else {
561 curr_grid_res_cnt[i] = spatial_resource_count[i]->GetAmount(cell_id);
562 }
563 }
564 return curr_grid_res_cnt;
565
566 }
567
GetFrozenResources(cAvidaContext & ctx,int cell_id) const568 const tArray<double> & cResourceCount::GetFrozenResources(cAvidaContext& ctx, int cell_id) const
569 // This differs from GetCellResources by leaving out DoUpdates which is
570 // useful inside methods that repeatedly call this before cells can change.
571 {
572 int num_resources = resource_count.GetSize();
573
574 for (int i = 0; i < num_resources; i++) {
575 if (geometry[i] == nGeometry::GLOBAL || geometry[i]==nGeometry::PARTIAL) {
576 curr_grid_res_cnt[i] = resource_count[i];
577 } else {
578 curr_grid_res_cnt[i] = spatial_resource_count[i]->GetAmount(cell_id);
579 }
580 }
581 return curr_grid_res_cnt;
582 }
583
GetResourcesGeometry() const584 const tArray<int> & cResourceCount::GetResourcesGeometry() const
585 {
586 return geometry;
587 }
588
GetSpatialRes(cAvidaContext & ctx)589 const tArray< tArray<double> > & cResourceCount::GetSpatialRes(cAvidaContext& ctx)
590 {
591 const int num_spatial_resources = spatial_resource_count.GetSize();
592 if (num_spatial_resources > 0) {
593 const int num_cells = spatial_resource_count[0]->GetSize();
594 DoUpdates(ctx);
595 for (int i = 0; i < num_spatial_resources; i++) {
596 for (int j = 0; j < num_cells; j++) {
597 curr_spatial_res_cnt[i][j] = spatial_resource_count[i]->GetAmount(j);
598 }
599 }
600 }
601 return curr_spatial_res_cnt;
602 }
603
Modify(cAvidaContext & ctx,const tArray<double> & res_change)604 void cResourceCount::Modify(cAvidaContext& ctx, const tArray<double> & res_change)
605 {
606 assert(resource_count.GetSize() == res_change.GetSize());
607
608 DoUpdates(ctx);
609 for (int i = 0; i < resource_count.GetSize(); i++) {
610 resource_count[i] += res_change[i];
611 assert(resource_count[i] >= 0.0);
612 }
613 }
614
615
Modify(cAvidaContext & ctx,int res_index,double change)616 void cResourceCount::Modify(cAvidaContext& ctx, int res_index, double change)
617 {
618 assert(res_index < resource_count.GetSize());
619
620 DoUpdates(ctx);
621 resource_count[res_index] += change;
622 assert(resource_count[res_index] >= 0.0);
623 }
624
ModifyCell(cAvidaContext & ctx,const tArray<double> & res_change,int cell_id)625 void cResourceCount::ModifyCell(cAvidaContext& ctx, const tArray<double> & res_change, int cell_id)
626 {
627 assert(resource_count.GetSize() == res_change.GetSize());
628
629 DoUpdates(ctx);
630 for (int i = 0; i < resource_count.GetSize(); i++) {
631 if (geometry[i] == nGeometry::GLOBAL || geometry[i]==nGeometry::PARTIAL) {
632 resource_count[i] += res_change[i];
633 assert(resource_count[i] >= 0.0);
634 } else {
635 double temp = spatial_resource_count[i]->Element(cell_id).GetAmount();
636 spatial_resource_count[i]->Rate(cell_id, res_change[i]);
637 /* Ideally the state of the cell's resource should not be set till
638 the end of the update so that all processes (inflow, outflow,
639 diffision, gravity and organism demand) have the same weight. However
640 waiting can cause problems with negative resources so we allow
641 the organism demand to work immediately on the state of the resource */
642
643 spatial_resource_count[i]->State(cell_id);
644 if(spatial_resource_count[i]->Element(cell_id).GetAmount() != temp){
645 spatial_resource_count[i]->SetModified(true);
646 }
647 assert(spatial_resource_count[i]->Element(cell_id).GetAmount() >= 0.0); //APW
648 }
649 }
650 }
651
Get(cAvidaContext & ctx,int res_index) const652 double cResourceCount::Get(cAvidaContext& ctx, int res_index) const
653 {
654 assert(res_index < resource_count.GetSize());
655 DoUpdates(ctx);
656 if (geometry[res_index] == nGeometry::GLOBAL || geometry[res_index]==nGeometry::PARTIAL) {
657 return resource_count[res_index];
658 } //else return spacial resource sum
659 return spatial_resource_count[res_index]->SumAll();
660 }
661
Set(cAvidaContext & ctx,int res_index,double new_level)662 void cResourceCount::Set(cAvidaContext& ctx, int res_index, double new_level)
663 {
664 assert(res_index < resource_count.GetSize());
665 DoUpdates(ctx);
666 if (geometry[res_index] == nGeometry::GLOBAL || geometry[res_index]==nGeometry::PARTIAL) {
667 resource_count[res_index] = new_level;
668 } else {
669 for(int i = 0; i < spatial_resource_count[res_index]->GetSize(); i++) {
670 spatial_resource_count[res_index]->SetCellAmount(i, new_level/spatial_resource_count[res_index]->GetSize());
671 }
672 }
673 }
674
ResizeSpatialGrids(int in_x,int in_y)675 void cResourceCount::ResizeSpatialGrids(int in_x, int in_y)
676 {
677 for (int i = 0; i < resource_count.GetSize(); i++) {
678 spatial_resource_count[i]->ResizeClear(in_x, in_y, geometry[i]);
679 curr_spatial_res_cnt[i].Resize(in_x * in_y);
680 }
681 }
682
GetCurrPeakX(cAvidaContext & ctx,int res_id) const683 int cResourceCount::GetCurrPeakX(cAvidaContext& ctx, int res_id) const
684 {
685 DoUpdates(ctx);
686 return spatial_resource_count[res_id]->GetCurrPeakX();
687 }
688
GetCurrPeakY(cAvidaContext & ctx,int res_id) const689 int cResourceCount::GetCurrPeakY(cAvidaContext& ctx, int res_id) const
690 {
691 DoUpdates(ctx);
692 return spatial_resource_count[res_id]->GetCurrPeakY();
693 }
694
GetFrozenPeakX(cAvidaContext & ctx,int res_id) const695 int cResourceCount::GetFrozenPeakX(cAvidaContext& ctx, int res_id) const
696 {
697 return spatial_resource_count[res_id]->GetCurrPeakX();
698 }
699
GetFrozenPeakY(cAvidaContext & ctx,int res_id) const700 int cResourceCount::GetFrozenPeakY(cAvidaContext& ctx, int res_id) const
701 {
702 return spatial_resource_count[res_id]->GetCurrPeakY();
703 }
704
GetWallCells(int res_id)705 tArray<int>* cResourceCount::GetWallCells(int res_id)
706 {
707 return spatial_resource_count[res_id]->GetWallCells();
708 }
709
GetMinUsedX(int res_id)710 int cResourceCount::GetMinUsedX(int res_id)
711 {
712 return spatial_resource_count[res_id]->GetMinUsedX();
713 }
714
GetMinUsedY(int res_id)715 int cResourceCount::GetMinUsedY(int res_id)
716 {
717 return spatial_resource_count[res_id]->GetMinUsedY();
718 }
719
GetMaxUsedX(int res_id)720 int cResourceCount::GetMaxUsedX(int res_id)
721 {
722 return spatial_resource_count[res_id]->GetMaxUsedX();
723 }
724
GetMaxUsedY(int res_id)725 int cResourceCount::GetMaxUsedY(int res_id)
726 {
727 return spatial_resource_count[res_id]->GetMaxUsedY();
728 }
729
730 ///// Private Methods /////////
DoUpdates(cAvidaContext & ctx,bool global_only) const731 void cResourceCount::DoUpdates(cAvidaContext& ctx, bool global_only) const
732 {
733 assert(update_time >= -EPSILON);
734
735 // Determine how many update steps have progressed
736 int num_steps = (int) (update_time / UPDATE_STEP);
737
738 // Preserve remainder of update_time
739 update_time -= num_steps * UPDATE_STEP;
740
741
742 while (num_steps > PRECALC_DISTANCE) {
743 for (int i = 0; i < resource_count.GetSize(); i++) {
744 if (geometry[i] == nGeometry::GLOBAL || geometry[i]==nGeometry::PARTIAL) {
745 resource_count[i] *= decay_precalc(i, PRECALC_DISTANCE);
746 resource_count[i] += inflow_precalc(i, PRECALC_DISTANCE);
747 }
748 }
749 num_steps -= PRECALC_DISTANCE;
750 }
751
752 for (int i = 0; i < resource_count.GetSize(); i++) {
753 if (geometry[i] == nGeometry::GLOBAL || geometry[i]==nGeometry::PARTIAL) {
754 resource_count[i] *= decay_precalc(i, num_steps);
755 resource_count[i] += inflow_precalc(i, num_steps);
756 }
757 }
758
759 if (global_only) return;
760
761 // If one (or more) complete update has occured update the spatial resources
762 while (m_spatial_update > m_last_updated) {
763 m_last_updated++;
764 for (int i = 0; i < resource_count.GetSize(); i++) {
765 if (geometry[i] != nGeometry::GLOBAL && geometry[i] != nGeometry::PARTIAL) {
766 spatial_resource_count[i]->UpdateCount(ctx);
767 spatial_resource_count[i]->Source(inflow_rate[i]);
768 spatial_resource_count[i]->Sink(decay_rate[i]);
769 if (spatial_resource_count[i]->GetCellListSize() > 0) {
770 spatial_resource_count[i]->CellInflow();
771 spatial_resource_count[i]->CellOutflow();
772 }
773 spatial_resource_count[i]->FlowAll();
774 spatial_resource_count[i]->StateAll();
775 // BDB: resource_count[i] = spatial_resource_count[i]->SumAll();
776 }
777 }
778 }
779 }
780
ReinitializeResources(cAvidaContext & ctx,double additional_resource)781 void cResourceCount::ReinitializeResources(cAvidaContext& ctx, double additional_resource)
782 {
783 for(int i = 0; i < resource_name.GetSize(); i++) {
784 Set(ctx, i, resource_initial[i] + additional_resource); //will cause problem if more than one resource is used. -- why? each resource is stored separately (BDC)
785
786 // Additionally, set any initial values given by the CELL command
787 spatial_resource_count[i]->ResetResourceCounts();
788 if (additional_resource != 0.0) {
789 spatial_resource_count[i]->RateAll(additional_resource);
790 spatial_resource_count[i]->StateAll();
791 }
792
793 } //End going through the resources
794 }
795
796 /*
797 * TODO: This is a duplicate of GetResourceCountID()
798 * Follow the same steps to get rid of it.
799 */
GetResourceByName(cString name) const800 int cResourceCount::GetResourceByName(cString name) const
801 {
802 int result = -1;
803
804 for(int i = 0; i < resource_name.GetSize(); i++)
805 {
806 if(resource_name[i] == name)
807 {
808 result = i;
809 }
810 }
811
812 return result;
813
814 }
815
816
817