1 /*
2 * cDeme.cc
3 * Avida
4 *
5 * Copyright 1999-2011 Michigan State University. All rights reserved.
6 *
7 *
8 * This file is part of Avida.
9 *
10 * Avida is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
12 *
13 * Avida is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License along with Avida.
17 * If not, see <http://www.gnu.org/licenses/>.
18 *
19 */
20
21 #include "cDeme.h"
22
23 #include "cBioGroup.h"
24 #include "cBioGroupManager.h"
25 #include "cClassificationManager.h"
26 #include "cEnvironment.h"
27 #include "cOrganism.h"
28 #include "cPhenotype.h"
29 #include "cPopulation.h"
30 #include "cPopulationCell.h"
31 #include "cResource.h"
32 #include "cStats.h"
33 #include "cWorld.h"
34 #include "cOrgMessagePredicate.h"
35 #include "cOrgMovementPredicate.h"
36 #include "cDemePredicate.h"
37 #include "cReactionResult.h" //@JJB**
38 #include "cTaskState.h" //@JJB**
39
40 #include <cmath>
41
42
43 /*! Constructor
44 */
cDeme()45 cDeme::cDeme()
46 : _id(0)
47 , width(0)
48 , replicateDeme(false)
49 , treatable(false)
50 , cur_birth_count(0)
51 , last_birth_count(0)
52 , cur_org_count(0)
53 , last_org_count(0)
54 , injected_count(0)
55 , birth_count_perslot(0)
56 , _age(0)
57 , generation(0)
58 , total_org_energy(0.0)
59 , time_used(0)
60 , gestation_time(0)
61 , cur_normalized_time_used(0.0)
62 , last_normalized_time_used(0.0)
63 , MSG_sendFailed(0)
64 , MSG_dropped(0)
65 , MSG_SuccessfullySent(0)
66 , energyInjectedIntoOrganisms(0.0)
67 , energyRemainingInDemeAtReplication(0.0)
68 , total_energy_testament(0.0)
69 , eventsTotal(0)
70 , eventsKilled(0)
71 , eventsKilledThisSlot(0)
72 , eventKillAttempts(0)
73 , eventKillAttemptsThisSlot(0)
74 , consecutiveSuccessfulEventPeriods(0)
75 , sleeping_count(0)
76 , avg_founder_generation(0.0)
77 , generations_per_lifetime(0.0)
78 , deme_resource_count(0)
79 , m_germline_genotype_id(0)
80 , points(0)
81 , migrations_out(0)
82 , migrations_in(0)
83 , suicides(0)
84 , m_network(0)
85 , m_num_reproductives(0)
86 , m_input_pointer(0) //@JJB**
87 , m_input_buf(0) //@JJB**
88 , m_output_buf(0) //@JJB**
89 , m_reaction_result(NULL) //@JJB**
90 {
91 }
92
~cDeme()93 cDeme::~cDeme()
94 {
95 if(m_network) delete m_network;
96
97 // Remove Task States @JJB**
98 tArray<cTaskState*> task_states(0);
99 m_task_states.GetValues(task_states);
100 for (int i = 0; i < task_states.GetSize(); i++) delete task_states[i];
101 delete m_reaction_result;
102 }
103
operator =(const cDeme & in_deme)104 cDeme& cDeme::operator=(const cDeme& in_deme) //@JJB**
105 {
106 m_world = in_deme.m_world;
107 _id = in_deme._id;
108 cell_ids = in_deme.cell_ids;
109 width = in_deme.width;
110 replicateDeme = in_deme.replicateDeme;
111 treatable = in_deme.treatable;
112 treatment_ages = in_deme.treatment_ages;
113 cur_birth_count = in_deme.cur_birth_count;
114 last_birth_count = in_deme.last_birth_count;
115 cur_org_count = in_deme.cur_org_count;
116 last_org_count = in_deme.last_org_count;
117 injected_count = in_deme.injected_count;
118 birth_count_perslot = in_deme.birth_count_perslot;
119 _age = in_deme._age;
120 generation = in_deme.generation;
121 total_org_energy = in_deme.total_org_energy;
122 time_used = in_deme.time_used;
123 gestation_time = in_deme.gestation_time;
124 cur_normalized_time_used = in_deme.cur_normalized_time_used;
125 last_normalized_time_used = in_deme.last_normalized_time_used;
126 MSG_sendFailed = in_deme.MSG_sendFailed;
127 MSG_dropped = in_deme.MSG_dropped;
128 MSG_SuccessfullySent = in_deme.MSG_SuccessfullySent;
129 energyInjectedIntoOrganisms = in_deme.energyInjectedIntoOrganisms;
130 energyRemainingInDemeAtReplication = in_deme.energyRemainingInDemeAtReplication;
131 total_energy_testament = in_deme.total_energy_testament;
132 eventsTotal = in_deme.eventsTotal;
133 eventsKilled = in_deme.eventsKilled;
134 eventsKilledThisSlot = in_deme.eventsKilledThisSlot;
135 eventKillAttempts = in_deme.eventKillAttempts;
136 eventKillAttemptsThisSlot = in_deme.eventKillAttemptsThisSlot;
137 consecutiveSuccessfulEventPeriods = in_deme.consecutiveSuccessfulEventPeriods;
138 sleeping_count = in_deme.sleeping_count;
139 energyUsage = in_deme.energyUsage;
140 total_energy_donated = in_deme.total_energy_donated;
141 total_energy_received = in_deme.total_energy_received;
142 total_energy_applied = in_deme.total_energy_applied;
143 cur_task_exe_count = in_deme.cur_task_exe_count;
144 cur_reaction_count = in_deme.cur_reaction_count;
145 last_task_exe_count = in_deme.last_task_exe_count;
146 last_reaction_count = in_deme.last_reaction_count;
147 cur_org_task_count = in_deme.cur_org_task_count;
148 cur_org_task_exe_count = in_deme.cur_org_task_exe_count;
149 cur_org_reaction_count = in_deme.cur_org_reaction_count;
150 last_org_task_count = in_deme.last_org_task_count;
151 last_org_task_exe_count = in_deme.last_org_task_exe_count;
152 last_org_reaction_count = in_deme.last_org_reaction_count;
153 avg_founder_generation = in_deme.avg_founder_generation;
154 generations_per_lifetime = in_deme.generations_per_lifetime;
155 _germline = in_deme._germline;
156 cell_events = in_deme.cell_events;
157 event_slot_end_points = in_deme.event_slot_end_points;
158 m_germline_genotype_id = in_deme.m_germline_genotype_id;
159 m_founder_genotype_ids = in_deme.m_founder_genotype_ids;
160 m_founder_phenotypes = in_deme.m_founder_phenotypes;
161 _current_merit = in_deme._current_merit;
162 _next_merit = in_deme._next_merit;
163 deme_pred_list = in_deme.deme_pred_list;
164 message_pred_list = in_deme.message_pred_list;
165 movement_pred_list = in_deme.movement_pred_list;
166 points = in_deme.points;
167 migrations_out = in_deme.migrations_out;
168 migrations_in = in_deme.migrations_in;
169 suicides = in_deme.suicides;
170 //m_network = in_deme.m_network;
171 m_inputs = in_deme.m_inputs;
172 m_input_buf = in_deme.m_input_buf;
173 m_output_buf = in_deme.m_output_buf;
174 //m_reaction_result = in_deme.m_reaction_result;
175 m_task_count = in_deme.m_task_count;
176 m_reaction_count = in_deme.m_reaction_count;
177 m_last_task_count = in_deme.m_last_task_count;
178 m_cur_reaction_add_reward = in_deme.m_cur_reaction_add_reward;
179 m_cur_bonus = in_deme.m_cur_bonus;
180 m_cur_merit = in_deme.m_cur_merit;
181 m_total_res_consumed = in_deme.m_total_res_consumed;
182 m_switch_penalties = in_deme.m_switch_penalties;
183 m_shannon_matrix = in_deme.m_shannon_matrix;
184 m_num_active = in_deme.m_num_active;
185 m_num_reproductives = in_deme.m_num_reproductives;
186
187 tList<cTaskState*> hash_values;
188 tList<void*> hash_keys;
189 in_deme.m_task_states.AsLists(hash_keys, hash_values);
190 tListIterator<cTaskState*> vit(hash_values);
191 tListIterator<void*> kit(hash_keys);
192 while(vit.Next() && kit.Next()) {
193 cTaskState* new_ts = new cTaskState(**(vit.Get()));
194 m_task_states.Set(*(kit.Get()), new_ts);
195 }
196
197 return *this;
198 }
199
Setup(int id,const tArray<int> & in_cells,int in_width,cWorld * world)200 void cDeme::Setup(int id, const tArray<int> & in_cells, int in_width, cWorld* world)
201 {
202 _id = id;
203 cell_ids = in_cells;
204 cur_birth_count = 0;
205 last_birth_count = 0;
206 cur_org_count = 0;
207 last_org_count = 0;
208 birth_count_perslot = 0;
209 m_world = world;
210
211 replicateDeme = false;
212
213 _current_merit = 1.0;
214 _next_merit = 1.0;
215
216 const int num_tasks = m_world->GetEnvironment().GetNumTasks();
217 const int num_reactions = m_world->GetEnvironment().GetNumReactions();
218
219 cur_task_exe_count.Resize(num_tasks);
220 cur_task_exe_count.SetAll(0);
221 cur_reaction_count.ResizeClear(num_reactions);
222 cur_reaction_count.SetAll(0);
223 last_task_exe_count.ResizeClear(num_tasks);
224 last_task_exe_count.SetAll(0);
225 last_reaction_count.ResizeClear(num_reactions);
226 last_reaction_count.SetAll(0);
227
228 cur_org_task_count.Resize(num_tasks);
229 cur_org_task_count.SetAll(0);
230 cur_org_task_exe_count.Resize(num_tasks);
231 cur_org_task_exe_count.SetAll(0);
232 cur_org_reaction_count.ResizeClear(num_reactions);
233 cur_org_reaction_count.SetAll(0);
234 last_org_task_count.ResizeClear(num_tasks);
235 last_org_task_count.SetAll(0);
236 last_org_task_exe_count.ResizeClear(num_tasks);
237 last_org_task_exe_count.SetAll(0);
238 last_org_reaction_count.ResizeClear(num_reactions);
239 last_org_reaction_count.SetAll(0);
240 m_total_res_consumed = 0;
241 m_switch_penalties = 0;
242 m_num_active = 0;
243 m_num_reproductives = 0;
244
245 m_input_buf = tBuffer<int>(m_world->GetEnvironment().GetInputSize()); //@JJB**
246 m_output_buf = tBuffer<int>(m_world->GetEnvironment().GetOutputSize()); //@JJB**
247 m_task_count.ResizeClear(num_tasks); //@JJB**
248 m_task_count.SetAll(0);
249 m_reaction_count.ResizeClear(m_world->GetEnvironment().GetReactionLib().GetSize()); //@JJB**
250 m_reaction_count.SetAll(0);
251 m_last_task_count.ResizeClear(num_tasks); //@JJB**
252 m_last_task_count.SetAll(0);
253 m_cur_reaction_add_reward.ResizeClear(m_world->GetEnvironment().GetReactionLib().GetSize()); //@JJB**
254 m_cur_reaction_add_reward.SetAll(0.0);
255 m_cur_bonus = m_world->GetConfig().DEFAULT_BONUS.Get(); //@JJB**
256 if (m_world->GetConfig().BASE_MERIT_METHOD.Get() == BASE_MERIT_CONST) {
257 m_cur_merit = m_world->GetConfig().BASE_CONST_MERIT.Get();
258 } else {
259 m_cur_merit = 1.0;
260 }
261
262 total_energy_donated = 0.0;
263 total_energy_received = 0.0;
264 total_energy_applied = 0.0;
265
266 // If width is negative, set it to the full number of cells.
267 width = in_width;
268 if (width < 1) width = cell_ids.GetSize();
269
270 // drain spacial energy resources and place energy in cells
271 }
272
273
GetCellID(int x,int y) const274 int cDeme::GetCellID(int x, int y) const
275 {
276 assert(x >= 0 && x < GetWidth());
277 assert(y >= 0 && y < GetHeight());
278
279 const int pos = y * width + x;
280 return cell_ids[pos];
281 }
282
283
284 /*! Note that for this method to work, we blatantly assume that IDs are in
285 monotonically increasing order!! */
GetCellPosition(int cellid) const286 std::pair<int, int> cDeme::GetCellPosition(int cellid) const
287 {
288 assert(cell_ids.GetSize()>0);
289 assert(GetWidth() > 0);
290 // cellid -= cell_ids[0];
291 // return std::make_pair(cellid % GetWidth(), cellid / GetWidth());
292 return std::make_pair(cellid % GetWidth(), ( cellid % cell_ids.GetSize() ) / GetWidth());
293 }
294
GetCell(int pos) const295 cPopulationCell& cDeme::GetCell(int pos) const
296 {
297 return m_world->GetPopulation().GetCell(cell_ids[pos]);
298 }
299
GetCell(int x,int y) const300 cPopulationCell& cDeme::GetCell(int x, int y) const
301 {
302 return m_world->GetPopulation().GetCell(GetCellID(x,y));
303 }
304
305
GetOrganism(int pos) const306 cOrganism* cDeme::GetOrganism(int pos) const
307 {
308 return GetCell(pos).GetOrganism();
309 }
310
GetGenotypeIDs()311 std::vector<int> cDeme::GetGenotypeIDs()
312 {
313 std::vector<int> genotype_ids;
314 for (int i = 0; i < GetSize(); i++) {
315 cPopulationCell& cell = GetCell(i);
316 if (cell.IsOccupied()) genotype_ids.push_back(cell.GetOrganism()->GetBioGroup("genotype")->GetID());
317 }
318
319 //assert(genotype_ids.size()>0); // How did we get to replication otherwise?
320 //@JEB some germline methods can result in empty source demes if they didn't produce a germ)
321
322 return genotype_ids;
323 }
324
325
326
GetNumOrgsWithOpinion() const327 int cDeme::GetNumOrgsWithOpinion() const
328 {
329 int demeSize = GetSize();
330 int count = 0;
331
332 for (int pos = 0; pos < demeSize; ++pos) {
333 cPopulationCell& cell = GetCell(pos);
334 if (cell.IsOccupied() && cell.GetOrganism()->HasOpinion())
335 ++count;
336 }
337
338 return count;
339 }
340
ProcessPreUpdate()341 void cDeme::ProcessPreUpdate()
342 {
343 deme_resource_count.SetSpatialUpdate(m_world->GetStats().GetUpdate());
344 }
345
ProcessUpdate(cAvidaContext & ctx)346 void cDeme::ProcessUpdate(cAvidaContext& ctx)
347 {
348 // test deme predicate
349 for (int i = 0; i < deme_pred_list.Size(); i++) {
350 if (deme_pred_list[i]->GetName() == "cDemeResourceThreshold") {
351 (*deme_pred_list[i])(ctx, &deme_resource_count);
352 }
353 }
354
355 energyUsage.Clear();
356
357 if(IsEmpty()) { // deme is not processed if no organisms are present
358 total_energy_testament = 0.0;
359 return;
360 }
361
362 if (m_world->GetConfig().ENERGY_ENABLED.Get()) {
363 for(int i = 0; i < GetSize(); i++) {
364 cPopulationCell& cell = GetCell(i);
365 if(cell.IsOccupied()) {
366 energyUsage.Add(cell.GetOrganism()->GetPhenotype().GetEnergyUsageRatio());
367 }
368 }
369 }
370
371 for(int i = 0; i < cell_events.Size(); i++) {
372 cDemeCellEvent& event = cell_events[i];
373
374 if(event.IsActive() && event.GetDelay() < _age && _age <= event.GetDelay()+event.GetDuration()) {
375 //remove energy from cells (should be done with outflow, but this will work for now)
376 int eventCell = event.GetNextEventCellID();
377 cResource* res = m_world->GetEnvironment().GetResourceLib().GetResource("CELL_ENERGY");
378
379 while(eventCell != -1) {
380 cPopulationCell& cell = m_world->GetPopulation().GetCell(GetCellID(eventCell));
381 if(event.GetEventID() == cell.GetCellData()){
382 if(cell.IsOccupied()) {
383 // remove energy from organism
384 cPhenotype& orgPhenotype = cell.GetOrganism()->GetPhenotype();
385 orgPhenotype.ReduceEnergy(orgPhenotype.GetStoredEnergy()*m_world->GetConfig().ATTACK_DECAY_RATE.Get());
386 }
387 //remove energy from cell... organism might not takeup all of a cell's energy
388 tArray<double> cell_resources = deme_resource_count.GetCellResources(eventCell, ctx); // uses global cell_id; is this a problem
389 cell_resources[res->GetID()] *= m_world->GetConfig().ATTACK_DECAY_RATE.Get();
390 deme_resource_count.ModifyCell(ctx, cell_resources, eventCell);
391 }
392 eventCell = event.GetNextEventCellID();
393 }
394 }
395
396
397
398 if(!event.IsActive() && event.GetDelay() == _age) {
399 eventsTotal++;
400 event.ActivateEvent(); //start event
401 int eventCell = event.GetNextEventCellID();
402 while(eventCell != -1) {
403 // place event ID in cells' data
404 if(event.IsDecayed()) {
405 m_world->GetPopulation().GetCell(GetCellID(eventCell)).SetCellData(event.GetEventIDDecay(GetCellPosition(eventCell)));
406 } else {
407 m_world->GetPopulation().GetCell(GetCellID(eventCell)).SetCellData(event.GetEventID());
408 }
409
410 // record activation of each cell in stats
411 std::pair<int, int> pos = GetCellPosition(eventCell);
412 m_world->GetStats().IncEventCount(pos.first, pos.second);
413
414
415 //TODO // increase outflow of energy from these cells if not event currently present
416
417
418 eventCell = event.GetNextEventCellID();
419 }
420 } else if (event.IsActive() && event.GetDelay()+event.GetDuration() == _age) {
421 int eventCell = event.GetNextEventCellID();
422 while (eventCell != -1) {
423 if (event.GetEventID() == m_world->GetPopulation().GetCell(GetCellID(eventCell)).GetCellData()) { // eventID == CellData
424 //set cell data to 0
425 m_world->GetPopulation().GetCell(GetCellID(eventCell)).SetCellData(0);
426
427 // TODO // remove energy outflow from these cells
428
429 }
430 eventCell = event.GetNextEventCellID();
431 }
432 event.DeactivateEvent(); //event over
433 }
434 }
435
436 for (vector<pair<int, int> >::iterator iter = event_slot_end_points.begin();
437 iter < event_slot_end_points.end();
438 iter++) {
439 if (_age == (*iter).first) {
440 // at end point
441 if (GetEventsKilledThisSlot() >=
442 m_world->GetConfig().DEMES_MIM_EVENTS_KILLED_RATIO.Get() * (*iter).second) {
443 consecutiveSuccessfulEventPeriods++;
444 } else {
445 consecutiveSuccessfulEventPeriods = 0;
446 }
447
448 // update stats.flow_rate_tuples
449 std::map<int, flow_rate_tuple>& flowRateTuples = m_world->GetStats().FlowRateTuples();
450
451 flowRateTuples[(*iter).second].orgCount.Add(GetOrgCount());
452 flowRateTuples[(*iter).second].eventsKilled.Add(GetEventsKilledThisSlot());
453 flowRateTuples[(*iter).second].attemptsToKillEvents.Add(GetEventKillAttemptsThisSlot());
454 flowRateTuples[(*iter).second].AvgEnergyUsageRatio.Add(energyUsage.Average());
455 flowRateTuples[(*iter).second].totalBirths.Add(birth_count_perslot);
456 flowRateTuples[(*iter).second].currentSleeping.Add(sleeping_count);
457 birth_count_perslot = 0;
458 eventsKilledThisSlot = 0;
459 eventKillAttemptsThisSlot = 0;
460 break;
461 }
462 }
463 ++_age;
464
465 // Let the network process the update too, if we have one.
466 if (IsNetworkInitialized()) m_network->ProcessUpdate();
467 }
468
469
470 /*! Called when an organism living in a cell in this deme is about to be killed.
471
472 This method is called from cPopulation::KillOrganism().
473 */
474
OrganismDeath(cPopulationCell & cell)475 void cDeme::OrganismDeath(cPopulationCell& cell)
476 {
477 // Clean up this deme's network, if we have one.
478 if (IsNetworkInitialized()) m_network->OrganismDeath(cell);
479
480 // Add information about the organisms tasks to the deme so that
481 // we can use it to compute shannon mutual information and measure
482 // division of labor.
483 if (m_world->GetConfig().DEMES_TRACK_SHANNON_INFO.Get()) UpdateShannon(cell);
484 }
485
486
487
Reset(cAvidaContext & ctx,bool resetResources,double deme_energy)488 void cDeme::Reset(cAvidaContext& ctx, bool resetResources, double deme_energy)
489 {
490 double additional_resource = 0.0;
491 // Handle energy model
492 if (m_world->GetConfig().ENERGY_ENABLED.Get())
493 {
494 total_energy_testament = 0.0;
495
496 if(m_world->GetConfig().ENERGY_PASSED_ON_DEME_REPLICATION_METHOD.Get() == 1) {
497 // split deme energy evenly between cell in deme
498 additional_resource = deme_energy; // spacial resource handles resource division
499 }
500 }
501
502 // Handle stat and resource resets
503 _age = 0;
504 time_used = 0;
505 cur_birth_count = 0;
506 cur_normalized_time_used = 0;
507 injected_count = 0;
508 birth_count_perslot = 0;
509 eventsTotal = 0;
510 eventsKilled = 0;
511 eventsKilledThisSlot = 0;
512 eventKillAttempts = 0;
513 eventKillAttemptsThisSlot = 0;
514 sleeping_count = 0;
515 MSG_sendFailed = 0;
516 MSG_dropped = 0;
517 MSG_SuccessfullySent = 0;
518 m_total_res_consumed = 0;
519 m_switch_penalties = 0;
520 m_num_active = 0;
521 m_shannon_matrix.clear();
522 m_num_reproductives = 0;
523
524 m_input_pointer = 0; //@JJB**
525 m_input_buf.Clear(); //@JJB**
526 m_output_buf.Clear(); //@JJB**
527 m_last_task_count = m_task_count; //@JJB**
528 m_task_count.SetAll(0); //@JJB**
529 m_reaction_count.SetAll(0); //@JJB**
530 m_cur_reaction_add_reward.SetAll(0.0); //@JJB**
531 m_cur_bonus = m_world->GetConfig().DEFAULT_BONUS.Get(); //@JJB**
532 if (m_world->GetConfig().BASE_MERIT_METHOD.Get() == BASE_MERIT_CONST) {
533 m_cur_merit = m_world->GetConfig().BASE_CONST_MERIT.Get();
534 } else {
535 m_cur_merit = 1.0;
536 }
537
538 // Reset Task States
539 tArray<cTaskState*> task_states(0);
540 m_task_states.GetValues(task_states);
541 for (int i = 0; i < task_states.GetSize(); i++) delete task_states[i];
542 m_task_states.ClearAll();
543
544 consecutiveSuccessfulEventPeriods = 0;
545
546 replicateDeme = false;
547
548 total_energy_donated = 0.0;
549 total_energy_received = 0.0;
550 total_energy_applied = 0.0;
551
552 cur_task_exe_count.SetAll(0);
553 cur_reaction_count.SetAll(0);
554
555 //reset remaining deme predicates
556 for (int i = 0; i < deme_pred_list.Size(); i++) {
557 deme_pred_list[i]->Reset();
558 }
559 //reset remaining message predicates
560 for (int i = 0; i < message_pred_list.Size(); i++) {
561 message_pred_list[i]->Reset();
562 }
563 //reset remaining message predicates
564 for (int i = 0; i < movement_pred_list.Size(); i++) {
565 movement_pred_list[i]->Reset();
566 }
567
568 if (resetResources) {
569 deme_resource_count.ReinitializeResources(ctx, additional_resource);
570 }
571
572 // Instead of polluting cDemeNetwork with Resets, we're just going to delete it,
573 // and go ahead and rely on lazy initialization to fill this back in.
574 if (m_network) {
575 delete m_network;
576 m_network = 0;
577 }
578 }
579
580
DivideReset(cAvidaContext & ctx,cDeme & parent_deme,bool resetResources,double deme_energy)581 void cDeme::DivideReset(cAvidaContext& ctx, cDeme& parent_deme, bool resetResources, double deme_energy)
582 {
583 // the parent might be us, so save this value...
584 double old_avg_founder_generation = parent_deme.GetAvgFounderGeneration();
585
586 // update our average founder generation
587 cDoubleSum gen;
588 for (int i=0; i< m_founder_phenotypes.GetSize(); i++) {
589 gen.Add( m_founder_phenotypes[i].GetGeneration() );
590 }
591 avg_founder_generation = gen.Average();
592
593 // update our generations per lifetime based on current founders and parents generation
594 generations_per_lifetime = avg_founder_generation - old_avg_founder_generation;
595
596 //Save statistics according to parent before reset.
597 generation = parent_deme.GetGeneration() + 1;
598 gestation_time = parent_deme.GetTimeUsed();
599 last_normalized_time_used = parent_deme.GetNormalizedTimeUsed();
600
601 last_task_exe_count = parent_deme.GetCurTaskExeCount();
602 last_reaction_count = parent_deme.GetCurReactionCount();
603
604 last_org_task_count = parent_deme.GetCurOrgTaskCount();
605 last_org_task_exe_count = parent_deme.GetCurOrgTaskExeCount();
606 last_org_reaction_count = parent_deme.GetCurOrgReactionCount();
607
608 last_org_count = parent_deme.GetLastOrgCount(); // Org count was updated upon KillAll()....
609 last_birth_count = parent_deme.GetBirthCount();
610
611 Reset(ctx, resetResources, deme_energy);
612 }
613
614
615 // Given the input deme founders and original ones,
616 // calculate how many generations this deme went through to divide.
UpdateGenerationsPerLifetime(double old_avg_founder_generation,tArray<cPhenotype> & new_founder_phenotypes)617 void cDeme::UpdateGenerationsPerLifetime(double old_avg_founder_generation, tArray<cPhenotype>& new_founder_phenotypes)
618 {
619 cDoubleSum gen;
620 for (int i=0; i< new_founder_phenotypes.GetSize(); i++) {
621 gen.Add( new_founder_phenotypes[i].GetGeneration() );
622 }
623 double new_avg_founder_generation = gen.Average();
624 generations_per_lifetime = new_avg_founder_generation - old_avg_founder_generation;
625 }
626
627 /*! Check every cell in this deme for a living organism. If found, kill it. */
KillAll(cAvidaContext & ctx)628 void cDeme::KillAll(cAvidaContext& ctx)
629 {
630 last_org_count = GetOrgCount();
631 for (int i=0; i<GetSize(); ++i) {
632 cPopulationCell& cell = m_world->GetPopulation().GetCell(cell_ids[i]);
633 if(cell.IsOccupied()) {
634 m_world->GetPopulation().KillOrganism(cell, ctx);
635 }
636 }
637
638 // HACK: organism are killed after DivideReset is called.
639 // need clear a deme before it is reset.
640 sleeping_count = 0;
641 }
642
UpdateStats()643 void cDeme::UpdateStats()
644 {
645 //save stats about what tasks our orgs were doing
646 //usually called before KillAll
647
648 for (int j = 0; j < cur_org_task_count.GetSize(); j++) {
649 int count = 0;
650 for (int k=0; k<GetSize(); k++) {
651 int cellid = GetCellID(k);
652 if(m_world->GetPopulation().GetCell(cellid).IsOccupied()) {
653 count += (m_world->GetPopulation().GetCell(cellid).GetOrganism()->GetPhenotype().GetLastTaskCount()[j] > 0);
654 }
655 cur_org_task_count[j] = count;
656 }
657 }
658
659 for(int j = 0; j < cur_org_task_exe_count.GetSize(); j++) {
660 int count = 0;
661 for(int k=0; k<GetSize(); k++) {
662 int cellid = GetCellID(k);
663 if(m_world->GetPopulation().GetCell(cellid).IsOccupied()) {
664 count += m_world->GetPopulation().GetCell(cellid).GetOrganism()->GetPhenotype().GetLastTaskCount()[j];
665 }
666 cur_org_task_exe_count[j] = count;
667 }
668 }
669
670 for(int j = 0; j < cur_org_reaction_count.GetSize(); j++) {
671 int count = 0;
672 for(int k=0; k<GetSize(); k++) {
673 int cellid = GetCellID(k);
674 if(m_world->GetPopulation().GetCell(cellid).IsOccupied()) {
675 count += m_world->GetPopulation().GetCell(cellid).GetOrganism()->GetPhenotype().GetLastReactionCount()[j];
676 }
677 cur_org_reaction_count[j] = count;
678 }
679 }
680 }
681
682
683 /*! Replacing this deme's germline has the effect of changing the deme's lineage.
684 There's still some work to do here; the lineage labels of the Genomes in the germline
685 are all messed up.
686 */
ReplaceGermline(const cGermline & germline)687 void cDeme::ReplaceGermline(const cGermline& germline) {
688 _germline = germline;
689 }
690
691
692 /*! If this method is called, this is the "parent" deme. As with individuals,
693 we need to rotate the heritable merit to the current merit.
694 */
UpdateDemeMerit()695 void cDeme::UpdateDemeMerit() {
696 assert(_next_merit.GetDouble()>=1.0);
697 _current_merit = _next_merit;
698 _next_merit = 1.0;
699 }
700
701
702 /*! Update this deme's merit from the given source.
703 */
UpdateDemeMerit(cDeme & source)704 void cDeme::UpdateDemeMerit(cDeme& source) {
705 _current_merit = source.GetHeritableDemeMerit();
706 _next_merit = 1.0;
707 assert(_current_merit.GetDouble()>=1.0);
708 }
709
710
ModifyDemeResCount(cAvidaContext & ctx,const tArray<double> & res_change,const int absolute_cell_id)711 void cDeme::ModifyDemeResCount(cAvidaContext& ctx, const tArray<double>& res_change, const int absolute_cell_id) {
712 // find relative cell_id in deme resource count
713 const int relative_cell_id = GetRelativeCellID(absolute_cell_id);
714 deme_resource_count.ModifyCell(ctx, res_change, relative_cell_id);
715 }
716
SetupDemeRes(int id,cResource * res,int verbosity,cWorld * world)717 void cDeme::SetupDemeRes(int id, cResource * res, int verbosity, cWorld* world) {
718 const double decay = 1.0 - res->GetOutflow();
719 //addjust the resources cell list pointer here if we want CELL env. commands to be replicated in each deme
720
721 int* temp = &id;
722
723 deme_resource_count.Setup(world, *temp, res->GetName(), res->GetInitial(),
724 res->GetInflow(), decay,
725 res->GetGeometry(), res->GetXDiffuse(),
726 res->GetXGravity(), res->GetYDiffuse(),
727 res->GetYGravity(), res->GetInflowX1(),
728 res->GetInflowX2(), res->GetInflowY1(),
729 res->GetInflowY2(), res->GetOutflowX1(),
730 res->GetOutflowX2(), res->GetOutflowY1(),
731 res->GetOutflowY2(), res->GetCellListPtr(),
732 res->GetCellIdListPtr(), verbosity,
733 res->GetDynamicResource(), res->GetPeaks(),
734 res->GetMinHeight(), res->GetMinRadius(), res->GetRadiusRange(),
735 res->GetAh(), res->GetAr(),
736 res->GetAcx(), res->GetAcy(),
737 res->GetHStepscale(), res->GetRStepscale(),
738 res->GetCStepscaleX(), res->GetCStepscaleY(),
739 res->GetHStep(), res->GetRStep(),
740 res->GetCStepX(), res->GetCStepY(),
741 res->GetUpdateDynamic(), res->GetPeakX(), res->GetPeakY(),
742 res->GetHeight(), res->GetSpread(), res->GetPlateau(), res->GetDecay(),
743 res->GetMaxX(), res->GetMaxY(), res->GetMinX(), res->GetMinY(),
744 res->GetAscaler(), res->GetUpdateStep(),
745 res->GetHalo(), res->GetHaloInnerRadius(), res->GetHaloWidth(),
746 res->GetHaloAnchorX(), res->GetHaloAnchorY(), res->GetMoveSpeed(),
747 res->GetPlateauInflow(), res->GetPlateauOutflow(), res->GetConeInflow(), res->GetConeOutflow(),
748 res->GetGradientInflow(), res->GetIsPlateauCommon(), res->GetFloor(), res->GetHabitat(),
749 res->GetMinSize(), res->GetMaxSize(), res->GetConfig(), res->GetCount(), res->GetResistance(),
750 res->GetInitialPlatVal(), res->GetThreshold(), res->GetRefuge(), res->GetGradient()
751 );
752
753 if(res->GetEnergyResource()) {
754 energy_res_ids.Push(id);
755 }
756 }
757
GetCellEnergy(int absolute_cell_id,cAvidaContext & ctx) const758 double cDeme::GetCellEnergy(int absolute_cell_id, cAvidaContext& ctx) const
759 {
760 assert(cell_ids[0] <= absolute_cell_id);
761 assert(absolute_cell_id <= cell_ids[cell_ids.GetSize()-1]);
762
763 double total_energy = 0.0;
764 int relative_cell_id = GetRelativeCellID(absolute_cell_id);
765 tArray<double> cell_resources = deme_resource_count.GetCellResources(relative_cell_id, ctx);
766
767 // sum all energy resources
768 for (int i = 0; i < energy_res_ids.GetSize(); i++) {
769 if (cell_resources[energy_res_ids[i]] > 0.0) {
770 total_energy += cell_resources[energy_res_ids[i]];
771 cell_resources[energy_res_ids[i]] *= -1.0;
772 }
773 }
774
775 return total_energy;
776 }
777
GetAndClearCellEnergy(int absolute_cell_id,cAvidaContext & ctx)778 double cDeme::GetAndClearCellEnergy(int absolute_cell_id, cAvidaContext& ctx)
779 {
780 assert(cell_ids[0] <= absolute_cell_id);
781 assert(absolute_cell_id <= cell_ids[cell_ids.GetSize()-1]);
782
783 double total_energy = 0.0;
784 int relative_cell_id = GetRelativeCellID(absolute_cell_id);
785 tArray<double> cell_resources = deme_resource_count.GetCellResources(relative_cell_id, ctx);
786
787 // sum all energy resources
788 for (int i = 0; i < energy_res_ids.GetSize(); i++) {
789 if (cell_resources[energy_res_ids[i]] > 0.0) {
790 total_energy += cell_resources[energy_res_ids[i]];
791 cell_resources[energy_res_ids[i]] *= -1.0;
792 }
793 }
794
795 // set energy resources to zero
796 deme_resource_count.ModifyCell(ctx, cell_resources, relative_cell_id);
797
798 return total_energy;
799 }
800
GiveBackCellEnergy(int absolute_cell_id,double value,cAvidaContext & ctx)801 void cDeme::GiveBackCellEnergy(int absolute_cell_id, double value, cAvidaContext& ctx)
802 {
803 assert(cell_ids[0] <= absolute_cell_id);
804 assert(absolute_cell_id <= cell_ids[cell_ids.GetSize()-1]);
805
806 int relative_cell_id = GetRelativeCellID(absolute_cell_id);
807 tArray<double> cell_resources = deme_resource_count.GetCellResources(relative_cell_id, ctx);
808
809 double amount_per_resource = value / energy_res_ids.GetSize();
810
811 // put back energy resources evenly
812 for(int i = 0; i < energy_res_ids.GetSize(); i++) {
813 cell_resources[energy_res_ids[i]] += amount_per_resource;
814 }
815 deme_resource_count.ModifyCell(ctx, cell_resources, relative_cell_id);
816 }
817
SetCellEvent(int x1,int y1,int x2,int y2,int delay,int duration,bool static_position,int total_events)818 void cDeme::SetCellEvent(int x1, int y1, int x2, int y2,
819 int delay, int duration, bool static_position, int total_events)
820 {
821 for (int i = 0; i < total_events; i++) {
822 cDemeCellEvent demeEvent = cDemeCellEvent(x1, y1, x2, y2, delay, duration, width, GetHeight(), static_position, this, m_world);
823 cell_events.Add(demeEvent);
824 }
825 }
826
827 /*void cDeme::SetCellEventGradient(int x1, int y1, int x2, int y2, int delay, int duration, bool static_pos, int time_to_live) {
828 cDemeCellEvent demeEvent = cDemeCellEvent(x1, y1, x2, y2, delay, duration, width, GetHeight(), static_pos, time_to_live, this);
829 demeEvent.DecayEventIDFromCenter();
830 cell_events.Add(demeEvent);
831 }*/
832
GetNumEvents()833 int cDeme::GetNumEvents()
834 {
835 return cell_events.Size();
836 }
837
SetCellEventSlots(int x1,int y1,int x2,int y2,int delay,int duration,bool static_position,int m_total_slots,int m_total_events_per_slot_max,int m_total_events_per_slot_min,int m_tolal_event_flow_levels)838 void cDeme::SetCellEventSlots(int x1, int y1, int x2, int y2, int delay, int duration,
839 bool static_position, int m_total_slots, int m_total_events_per_slot_max,
840 int m_total_events_per_slot_min, int m_tolal_event_flow_levels) {
841 assert(cell_events.Size() == 0); // not designed to be used with other cell events
842 assert(m_world->GetConfig().DEMES_MAX_AGE.Get() >= m_total_slots);
843
844 int flow_level_increment = (m_total_events_per_slot_max - m_total_events_per_slot_min) / (m_tolal_event_flow_levels-1);
845 int slot_length = m_world->GetConfig().DEMES_MAX_AGE.Get() / m_total_slots;
846
847 // setup stats tuples
848
849 for (int i = 0; i < m_total_slots; i++) {
850 int slot_flow_level = flow_level_increment * m_world->GetRandom().GetInt(m_tolal_event_flow_levels) + m_total_events_per_slot_min; // number of event during this slot
851 int slot_delay = i * slot_length;
852 event_slot_end_points.push_back(make_pair(slot_delay+slot_length, slot_flow_level)); // last slot is never reached it is == to MAX_AGE
853
854 for (int k = 0; k < slot_flow_level; k++) {
855 cDemeCellEvent demeEvent = cDemeCellEvent(x1, y1, x2, y2, delay, duration, width, GetHeight(), static_position, this, m_world);
856 demeEvent.ConfineToTimeSlot(slot_delay, slot_delay+slot_length);
857 cell_events.Add(demeEvent);
858 }
859 }
860
861 // setup stats.flow_rate_tuples
862 std::map<int, flow_rate_tuple>& flowRateTuples = m_world->GetStats().FlowRateTuples();
863
864 for (int i = m_total_events_per_slot_min; i <= m_total_events_per_slot_max; i+=flow_level_increment) {
865 flowRateTuples[i].orgCount.Clear();
866 flowRateTuples[i].eventsKilled.Clear();
867 flowRateTuples[i].attemptsToKillEvents.Clear();
868 flowRateTuples[i].AvgEnergyUsageRatio.Clear();
869 flowRateTuples[i].totalBirths.Clear();
870 flowRateTuples[i].currentSleeping.Clear();
871 }
872 }
873
KillCellEvent(const int eventID)874 bool cDeme::KillCellEvent(const int eventID)
875 {
876 eventKillAttemptsThisSlot++;
877
878 if (eventID <= 0) return false;
879 for( int i = 0; i < cell_events.Size(); i++) {
880 cDemeCellEvent& event = cell_events[i];
881 if(event.IsActive() && event.GetEventID() == eventID) {
882 // remove event ID from all cells
883 int eventCell = event.GetNextEventCellID();
884 while(eventCell != -1) {
885 if(event.GetEventID() == m_world->GetPopulation().GetCell(GetCellID(eventCell)).GetCellData()) { // eventID == CellData
886 //set cell data to 0
887 m_world->GetPopulation().GetCell(GetCellID(eventCell)).SetCellData(0);
888
889 // TODO // remove energy outflow from these cells
890
891 }
892 eventCell = event.GetNextEventCellID();
893 }
894 event.DeactivateEvent(); //event over
895 eventsKilled++;
896 eventsKilledThisSlot++;
897 eventKillAttempts++;
898 return true;
899 }
900 }
901 return false;
902 }
903
CalculateTotalEnergy(cAvidaContext & ctx) const904 double cDeme::CalculateTotalEnergy(cAvidaContext& ctx) const
905 {
906 assert(m_world->GetConfig().ENERGY_ENABLED.Get());
907
908 double energy_sum = 0.0;
909 for (int i=0; i<GetSize(); i++) {
910 int cellid = GetCellID(i);
911 cPopulationCell& cell = m_world->GetPopulation().GetCell(cellid);
912 if(cell.IsOccupied()) {
913 cOrganism* organism = cell.GetOrganism();
914 cPhenotype& phenotype = organism->GetPhenotype();
915 energy_sum += phenotype.GetStoredEnergy();
916 } else {
917 double energy_in_cell = GetCellEnergy(cellid, ctx);
918 energy_sum += energy_in_cell * m_world->GetConfig().FRAC_ENERGY_TRANSFER.Get();
919 }
920 }
921 return energy_sum;
922 }
923
CalculateTotalInitialEnergyResources() const924 double cDeme::CalculateTotalInitialEnergyResources() const
925 {
926 assert(m_world->GetConfig().ENERGY_ENABLED.Get());
927
928 double energy_sum = 0.0;
929 for(int i = 0; i < energy_res_ids.GetSize(); i++) {
930 energy_sum += deme_resource_count.GetInitialResourceValue(i);
931 }
932 return energy_sum;
933 }
934
935
936 // --- Founder list management --- //
937
AddFounder(cBioGroup * bg,cPhenotype * _in_phenotype)938 void cDeme::AddFounder(cBioGroup* bg, cPhenotype * _in_phenotype)
939 {
940 // save genotype id
941 m_founder_genotype_ids.Push( bg->GetID() );
942 cPhenotype phenotype;
943 if (_in_phenotype) phenotype = *_in_phenotype;
944 m_founder_phenotypes.Push( phenotype );
945
946 bg->AddPassiveReference();
947 }
948
ClearFounders()949 void cDeme::ClearFounders()
950 {
951 // check for unused genotypes, now that we're done with these
952 for (int i=0; i<m_founder_genotype_ids.GetSize(); i++) {
953
954 cBioGroup* bg = m_world->GetClassificationManager().GetBioGroupManager("genotype")->GetBioGroup(m_founder_genotype_ids[i]);
955 assert(bg);
956 bg->RemovePassiveReference();
957 }
958
959 // empty our list
960 m_founder_genotype_ids.ResizeClear(0);
961 m_founder_phenotypes.ResizeClear(0);
962 }
963
ReplaceGermline(cBioGroup * bg)964 void cDeme::ReplaceGermline(cBioGroup* bg)
965 {
966 // same genotype, no changes
967 if (m_germline_genotype_id == bg->GetID()) return;
968
969 // first, save and put a hold on new germline genotype
970 int prev_germline_genotype_id = m_germline_genotype_id;
971 m_germline_genotype_id = bg->GetID();
972 bg->AddPassiveReference();
973
974
975 // next, if we previously were saving a germline genotype, free it
976 cBioGroup* pbg = m_world->GetClassificationManager().GetBioGroupManager("genotype")->GetBioGroup(prev_germline_genotype_id);
977 if (pbg) pbg->RemovePassiveReference();
978 }
979
DemePredSatisfiedPreviously()980 bool cDeme::DemePredSatisfiedPreviously()
981 {
982 for(int i = 0; i < deme_pred_list.Size(); i++) {
983 if(deme_pred_list[i]->PreviouslySatisfied()) {
984 deme_pred_list[i]->UpdateStats(m_world->GetStats());
985 return true;
986 }
987 }
988 return false;
989 }
990
MsgPredSatisfiedPreviously()991 bool cDeme::MsgPredSatisfiedPreviously()
992 {
993 for(int i = 0; i < message_pred_list.Size(); i++) {
994 if(message_pred_list[i]->PreviouslySatisfied()) {
995 message_pred_list[i]->UpdateStats(m_world->GetStats());
996 return true;
997 }
998 }
999 return false;
1000 }
1001
MovPredSatisfiedPreviously()1002 bool cDeme::MovPredSatisfiedPreviously()
1003 {
1004 for(int i = 0; i < movement_pred_list.Size(); i++) {
1005 if(movement_pred_list[i]->PreviouslySatisfied()) {
1006 movement_pred_list[i]->UpdateStats(m_world->GetStats());
1007 return true;
1008 }
1009 }
1010 return false;
1011 }
1012
GetNumDemePredicates()1013 int cDeme::GetNumDemePredicates()
1014 {
1015 return deme_pred_list.Size();
1016 }
1017
GetNumMessagePredicates()1018 int cDeme::GetNumMessagePredicates()
1019 {
1020 return message_pred_list.Size();
1021 }
1022
GetNumMovementPredicates()1023 int cDeme::GetNumMovementPredicates()
1024 {
1025 return movement_pred_list.Size();
1026 }
1027
GetDemePredicate(int i)1028 cDemePredicate* cDeme::GetDemePredicate(int i)
1029 {
1030 assert(i < deme_pred_list.Size());
1031 return deme_pred_list[i];
1032 }
1033
GetMsgPredicate(int i)1034 cOrgMessagePredicate* cDeme::GetMsgPredicate(int i)
1035 {
1036 assert(i < message_pred_list.Size());
1037 return message_pred_list[i];
1038 }
1039
GetMovPredicate(int i)1040 cOrgMovementPredicate* cDeme::GetMovPredicate(int i)
1041 {
1042 assert(i < movement_pred_list.Size());
1043 return movement_pred_list[i];
1044 }
1045
AddDemeResourceThresholdPredicate(cString resourceName,cString comparisonOperator,double threasholdValue)1046 void cDeme::AddDemeResourceThresholdPredicate(cString resourceName, cString comparisonOperator, double threasholdValue)
1047 {
1048 cDemeResourceThresholdPredicate* pred =
1049 new cDemeResourceThresholdPredicate(resourceName, comparisonOperator, threasholdValue);
1050 deme_pred_list.Add(pred);
1051
1052 cString name = resourceName + " " + comparisonOperator +
1053 cStringUtil::Stringf(" %f", threasholdValue);
1054 m_world->GetStats().AddDemeResourceThresholdPredicate(name);
1055 }
1056
AddEventReceivedCenterPred(int times)1057 void cDeme::AddEventReceivedCenterPred(int times)
1058 {
1059 if (cell_events.Size() == 0) {
1060 cerr<<"Error: An EventReceivedCenterPred cannot be created until a CellEvent is added.\n";
1061 exit(1);
1062 }
1063
1064 for (int i = 0; i < cell_events.Size(); i++) {
1065 if (!cell_events[i].IsDead()) {
1066 int sink_cell = GetCellID(GetSize()/2);
1067 cOrgMessagePred_EventReceivedCenter* pred =
1068 new cOrgMessagePred_EventReceivedCenter(&cell_events[i], sink_cell, times);
1069 m_world->GetStats().AddMessagePredicate(pred);
1070 message_pred_list.Add(pred);
1071 }
1072 }
1073 }
1074
AddEventReceivedLeftSidePred(int times)1075 void cDeme::AddEventReceivedLeftSidePred(int times)
1076 {
1077 if (cell_events.Size() == 0) {
1078 cerr<<"Error: An EventReceivedLeftSidePred cannot be created until a CellEvent is added.\n";
1079 exit(1);
1080 }
1081
1082 for (int i = 0; i < cell_events.Size(); i++) {
1083 if (!cell_events[i].IsDead()) {
1084 cOrgMessagePred_EventReceivedLeftSide* pred =
1085 new cOrgMessagePred_EventReceivedLeftSide(&cell_events[i], m_world->GetPopulation(), times);
1086 m_world->GetStats().AddMessagePredicate(pred);
1087 message_pred_list.Add(pred);
1088 }
1089 }
1090 }
1091
AddEventMoveCenterPred(int times)1092 void cDeme::AddEventMoveCenterPred(int times)
1093 {
1094 if (cell_events.Size() == 0) {
1095 cerr<<"Error: An EventMovedIntoCenter cannot be created until a CellEvent is added.\n";
1096 exit(1);
1097 }
1098 for (int i = 0; i < cell_events.Size(); i++) {
1099 if (!cell_events[i].IsDead()) {
1100 cOrgMovementPred_EventMovedIntoCenter* pred = new cOrgMovementPred_EventMovedIntoCenter(&cell_events[i], m_world->GetPopulation(), times);
1101 m_world->GetStats().AddMovementPredicate(pred);
1102 movement_pred_list.Add(pred);
1103 }
1104 }
1105 }
1106
1107
AddEventMoveBetweenTargetsPred(int times)1108 void cDeme::AddEventMoveBetweenTargetsPred(int times)
1109 {
1110 if (cell_events.Size() == 0) {
1111 cerr << "Error: An EventMoveBetweenTargets cannot be created until at least one CellEvent is added.\n";
1112 exit(1);
1113 }
1114
1115 tVector<cDemeCellEvent *> alive_events;
1116
1117 for (int i = 0; i < cell_events.Size(); i++) {
1118 if (!cell_events[i].IsDead()) {
1119 alive_events.Add(&cell_events[i]);
1120 }
1121 }
1122
1123 cOrgMovementPred_EventMovedBetweenTargets* pred = new cOrgMovementPred_EventMovedBetweenTargets(alive_events, m_world->GetPopulation(), times);
1124 m_world->GetStats().AddMovementPredicate(pred);
1125 movement_pred_list.Add(pred);
1126 }
1127
1128
AddEventEventNUniqueIndividualsMovedIntoTargetPred(int times)1129 void cDeme::AddEventEventNUniqueIndividualsMovedIntoTargetPred(int times)
1130 {
1131 if (cell_events.Size() == 0) {
1132 cerr << "Error: An EventMovedIntoCenter cannot be created until a CellEvent is added.\n";
1133 exit(1);
1134 }
1135 for (int i = 0; i < cell_events.Size(); i++) {
1136 if (!cell_events[i].IsDead()) {
1137 cOrgMovementPred_EventNUniqueIndividualsMovedIntoTarget* pred =
1138 new cOrgMovementPred_EventNUniqueIndividualsMovedIntoTarget(&cell_events[i], m_world->GetPopulation(), times);
1139 m_world->GetStats().AddMovementPredicate(pred);
1140 movement_pred_list.Add(pred);
1141 }
1142 }
1143 }
1144
AddPheromone(int absolute_cell_id,double value,cAvidaContext & ctx)1145 void cDeme::AddPheromone(int absolute_cell_id, double value, cAvidaContext& ctx)
1146 {
1147 assert(cell_ids[0] <= absolute_cell_id);
1148 assert(absolute_cell_id <= cell_ids[cell_ids.GetSize()-1]);
1149
1150 // cPopulation& pop = m_world->GetPopulation();
1151
1152 int relative_cell_id = GetRelativeCellID(absolute_cell_id);
1153 tArray<double> cell_resources = deme_resource_count.GetCellResources(relative_cell_id, ctx);
1154
1155 for (int i = 0; i < deme_resource_count.GetSize(); i++) {
1156 if (strcmp(deme_resource_count.GetResName(i), "pheromone") == 0) {
1157 // There should only be one "pheromone" resource, so no need to divvy value up
1158 cell_resources[i] = value;
1159 }
1160 else {
1161 cell_resources[i] = 0;
1162 }
1163 }
1164
1165 //It appears that ModifyCell adds the amount of resources specified in the cell_resources array, so I'm just
1166 //settign the element to the value I want to add instead of setting the element to the current value plus the amount to add
1167 // Ask Ben why he does it differently in GiveBackCellEnergy()
1168
1169 deme_resource_count.ModifyCell(ctx, cell_resources, relative_cell_id);
1170
1171 // CellData-based version
1172 //const int newval = pop.GetCell(absolute_cell_id).GetCellData() + (int) round(value);
1173 //pop.GetCell(absolute_cell_id).SetCellData(newval);
1174
1175 } //End AddPheromone()
1176
GetSpatialResource(int rel_cellid,int resource_id,cAvidaContext & ctx) const1177 double cDeme::GetSpatialResource(int rel_cellid, int resource_id, cAvidaContext& ctx) const
1178 {
1179 assert(rel_cellid >= 0);
1180 assert(rel_cellid < GetSize());
1181 assert(resource_id >= 0);
1182 assert(resource_id < deme_resource_count.GetSize());
1183
1184 tArray<double> cell_resources = deme_resource_count.GetCellResources(rel_cellid, ctx);
1185 return cell_resources[resource_id];
1186 }
1187
AdjustSpatialResource(cAvidaContext & ctx,int rel_cellid,int resource_id,double amount)1188 void cDeme::AdjustSpatialResource(cAvidaContext& ctx, int rel_cellid, int resource_id, double amount)
1189 {
1190 assert(rel_cellid >= 0);
1191 assert(rel_cellid < GetSize());
1192 assert(resource_id >= 0);
1193 assert(resource_id < deme_resource_count.GetSize());
1194
1195 tArray<double> res_change;
1196 res_change.Resize(deme_resource_count.GetSize(), 0);
1197 res_change[resource_id] = amount;
1198
1199 deme_resource_count.ModifyCell(ctx, res_change, rel_cellid);
1200 }
1201
AdjustResource(cAvidaContext & ctx,int resource_id,double amount)1202 void cDeme::AdjustResource(cAvidaContext& ctx, int resource_id, double amount)
1203 {
1204 double new_amount = deme_resource_count.Get(ctx, resource_id) + amount;
1205 deme_resource_count.Set(ctx, resource_id, new_amount);
1206 }
1207
GetSlotFlowRate() const1208 int cDeme::GetSlotFlowRate() const
1209 {
1210 vector<pair<int, int> >::const_iterator iter = event_slot_end_points.begin();
1211 while (iter != event_slot_end_points.end()) {
1212 if (GetAge() <= (*iter).first) {
1213 return (*iter).second;
1214 }
1215 iter++;
1216 }
1217
1218 // assert(false); // slots must be of equal size and fit perfectally in deme lifetime
1219 return 0;
1220 }
1221
1222
1223 // Return whether or not this deme is treatable at the given age (updates).
1224 // If a deme is not treatable, this function will always return false.
IsTreatableAtAge(const int age)1225 bool cDeme::IsTreatableAtAge(const int age)
1226 {
1227 if (isTreatable()) {
1228 if (treatment_ages.size() == 0) { return false; } // implies treatable every update
1229 set<int>::iterator it;
1230 it = treatment_ages.find(age);
1231 if(it != treatment_ages.end()) return true;
1232 }
1233
1234 return false;
1235
1236 } //End cDeme::IsTreatableAtAge()
1237
1238
1239 /*! Retrieve the network object, initializing it as needed.
1240 */
GetNetwork()1241 cDemeNetwork& cDeme::GetNetwork()
1242 {
1243 InitNetworkCreation();
1244 return *m_network;
1245 }
1246
ResetInputs(cAvidaContext & ctx)1247 void cDeme::ResetInputs(cAvidaContext& ctx)
1248 {
1249 m_world->GetEnvironment().SetupInputs(ctx, m_inputs);
1250 }
1251
1252 //@JJB**
GetNextDemeInput(cAvidaContext & ctx)1253 int cDeme::GetNextDemeInput(cAvidaContext& ctx)
1254 {
1255 // Hack protection
1256 if (!m_inputs.GetSize()) ResetInputs(ctx);
1257 m_input_pointer %= m_inputs.GetSize();
1258 return m_inputs[m_input_pointer++];
1259 }
1260
1261 //@JJB**
DoDemeInput(int value)1262 void cDeme::DoDemeInput(int value)
1263 {
1264 m_input_buf.Add(value);
1265 }
1266
1267 //@JJB**
DoDemeOutput(cAvidaContext & ctx,int value)1268 void cDeme::DoDemeOutput(cAvidaContext& ctx, int value)
1269 {
1270 // Add value to the output buffer
1271 m_output_buf.Add(value);
1272
1273 // Needed to setup taskctx, but will not actually be used
1274 tList<tBuffer<int> > other_input_list;
1275 tList<tBuffer<int> > other_output_list;
1276 tSmartArray<int> ext_mem;
1277
1278 // Setup the task context
1279 cTaskContext taskctx(NULL, m_input_buf, m_output_buf, other_input_list, other_output_list,
1280 ext_mem, false, NULL, this);
1281 taskctx.SetTaskStates(&m_task_states);
1282
1283 const cEnvironment& env = m_world->GetEnvironment();
1284 const int num_resources = env.GetResourceLib().GetSize();
1285 const int num_tasks = env.GetNumTasks();
1286 const int num_reactions = env.GetReactionLib().GetSize();
1287
1288 // Create a new reaction result
1289 if (!m_reaction_result) m_reaction_result = new cReactionResult(num_resources, num_tasks, num_reactions);
1290 cReactionResult& result = *m_reaction_result;
1291
1292 // Not actually set up for deme's using resources during reactions
1293 tArray<double> res_in;
1294 tArray<double> rbins_in;
1295
1296 // The environment, evaluates if a task and if a resulting reaction were completed
1297 bool found = env.TestOutput(ctx, result, taskctx, m_task_count, m_reaction_count, res_in, rbins_in);
1298
1299 // No task completed, end here
1300 if (found == false) {
1301 result.Invalidate();
1302 return;
1303 }
1304
1305 // Update records with the results..
1306
1307 // Update deme's task and reaction counters
1308 for (int i = 0; i < num_tasks; i++) {
1309 if (result.TaskDone(i) == true) {
1310 m_task_count[i]++;
1311 }
1312 }
1313 for (int i = 0; i < num_reactions; i++) {
1314 if (result.ReactionTriggered(i) == true) {
1315 m_reaction_count[i]++;
1316 }
1317 }
1318
1319 // Update stats task totals
1320 for (int i = 0; i < num_tasks; i++) {
1321 if (result.TaskDone(i) && !m_last_task_count[i]) {
1322 m_world->GetStats().AddNewTaskCount(i);
1323 int prev_num_tasks = 0;
1324 int cur_num_tasks = 0;
1325 for (int j = 0; j < num_tasks; j++) {
1326 if (m_last_task_count[j] > 0) prev_num_tasks++;
1327 if (m_task_count[j] > 0) cur_num_tasks++;
1328 }
1329 m_world->GetStats().AddOtherTaskCounts(i, prev_num_tasks, cur_num_tasks);
1330 }
1331 }
1332
1333 // Update stats reaction totals
1334 for (int i = 0; i < num_reactions; i++) {
1335 m_cur_reaction_add_reward[i] += result.GetReactionAddBonus(i);
1336 if (result.ReactionTriggered(i) && last_reaction_count[i] == 0) {
1337 m_world->GetStats().AddNewReactionCount(i);
1338 }
1339 }
1340
1341 // Update deme's merit bonus from reaction
1342 m_cur_bonus *= result.GetMultBonus(); //**
1343 m_cur_bonus += result.GetAddBonus(); //**
1344
1345 //**
1346 //if (result.GetActiveDeme()) {
1347 // double deme_bonus = GetHeritableDemeMerit().GetDouble();
1348 // deme_bonus *= result.GetMultDemeBonus();
1349 // deme_bonus += result.GetAddDemeBonus();
1350 // UpdateHeritableDemeMerit(deme_bonus);
1351 //}
1352
1353 // If applying merit changes immediately, update the deme's merit merit
1354 if (m_world->GetConfig().MERIT_INC_APPLY_IMMEDIATE.Get()) {
1355 UpdateCurMerit();
1356 }
1357
1358 //**
1359 //for (int i = 0; i < num_tasks; i++) {
1360 // if (result.TaskDone(i) == true) AddCurTask(i);
1361 //}
1362 //for (int i = 0; i < num_reactions; i++) {
1363 // if (result.ReactionTriggered(i) == true) AddCurReaction(i);
1364 //}
1365
1366 result.Invalidate();
1367 }
1368
1369 //@JJB**
UpdateCurMerit()1370 void cDeme::UpdateCurMerit()
1371 {
1372 double merit_base;
1373 if (m_world->GetConfig().BASE_MERIT_METHOD.Get() == BASE_MERIT_CONST) {
1374 merit_base = m_world->GetConfig().BASE_CONST_MERIT.Get();
1375 } else {
1376 merit_base = 1.0;
1377 }
1378 m_cur_merit = merit_base * m_cur_bonus;
1379 }
1380
1381 //@JJB**
CalcCurMerit()1382 cMerit cDeme::CalcCurMerit()
1383 {
1384 cMerit cur_merit;
1385 double merit_base;
1386 if (m_world->GetConfig().BASE_MERIT_METHOD.Get() == BASE_MERIT_CONST) {
1387 merit_base = m_world->GetConfig().BASE_CONST_MERIT.Get();
1388 } else {
1389 merit_base = 1.0;
1390 }
1391 cur_merit = merit_base * m_cur_bonus;
1392 return cur_merit;
1393 }
1394
1395 // Returns the minimum number of times any of the reactions were performed
MinNumTimesReactionPerformed()1396 int cDeme::MinNumTimesReactionPerformed()
1397 {
1398 int min = 100000000;
1399 for (int i = 0; i < cur_reaction_count.GetSize(); i++ ) {
1400 if (cur_reaction_count[i] < min) { min = cur_reaction_count[i]; }
1401 }
1402 return min;
1403 }
1404
1405
1406 // Track the number of resources used.
GetTotalResourceAmountConsumed() const1407 double cDeme::GetTotalResourceAmountConsumed() const
1408 {
1409 double total = m_total_res_consumed;
1410 return total;
1411 }
1412
1413
1414 /** calculate shannon mutual entropy for the relationship between individuals and tasks. A few notes:
1415 - pi - probability of an individual - 1/# orgs that do something
1416 - pij - probability an individual performs a task:
1417 (num times ind performs this task / num tasks performed by ind) / num orgs that did stuff
1418
1419
1420 Because we are only including organisms that do a task in this calculation, we also need a separate
1421 stat that tracks the percentage of a deme that does stuff. */
1422
GetShannonMutualInformation()1423 double cDeme::GetShannonMutualInformation()
1424 {
1425 // exit if 0 or 1 organism did stuff.
1426 if ((m_num_active == 0) || (m_num_active ==1)) return 0.0;
1427
1428 int num_org = m_shannon_matrix.size();
1429 int num_task = m_shannon_matrix[0].size();
1430 double* ptask_array = new double[num_task];
1431 // create ptasks
1432
1433 for (int j=0; j<num_task; j++){
1434 ptask_array[j] =0;
1435
1436 for (int i=0; i<num_org; i++) {
1437 ptask_array[j] += m_shannon_matrix[i][j];
1438 }
1439 ptask_array[j] /= m_num_active;
1440 }
1441
1442 double shannon_sum = 0.0;
1443 double shannon_change = 0.0;
1444 double pij = 0.0;
1445 double pi = 1.0/m_num_active;
1446 double pj = 0;
1447 double pij_sum = 0.0;
1448 // calculate shannon mutual information
1449 for (int i=0; i<num_org; i++) {
1450 for (int j=0; j<num_task; j++) {
1451 pij = m_shannon_matrix[i][j]/m_num_active;
1452 pj = ptask_array[j];
1453 pij_sum += pij;
1454 if (pi && pj && pij) {
1455 shannon_change= (pij * log(pij / (pi * pj)));
1456 shannon_sum += shannon_change;
1457 }
1458 }
1459 }
1460
1461 // shannon_sum /= log((double)m_num_active);
1462 delete ptask_array;
1463 return shannon_sum;
1464 }
1465
GetNumOrgsPerformedReaction()1466 double cDeme::GetNumOrgsPerformedReaction()
1467 {
1468 return m_num_active;
1469 }
1470
1471
1472 /* Update the task counts for an organism. if we are tracking repro, then the repro count replaces the first task count */
UpdateShannon(cPopulationCell & cell)1473 void cDeme::UpdateShannon(cPopulationCell& cell)
1474 {
1475 int org_react_count = 0;
1476 vector<double> org_row;
1477 if (cell.IsOccupied()) {
1478 cOrganism* organism = cell.GetOrganism();
1479 cPhenotype& phenotype = organism->GetPhenotype();
1480 const tArray<int> curr_react = phenotype.GetCumulativeReactionCount();
1481 org_row.resize(curr_react.GetSize(), 0.0);
1482 // track number of reproductives
1483 if ((phenotype.GetNumDivides() - phenotype.GetNumDivideFailed()) > 0) {
1484 ++m_num_reproductives;
1485 }
1486
1487 for (int j=0; j<curr_react.GetSize(); j++) {
1488
1489 // we are tracking repro as a task
1490 if ((m_world->GetConfig().DEMES_TRACK_SHANNON_INFO.Get() == 2) && (j==0)) {
1491 org_react_count += (phenotype.GetNumDivides() - phenotype.GetNumDivideFailed());
1492 org_row[j] = (phenotype.GetNumDivides() - phenotype.GetNumDivideFailed());
1493 } else {
1494 org_react_count += curr_react[j];
1495 org_row[j] = curr_react[j];
1496 }
1497 }
1498
1499 if (org_react_count > 0) {
1500 // normalize the data for the current organism.
1501 for (int j=0; j<curr_react.GetSize(); j++){
1502 if (org_row[j]) org_row[j] /= org_react_count;
1503 }
1504 m_num_active++;
1505 m_shannon_matrix.push_back(org_row);
1506 }
1507 }
1508 }
1509
UpdateShannonAll()1510 void cDeme::UpdateShannonAll()
1511 {
1512 for (int i=0; i<GetSize(); ++i) {
1513 cPopulationCell& cell = m_world->GetPopulation().GetCell(GetCellID(i));
1514 UpdateShannon(cell);
1515 }
1516 }
1517
GetPercentReproductives()1518 double cDeme::GetPercentReproductives()
1519 {
1520 double per = (m_num_reproductives/((double)injected_count + (double)cur_birth_count))*100;
1521 return per;
1522 }
1523
ClearShannonInformationStats()1524 void cDeme::ClearShannonInformationStats()
1525 {
1526 m_total_res_consumed = 0;
1527 m_switch_penalties = 0;
1528 m_num_active = 0;
1529 m_shannon_matrix.clear();
1530 m_num_reproductives = 0;
1531 }
1532
1533
1534
1535 /* Returns the average number of mutations that have occured to the deme's germline as a the result damage accrued by performing tasks. */
GetAveVarGermMut()1536 std::pair<double, double> cDeme::GetAveVarGermMut()
1537 {
1538
1539 if ((m_world->GetConfig().DEMES_ORGANISM_SELECTION.Get() != 7) &&
1540 (m_world->GetConfig().DEMES_ORGANISM_SELECTION.Get() != 8)) {
1541 cPopulationCell& c = GetCell(0);
1542 if (c.IsOccupied()) { c.GetOrganism()->JoinGermline(); }
1543 }
1544 cDoubleSum mut_count;
1545
1546 for (int i=0; i<GetSize(); ++i) {
1547 cPopulationCell& cell = GetCell(i);
1548 if (cell.IsOccupied()) {
1549 cOrganism* o = cell.GetOrganism();
1550 if (o->IsGermline()) {
1551 mut_count.Add(o->GetNumOfPointMutationsApplied());
1552 }
1553 }
1554 }
1555
1556 return (make_pair(mut_count.Average(), mut_count.Variance()));
1557 }
1558
GetAveVarSomaMut()1559 std::pair<double, double> cDeme::GetAveVarSomaMut()
1560 {
1561
1562 if ((m_world->GetConfig().DEMES_ORGANISM_SELECTION.Get() != 7) &&
1563 (m_world->GetConfig().DEMES_ORGANISM_SELECTION.Get() != 8)) {
1564 cPopulationCell& c = GetCell(0);
1565 if (c.IsOccupied()) { c.GetOrganism()->JoinGermline(); }
1566 }
1567
1568 cDoubleSum mut_count;
1569
1570 for (int i=0; i<GetSize(); ++i) {
1571
1572 cPopulationCell& cell = GetCell(i);
1573 if (cell.IsOccupied()) {
1574 cOrganism* o = cell.GetOrganism();
1575 if (!o->IsGermline()) {
1576 mut_count.Add(o->GetNumOfPointMutationsApplied());
1577 }
1578 }
1579 }
1580
1581 return (make_pair(mut_count.Average(), mut_count.Variance()));
1582 }
1583
GetGermlineNumPercent()1584 std::pair<double, double> cDeme::GetGermlineNumPercent() {
1585 double count = 0;
1586 double total_count = 0;
1587 for (int i=0; i<GetSize(); ++i) {
1588
1589 cPopulationCell& cell = GetCell(i);
1590 if (cell.IsOccupied()) {
1591 cOrganism* o = cell.GetOrganism();
1592 if (o->IsGermline()) {
1593 ++count;
1594 }
1595 ++total_count;
1596 }
1597 }
1598 total_count = (count/total_count) * 100;
1599 return (make_pair(count, total_count));
1600 }
1601
GetAveVarGermWorkLoad()1602 std::pair<double, double> cDeme::GetAveVarGermWorkLoad() {
1603 cDoubleSum work_load;
1604
1605 for (int i=0; i<GetSize(); ++i) {
1606
1607 cPopulationCell& cell = GetCell(i);
1608 if (cell.IsOccupied()) {
1609 cOrganism* o = cell.GetOrganism();
1610 double cur_org_w = 0;
1611 if (o->IsGermline()) {
1612 // Compute workload...
1613 const tArray<int> curr_react = o->GetPhenotype().GetCumulativeReactionCount();
1614 for (int j=0; j<curr_react.GetSize(); j++) {
1615 double weight = 1.0;
1616 // weight each reaction according to.
1617 if (m_world->GetConfig().INST_POINT_MUT_SLOPE.Get() > 0.0) {
1618 weight = m_world->GetConfig().INST_POINT_MUT_SLOPE.Get() * j;
1619 }
1620 // The first task never has any work associated with it.
1621 if (j > 0) {
1622 cur_org_w += weight * curr_react[j];
1623 }
1624 }
1625 work_load.Add(cur_org_w);
1626 }
1627 }
1628 }
1629
1630 return (make_pair(work_load.Average(), work_load.Variance()));
1631
1632 }
1633
GetAveVarSomaWorkLoad()1634 std::pair<double, double> cDeme::GetAveVarSomaWorkLoad() {
1635 cDoubleSum work_load;
1636
1637 for (int i=0; i<GetSize(); ++i) {
1638
1639 cPopulationCell& cell = GetCell(i);
1640 if (cell.IsOccupied()) {
1641 cOrganism* o = cell.GetOrganism();
1642 double cur_org_w = 0;
1643 if (!o->IsGermline()) {
1644 // Compute workload...
1645 const tArray<int> curr_react = o->GetPhenotype().GetCumulativeReactionCount();
1646 for (int j=0; j<curr_react.GetSize(); j++) {
1647 double weight = 1.0;
1648 // weight each reaction according to.
1649 if (m_world->GetConfig().INST_POINT_MUT_SLOPE.Get() > 0.0) {
1650 weight = m_world->GetConfig().INST_POINT_MUT_SLOPE.Get() * j;
1651 }
1652 // The first task never has any work associated with it.
1653 if (j > 0) {
1654 cur_org_w += weight * curr_react[j];
1655 }
1656 }
1657 work_load.Add(cur_org_w);
1658 }
1659 }
1660 }
1661
1662 return (make_pair(work_load.Average(), work_load.Variance()));
1663
1664 }
1665
1666
GetAveVarWorkLoad()1667 std::pair<double, double> cDeme::GetAveVarWorkLoad() {
1668 cDoubleSum work_load;
1669
1670 for (int i=0; i<GetSize(); ++i) {
1671
1672 cPopulationCell& cell = GetCell(i);
1673 if (cell.IsOccupied()) {
1674 cOrganism* o = cell.GetOrganism();
1675 double cur_org_w = 0;
1676 // Compute workload...
1677 const tArray<int> curr_react = o->GetPhenotype().GetCumulativeReactionCount();
1678 for (int j=0; j<curr_react.GetSize(); j++) {
1679 double weight = 1.0;
1680 // weight each reaction according to.
1681 if (m_world->GetConfig().INST_POINT_MUT_SLOPE.Get() > 0.0) {
1682 weight = m_world->GetConfig().INST_POINT_MUT_SLOPE.Get() * j;
1683 }
1684 // The first task never has any work associated with it.
1685 if (j > 0) {
1686 cur_org_w += weight * curr_react[j];
1687 }
1688 work_load.Add(cur_org_w);
1689 }
1690 }
1691 }
1692
1693 return (make_pair(work_load.Average(), work_load.Variance()));
1694
1695 }
1696
1697