diff --git a/source/Evolve/NK.h b/source/Evolve/NK.h index 4e31a0e480..d0a534d6f8 100644 --- a/source/Evolve/NK.h +++ b/source/Evolve/NK.h @@ -141,6 +141,21 @@ namespace emp { return total; } + /// Get the fitness of a site in a bitstring + // (pass by value so can be modified.) + double GetSiteFitness(size_t n, BitVector genome) const { + emp_assert(genome.GetSize() == N, genome.GetSize(), N); + + // Use a double-length genome to easily handle wrap-around. + genome.Resize(N*2); + genome |= (genome << N); + + size_t mask = emp::MaskLow(K+1); + const size_t cur_val = (genome >> n).GetUInt(0) & mask; + return GetFitness(n, cur_val); + } + + void SetState(size_t n, size_t state, double in_fit) { landscape[n][state] = in_fit; } void RandomizeStates(Random & random, size_t num_states=1) { diff --git a/source/Evolve/OEE.h b/source/Evolve/OEE.h new file mode 100644 index 0000000000..6a5bcba486 --- /dev/null +++ b/source/Evolve/OEE.h @@ -0,0 +1,235 @@ +#ifndef EMP_OEE_STATS_H +#define EMP_OEE_STATS_H + +#include "Systematics.h" +#include "bloom_filter.hpp" +#include "base/vector.h" +#include "base/Ptr.h" +#include "tools/set_utils.h" +#include "tools/vector_utils.h" + +#include + +namespace emp { + + // Setup possible types for keeping track of what we've seen for novelty + + template + class SeenSet { + std::set s; + public: + using skel_t = SKEL_TYPE; + // Placeholders to ensure that constructor signature is same as bloom filter + SeenSet(int placeholder_1 = 200000, double placeholder_2 = 0.0001) { ; } + void insert(const skel_t & val) {s.insert(val);} + bool contains(const skel_t & val) {return Has(s, val);} + }; + + class SeenBloomFilter { + bloom_filter b; + + public: + + using skel_t = std::string; + SeenBloomFilter(int bloom_count = 200000, double false_positive = 0.0001) { + bloom_parameters parameters; + + // How many elements roughly do we expect to insert? + parameters.projected_element_count = bloom_count; + + // Maximum tolerable false positive probability? (0,1) + parameters.false_positive_probability = false_positive; + + if (!parameters) + { + std::cout << "Error - Invalid set of bloom filter parameters!" << std::endl; + } + + parameters.compute_optimal_parameters(); + b = bloom_filter(parameters); + } + + void insert(const skel_t & val) {b.insert(val);} + bool contains(const skel_t & val) {return b.contains(val);} + }; + + + template > + class OEETracker { + private: + using systematics_t = SYSTEMATICS_TYPE; + using taxon_t = typename systematics_t::taxon_t; + using info_t = typename systematics_t::info_t; + using hash_t = typename Ptr::hash_t; + using fun_calc_complexity_t = std::function; + using fun_calc_data_t = std::function; // TODO: Allow other skeleton types + + struct snapshot_info_t { + Ptr taxon = nullptr; // This is what the systematics manager has + Ptr skel = nullptr; + int count = 0; // Count of this taxon at time of snapshot + + ~snapshot_info_t() {if (skel){skel.Delete();}} + // bool operator==(const snapshot_info_t & other) const {return other.taxon == taxon;} + }; + + std::deque> snapshots; + std::deque snapshot_times; + Ptr systematics_manager; + + std::map prev_coal_set; + // std::unordered_set seen; + + fun_calc_data_t skeleton_fun; + fun_calc_complexity_t complexity_fun; + int generation_interval = 10; + int resolution = 10; + + DataManager data_nodes; + SEEN_TYPE seen; + bool prune_top; + + public: + OEETracker(Ptr s, fun_calc_data_t d, fun_calc_complexity_t c, bool remove_top = false, int bloom_count = 200000, double bloom_false_positive = .0001) : + systematics_manager(s), skeleton_fun(d), complexity_fun(c), seen(bloom_count, bloom_false_positive), prune_top(remove_top) { + + emp_assert(s->GetStoreAncestors(), "OEE tracker only works with systematics manager where store_ancestor is set to true"); + + data_nodes.New("change"); + data_nodes.New("novelty"); + data_nodes.New("diversity"); + data_nodes.New("complexity"); + + } + + ~OEETracker() {} + + int GetResolution() const {return resolution;} + int GetGenerationInterval() const {return generation_interval;} + + void SetResolution(int r) {resolution = r;} + void SetGenerationInterval(int g) {generation_interval = g;} + + void Update(size_t gen, int ud = -1) { + if (Mod((int)gen, resolution) == 0) { + if (ud == -1) { + ud = gen; + } + auto & sys_active = systematics_manager->GetActive(); + + snapshots.emplace_back(sys_active.size()); + int i = 0; + for (auto tax : sys_active) { + snapshots.back()[i].taxon = tax; + info_t info = tax->GetInfo(); + snapshots.back()[i].skel.New(skeleton_fun(info)); + snapshots.back()[i].count = tax->GetNumOrgs(); + i++; + } + + snapshot_times.push_back(ud); + if ((int)snapshots.size() > generation_interval/resolution + 1) { + if (prune_top) { + systematics_manager->RemoveBefore(snapshot_times.front() - 1); + } + snapshot_times.pop_front(); + + snapshots.pop_front(); + } + CalcStats(ud); + } + } + + void CalcStats(size_t ud) { + std::map coal_set = CoalescenceFilter(ud); + int change = 0; + int novelty = 0; + double most_complex = 0; + double diversity = 0; + if (coal_set.size() > 0) { + diversity = Entropy(coal_set, [](std::pair entry){return entry.second;}); + } + + for (auto & tax : coal_set) { + if (!Has(prev_coal_set, tax.first)) { + change++; + } + if (!seen.contains(tax.first)) { + novelty++; + seen.insert(tax.first); + } + double complexity = complexity_fun(tax.first); + if (complexity > most_complex) { + most_complex = complexity; + } + } + + data_nodes.Get("change").Add(change); + data_nodes.Get("novelty").Add(novelty); + data_nodes.Get("diversity").Add(diversity); + data_nodes.Get("complexity").Add(most_complex); + + std::swap(prev_coal_set, coal_set); + } + + std::map CoalescenceFilter(size_t ud) { + + emp_assert(emp::Mod(generation_interval, resolution) == 0, "Generation interval must be a multiple of resolution", generation_interval, resolution); + + std::map res; + + if ((int)snapshots.size() <= generation_interval/resolution) { + return res; + } + + std::set> extant_canopy_roots = systematics_manager->GetCanopyExtantRoots(snapshot_times.front()); + for ( snapshot_info_t & t : snapshots.front()) { + if (Has(extant_canopy_roots, t.taxon)) { + if (Has(res, *(t.skel))) { + res[*(t.skel)] += t.count; + } else { + res[*(t.skel)] = t.count; + } + } + } + + return res; + } + + + Ptr> GetDataNode(const std::string & name) { + return &(data_nodes.Get(name)); + } + + }; + + // Helper function for skeletonization when organism is a sequence of + + // Assumes org is sequence of inst_type + template + emp::vector Skeletonize(ORG_TYPE & org, const INST_TYPE null_value, std::function fit_fun) { + emp_assert(org.size() > 0, "Empty org passed to skeletonize"); + + emp::vector skeleton; + // Some fitness functions may require the org to be const and smoe may require it to not be + // We can let the compiler deducce whetehr ORG_TYPE is const or not. + // But the test org really needs to not be const + typename std::remove_const::type test_org = ORG_TYPE(org); + double fitness = fit_fun(test_org); + + for (int i = 0; i < (int)org.size(); i++) { + test_org[i] = null_value; + double new_fitness = fit_fun(test_org); + if (new_fitness < fitness) { + skeleton.push_back(org[i]); + } + test_org[i] = org[i]; + } + + return skeleton; + } + + +} + +#endif \ No newline at end of file diff --git a/source/Evolve/Resource.h b/source/Evolve/Resource.h index 17ab22fdc7..96c9ea8a58 100644 --- a/source/Evolve/Resource.h +++ b/source/Evolve/Resource.h @@ -60,8 +60,8 @@ namespace emp { }; template - void ResourceSelect(World & world, const emp::vector< std::function > & extra_funs, - emp::vector & pools, size_t t_size, size_t tourny_count=1, double frac = .0025, double max_bonus = 5, double cost = 0, bool use_base = true) { + void ResourceSelect(World & world, emp::vector< std::function > & extra_funs, + emp::vector & pools, size_t t_size, size_t tourny_count=1, double frac = .0025, double max_bonus = 5, double cost = 0, bool use_base = true, double min_score = 0) { emp_assert(world.GetFitFun(), "Must define a base fitness function"); emp_assert(world.GetSize() > 0); @@ -98,7 +98,7 @@ namespace emp { cur_fit = emp::Pow(cur_fit, 2.0); // if (org_id==0) {std::cout << "Allele: " << world[org_id][ex_id] <<" Curr fit: " << extra_funs[ex_id](world[org_id]) << " Curr fit squared: " << cur_fit << " Amount: " << pools[ex_id].GetAmount() << " Frac: " << frac;} cur_fit *= frac*(pools[ex_id].GetAmount()-cost); - if (cur_fit > 0) { + if (cur_fit > min_score) { cur_fit -= cost; } else { cur_fit = 0; diff --git a/source/Evolve/Systematics.h b/source/Evolve/Systematics.h index 51f5a589f4..d353d49626 100644 --- a/source/Evolve/Systematics.h +++ b/source/Evolve/Systematics.h @@ -26,9 +26,11 @@ #include #include #include +#include #include "../base/Ptr.h" #include "../control/Signal.h" +#include "../data/DataFile.h" #include "../data/DataManager.h" #include "../data/DataNode.h" #include "../tools/info_theory.h" @@ -37,6 +39,9 @@ #include "../tools/stats.h" #include "../tools/string_utils.h" +#include "SystematicsAnalysis.h" +#include "World_structure.h" + namespace emp { /// The systematics manager allows an optional second template type that @@ -52,6 +57,21 @@ namespace emp { using has_phen_t = std::false_type; }; /// The default - an empty struct + struct fitness { + using has_fitness_t = std::true_type; + using has_mutations_t = std::false_type; + using has_phen_t = std::false_type; + + DataNode fitness; /// This taxon's fitness (for assessing deleterious mutational steps) + void RecordFitness(double fit) { + fitness.Add(fit); + } + + const double GetFitness() const { + return fitness.GetMean(); + } + }; + template struct mut_landscape_info { /// Track information related to the mutational landscape /// Maps a string representing a type of mutation to a count representing @@ -110,12 +130,14 @@ namespace emp { size_t id; ///< ID for this Taxon (Unique within this Systematics) const info_t info; ///< Details for the organims associated within this taxanomic group. Ptr parent; ///< Pointer to parent group (nullptr if injected) + std::set > offspring; ///< Pointers to all immediate offspring taxa size_t num_orgs; ///< How many organisms currently exist of this group? size_t tot_orgs; ///< How many organisms have ever existed of this group? size_t num_offspring; ///< How many direct offspring groups exist from this one. size_t total_offspring; ///< How many total extant offspring taxa exist from this one (i.e. including indirect) size_t depth; ///< How deep in tree is this node? (Root is 0) double origination_time; ///< When did this taxon first appear in the population? + double destruction_time; ///< When did this taxon leave the population? DATA_STRUCT data; ///< A struct for storing additional information about this taxon @@ -125,12 +147,18 @@ namespace emp { Taxon(size_t _id, const info_t & _info, Ptr _parent=nullptr) : id (_id), info(_info), parent(_parent) , num_orgs(0), tot_orgs(0), num_offspring(0), total_offspring(0) - , depth(parent ? (parent->depth+1) : 0) { ; } - Taxon(const Taxon &) = delete; + , depth(parent ? (parent->depth+1) : 0) + , destruction_time(std::numeric_limits::infinity()) { ; } + // Taxon(const Taxon &) = delete; + Taxon(const Taxon &) = default; // TODO: Check with Charles about this Taxon(Taxon &&) = default; Taxon & operator=(const Taxon &) = delete; Taxon & operator=(Taxon &&) = default; + bool operator<(const Taxon & other) const { + return id < other.GetID(); + } + /// Get a unique ID for this taxon; IDs are assigned sequentially, so newer taxa have higher IDs. size_t GetID() const { return id; } @@ -139,6 +167,7 @@ namespace emp { /// Retrieve a pointer to the parent Taxon. Ptr GetParent() const { return parent; } + void NullifyParent() {parent = nullptr;} /// Get the number of living organisms currently associated with this Taxon. size_t GetNumOrgs() const { return num_orgs; } @@ -155,14 +184,25 @@ namespace emp { data_t & GetData() {return data;} const data_t & GetData() const {return data;} + std::set > GetOffspring() {return offspring;} + + void SetData(data_t d) {data = d;} + double GetOriginationTime() const {return origination_time;} void SetOriginationTime(double time) {origination_time = time;} + double GetDestructionTime() const {return destruction_time;} + void SetDestructionTime(double time) {destruction_time = time;} + /// Add a new organism to this Taxon. void AddOrg() { ++num_orgs; ++tot_orgs; } /// Add a new offspring Taxon to this one. - void AddOffspring() { ++num_offspring; AddTotalOffspring();} + void AddOffspring(Ptr offspring_tax) { + ++num_offspring; + offspring.insert(offspring_tax); + AddTotalOffspring(); + } /// Recursively increment total offspring count for this and all ancestors // Should this be protected or private or something? @@ -187,10 +227,15 @@ namespace emp { return num_orgs; } + void RemoveFromOffspring(Ptr offspring_tax) { + offspring.erase(offspring_tax); + } + /// Remove and offspring taxa after its entire sub-tree has died out (pruning) - bool RemoveOffspring() { - emp_assert(num_offspring > 0); + bool RemoveOffspring(Ptr offspring_tax) { + emp_assert(num_offspring > 0, num_offspring, id); --num_offspring; + RemoveFromOffspring(offspring_tax); // If we are out of BOTH offspring and organisms, this Taxon should deactivate. return num_orgs || num_offspring; @@ -224,6 +269,7 @@ namespace emp { size_t org_count; ///< How many organisms are currently active? size_t total_depth; ///< Sum of taxa depths for calculating average. size_t num_roots; ///< How many distint injected ancestors are currently in population? + int max_depth; ///< Depth of deepest taxon. -1 means needs to be recalculated size_t next_id; ///< What ID value should the next new taxon have? size_t curr_update; @@ -234,7 +280,7 @@ namespace emp { SystematicsBase(bool _active=true, bool _anc=true, bool _all=false, bool _pos=true) : store_active(_active), store_ancestors(_anc), store_outside(_all) , archive(store_ancestors || store_outside), store_position(_pos), track_synchronous(false) - , org_count(0), total_depth(0), num_roots(0), next_id(0), curr_update(0) { ; } + , org_count(0), total_depth(0), num_roots(0), max_depth(0), next_id(0), curr_update(0) { ; } virtual ~SystematicsBase(){;} @@ -265,6 +311,9 @@ namespace emp { /// How many independent trees are being tracked? size_t GetNumRoots() const { return num_roots; } + /// What ID will the next taxon have? + size_t GetNextID() const {return next_id;} + /// What is the average phylogenetic depth of organisms in the population? double GetAveDepth() const { return ((double) total_depth) / (double) org_count; } @@ -326,21 +375,25 @@ namespace emp { virtual size_t GetNumOutside() const = 0; virtual size_t GetTreeSize() const = 0; virtual size_t GetNumTaxa() const = 0; + virtual int GetMaxDepth() = 0; virtual int GetPhylogeneticDiversity() const = 0; virtual double GetMeanPairwiseDistance(bool branch_only) const = 0; virtual double GetSumPairwiseDistance(bool branch_only) const = 0; virtual double GetVariancePairwiseDistance(bool branch_only) const = 0; virtual emp::vector GetPairwiseDistances(bool branch_only) const = 0; + virtual int SackinIndex() const = 0; + virtual double CollessLikeIndex() const = 0; virtual int GetMRCADepth() const = 0; - virtual void AddOrg(ORG && org, int pos, int update, bool next) = 0; - virtual void AddOrg(ORG & org, int pos, int update, bool next) = 0; - virtual bool RemoveOrg(int pos) = 0; - virtual bool RemoveNextOrg(int pos) = 0; + virtual void AddOrg(ORG && org, WorldPosition pos, int update) = 0; + virtual void AddOrg(ORG & org, WorldPosition pos, int update) = 0; + virtual bool RemoveOrg(WorldPosition pos, int time=-1) = 0; + virtual void RemoveOrgAfterRepro(WorldPosition pos, int time=-1) = 0; + // virtual bool RemoveNextOrg(WorldPosition pos, int time=-1) = 0; virtual void PrintStatus(std::ostream & os) const = 0; virtual double CalcDiversity() const = 0; virtual void Update() = 0; virtual void SetNextParent(int pos) = 0; - + virtual void SetNextParent(WorldPosition & pos) = 0; }; /// @brief A tool to track phylogenetic relationships among organisms. @@ -355,12 +408,16 @@ namespace emp { class Systematics : public SystematicsBase { private: using parent_t = SystematicsBase; + public: using taxon_t = Taxon; + using info_t = ORG_INFO; + private: using hash_t = typename Ptr::hash_t; using fun_calc_info_t = std::function; fun_calc_info_t calc_info_fun; Ptr next_parent; + Ptr most_recent; using parent_t::store_active; using parent_t::store_ancestors; @@ -373,7 +430,10 @@ namespace emp { using parent_t::num_roots; using parent_t::next_id; using parent_t::curr_update; + using parent_t::max_depth; + + public: using typename parent_t::data_ptr_t; using parent_t::GetNumActive; using parent_t::GetNumAncestors; @@ -395,7 +455,8 @@ namespace emp { using parent_t::GetMRCADepth; using parent_t::AddOrg; using parent_t::RemoveOrg; - using parent_t::RemoveNextOrg; + using parent_t::RemoveOrgAfterRepro; + // using parent_t::RemoveNextOrg; // using parent_t::Parent; using parent_t::PrintStatus; // using parent_t::PrintLineage; @@ -412,15 +473,35 @@ namespace emp { using parent_t::AddVolatilityDataNode; using parent_t::AddUniqueTaxaDataNode; using parent_t::AddMutationCountDataNode; + using parent_t::GetMaxDepth; + + struct SnapshotInfo { + using snapshot_fun_t = std::function; + snapshot_fun_t fun; + std::string key; + std::string desc; + + SnapshotInfo(const snapshot_fun_t & _fun, const std::string & _key, const std::string & _desc="") + : fun(_fun), + key(_key), + desc(_desc) + { ; } + }; + + emp::vector user_snapshot_funs; std::unordered_set< Ptr, hash_t > active_taxa; ///< A set of all living taxa. std::unordered_set< Ptr, hash_t > ancestor_taxa; ///< A set of all dead, ancestral taxa. std::unordered_set< Ptr, hash_t > outside_taxa; ///< A set of all dead taxa w/o descendants. + Ptr to_be_removed = nullptr; + int removal_time = -1; + int removal_pos = -1; + emp::vector > taxon_locations; emp::vector > next_taxon_locations; - Signal)> on_new_sig; ///< Trigger when any organism is pruned from tree + Signal, ORG & org)> on_new_sig; ///< Trigger when any organism is pruned from tree Signal)> on_prune_sig; ///< Trigger when any organism is pruned from tree mutable Ptr mrca; ///< Most recent common ancestor in the population. @@ -429,11 +510,12 @@ namespace emp { void Prune(Ptr taxon); /// Called when an offspring taxa has been deleted. - void RemoveOffspring(Ptr taxon); + void RemoveOffspring(Ptr offspring, Ptr taxon); /// Called when there are no more living members of a taxon. There may be descendants. - void MarkExtinct(Ptr taxon); + void MarkExtinct(Ptr taxon, int time=-1); + public: @@ -464,6 +546,15 @@ namespace emp { void Update() { ++curr_update; if (track_synchronous) { + + // Clear pending removal + if (to_be_removed != nullptr) { + RemoveOrg(to_be_removed, removal_time); + taxon_locations[removal_pos] = nullptr; + to_be_removed = nullptr; + removal_pos = -1; + } + std::swap(taxon_locations, next_taxon_locations); next_taxon_locations.resize(0); } @@ -475,7 +566,8 @@ namespace emp { std::unordered_set< Ptr, hash_t > * GetActivePtr() { return &active_taxa; } const std::unordered_set< Ptr, hash_t > & GetActive() const { return active_taxa; } const std::unordered_set< Ptr, hash_t > & GetAncestors() const { return ancestor_taxa; } - + const std::unordered_set< Ptr, hash_t > & GetOutside() const { return outside_taxa; } + /// How many taxa are still active in the population? size_t GetNumActive() const { return active_taxa.size(); } @@ -491,12 +583,36 @@ namespace emp { /// How many taxa are stored in total? size_t GetNumTaxa() const { return GetTreeSize() + GetNumOutside(); } + int GetMaxDepth() { + if (max_depth != -1) { + return max_depth; + } + + for (auto tax : active_taxa) { + int depth = tax->GetDepth(); + if (depth > max_depth) { + max_depth = depth; + } + } + return max_depth; + } + + void SetNextParent(WorldPosition & pos) { + emp_assert(pos.IsActive() || !pos.IsValid()); + if (!pos.IsValid()) { + next_parent = nullptr; + } else { + next_parent = taxon_locations[pos.GetIndex()]; + } + } + void SetNextParent(int pos) { emp_assert(pos < (int)taxon_locations.size(), "Invalid parent", pos, taxon_locations.size()); if (pos == -1) { next_parent = nullptr; } else { emp_assert(pos >= 0, "Invalid parent", pos); + emp_assert(taxon_locations[pos], pos); next_parent = taxon_locations[pos]; } } @@ -505,7 +621,15 @@ namespace emp { next_parent = p; } - SignalKey OnNew(std::function)> & fun) { return on_new_sig.AddAction(fun); } + Ptr GetNextParent() { + return next_parent; + } + + Ptr GetMostRecent() { + return most_recent; + } + + SignalKey OnNew(std::function, ORG & org)> & fun) { return on_new_sig.AddAction(fun); } /// Privide a function for Systematics to call each time a taxon is about to be pruned. /// Trigger: Taxon is about to be killed @@ -526,7 +650,7 @@ namespace emp { return node; } - virtual data_ptr_t AddPairwiseDistanceDataNode(const std::string & name = "pairwise_distances") { + virtual data_ptr_t AddPairwiseDistanceDataNode(const std::string & name = "pairwise_distance") { auto node = AddDataNode(name); node->AddPullSet([this](){ return GetPairwiseDistances(); @@ -643,6 +767,21 @@ namespace emp { return node; } + /// Add a new snapshot function. + /// When a snapshot of the systematics is taken, in addition to the default + /// set of functions, all user-added snapshot functions are run. Functions + /// take a reference to a taxon as input and return the string to be dumped + /// in the file at the given key. + void AddSnapshotFun(const std::function & fun, + const std::string & key, const std::string & desc="") { + user_snapshot_funs.emplace_back(fun, key, desc); + } + + bool IsTaxonAt(int id) { + emp_assert(id < (int) taxon_locations.size(), "Invalid taxon location", id, taxon_locations.size()); + return taxon_locations[id]; + } + Ptr GetTaxonAt(int id) { emp_assert(id < (int) taxon_locations.size(), "Invalid taxon location", id, taxon_locations.size()); emp_assert(taxon_locations[id], "No taxon at specified location"); @@ -709,12 +848,12 @@ namespace emp { emp_assert(time != -1 && "Invalid time - are you passing time to rg?", time); emp_assert(time >= tax->GetOriginationTime() - && "GetEvolutionaryDistinctiveness recieved a time that is earlier than the taxon's origination time."); + && "GetEvolutionaryDistinctiveness recieved a time that is earlier than the taxon's origination time.", tax->GetOriginationTime(), time); while (test_taxon) { - emp_assert(test_taxon->GetOriginationTime() != -1 && - "Invalid time - are you passing time to rg?"); + // emp_assert(test_taxon->GetOriginationTime() != -1 && + // "Invalid time - are you passing time to rg?", time); depth += time - test_taxon->GetOriginationTime(); // std::cout << "Tax: " << test_taxon->GetID() << " depth: " << depth << " time: " << time << " Orig: " << test_taxon->GetOriginationTime() << " divisor: " << divisor << std::endl; @@ -865,7 +1004,7 @@ namespace emp { next_pointers.erase(tax.first); Ptr test_taxon = tax.first->GetParent(); - while (test_taxon && test_taxon->GetNumOff() == 1 && test_taxon->GetNumOrgs() < 0) { + while (test_taxon && test_taxon->GetNumOff() == 1 && test_taxon->GetNumOrgs() == 0) { if (!branch_only) { for (size_t i = 0; i < new_dist_vec.size(); i++){ new_dist_vec[i]++; @@ -900,6 +1039,35 @@ namespace emp { } + /** + * Returns a vector containing all taxa from @param time_point that were + * + * */ + std::set> GetCanopyExtantRoots(int time_point = 0) const { + // NOTE: This could be made faster by doing something similar to the pairwise distance + // function + + std::set< Ptr> result; + // std::cout << "starting " << time_point << std::endl; + for (Ptr tax : active_taxa) { + // std::cout << tax->GetInfo() << std::endl; + while (tax) { + // std::cout << tax->GetInfo() << " " << tax->GetOriginationTime() << " " << tax->GetDestructionTime() << std::endl; + if (tax->GetOriginationTime() <= time_point && tax->GetDestructionTime() > time_point ) { + result.insert(tax); + // std::cout << "inserting " << tax->GetInfo() << std::endl; + break; + } + tax = tax->GetParent(); + } + } + + return result; + + } + + + /** Counts the total number of ancestors between @param tax and MRCA, if there is one. If * there is no common ancestor, distance to the root of this tree is calculated instead. */ @@ -941,7 +1109,255 @@ namespace emp { return depth; } + /** Calculate Sackin Index of this tree (Sackin, 1972; reviewed in Shao, 1990). + * Measures tree balance + */ + int SackinIndex() const { + int sackin = 0; + + for (auto taxon : active_taxa) { + sackin += GetBranchesToRoot(taxon) + 1; // Sackin index counts root as branch + } + + return sackin; + } + + + // Graph ToGraph() const { + + // std::map, int> ids; + // int next_id = 0; + + // for (Ptr tax : active_taxa) { + // ids[tax] = next_id; + // next_id++; + // } + + // for (Ptr tax : ancestor_taxa) { + // ids[tax] = next_id; + // next_id++; + // } + + // for (Ptr tax : outside_taxa) { + // ids[tax] = next_id; + // next_id++; + // } + + // Graph g(next_id); + + // for (Ptr tax : active_taxa) { + // if (tax->GetParent()) { + // g.AddEdge(ids[tax->GetParent()], ids[tax]); + // } + // } + + // for (Ptr tax : ancestor_taxa) { + // if (tax->GetParent()) { + // g.AddEdge(ids[tax->GetParent()], ids[tax]); + // } + // } + + // for (Ptr tax : outside_taxa) { + // if (tax->GetParent()) { + // g.AddEdge(ids[tax->GetParent()], ids[tax]); + // } + // } + + // return g; + // } + + // Graph ToMinimalGraph() const { + // std::map, int> ids; + // int next_id = 0; + + // for (Ptr tax : active_taxa) { + // if (tax->GetNumOff() == 1) { + // continue; + // } + // ids[tax] = next_id; + // next_id++; + // } + + // for (Ptr tax : ancestor_taxa) { + // if (tax->GetNumOff() == 1) { + // continue; + // } + // ids[tax] = next_id; + // next_id++; + // } + + // for (Ptr tax : outside_taxa) { + // if (tax->GetNumOff() == 1) { + // continue; + // } + // ids[tax] = next_id; + // next_id++; + // } + + // Graph g(next_id); + + // for (Ptr tax : active_taxa) { + // if (tax->GetNumOff() == 1) { + // continue; + // } + + // Ptr parent = tax->GetParent(); + // while (parent) { + // if (parent->GetNumOff() == 1) { + // parent = parent->GetParent(); + // } else { + // g.AddEdge(ids[parent], ids[tax]); + // } + // } + // } + + // for (Ptr tax : ancestor_taxa) { + // if (tax->GetNumOff() == 1) { + // continue; + // } + + // Ptr parent = tax->GetParent(); + // while (parent) { + // if (parent->GetNumOff() == 1) { + // parent = parent->GetParent(); + // } else { + // g.AddEdge(ids[parent], ids[tax]); + // } + // } + // } + + // for (Ptr tax : outside_taxa) { + // if (tax->GetNumOff() == 1) { + // continue; + // } + + // Ptr parent = tax->GetParent(); + // while (parent) { + // if (parent->GetNumOff() == 1) { + // parent = parent->GetParent(); + // } else { + // g.AddEdge(ids[parent], ids[tax]); + // } + // } + // } + + // return g; + // } + + struct CollessStruct { + double total = 0; + emp::vector ns; + }; + + CollessStruct RecursiveCollessStep(Ptr curr) const { + CollessStruct result; + + while (curr->GetNumOff() == 1) { + curr = *(curr->GetOffspring().begin()); + } + + if (curr->GetNumOff() == 0) { + result.ns.push_back(0); // Node itself is calculated at level above + return result; + } + + for (Ptr off : curr->GetOffspring()) { + // std::cout << "Recursing on ID: " << off->GetID() << " Offspring: " << off->GetTotalOffspring() << std::endl; + + CollessStruct new_result = RecursiveCollessStep(off); + result.ns.push_back(Sum(new_result.ns) + log(off->GetOffspring().size() + exp(1))); + result.total += new_result.total; + } + + // std::cout << "Evaluating: " << curr->GetID() << std::endl; + + double med = Median(result.ns); + double sum_diffs = 0; + // std::cout << "Median: " << med << std::endl; + for (double n : result.ns) { + // std::cout << n << std::endl; + sum_diffs += std::abs(n-med); + } + // std::cout << "Sumdiffs: " << sum_diffs << " n: " << result.ns.size() << " average: " << sum_diffs/result.ns.size() << std::endl; + result.total += sum_diffs/result.ns.size(); + return result; + } + + /** Calculate Colless Index of this tree (Colless, 1982; reviewed in Shao, 1990). + * Measures tree balance. The standard Colless index only works for bifurcating trees, + * so this will be a Colless-like Index, as suggested in + * "Sound Colless-like balance indices for multifurcating trees" (Mir, 2018, PLoS One) + */ + double CollessLikeIndex() const { + GetMRCA(); + + return RecursiveCollessStep(mrca).total; + } + + + + void RemoveBefore(int ud) { + + // @ELD: This would be such a nice way to do it + // but we can't because we need to notify offspring + // when their parents are un-tracked + // std::set> to_remove; + // for (Ptr tax : ancestor_taxa) { + // if (tax->GetDestructionTime() < ud) { + // to_remove.insert(tax); + // } + // } + // for (Ptr tax : to_remove) { + // ancestor_taxa.erase(tax); + // tax.Delete(); + // } + + std::map, std::set>> to_remove; + + for (Ptr tax : active_taxa) { + Ptr curr = tax; + + while (curr && !CanRemove(curr->GetParent(), ud)) { + curr = curr->GetParent(); + } + + if (curr) { + Ptr next = curr->GetParent(); + while (next) { + to_remove[next].insert(curr); + curr = next; + next = next->GetParent(); + } + } + } + // std::cout << "About to remove " << to_remove.size() << " orgs" << std::endl; + for (std::pair, std::set>> el : to_remove) { + emp_assert(el.first->GetDestructionTime() < ud, el.first->GetDestructionTime(), ud); + if (el.first->GetNumOff() == el.second.size()) { + // Everything is account for + for (auto tax : el.second) { + tax->NullifyParent(); + } + ancestor_taxa.erase(el.first); + el.first.Delete(); + } + } + + } + + bool CanRemove(Ptr t, int ud) { + if (!t) { + return false; + } + while (t) { + if (t->GetNumOrgs() > 0 || t->GetDestructionTime() >= ud) { + return false; + } + t = t->GetParent(); + } + return true; + } /// Request a pointer to the Most-Recent Common Ancestor for the population. Ptr GetMRCA() const; @@ -953,22 +1369,25 @@ namespace emp { /// If you would like the systematics manager to track taxon age, you can also supply /// the update at which the taxon is being added. /// return a pointer for the associated taxon. - void AddOrg(ORG && org, int pos, int update=-1, bool next=false); - Ptr AddOrg(ORG && org, int pos, Ptr parent=nullptr, int update=-1, bool next=false); - Ptr AddOrg(ORG && org, Ptr parent=nullptr, int update=-1, bool next=false); + void AddOrg(ORG && org, WorldPosition pos, int update=-1); + Ptr AddOrg(ORG && org, WorldPosition pos, Ptr parent=nullptr, int update=-1); + Ptr AddOrg(ORG && org, Ptr parent=nullptr, int update=-1); - void AddOrg(ORG & org, int pos, int update=-1, bool next=false); - Ptr AddOrg(ORG & org, int pos, Ptr parent=nullptr, int update=-1, bool next=false); - Ptr AddOrg(ORG & org, Ptr parent=nullptr, int update=-1, bool next=false); + void AddOrg(ORG & org, WorldPosition pos, int update=-1); + Ptr AddOrg(ORG & org, WorldPosition pos, Ptr parent=nullptr, int update=-1); + Ptr AddOrg(ORG & org, Ptr parent=nullptr, int update=-1); /// Remove an instance of an organism; track when it's gone. - bool RemoveOrg(int pos); - bool RemoveOrg(Ptr taxon); + bool RemoveOrg(WorldPosition pos, int time=-1); + bool RemoveOrg(Ptr taxon, int time=-1); + + void RemoveOrgAfterRepro(WorldPosition pos, int time=-1); + void RemoveOrgAfterRepro(Ptr taxon, int time=-1); /// Remove org from next population (for use with synchronous generations) - bool RemoveNextOrg(int pos); - bool RemoveNextOrg(Ptr taxon); + // bool RemoveNextOrg(WorldPosition pos, int time=-1); + // bool RemoveNextOrg(Ptr taxon, int time=-1); /// Climb up a lineage... Ptr Parent(Ptr taxon) const; @@ -979,6 +1398,8 @@ namespace emp { /// Print whole lineage. void PrintLineage(Ptr taxon, std::ostream & os=std::cout) const; + void Snapshot(const std::string & file_path) const; + /// Calculate the genetic diversity of the population. double CalcDiversity() const; @@ -994,16 +1415,21 @@ namespace emp { template void Systematics::Prune(Ptr taxon) { on_prune_sig.Trigger(taxon); - RemoveOffspring( taxon->GetParent() ); // Notify parent of the pruning. + RemoveOffspring( taxon, taxon->GetParent() ); // Notify parent of the pruning. if (store_ancestors) ancestor_taxa.erase(taxon); // Clear from ancestors set (if there) if (store_outside) outside_taxa.insert(taxon); // Add to outside set (if tracked) - else taxon.Delete(); // ...or else get rid of it. + else { + if (taxon == mrca) { + mrca = nullptr; + } + taxon.Delete(); // ...or else get rid of it. + } } template - void Systematics::RemoveOffspring(Ptr taxon) { + void Systematics::RemoveOffspring(Ptr offspring, Ptr taxon) { if (!taxon) { num_roots--; return; } // Offspring was root; remove and return. - bool still_active = taxon->RemoveOffspring(); // Taxon still active w/ 1 fewer offspring? + bool still_active = taxon->RemoveOffspring(offspring); // Taxon still active w/ 1 fewer offspring? if (!still_active) Prune(taxon); // If out of offspring, remove from tree. // If the taxon is still active AND the is the current mrca AND now has only one offspring, @@ -1013,10 +1439,15 @@ namespace emp { // Mark a taxon extinct if there are no more living members. There may be descendants. template - void Systematics::MarkExtinct(Ptr taxon) { + void Systematics::MarkExtinct(Ptr taxon, int time) { emp_assert(taxon); emp_assert(taxon->GetNumOrgs() == 0); + if (max_depth == (int)taxon->GetDepth()) { + // We no longer know the max depth + max_depth = -1; + } + if (taxon->GetParent()) { // Update extant descendant count for all ancestors taxon->GetParent()->RemoveTotalOffspring(); @@ -1024,11 +1455,20 @@ namespace emp { if (store_active) active_taxa.erase(taxon); if (!archive) { // If we don't archive taxa, delete them. + for (Ptr off_tax : taxon->GetOffspring()) { + off_tax->NullifyParent(); + } + taxon.Delete(); return; } + // std::cout << "About to set destruction time " << time << std::endl; + // Only need to track destruction time if we're archiving taxa + taxon->SetDestructionTime(time); - if (store_ancestors) ancestor_taxa.insert(taxon); // Move taxon to ancestors... + if (store_ancestors) { + ancestor_taxa.insert(taxon); // Move taxon to ancestors... + } if (taxon->GetNumOff() == 0) Prune(taxon); // ...and prune from there if needed. } @@ -1049,7 +1489,9 @@ namespace emp { Ptr test_taxon = candidate->GetParent(); while (test_taxon) { emp_assert(test_taxon->GetNumOff() >= 1); - if (test_taxon->GetNumOff() > 1) candidate = test_taxon; + // If the test_taxon is dead, we only want to update candidate when we hit a new branch point + // If test_taxon is still alive, though, we always need to update it + if (test_taxon->GetNumOff() > 1 || test_taxon->GetNumOrgs() > 0) candidate = test_taxon; test_taxon = test_taxon->GetParent(); } mrca = candidate; @@ -1071,10 +1513,10 @@ namespace emp { // Can't return a pointer for the associated taxon because of obnoxious inheritance problems template // Ptr::taxon_t> - void Systematics::AddOrg(ORG & org, int pos, int update, bool next) { + void Systematics::AddOrg(ORG & org, WorldPosition pos, int update) { emp_assert(store_position, "Trying to pass position to a systematics manager that can't use it"); // emp_assert(next_parent, "Adding organism with no parent specified and no next_parent set"); - AddOrg(org, pos, next_parent, update, next); + AddOrg(org, pos, next_parent, update); next_parent = nullptr; } @@ -1082,10 +1524,10 @@ namespace emp { // Can't return a pointer for the associated taxon because of obnoxious inheritance problems template // Ptr::taxon_t> - void Systematics::AddOrg(ORG && org, int pos, int update, bool next) { + void Systematics::AddOrg(ORG && org, WorldPosition pos, int update) { emp_assert(store_position, "Trying to pass position to a systematics manager that can't use it"); // emp_assert(next_parent, "Adding organism with no parent specified and no next_parent set"); - AddOrg(org, pos, next_parent, update, next); + AddOrg(org, pos, next_parent, update); next_parent = nullptr; } @@ -1093,30 +1535,32 @@ namespace emp { // Version for if you aren't tracking positions template Ptr::taxon_t> - Systematics::AddOrg(ORG & org, Ptr parent, int update, bool next) { - return AddOrg(org, -1, parent, update, next); + Systematics::AddOrg(ORG & org, Ptr parent, int update) { + return AddOrg(org, -1, parent, update); } // Version for if you aren't tracking positions template Ptr::taxon_t> - Systematics::AddOrg(ORG && org, Ptr parent, int update, bool next) { - return AddOrg(org, -1, parent, update, next); + Systematics::AddOrg(ORG && org, Ptr parent, int update) { + emp_assert(!store_position && + "Trying to add org to position-tracking systematics manager without position. Either specify a valid position or turn of position tracking for systematic manager.", store_position); + return AddOrg(org, WorldPosition::invalid_id, parent, update); } // Add information about a new organism, including its stored info and parent's taxon; // return a pointer for the associated taxon. template Ptr::taxon_t> - Systematics::AddOrg(ORG && org, int pos, Ptr parent, int update, bool next) { - return AddOrg(org, pos, parent, update, next); + Systematics::AddOrg(ORG && org, WorldPosition pos, Ptr parent, int update) { + return AddOrg(org, pos, parent, update); } // Add information about a new organism, including its stored info and parent's taxon; // return a pointer for the associated taxon. template Ptr::taxon_t> - Systematics::AddOrg(ORG & org, int pos, Ptr parent, int update, bool next) { + Systematics::AddOrg(ORG & org, WorldPosition pos, Ptr parent, int update) { org_count++; // Keep count of how many organisms are being tracked. ORG_INFO info = calc_info_fun(org); @@ -1131,75 +1575,95 @@ namespace emp { } cur_taxon = NewPtr(++next_id, info, parent); // Build new taxon. - on_new_sig.Trigger(cur_taxon); + if (max_depth != -1 && (int)cur_taxon->GetDepth() > max_depth) { + max_depth = cur_taxon->GetDepth(); + } + on_new_sig.Trigger(cur_taxon, org); if (store_active) active_taxa.insert(cur_taxon); // Store new taxon. - if (parent) parent->AddOffspring(); // Track tree info. + if (parent) parent->AddOffspring(cur_taxon); // Track tree info. cur_taxon->SetOriginationTime(update); } - - if (store_position && pos >= 0) { - if (next) { - if (pos >= (int)next_taxon_locations.size()) { - next_taxon_locations.resize(pos+1); + // std::cout << "about to store poisiton" << std::endl; + if (store_position && pos.GetIndex() >= 0) { + if (pos.GetPopID()) { + if (pos.GetIndex() >= next_taxon_locations.size()) { + next_taxon_locations.resize(pos.GetIndex()+1); } - next_taxon_locations[pos] = cur_taxon; + next_taxon_locations[pos.GetIndex()] = cur_taxon; } else { - if (pos >= (int)taxon_locations.size()) { - taxon_locations.resize(pos+1); + if (pos.GetIndex() >= taxon_locations.size()) { + taxon_locations.resize(pos.GetIndex()+1); } - taxon_locations[pos] = cur_taxon; + taxon_locations[pos.GetIndex()] = cur_taxon; } } cur_taxon->AddOrg(); // Record the current organism in its taxon. total_depth += cur_taxon->GetDepth(); // Track the total depth (for averaging) + + if (to_be_removed) { + RemoveOrg(to_be_removed, removal_time); + to_be_removed = nullptr; + } + + most_recent = cur_taxon; return cur_taxon; // Return the taxon used. } - // Remove an instance of an organism; track when it's gone. template - bool Systematics::RemoveOrg(int pos) { + void Systematics::RemoveOrgAfterRepro(WorldPosition pos, int time) { emp_assert(store_position, "Trying to remove org based on position from systematics manager that doesn't track it."); - emp_assert(pos < (int)taxon_locations.size(), "Invalid position requested for removal", pos, taxon_locations.size()); - bool active = RemoveOrg(taxon_locations[pos]); - taxon_locations[pos] = nullptr; - return active; - } + + if (pos.GetIndex() >= taxon_locations.size() || !taxon_locations[pos.GetIndex()]) { + // There's not actually a taxon here + return; + } - // Remove an instance of an organism; track when it's gone. + RemoveOrgAfterRepro(taxon_locations[pos.GetIndex()], time); + removal_pos = pos.GetIndex(); + } + template - bool Systematics::RemoveOrg(Ptr taxon) { - emp_assert(taxon); - - // Update stats - org_count--; - total_depth -= taxon->GetDepth(); - - // emp_assert(Has(active_taxa, taxon)); - const bool active = taxon->RemoveOrg(); - if (!active) MarkExtinct(taxon); - - return active; + void Systematics::RemoveOrgAfterRepro(Ptr taxon, int time) { + if (to_be_removed != nullptr) { + RemoveOrg(to_be_removed, removal_time); + taxon_locations[removal_pos] = nullptr; + to_be_removed = nullptr; + removal_pos = -1; + } + to_be_removed = taxon; + // std::cout << "Setting remove time to " << time << std::endl; + removal_time = time; } + // Remove an instance of an organism; track when it's gone. template - bool Systematics::RemoveNextOrg(int pos) { - emp_assert(track_synchronous, "Calling RemoveNextOrg on non-synchronous population. Did you mean to use RemoveOrg?"); + bool Systematics::RemoveOrg(WorldPosition pos, int time) { emp_assert(store_position, "Trying to remove org based on position from systematics manager that doesn't track it."); - emp_assert(pos < (int)next_taxon_locations.size(), "Invalid position requested for removal", pos, taxon_locations.size()); - - bool active = RemoveOrg(next_taxon_locations[pos]); - next_taxon_locations[pos] = nullptr; - return active; + + if (pos.GetPopID() == 0) { + emp_assert(pos.GetIndex() < taxon_locations.size(), "Invalid position requested for removal", pos.GetIndex(), taxon_locations.size()); + bool active = false; + if (taxon_locations[pos.GetIndex()]) { + //TODO: Figure out how this can ever not be true + active = RemoveOrg(taxon_locations[pos.GetIndex()], time); + } + taxon_locations[pos.GetIndex()] = nullptr; + return active; + } else { + emp_assert(pos.GetIndex() < next_taxon_locations.size(), "Invalid position requested for removal", pos.GetIndex(), taxon_locations.size()); + bool active = RemoveOrg(next_taxon_locations[pos.GetIndex()], time); + next_taxon_locations[pos.GetIndex()] = nullptr; + return active; + } } // Remove an instance of an organism; track when it's gone. template - bool Systematics::RemoveNextOrg(Ptr taxon) { - emp_assert(track_synchronous, "Calling RemoveNextOrg on non-synchronous population. Did you mean to use RemoveOrg?"); + bool Systematics::RemoveOrg(Ptr taxon, int time) { emp_assert(taxon); // Update stats @@ -1208,12 +1672,11 @@ namespace emp { // emp_assert(Has(active_taxa, taxon)); const bool active = taxon->RemoveOrg(); - if (!active) MarkExtinct(taxon); + if (!active) MarkExtinct(taxon, time); return active; } - // Climb up a lineage... template Ptr::taxon_t> Systematics::Parent(Ptr taxon) const { @@ -1231,25 +1694,26 @@ namespace emp { << " store_outside=" << store_outside << " archive=" << archive << " next_id=" << next_id + << " synchronous=" << track_synchronous << std::endl; os << "Active count: " << active_taxa.size(); for (const auto & x : active_taxa) { os << " [" << x->GetID() << "|" << x->GetNumOrgs() << "," << x->GetNumOff() << "|" - << ((bool) x->GetParent()) << "]"; + << (x->GetParent() ? emp::to_string(x->GetParent()->GetID()) : "null") << "]"; } os << std::endl; os << "Ancestor count: " << ancestor_taxa.size(); for (const auto & x : ancestor_taxa) { os << " [" << x->GetID() << "|" << x->GetNumOrgs() << "," << x->GetNumOff() << "|" - << ((bool) x->GetParent()) << "]"; + << (x->GetParent() ? emp::to_string(x->GetParent()->GetID()) : "null") << "]"; } os << std::endl; os << "Outside count: " << outside_taxa.size(); for (const auto & x : outside_taxa) { os << " [" << x->GetID() << "|" << x->GetNumOrgs() << "," << x->GetNumOff() << "|" - << ((bool) x->GetParent()) << "]"; + << (x->GetParent() ? emp::to_string(x->GetParent()->GetID()) : "null") << "]"; } os << std::endl; } @@ -1264,12 +1728,110 @@ namespace emp { } } + /// Take a snapshot of current state of taxon phylogeny. + /// WARNING: Current, this function assumes one parent taxon per-taxon. + template + void Systematics::Snapshot(const std::string & file_path) const { + emp::DataFile file(file_path); + Ptr cur_taxon; + emp::vector> wrapped_user_funs; + // Add default functions to file. + // - id: systematics ID for taxon + std::function get_id = [&cur_taxon]() { + return cur_taxon->GetID(); + }; + file.AddFun(get_id, "id", "Systematics ID for this taxon."); + + // - ancestor_list: ancestor list for taxon + std::function get_ancestor_list = [&cur_taxon]() -> std::string { + if (cur_taxon->GetParent() == nullptr) { return "[NONE]"; } + return "[" + to_string(cur_taxon->GetParent()->GetID()) + "]"; + }; + file.AddFun(get_ancestor_list, "ancestor_list", "Ancestor list for this taxon."); + + // - origin_time: When did this taxon first appear in the population? + std::function get_origin_time = [&cur_taxon]() { + return cur_taxon->GetOriginationTime(); + }; + file.AddFun(get_origin_time, "origin_time", "When did this taxon first appear in the population?"); + + // - destruction_time: When did this taxon leave the population? + std::function get_destruction_time = [&cur_taxon]() { + return cur_taxon->GetDestructionTime(); + }; + file.AddFun(get_destruction_time, "destruction_time", "When did this taxon leave the population?"); + + // - num_orgs: How many organisms currently exist of this group? + std::function get_num_orgs = [&cur_taxon]() { + return cur_taxon->GetNumOrgs(); + }; + file.AddFun(get_num_orgs, "num_orgs", "How many organisms currently exist of this group?"); + + // - tot_orgs: How many organisms have ever existed of this group? + std::function get_tot_orgs = [&cur_taxon]() { + return cur_taxon->GetTotOrgs(); + }; + file.AddFun(get_tot_orgs, "tot_orgs", "How many organisms have ever existed of this group?"); + + // - num_offspring: How many direct offspring groups exist from this one. + std::function get_num_offspring = [&cur_taxon]() { + return cur_taxon->GetNumOff(); + }; + file.AddFun(get_num_offspring, "num_offspring", "How many direct offspring groups exist from this one."); + + // - total_offspring: How many total extant offspring taxa exist from this one (i.e. including indirect) + std::function get_total_offspring = [&cur_taxon]() { + return cur_taxon->GetTotalOffspring(); + }; + file.AddFun(get_total_offspring, "total_offspring", "How many total extant offspring taxa exist from this one (i.e. including indirect)"); + + // - depth: How deep in tree is this node? (Root is 0) + std::function get_depth = [&cur_taxon]() { + return cur_taxon->GetDepth(); + }; + file.AddFun(get_depth, "depth", "How deep in tree is this node? (Root is 0)"); + + // Add user-added functions to file. + for (size_t i = 0; i < user_snapshot_funs.size(); ++i) { + wrapped_user_funs.emplace_back([this, i, &cur_taxon]() -> std::string { + return user_snapshot_funs[i].fun(*cur_taxon); + }); + } + + // Need to add file functions after wrapping to preserve integrity of + // function reference being passed to the data file object. + for (size_t i = 0; i < user_snapshot_funs.size(); ++i) { + file.AddFun(wrapped_user_funs[i], user_snapshot_funs[i].key, user_snapshot_funs[i].desc); + } + + // Output header information. + file.PrintHeaderKeys(); + + // Update file w/active taxa information + for (auto tax : active_taxa) { + cur_taxon = tax; + file.Update(); + } + + // Update file w/ancestor taxa information + for (auto tax : ancestor_taxa) { + cur_taxon = tax; + file.Update(); + } + + // Update file w/outside taxa information + for (auto tax : outside_taxa) { + cur_taxon = tax; + file.Update(); + } + + } + // Calculate the genetic diversity of the population. template double Systematics::CalcDiversity() const { return emp::Entropy(active_taxa, [](Ptr x){ return x->GetNumOrgs(); }, (double) org_count); } - } #endif diff --git a/source/Evolve/SystematicsAnalysis.h b/source/Evolve/SystematicsAnalysis.h index 9779f0e401..3ffe09c8f9 100644 --- a/source/Evolve/SystematicsAnalysis.h +++ b/source/Evolve/SystematicsAnalysis.h @@ -12,7 +12,6 @@ #include "Systematics.h" - // Mutation info functions. Assumes each taxon has a struct containing an unordered map // with keys that are strings indicating types of mutations and keys that are numbers // indicating the number of that type of mutation that occurred to make this taxon from @@ -20,6 +19,20 @@ namespace emp { + template + Ptr FindDominant(systematics_t & s) { + double best = -999999; + Ptr best_tax = nullptr; + for (Ptr tax : s.GetActive()) { + double f = tax->GetData().GetFitness(); + if (f > best) { + best = f; + best_tax = tax; + } + } + return best_tax; + } + /// Returns the total number of times a mutation of type @param type /// that along @param taxon 's lineage. (Different from CountMuts in /// that CountMuts sums them whereas CountMutSteps would count two @@ -142,18 +155,14 @@ namespace emp { /// along @param taxon 's lineage. template int CountUniquePhenotypes(Ptr taxon) { - int count = 0; std::setGetData().phenotype)> seen; while (taxon) { - if (!Has(seen, taxon->GetData().phenotype)) { - count++; - seen.insert(taxon->GetData().phenotype); - } + seen.insert(taxon->GetData().phenotype); taxon = taxon->GetParent(); } - return count; + return seen.size(); } }; diff --git a/source/Evolve/World.h b/source/Evolve/World.h index 9955d18dc3..3ed08a652e 100644 --- a/source/Evolve/World.h +++ b/source/Evolve/World.h @@ -411,6 +411,14 @@ namespace emp { systematics_labels.erase(label) ; } + template + Ptr> AddSystematics(std::function calc_taxon, bool active=true, bool anc=true, bool all=true, bool pos=true, std::string label="systematics" ) { + Ptr> sys_ptr; + sys_ptr.New(calc_taxon, active, anc, all, pos); + AddSystematics(sys_ptr, label); + return sys_ptr; + } + template void AddSystematics(Ptr > s, std::string label="systematics") { if (Has(systematics_labels, label)) { @@ -897,7 +905,7 @@ namespace emp { if (pos.IsActive()) { before_placement_sig.Trigger(*new_org, pos.GetIndex()); } for (Ptr > s : systematics) { - s->SetNextParent((int) p_pos.GetIndex()); + s->SetNextParent(p_pos); } // Clear out any old organism at this position. @@ -911,7 +919,7 @@ namespace emp { // Track the new systematics info for (Ptr > s : systematics) { - s->AddOrg(*new_org, (int) pos.GetIndex(), (int) update, !pos.IsActive()); + s->AddOrg(*new_org, pos, (int) update); } SetupOrg(*new_org, pos, *random_ptr); @@ -932,14 +940,12 @@ namespace emp { if (pos.IsActive()) { --num_orgs; // Track one fewer organisms in the population if (cache_on) ClearCache(id); // Delete any cached info about this organism - for (Ptr > s : systematics) { - s->RemoveOrg((int) pos.GetIndex()); // Notify systematics about organism removal - } - } else { - for (Ptr > s : systematics) { - s->RemoveNextOrg((int) pos.GetIndex()); // Notify systematics about organism removal - } + } + + for (Ptr > s : systematics) { + s->RemoveOrgAfterRepro(pos, update); // Notify systematics about organism removal } + } template diff --git a/source/Evolve/World_output.h b/source/Evolve/World_output.h index 2e52050cf2..a3c8820257 100644 --- a/source/Evolve/World_output.h +++ b/source/Evolve/World_output.h @@ -16,6 +16,20 @@ namespace emp { + template + DataFile & AddOEEFile(WORLD_TYPE & world, OEE_TYPE & oee_tracker, const std::string fpath = "oee_data.csv") { + auto & file = world.SetupFile(fpath); + std::function get_update = [&world](){return world.GetUpdate();}; + + file.AddFun(get_update, "update", "Update"); + file.AddCurrent(*oee_tracker.GetDataNode("change"), "change", "change potential"); + file.AddCurrent(*oee_tracker.GetDataNode("novelty"), "novelty", "novelty potential"); + file.AddCurrent(*oee_tracker.GetDataNode("diversity"), "ecology", "ecology potential"); + file.AddCurrent(*oee_tracker.GetDataNode("complexity"), "complexity", "complexity potential"); + file.PrintHeaderKeys(); + return file; + } + template DataFile & AddPhylodiversityFile(WORLD_TYPE & world, int systematics_id=0, const std::string & fpath="phylodiversity.csv"){ auto & file = world.SetupFile(fpath); diff --git a/source/Evolve/World_structure.h b/source/Evolve/World_structure.h index 6332f11286..baf119a0dd 100644 --- a/source/Evolve/World_structure.h +++ b/source/Evolve/World_structure.h @@ -227,7 +227,7 @@ namespace emp { } const size_t num_traits = traits.GetSize(); size_t trait_size = 1; - while (Pow(trait_size+1, num_traits) < world.GetSize()) trait_size++; + while (std::pow(trait_size+1, num_traits) < world.GetSize()) trait_size++; trait_counts.resize(num_traits, trait_size); SetMapElites(world, traits, trait_counts); } diff --git a/source/Evolve/bloom_filter.hpp b/source/Evolve/bloom_filter.hpp new file mode 100644 index 0000000000..975467c813 --- /dev/null +++ b/source/Evolve/bloom_filter.hpp @@ -0,0 +1,735 @@ +/* + ********************************************************************* + * * + * Open Bloom Filter * + * * + * Author: Arash Partow - 2000 * + * URL: http://www.partow.net * + * URL: http://www.partow.net/programming/hashfunctions/index.html * + * * + * Copyright notice: * + * Free use of the Open Bloom Filter Library is permitted under the * + * guidelines and in accordance with the MIT License. * + * http://www.opensource.org/licenses/MIT * + * * + ********************************************************************* +*/ + + +#ifndef INCLUDE_BLOOM_FILTER_HPP +#define INCLUDE_BLOOM_FILTER_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + + +static const std::size_t bits_per_char = 0x08; // 8 bits in 1 char(unsigned) + +static const unsigned char bit_mask[bits_per_char] = { + 0x01, //00000001 + 0x02, //00000010 + 0x04, //00000100 + 0x08, //00001000 + 0x10, //00010000 + 0x20, //00100000 + 0x40, //01000000 + 0x80 //10000000 + }; + +class bloom_parameters +{ +public: + + bloom_parameters() + : minimum_size(1), + maximum_size(std::numeric_limits::max()), + minimum_number_of_hashes(1), + maximum_number_of_hashes(std::numeric_limits::max()), + projected_element_count(10000), + false_positive_probability(1.0 / projected_element_count), + random_seed(0xA5A5A5A55A5A5A5AULL) + {} + + virtual ~bloom_parameters() + {} + + inline bool operator!() + { + return (minimum_size > maximum_size) || + (minimum_number_of_hashes > maximum_number_of_hashes) || + (minimum_number_of_hashes < 1) || + (0 == maximum_number_of_hashes) || + (0 == projected_element_count) || + (false_positive_probability < 0.0) || + (std::numeric_limits::infinity() == std::abs(false_positive_probability)) || + (0 == random_seed) || + (0xFFFFFFFFFFFFFFFFULL == random_seed); + } + + // Allowable min/max size of the bloom filter in bits + unsigned long long int minimum_size; + unsigned long long int maximum_size; + + // Allowable min/max number of hash functions + unsigned int minimum_number_of_hashes; + unsigned int maximum_number_of_hashes; + + // The approximate number of elements to be inserted + // into the bloom filter, should be within one order + // of magnitude. The default is 10000. + unsigned long long int projected_element_count; + + // The approximate false positive probability expected + // from the bloom filter. The default is assumed to be + // the reciprocal of the projected_element_count. + double false_positive_probability; + + unsigned long long int random_seed; + + struct optimal_parameters_t + { + optimal_parameters_t() + : number_of_hashes(0), + table_size(0) + {} + + unsigned int number_of_hashes; + unsigned long long int table_size; + }; + + optimal_parameters_t optimal_parameters; + + virtual bool compute_optimal_parameters() + { + /* + Note: + The following will attempt to find the number of hash functions + and minimum amount of storage bits required to construct a bloom + filter consistent with the user defined false positive probability + and estimated element insertion count. + */ + + if (!(*this)) + return false; + + double min_m = std::numeric_limits::infinity(); + double min_k = 0.0; + double k = 1.0; + + while (k < 1000.0) + { + const double numerator = (- k * projected_element_count); + const double denominator = std::log(1.0 - std::pow(false_positive_probability, 1.0 / k)); + + const double curr_m = numerator / denominator; + + if (curr_m < min_m) + { + min_m = curr_m; + min_k = k; + } + + k += 1.0; + } + + optimal_parameters_t& optp = optimal_parameters; + + optp.number_of_hashes = static_cast(min_k); + + optp.table_size = static_cast(min_m); + + optp.table_size += (((optp.table_size % bits_per_char) != 0) ? (bits_per_char - (optp.table_size % bits_per_char)) : 0); + + if (optp.number_of_hashes < minimum_number_of_hashes) + optp.number_of_hashes = minimum_number_of_hashes; + else if (optp.number_of_hashes > maximum_number_of_hashes) + optp.number_of_hashes = maximum_number_of_hashes; + + if (optp.table_size < minimum_size) + optp.table_size = minimum_size; + else if (optp.table_size > maximum_size) + optp.table_size = maximum_size; + + return true; + } + +}; + +class bloom_filter +{ +protected: + + typedef unsigned int bloom_type; + typedef unsigned char cell_type; + typedef std::vector table_type; + +public: + + bloom_filter() + : salt_count_(0), + table_size_(0), + projected_element_count_(0), + inserted_element_count_ (0), + random_seed_(0), + desired_false_positive_probability_(0.0) + {} + + bloom_filter(const bloom_parameters& p) + : projected_element_count_(p.projected_element_count), + inserted_element_count_(0), + random_seed_((p.random_seed * 0xA5A5A5A5) + 1), + desired_false_positive_probability_(p.false_positive_probability) + { + salt_count_ = p.optimal_parameters.number_of_hashes; + table_size_ = p.optimal_parameters.table_size; + + generate_unique_salt(); + + bit_table_.resize(table_size_ / bits_per_char, static_cast(0x00)); + } + + bloom_filter(const bloom_filter& filter) + { + this->operator=(filter); + } + + inline bool operator == (const bloom_filter& f) const + { + if (this != &f) + { + return + (salt_count_ == f.salt_count_ ) && + (table_size_ == f.table_size_ ) && + (bit_table_.size() == f.bit_table_.size() ) && + (projected_element_count_ == f.projected_element_count_ ) && + (inserted_element_count_ == f.inserted_element_count_ ) && + (random_seed_ == f.random_seed_ ) && + (desired_false_positive_probability_ == f.desired_false_positive_probability_) && + (salt_ == f.salt_ ) && + (bit_table_ == f.bit_table_ ) ; + } + else + return true; + } + + inline bool operator != (const bloom_filter& f) const + { + return !operator==(f); + } + + inline bloom_filter& operator = (const bloom_filter& f) + { + if (this != &f) + { + salt_count_ = f.salt_count_; + table_size_ = f.table_size_; + bit_table_ = f.bit_table_; + salt_ = f.salt_; + + projected_element_count_ = f.projected_element_count_; + inserted_element_count_ = f.inserted_element_count_; + + random_seed_ = f.random_seed_; + + desired_false_positive_probability_ = f.desired_false_positive_probability_; + } + + return *this; + } + + virtual ~bloom_filter() + {} + + inline bool operator!() const + { + return (0 == table_size_); + } + + inline void clear() + { + std::fill(bit_table_.begin(), bit_table_.end(), static_cast(0x00)); + inserted_element_count_ = 0; + } + + inline void insert(const unsigned char* key_begin, const std::size_t& length) + { + std::size_t bit_index = 0; + std::size_t bit = 0; + + for (std::size_t i = 0; i < salt_.size(); ++i) + { + compute_indices(hash_ap(key_begin, length, salt_[i]), bit_index, bit); + + bit_table_[bit_index / bits_per_char] |= bit_mask[bit]; + } + + ++inserted_element_count_; + } + + template + inline void insert(const T& t) + { + // Note: T must be a C++ POD type. + insert(reinterpret_cast(&t),sizeof(T)); + } + + inline void insert(const std::string& key) + { + insert(reinterpret_cast(key.data()),key.size()); + } + + inline void insert(const char* data, const std::size_t& length) + { + insert(reinterpret_cast(data),length); + } + + template + inline void insert(const InputIterator begin, const InputIterator end) + { + InputIterator itr = begin; + + while (end != itr) + { + insert(*(itr++)); + } + } + + inline virtual bool contains(const unsigned char* key_begin, const std::size_t length) const + { + std::size_t bit_index = 0; + std::size_t bit = 0; + + for (std::size_t i = 0; i < salt_.size(); ++i) + { + compute_indices(hash_ap(key_begin, length, salt_[i]), bit_index, bit); + + if ((bit_table_[bit_index / bits_per_char] & bit_mask[bit]) != bit_mask[bit]) + { + return false; + } + } + + return true; + } + + template + inline bool contains(const T& t) const + { + return contains(reinterpret_cast(&t),static_cast(sizeof(T))); + } + + inline bool contains(const std::string& key) const + { + return contains(reinterpret_cast(key.c_str()),key.size()); + } + + inline bool contains(const char* data, const std::size_t& length) const + { + return contains(reinterpret_cast(data),length); + } + + template + inline InputIterator contains_all(const InputIterator begin, const InputIterator end) const + { + InputIterator itr = begin; + + while (end != itr) + { + if (!contains(*itr)) + { + return itr; + } + + ++itr; + } + + return end; + } + + template + inline InputIterator contains_none(const InputIterator begin, const InputIterator end) const + { + InputIterator itr = begin; + + while (end != itr) + { + if (contains(*itr)) + { + return itr; + } + + ++itr; + } + + return end; + } + + inline virtual unsigned long long int size() const + { + return table_size_; + } + + inline unsigned long long int element_count() const + { + return inserted_element_count_; + } + + inline double effective_fpp() const + { + /* + Note: + The effective false positive probability is calculated using the + designated table size and hash function count in conjunction with + the current number of inserted elements - not the user defined + predicated/expected number of inserted elements. + */ + return std::pow(1.0 - std::exp(-1.0 * salt_.size() * inserted_element_count_ / size()), 1.0 * salt_.size()); + } + + inline bloom_filter& operator &= (const bloom_filter& f) + { + /* intersection */ + if ( + (salt_count_ == f.salt_count_ ) && + (table_size_ == f.table_size_ ) && + (random_seed_ == f.random_seed_) + ) + { + for (std::size_t i = 0; i < bit_table_.size(); ++i) + { + bit_table_[i] &= f.bit_table_[i]; + } + } + + return *this; + } + + inline bloom_filter& operator |= (const bloom_filter& f) + { + /* union */ + if ( + (salt_count_ == f.salt_count_ ) && + (table_size_ == f.table_size_ ) && + (random_seed_ == f.random_seed_) + ) + { + for (std::size_t i = 0; i < bit_table_.size(); ++i) + { + bit_table_[i] |= f.bit_table_[i]; + } + } + + return *this; + } + + inline bloom_filter& operator ^= (const bloom_filter& f) + { + /* difference */ + if ( + (salt_count_ == f.salt_count_ ) && + (table_size_ == f.table_size_ ) && + (random_seed_ == f.random_seed_) + ) + { + for (std::size_t i = 0; i < bit_table_.size(); ++i) + { + bit_table_[i] ^= f.bit_table_[i]; + } + } + + return *this; + } + + inline const cell_type* table() const + { + return bit_table_.data(); + } + + inline std::size_t hash_count() + { + return salt_.size(); + } + +protected: + + inline virtual void compute_indices(const bloom_type& hash, std::size_t& bit_index, std::size_t& bit) const + { + bit_index = hash % table_size_; + bit = bit_index % bits_per_char; + } + + void generate_unique_salt() + { + /* + Note: + A distinct hash function need not be implementation-wise + distinct. In the current implementation "seeding" a common + hash function with different values seems to be adequate. + */ + const unsigned int predef_salt_count = 128; + + static const bloom_type predef_salt[predef_salt_count] = + { + 0xAAAAAAAA, 0x55555555, 0x33333333, 0xCCCCCCCC, + 0x66666666, 0x99999999, 0xB5B5B5B5, 0x4B4B4B4B, + 0xAA55AA55, 0x55335533, 0x33CC33CC, 0xCC66CC66, + 0x66996699, 0x99B599B5, 0xB54BB54B, 0x4BAA4BAA, + 0xAA33AA33, 0x55CC55CC, 0x33663366, 0xCC99CC99, + 0x66B566B5, 0x994B994B, 0xB5AAB5AA, 0xAAAAAA33, + 0x555555CC, 0x33333366, 0xCCCCCC99, 0x666666B5, + 0x9999994B, 0xB5B5B5AA, 0xFFFFFFFF, 0xFFFF0000, + 0xB823D5EB, 0xC1191CDF, 0xF623AEB3, 0xDB58499F, + 0xC8D42E70, 0xB173F616, 0xA91A5967, 0xDA427D63, + 0xB1E8A2EA, 0xF6C0D155, 0x4909FEA3, 0xA68CC6A7, + 0xC395E782, 0xA26057EB, 0x0CD5DA28, 0x467C5492, + 0xF15E6982, 0x61C6FAD3, 0x9615E352, 0x6E9E355A, + 0x689B563E, 0x0C9831A8, 0x6753C18B, 0xA622689B, + 0x8CA63C47, 0x42CC2884, 0x8E89919B, 0x6EDBD7D3, + 0x15B6796C, 0x1D6FDFE4, 0x63FF9092, 0xE7401432, + 0xEFFE9412, 0xAEAEDF79, 0x9F245A31, 0x83C136FC, + 0xC3DA4A8C, 0xA5112C8C, 0x5271F491, 0x9A948DAB, + 0xCEE59A8D, 0xB5F525AB, 0x59D13217, 0x24E7C331, + 0x697C2103, 0x84B0A460, 0x86156DA9, 0xAEF2AC68, + 0x23243DA5, 0x3F649643, 0x5FA495A8, 0x67710DF8, + 0x9A6C499E, 0xDCFB0227, 0x46A43433, 0x1832B07A, + 0xC46AFF3C, 0xB9C8FFF0, 0xC9500467, 0x34431BDF, + 0xB652432B, 0xE367F12B, 0x427F4C1B, 0x224C006E, + 0x2E7E5A89, 0x96F99AA5, 0x0BEB452A, 0x2FD87C39, + 0x74B2E1FB, 0x222EFD24, 0xF357F60C, 0x440FCB1E, + 0x8BBE030F, 0x6704DC29, 0x1144D12F, 0x948B1355, + 0x6D8FD7E9, 0x1C11A014, 0xADD1592F, 0xFB3C712E, + 0xFC77642F, 0xF9C4CE8C, 0x31312FB9, 0x08B0DD79, + 0x318FA6E7, 0xC040D23D, 0xC0589AA7, 0x0CA5C075, + 0xF874B172, 0x0CF914D5, 0x784D3280, 0x4E8CFEBC, + 0xC569F575, 0xCDB2A091, 0x2CC016B4, 0x5C5F4421 + }; + + if (salt_count_ <= predef_salt_count) + { + std::copy(predef_salt, + predef_salt + salt_count_, + std::back_inserter(salt_)); + + for (std::size_t i = 0; i < salt_.size(); ++i) + { + /* + Note: + This is done to integrate the user defined random seed, + so as to allow for the generation of unique bloom filter + instances. + */ + salt_[i] = salt_[i] * salt_[(i + 3) % salt_.size()] + static_cast(random_seed_); + } + } + else + { + std::copy(predef_salt, predef_salt + predef_salt_count, std::back_inserter(salt_)); + + srand(static_cast(random_seed_)); + + while (salt_.size() < salt_count_) + { + bloom_type current_salt = static_cast(rand()) * static_cast(rand()); + + if (0 == current_salt) + continue; + + if (salt_.end() == std::find(salt_.begin(), salt_.end(), current_salt)) + { + salt_.push_back(current_salt); + } + } + } + } + + inline bloom_type hash_ap(const unsigned char* begin, std::size_t remaining_length, bloom_type hash) const + { + const unsigned char* itr = begin; + unsigned int loop = 0; + + while (remaining_length >= 8) + { + const unsigned int& i1 = *(reinterpret_cast(itr)); itr += sizeof(unsigned int); + const unsigned int& i2 = *(reinterpret_cast(itr)); itr += sizeof(unsigned int); + + hash ^= (hash << 7) ^ i1 * (hash >> 3) ^ + (~((hash << 11) + (i2 ^ (hash >> 5)))); + + remaining_length -= 8; + } + + if (remaining_length) + { + if (remaining_length >= 4) + { + const unsigned int& i = *(reinterpret_cast(itr)); + + if (loop & 0x01) + hash ^= (hash << 7) ^ i * (hash >> 3); + else + hash ^= (~((hash << 11) + (i ^ (hash >> 5)))); + + ++loop; + + remaining_length -= 4; + + itr += sizeof(unsigned int); + } + + if (remaining_length >= 2) + { + const unsigned short& i = *(reinterpret_cast(itr)); + + if (loop & 0x01) + hash ^= (hash << 7) ^ i * (hash >> 3); + else + hash ^= (~((hash << 11) + (i ^ (hash >> 5)))); + + ++loop; + + remaining_length -= 2; + + itr += sizeof(unsigned short); + } + + if (remaining_length) + { + hash += ((*itr) ^ (hash * 0xA5A5A5A5)) + loop; + } + } + + return hash; + } + + std::vector salt_; + std::vector bit_table_; + unsigned int salt_count_; + unsigned long long int table_size_; + unsigned long long int projected_element_count_; + unsigned long long int inserted_element_count_; + unsigned long long int random_seed_; + double desired_false_positive_probability_; +}; + +inline bloom_filter operator & (const bloom_filter& a, const bloom_filter& b) +{ + bloom_filter result = a; + result &= b; + return result; +} + +inline bloom_filter operator | (const bloom_filter& a, const bloom_filter& b) +{ + bloom_filter result = a; + result |= b; + return result; +} + +inline bloom_filter operator ^ (const bloom_filter& a, const bloom_filter& b) +{ + bloom_filter result = a; + result ^= b; + return result; +} + +class compressible_bloom_filter : public bloom_filter +{ +public: + + compressible_bloom_filter(const bloom_parameters& p) + : bloom_filter(p) + { + size_list.push_back(table_size_); + } + + inline unsigned long long int size() const + { + return size_list.back(); + } + + inline bool compress(const double& percentage) + { + if ( + (percentage < 0.0) || + (percentage >= 100.0) + ) + { + return false; + } + + unsigned long long int original_table_size = size_list.back(); + unsigned long long int new_table_size = static_cast((size_list.back() * (1.0 - (percentage / 100.0)))); + + new_table_size -= new_table_size % bits_per_char; + + if ( + (bits_per_char > new_table_size) || + (new_table_size >= original_table_size) + ) + { + return false; + } + + desired_false_positive_probability_ = effective_fpp(); + + const unsigned long long int new_tbl_raw_size = new_table_size / bits_per_char; + + table_type tmp(new_tbl_raw_size); + + std::copy(bit_table_.begin(), bit_table_.begin() + new_tbl_raw_size, tmp.begin()); + + typedef table_type::iterator itr_t; + + itr_t itr = bit_table_.begin() + (new_table_size / bits_per_char); + itr_t end = bit_table_.begin() + (original_table_size / bits_per_char); + itr_t itr_tmp = tmp.begin(); + + while (end != itr) + { + *(itr_tmp++) |= (*itr++); + } + + std::swap(bit_table_, tmp); + + size_list.push_back(new_table_size); + + return true; + } + +private: + + inline void compute_indices(const bloom_type& hash, std::size_t& bit_index, std::size_t& bit) const + { + bit_index = hash; + + for (std::size_t i = 0; i < size_list.size(); ++i) + { + bit_index %= size_list[i]; + } + + bit = bit_index % bits_per_char; + } + + std::vector size_list; +}; + +#endif + + +/* + Note 1: + If it can be guaranteed that bits_per_char will be of the form 2^n then + the following optimization can be used: + + bit_table_[bit_index >> n] |= bit_mask[bit_index & (bits_per_char - 1)]; + + Note 2: + For performance reasons where possible when allocating memory it should + be aligned (aligned_alloc) according to the architecture being used. +*/ \ No newline at end of file diff --git a/source/config/config.h b/source/config/config.h index 99c83903ec..f2465c1e3b 100644 --- a/source/config/config.h +++ b/source/config/config.h @@ -52,6 +52,7 @@ #include "../tools/string_utils.h" #include "ConfigManager.h" + namespace emp { using namespace std::placeholders; @@ -186,6 +187,8 @@ namespace emp { ~ConfigGroup() { ; } size_t GetSize() const { return entry_set.size(); } + std::string GetName() const {return name;} + std::string GetDesc() const {return desc;} ConfigEntry * GetEntry(size_t id) { return entry_set[id]; } ConfigEntry * GetLastEntry() { emp_assert(GetSize() > 0); return entry_set.back(); } @@ -377,6 +380,8 @@ namespace emp { for (auto & x : type_manager_map) delete x.second; } + friend class ConfigWebUI; + ConfigEntry * operator[](const std::string & name) { return var_map[name]; } auto begin() -> decltype(var_map.begin()) { return var_map.begin(); } auto end() -> decltype(var_map.end()) { return var_map.end(); } diff --git a/source/config/config_web_interface.h b/source/config/config_web_interface.h new file mode 100644 index 0000000000..4abeba669e --- /dev/null +++ b/source/config/config_web_interface.h @@ -0,0 +1,146 @@ +#ifndef EMP_CONFIG_WEB_INTERFACE_H +#define EMP_CONFIG_WEB_INTERFACE_H + +#include "config.h" +#include "../web/Div.h" +#include "../web/Element.h" +#include "../web/Input.h" + +#include +#include +#include "../tools/set_utils.h" +#include "../tools/string_utils.h" + +namespace emp { + + class ConfigWebUI { + private: + inline static std::set numeric_types = {"int", "double", "float", "uint32_t", "uint64_t", "size_t"}; + Config & config; + web::Div settings_div; + std::set exclude; + std::map group_divs; + std::map input_map; + std::function on_change_fun = [](const std::string & val){;}; + std::function format_label_fun = [](std::string name){ + emp::vector sliced = slice(name, '_'); + return to_titlecase(join(sliced, " ")); + }; + public: + ConfigWebUI(Config & c, const std::string & div_name = "settings_div") + : config(c), settings_div(div_name) {;} + + void SetOnChangeFun(std::function fun) {on_change_fun = fun;} + + template + void SetDefaultRangeFloatingPoint(web::Input & input, T val) { + if (val > 0 && val < 1) { + // This is a common range for numbers to be in + input.Min(0); + if (val > .1) { + input.Max(1); + } else { + input.Max(val * 100); + } + input.Step(val/10.0); + } else if (val > 0) { + // Assume this is a positive number + input.Min(0); + input.Max(val * 10); + input.Step(val/10.0); + } else if (val < 0) { + input.Min(val * 10); // since val is negative + input.Max(val * -10); + input.Step(val/-10.0); // A negative step would be confusing + } + + // Otherwise val is 0 and we have nothing to go on + } + + void SetDefaultRangeFixedPoint(web::Input & input, int val) { + // Default step is 1, which should be fine for fixed point + + if (val > 0) { + // Assume this is a positive number + input.Min(0); + input.Max(val * 10); + } else if (val < 0) { + input.Min(val * 10); // since val is negative + input.Max(val * -10); + } + + // Otherwise val is 0 and we have nothing to go on + } + + void ExcludeConfig(std::string setting) { + exclude.insert(setting); + } + + void Setup(const std::string & id_prefix = "settings_") { + + for (auto group : config.group_set) { + // std::cout << "GROUP: " << group->GetName() << std::endl; + std::string group_name = group->GetName(); + group_divs[group_name] = web::Div(id_prefix + group_name); + group_divs[group_name] << "

" << group->GetDesc() << "

"; + for (size_t i = 0; i < group->GetSize(); i++) { + // std::cout << group->GetEntry(i)->GetType() << std::endl; + std::string name = group->GetEntry(i)->GetName(); + if (Has(exclude, name)) { + continue; + } + std::string type = group->GetEntry(i)->GetType(); + std::string value = group->GetEntry(i)->GetValue(); + + if (Has(numeric_types, type)) { + input_map[name] = emp::web::Input( + [this, name](std::string val){ + std::cout << name << " " << val << " " <(value)); + } else if (type == "float") { + SetDefaultRangeFloatingPoint(input_map[name], emp::from_string(value)); + } else { + // TODO: Correctly handle all types (although I'm not sure it actually matters?) + SetDefaultRangeFixedPoint(input_map[name], emp::from_string(value)); + } + + } else if (type == "bool") { + input_map[name] = emp::web::Input( + [this, name](std::string val){config.Set(name, val); + on_change_fun(val);}, + "checkbox", format_label_fun(name), name + "_input_checkbox" + ); + } else { + input_map[name] = emp::web::Input( + [this, name](std::string val){config.Set(name, val); + on_change_fun(val);}, + "text", format_label_fun(name), name + "_input_textbox" + ); + + } + + input_map[name].Value(value); + + group_divs[group_name] << web::Element("form") << input_map[name]; + + } + settings_div << group_divs[group_name]; + } + + } + + web::Div & GetDiv() {return settings_div;} + + }; + +} + +#endif \ No newline at end of file diff --git a/source/data/DataFile.h b/source/data/DataFile.h index b3e1fcbdac..88f6176075 100644 --- a/source/data/DataFile.h +++ b/source/data/DataFile.h @@ -34,10 +34,12 @@ namespace emp { protected: using fun_t = void(std::ostream &); using time_fun_t = std::function; + using pre_fun_t = std::function; std::string filename; ///< Name of the file that we are printing to (if one exists) std::ostream * os; ///< Stream to print to. FunctionSet funs; ///< Set of functions to call, one per column in the file. + FunctionSet pre_funs; ///< Set of functions to call before calculating data. emp::vector keys; ///< Keywords associated with each column. emp::vector descs; ///< Full description for each column. time_fun_t timing_fun; ///< Function to determine updates to print on (default: all) @@ -55,12 +57,12 @@ namespace emp { #else new std::ofstream(in_filename) #endif - ), funs(), keys(), descs() + ), funs(), pre_funs(), keys(), descs() , timing_fun([](size_t){return true;}) , line_begin(b), line_spacer(s), line_end(e) { ; } DataFile(std::ostream & in_os, const std::string & b="", const std::string & s=",", const std::string & e="\n") - : filename(), os(&in_os), funs(), keys(), descs(), timing_fun([](size_t){return true;}) + : filename(), os(&in_os), funs(), pre_funs(), keys(), descs(), timing_fun([](size_t){return true;}) , line_begin(b), line_spacer(s), line_end(e) { ; } DataFile(const DataFile &) = default; @@ -140,6 +142,7 @@ namespace emp { /// Run all of the functions and print the results as a new line in the file virtual void Update() { + pre_funs.Run(); *os << line_begin; for (size_t i = 0; i < funs.size(); i++) { if (i > 0) *os << line_spacer; @@ -154,6 +157,10 @@ namespace emp { virtual void Update(size_t update) { if (timing_fun(update)) Update(); } + void AddPreFun(pre_fun_t fun) { + pre_funs.Add(fun); + } + /// If a function takes an ostream, pass in the correct one. /// Generic function for adding a column to the DataFile. In practice, you probably /// want to call one of the more specific ones. @@ -331,12 +338,17 @@ namespace emp { /// Requires that @param node have the data::Stats or data::FullStats modifier. /// @param key and @param desc will have the name of the stat appended to the beginning. /// Note: excludes standard deviation, because it is easily calculated from variance + /// Note: Setting @param pull and/or @param reset to true only pulls on first statistic + /// calculated and only resets on the last. Otherwise there would be a risk of data loss or + /// at least massive replication of computational effort. Even still, think carefully + /// before setting either of these to true when you're drawing varied information from the + /// same node. template void AddStats(DataNode & node, const std::string & key="", const std::string & desc="", const bool & reset=false, const bool & pull=false) { - AddMean(node, "mean_" + key, "mean of " + desc, reset, pull); - AddMin(node, "min_" + key, "min of " + desc, reset, pull); - AddMax(node, "max_" + key, "max of " + desc, reset, pull); - AddVariance(node, "variance_" + key, "variance of " + desc, reset, pull); + AddMean(node, "mean_" + key, "mean of " + desc, false, pull); + AddMin(node, "min_" + key, "min of " + desc); + AddMax(node, "max_" + key, "max of " + desc); + AddVariance(node, "variance_" + key, "variance of " + desc, reset); } @@ -345,11 +357,16 @@ namespace emp { /// Requires that @param node have the data::Stats or data::FullStats modifier. /// @param key and @param desc will have the name of the stat appended to the beginning. /// Note: excludes standard deviation, because it is easily calculated from variance + /// Note: Setting @param pull and/or @param reset to true only pulls on first statistic + /// calculated and only resets on the last. Otherwise there would be a risk of data loss or + /// at least massive replication of computational effort. Even still, think carefully + /// before setting either of these to true when you're drawing varied information from the + /// same node. template void AddAllStats(DataNode & node, const std::string & key="", const std::string & desc="", const bool & reset=false, const bool & pull=false) { - AddStats(node, key, desc, reset, pull); - AddSkew(node, "skew_" + key, "skew of " + desc, reset, pull); - AddKurtosis(node, "kurtosis_" + key, "kurtosis of " + desc, reset, pull); + AddStats(node, key, desc, false, pull); + AddSkew(node, "skew_" + key, "skew of " + desc); + AddKurtosis(node, "kurtosis_" + key, "kurtosis of " + desc, reset); } /// Add a function that always pulls the count of the @param bin_id 'th bin of the histogram @@ -368,6 +385,28 @@ namespace emp { return Add(in_fun, key, desc); } + /// Add a set of functions that always pull the count of each bin of the histogram + /// from @param node. Requires that @param node have the data::Histogram modifier. + /// If @param reset is set true, we will call Reset on that DataNode after pulling the + /// current value from the node + /// Note: Setting @param pull and/or @param reset to true only pulls on first statistic + /// calculated and only resets on the last. Otherwise there would be a risk of data loss or + /// at least massive replication of computational effort. Even still, think carefully + /// before setting either of these to true when you're drawing varied information from the + /// same node. + template + void AddAllHistBins(DataNode & node, const std::string & key="", const std::string & desc="", const bool & reset=false, const bool & pull=false) { + bool actual_reset = false; + bool actual_pull = pull; + for (size_t i = 0; i < node.GetHistCounts().size(); i++) { + if (i == node.GetHistCounts().size() - 1) { + actual_reset = reset; + } + AddHistBin(node, i, key + "_bin_" + emp::to_string(i), desc + " bin " + emp::to_string(i), actual_reset, actual_pull); + actual_pull = false; + } + } + /// Add a function that always pulls the inferiority (mean divided by max) from the DataNode @param node. /// Requires that @param node have the data::Range or data::FullRange modifier. /// If @param reset is set true, we will call Reset on that DataNode after pulling the diff --git a/source/data/DataManager.h b/source/data/DataManager.h index ef8a4618c5..5fb20e2779 100644 --- a/source/data/DataManager.h +++ b/source/data/DataManager.h @@ -86,7 +86,7 @@ namespace emp { /// Returns a reference to the node named @param name. /// Throws an error if there is no node with that name in this manager node_t & Get(const std::string & name) { - emp_assert(Has(node_map, name), name); + emp_assert(Has(node_map, name), name, emp::to_string(Keys(node_map))); return *(node_map[name]); } @@ -104,7 +104,7 @@ namespace emp { * my_data_manager.AddData("my_node_name", 1, 2, 3, 4, 5); */ template void AddData(const std::string & name, Ts... extra) { - emp_assert(Has(node_map, name), name); + emp_assert(Has(node_map, name), name, emp::to_string(Keys(node_map))); node_map[name]->Add(extra...); } diff --git a/source/data/DataNode.h b/source/data/DataNode.h index 8b42e1bed9..a20fd6dda7 100644 --- a/source/data/DataNode.h +++ b/source/data/DataNode.h @@ -547,6 +547,8 @@ namespace emp { VAL_TYPE width; ///< How wide is the overall histogram? IndexMap bins; ///< Map of values to which bin they fall in. emp::vector counts; ///< Counts in each bin. + size_t overflow; + size_t underflow; using this_t = DataNodeModule; using parent_t = DataNodeModule; @@ -555,7 +557,7 @@ namespace emp { using base_t::val_count; public: - DataNodeModule() : offset(0.0), width(0), bins(), counts() { ; } + DataNodeModule() : offset(0.0), width(0), bins(), counts(), overflow(0), underflow(0) { ; } /// Returns the minimum value this histogram is capable of containing /// (i.e. the minimum value for the first bin) @@ -574,6 +576,14 @@ namespace emp { /// Return a vector containing the count of items in each bin of the histogram const emp::vector & GetHistCounts() const { return counts; } + /// Return the count of numbers added to this histogram that were above the + /// upper bound on the histogram + int GetOverflow() const {return overflow;} + + /// Return the count of numbers added to this histogram that were belowed the + /// allowed lower bound + int GetUnderflow() const {return underflow;} + /// Return a vector containing the lowest value allowed in each bin. emp::vector GetBinMins() const { emp::vector bin_mins(counts.size()); @@ -603,9 +613,15 @@ namespace emp { /// Add @param val to the DataNode void AddDatum(const VAL_TYPE & val) { - size_t bin_id = bins.Index((double) (val - offset)); - // size_t bin_id = counts.size() * ((double) (val - offset)) / (double) width; - counts[bin_id]++; + if ((double)(val - offset) >= width) { + overflow++; + } else if (val < offset) { + underflow++; + } else { + size_t bin_id = bins.Index((double) (val - offset)); + // size_t bin_id = counts.size() * ((double) (val - offset)) / (double) width; + counts[bin_id]++; + } parent_t::AddDatum(val); } diff --git a/source/hardware/AvidaGP.h b/source/hardware/AvidaGP.h index 394854449d..f083e7ebff 100644 --- a/source/hardware/AvidaGP.h +++ b/source/hardware/AvidaGP.h @@ -255,6 +255,10 @@ namespace emp { return genome < other.genome; } + bool operator!=(const this_t & other) const { + return genome != other.genome; + } + /// Reset the entire CPU to a starting state, without a genome. void Reset() { genome.sequence.resize(0); // Clear out genome @@ -294,8 +298,10 @@ namespace emp { // Accessors Ptr GetInstLib() const { return genome.inst_lib; } inst_t GetInst(size_t pos) const { return genome.sequence[pos]; } + inst_t& operator[](size_t pos) {return genome.sequence[pos]; } // Alias for compatability with tools const genome_t & GetGenome() const { return genome; } const size_t GetSize() const { return genome.sequence.size(); } + const size_t size() const { return GetSize(); } // Alias for compatability with tools double GetReg(size_t id) const { return regs[id]; } double GetInput(int id) const { return Find(inputs, id, 0.0); } const std::unordered_map & GetInputs() const { return inputs; } @@ -567,5 +573,13 @@ namespace emp { }; } +namespace std { + + /// operator<< to work with ostream (must be in std to work) + inline std::ostream & operator<<(std::ostream & out, const emp::AvidaGP & org) { + org.PrintGenome(out); + return out; + } +} #endif diff --git a/source/hardware/BitSorter.h b/source/hardware/BitSorter.h index 16972b6fc9..bc9dfb61b6 100644 --- a/source/hardware/BitSorter.h +++ b/source/hardware/BitSorter.h @@ -19,8 +19,10 @@ namespace emp { class BitSorter { - private: - using bits_t = uint32_t; ///< Type used to represent pairs if posisions as bit masks. + public: + using bits_t = uint32_t; ///< Type used to represent pairs if posisions as bit masks. + protected: + emp::vector compare_set; ///< Comparators, in order (pairs of 1's in bitstring) public: @@ -32,8 +34,36 @@ namespace emp { BitSorter & operator=(const BitSorter &) = default; BitSorter & operator=(BitSorter &&) = default; + bool operator!=(const BitSorter & other) const { + return compare_set != other.compare_set; + } + + bool operator<(const BitSorter & other) const { + return compare_set < other.compare_set; + } + + /// How many comparators are in this sorting network. size_t GetSize() const { return compare_set.size(); } + size_t size() const { return GetSize(); } + + std::pair GetComparator(size_t idx) { + emp_assert(idx < compare_set.size(), idx, compare_set.size()); + bits_t curr = compare_set[idx]; + bits_t pos1 = pop_bit(curr); + bits_t pos2 = find_bit(curr); + return std::make_pair(pos1, pos2); + } + bits_t & operator[](size_t idx) { + return compare_set[idx]; + } + bits_t GetBits(size_t idx) { + return compare_set[idx]; + } + + void Clear() { + compare_set.clear(); + } /// If this network is compressed as far as possible, what will the max depth of each position be? void CalcDepth(size_t num_bits, emp::vector & depth_vals) const { @@ -153,4 +183,13 @@ namespace emp { } +namespace std { + + /// operator<< to work with ostream (must be in std to work) + inline std::ostream & operator<<(std::ostream & out, const emp::BitSorter & bitsort) { + out << bitsort.AsString(); + return out; + } +} + #endif \ No newline at end of file diff --git a/source/tools/File.h b/source/tools/File.h index ae2eb4bf84..de75a371d2 100644 --- a/source/tools/File.h +++ b/source/tools/File.h @@ -62,6 +62,9 @@ namespace emp { /// Compatibility with size() size_t size() const { return lines.size(); } + /// Return entire text of the file + emp::vector GetAllLines() {return lines;} + /// Index into a specific line in this file. std::string & operator[](size_t pos) { return lines[pos]; } diff --git a/source/tools/Graph.h b/source/tools/Graph.h index 739c914233..16e5a51176 100644 --- a/source/tools/Graph.h +++ b/source/tools/Graph.h @@ -27,9 +27,10 @@ namespace emp { class Node { private: BitVector edge_set; /// What other node IDs is this one connected to? + std::string label; public: - Node(size_t num_nodes) : edge_set(num_nodes) { ; } - Node(const Node & in_node) : edge_set(in_node.edge_set) { ; } + Node(size_t num_nodes) : edge_set(num_nodes), label("") { ; } + Node(const Node & in_node) : edge_set(in_node.edge_set), label(in_node.label) { ; } ~Node() { ; } /// Set this node to have the same connections as another node. @@ -64,6 +65,15 @@ namespace emp { /// Identify how many other nodes from a provided set (a BitVector) this one is connected to. size_t GetMaskedDegree(const BitVector & mask) const { return (mask & edge_set).CountOnes(); } + + void SetLabel(std::string lab) { + label = lab; + } + + std::string GetLabel() { + return label; + } + }; protected: @@ -90,6 +100,9 @@ namespace emp { return edge_count; } + Node GetNode(int i) {return nodes[i];} + emp::vector GetNodes(){return nodes;} + /// Change the number of vertices in this graph. void Resize(size_t new_size) { nodes.resize(new_size, new_size); @@ -106,17 +119,44 @@ namespace emp { } /// Get the degree of a specified node. + /// For directed graphs, this is the out-degree size_t GetDegree(size_t id) const { emp_assert(id < nodes.size()); return nodes[id].GetDegree(); } + /// Get the in-degree (number of incoming edges) + /// of the node @param id. + /// This should only be used for directed graphs (for + /// undirected graphs, GetDegree() is equivalent and faster) + size_t GetInDegree(size_t id) const { + size_t count = 0; + for (auto & node : nodes) { + // Node is allowed to to have edge to itself so it's + // okay that we don't exclude it + if (node.HasEdge(id)) { + count++; + } + } + return count; + } + /// Get how many of a set of nodes that a specified node is connected to. size_t GetMaskedDegree(size_t id, const BitVector & mask) const { emp_assert(id < nodes.size()); return nodes[id].GetMaskedDegree(mask); } + /// Set label of node @param id + void SetLabel(size_t id, std::string lab) { + nodes[id].SetLabel(lab); + } + + /// Get label of node @param id + std::string GetLabel(size_t id) { + return nodes[id].GetLabel(); + } + /// Determine if a specific edge is included in this graph. bool HasEdge(size_t from, size_t to) const { emp_assert(from < nodes.size() && to < nodes.size()); @@ -289,6 +329,8 @@ namespace emp { } } + emp::vector > GetWeights(){return weights;} + }; } diff --git a/source/tools/Random.h b/source/tools/Random.h index 6291799239..7fbe1f82c7 100644 --- a/source/tools/Random.h +++ b/source/tools/Random.h @@ -26,7 +26,7 @@ namespace emp { /// Middle Square Weyl Sequence: A versatile and non-patterned pseudo-random-number /// generator. - /// Based on: https://en.wikipedia.org/wiki/Middle-square_method + /// Based on: https://en.wikipedia.org/wiki/Middle-square_method class Random { protected: @@ -248,6 +248,19 @@ namespace emp { return k; } + inline uint32_t GetRandGeometric(double p){ + emp_assert(p >= 0 && p <= 1, "Pobabilities must be between 0 and 1"); + // TODO: When we have warnings, add one for passing a really small number to + // this function. Alternatively, make this function not ludicrously slow with small numbers. + // Looks like return floor(ln(GetDouble())/ln(1-p)) might be sufficient? + if (p == 0) { + return std::numeric_limits::infinity(); + } + uint32_t result = 1; + while (!P(p)) { result++;} + return result; + } + }; diff --git a/source/tools/distances.h b/source/tools/distances.h new file mode 100644 index 0000000000..18a7c2226d --- /dev/null +++ b/source/tools/distances.h @@ -0,0 +1,53 @@ +/** + * @note This file is part of Empirical, https://github.com/devosoft/Empirical + * @copyright Copyright (C) Michigan State University, MIT Software license; see doc/LICENSE.md + * @date 2017-2018 + * + * @file distances.h + * @brief Library of commonly used distance functions + * @note Status: BETA + */ + +#ifndef EMP_DISTANCES_H +#define EMP_DISTANCES_H + +#include "math.h" +#include "../meta/type_traits.h" + +namespace emp { + + /// Calculate Euclidean distance between two containers + template + typename std::enable_if::value, double>::type + EuclideanDistance(C & p1, C & p2) { + emp_assert(p1.size() == p2.size() + && "Cannot calculate euclidean distance between two containers of different lengths."); + + double dist = 0; + for (size_t i = 0; i < p1.size(); ++i) { + dist += pow(p1[i] - p2[i], 2); + } + + return sqrt(dist); + } + + /// Calculate Euclidean distance between two containers of pointers (de-referencing the pointers) + template + typename std::enable_if::value, double>::type + EuclideanDistance(C & p1, C & p2) { + + emp_assert(p1.size() == p2.size() + && "Cannot calculate euclidean distance between two containers of different lengths."); + + double dist = 0; + for (size_t i = 0; i < p1.size(); ++i) { + dist += emp::Pow(*p1[i] - *p2[i], 2); + } + + return sqrt(dist); + } + + +} + +#endif \ No newline at end of file diff --git a/source/tools/map_utils.h b/source/tools/map_utils.h index b6287e9741..f093fe7ea2 100644 --- a/source/tools/map_utils.h +++ b/source/tools/map_utils.h @@ -14,6 +14,8 @@ #include #include +#include "../base/vector.h" + namespace emp { /// Take any map type, and run find to determine if a key is present. @@ -23,6 +25,19 @@ namespace emp { } + template + inline auto Keys( const MAP_T & in_map) -> emp::vectorfirst)>::type> { + using KEY_T = typename std::remove_constfirst)>::type; + emp::vector keys; + for (auto it : in_map) { + keys.push_back(it.first); + } + + return keys; + + } + + /// Take any map, run find() member function, and return the result found /// (or default value if no results found). template diff --git a/source/tools/math.h b/source/tools/math.h index 3d9c2b28b3..a61dee4197 100644 --- a/source/tools/math.h +++ b/source/tools/math.h @@ -315,6 +315,15 @@ namespace emp { return *max_found; } + inline constexpr int Factorial(int i) { + int result = 1; + while (i > 0) { + result *= i; + i--; + } + return result; + } + } #endif diff --git a/source/tools/set_utils.h b/source/tools/set_utils.h index 7f4bc7f162..0699f23fd2 100644 --- a/source/tools/set_utils.h +++ b/source/tools/set_utils.h @@ -58,7 +58,7 @@ namespace emp { /// Compute the set difference of @param s1 and @param s2 (elements that are in S1 but no S2) template - std::set difference(emp::vector s1, emp::vector & s2) { + std::set difference(emp::vector s1, emp::vector s2) { // Based on PierreBdR's answer to https://stackoverflow.com/questions/283977/c-stl-set-difference std::sort(s1.begin(), s1.end()); // set_difference expects sorted things std::sort(s2.begin(), s2.end()); // set_difference expects sorted things @@ -93,7 +93,7 @@ namespace emp { /// Compute the set intersection of @param s1 and @param s2 (elements that are in both S1 and S2) template - std::set intersection(std::set s1, std::set s2) { + std::set intersection(std::set & s1, std::set & s2) { std::set result; std::set_intersection(s1.begin(), s1.end(), s2.begin(), s2.end(), std::inserter(result, result.end())); @@ -114,7 +114,7 @@ namespace emp { /// Compute the set intersection of @param s1 and @param s2 (elements that are in both S1 and S2) template - std::set intersection(std::set s1, emp::vector s2) { + std::set intersection(std::set & s1, emp::vector s2) { std::sort(s2.begin(), s2.end()); // set_intersection expects sorted things std::set result; @@ -125,7 +125,7 @@ namespace emp { /// Compute the set intersection of @param s1 and @param s2 (elements that are in both S1 and S2) template - std::set intersection(emp::vector s1, std::set s2) { + std::set intersection(emp::vector s1, std::set & s2) { std::sort(s1.begin(), s1.end()); // set_intersection expects sorted things std::set result; @@ -136,7 +136,7 @@ namespace emp { /// Compute the set union of @param s1 and @param s2 (elements that are in either S1 or S2) template - std::set set_union(std::set s1, std::set s2) { + std::set set_union(std::set & s1, std::set & s2) { std::set result; std::set_union(s1.begin(), s1.end(), s2.begin(), s2.end(), std::inserter(result, result.end())); @@ -157,7 +157,7 @@ namespace emp { /// Compute the set union of @param s1 and @param s2 (elements that are in either S1 or S2) template - std::set set_union(std::set s1, emp::vector s2) { + std::set set_union(std::set & s1, emp::vector s2) { std::sort(s2.begin(), s2.end()); // set_union expects sorted things std::set result; @@ -168,7 +168,7 @@ namespace emp { /// Compute the set union of @param s1 and @param s2 (elements that are in either S1 or S2) template - std::set set_union(emp::vector s1, std::set s2) { + std::set set_union(emp::vector s1, std::set & s2) { std::sort(s1.begin(), s1.end()); // set_union expects sorted things std::set result; @@ -179,7 +179,7 @@ namespace emp { /// Compute the set symmetric_difference of @param s1 and @param s2 (elements that are in either S1 or S2 but not both) template - std::set symmetric_difference(std::set s1, std::set s2) { + std::set symmetric_difference(std::set & s1, std::set & s2) { std::set result; std::set_symmetric_difference(s1.begin(), s1.end(), s2.begin(), s2.end(), std::inserter(result, result.end())); @@ -200,7 +200,7 @@ namespace emp { /// Compute the set symmetric_difference of @param s1 and @param s2 (elements that are in either S1 or S2 but not both) template - std::set symmetric_difference(std::set s1, emp::vector s2) { + std::set symmetric_difference(std::set & s1, emp::vector s2) { std::sort(s2.begin(), s2.end()); // set_symmetric_difference expects sorted things std::set result; @@ -211,7 +211,7 @@ namespace emp { /// Compute the set symmetric_difference of @param s1 and @param s2 (elements that are in either S1 or S2 but not both) template - std::set symmetric_difference(emp::vector s1, std::set s2) { + std::set symmetric_difference(emp::vector s1, std::set & s2) { std::sort(s1.begin(), s1.end()); // set_symmetric_difference expects sorted things std::set result; diff --git a/source/tools/spatial_stats.h b/source/tools/spatial_stats.h new file mode 100644 index 0000000000..9ae5830767 --- /dev/null +++ b/source/tools/spatial_stats.h @@ -0,0 +1,121 @@ +/** + * @note This file is part of Empirical, https://github.com/devosoft/Empirical + * @copyright Copyright (C) Michigan State University, MIT Software license; see doc/LICENSE.md + * @date 2016-2019 + * + * @file spatial_stats.h + * @brief Functions for calculating various spatial statistics. + * @note Status: BETA + */ + + +#ifndef EMP_SPATIAL_STATS_H +#define EMP_SPATIAL_STATS_H + +#include "../base/vector.h" +#include "../Evolve/World.h" +#include "stats.h" + +#include + +namespace emp { + + + // Point on grid stats + + template + double GridPointDensity(emp::World & w, size_t xid, size_t yid, int neighborhood_size = 2) { + double density = 0; + int x_modifier = 0; + int y_modifier = 0; + + if ((int)xid < neighborhood_size) { + x_modifier = neighborhood_size - xid; + } else if ((int)xid + neighborhood_size >= (int)w.GetWidth()) { + x_modifier = xid + neighborhood_size - w.GetWidth() + 1; + } + + if ((int)yid < neighborhood_size) { + y_modifier = neighborhood_size - yid; + } else if ((int)yid + neighborhood_size >= (int)w.GetHeight()) { + y_modifier = yid + neighborhood_size - w.GetHeight() + 1; + } + + double denominator = (neighborhood_size*2 + 1 - x_modifier) * (neighborhood_size * 2 + 1 - y_modifier); + + for (size_t x = std::max((int)xid - neighborhood_size, 0); (int)x <= std::min((int)xid + neighborhood_size, (int)w.GetWidth()-1); x++) { + for (size_t y = std::max((int)yid - neighborhood_size, 0); (int)y <= std::min((int)yid + neighborhood_size, (int)w.GetHeight()-1); y++) { + if (w.IsOccupied(x+y*w.GetWidth())) { + density++; + } + } + } + + return density/denominator; + } + + template + double GridPointShannonEntropy(emp::World & w, size_t xid, size_t yid, int neighborhood_size = 2) { + // double density = 0; + // int x_modifier = 0; + // int y_modifier = 0; + // if ((int)xid < neighborhood_size) { + // x_modifier = neighborhood_size - xid; + // } else if (xid + neighborhood_size >= w.GetWidth()) { + // x_modifier = xid + neighborhood_size - w.GetWidth() + 1; + // } + + // if ((int)yid < neighborhood_size) { + // y_modifier = neighborhood_size - yid; + // } else if (yid + neighborhood_size >= w.GetHeight()) { + // y_modifier = yid + neighborhood_size - w.GetHeight() + 1; + // } + + // double denominator = (neighborhood_size*2 + 1 - x_modifier) * (neighborhood_size * 2 + 1 - y_modifier); + + emp::vector orgs; + + for (size_t x = std::max((int)xid - neighborhood_size, 0); x <= std::min((int)xid + neighborhood_size, (int)w.GetWidth()-1); x++) { + for (size_t y = std::max((int)yid - neighborhood_size, 0); y <= std::min((int)yid + neighborhood_size, (int)w.GetHeight()-1); y++) { + if (w.IsOccupied(x+y*w.GetWidth())) { + orgs.push_back(w.GetOrg(x,y)); + } + } + } + + return ShannonEntropy(orgs); + } + + // Grid stats + template + emp::vector> GridDensity(emp::World & w, int neighborhood_size = 2) { + emp::vector> densities; + emp_assert(w.GetAttribute("PopStruct") == "Grid", "Trying to calculate grid statistics on non-grid world. Did you forget to call SetPopStruct_Grid() on this world?"); + + densities.resize(w.GetHeight()); + for (size_t y = 0; y < w.GetHeight(); y++) { + densities[y].resize(w.GetWidth()); + for (size_t x = 0; x < w.GetWidth(); x++) { + densities[y][x] = GridPointDensity(w, x, y, neighborhood_size); + } + } + return densities; + } + + template + emp::vector> GridShannonEntropy(emp::World & w, int neighborhood_size = 2) { + emp::vector> diversities; + emp_assert(w.GetAttribute("PopStruct") == "Grid", "Trying to calculate grid statistics on non-grid world. Did you forget to call SetPopStruct_Grid() on this world?"); + + diversities.resize(w.GetHeight()); + for (size_t y = 0; y < w.GetHeight(); y++) { + diversities[y].resize(w.GetWidth()); + for (size_t x = 0; x < w.GetWidth(); x++) { + diversities[y][x] = GridPointShannonEntropy(w, x, y, neighborhood_size); + } + } + return diversities; + } +} + +#endif \ No newline at end of file diff --git a/source/tools/stats.h b/source/tools/stats.h index e865142fd5..b5d4870230 100644 --- a/source/tools/stats.h +++ b/source/tools/stats.h @@ -23,6 +23,7 @@ #include "../base/map.h" #include "../base/vector.h" #include "../meta/type_traits.h" +#include "vector_utils.h" #include "math.h" namespace emp { @@ -72,7 +73,7 @@ namespace emp { // Shannon entropy calculation double result = 0; - for (auto element : counts) { + for (auto & element : counts) { double p = double(element.second)/elements.size(); result += p * Log2(p); } @@ -90,6 +91,17 @@ namespace emp { return (double)Sum(elements)/ (double) elements.size(); } + template + emp::sfinae_decoy + Median(C elements) { + emp::Sort(elements); + if (elements.size() % 2 == 1) { + return elements[elements.size() / 2]; + } else { + return (elements[elements.size() / 2 - 1] + elements[elements.size() / 2])/2.0; + } + } + /// Calculate variance of the members of the container passed /// Only works on containers with a scalar member type diff --git a/source/tools/string_utils.h b/source/tools/string_utils.h index 36c96147b1..5650cec2dc 100644 --- a/source/tools/string_utils.h +++ b/source/tools/string_utils.h @@ -19,6 +19,7 @@ #include #include #include +#include #include #include @@ -287,7 +288,28 @@ namespace emp { return value; } - /// Convert an integer to a roman numeral string. + /// Make first letter of each word upper case + static inline std::string to_titlecase(std::string value) { + constexpr int char_shift = 'a' - 'A'; + bool next_upper = true; + for (size_t i = 0; i < value.size(); i++) { + if (next_upper && value[i] >= 'a' && value[i] <= 'z') { + value[i] = (char) (value[i] - char_shift); + } else if (!next_upper && value[i] >= 'A' && value[i] <= 'Z') { + value[i] = (char) (value[i] + char_shift); + } + + if (value[i] == ' ') { + next_upper = true; + } else { + next_upper = false; + } + } + return value; + } + + + // Convert an integer to a roman numeral string. static inline std::string to_roman_numeral(int val, const std::string & prefix="") { std::string ret_string(prefix); if (val < 0) ret_string += to_roman_numeral(-val, "-"); @@ -865,6 +887,24 @@ namespace emp { ss >> out_val; return out_val; } + + template + inline std::string join(const emp::vector & v, std::string join_str) { + + if (v.size() == 0) { + return ""; + } else if (v.size() == 1) { + return to_string(v[0]); + } else { + std::stringstream res; + res << v[0]; + for (size_t i = 1; i < v.size(); i++) { + res << join_str; + res << to_string(v[i]); + } + return res.str(); + } + } // -------- Functions that operate on VECTORS of strings -------- diff --git a/source/tools/vector_utils.h b/source/tools/vector_utils.h index 17dac66153..c0f5d00c73 100644 --- a/source/tools/vector_utils.h +++ b/source/tools/vector_utils.h @@ -14,6 +14,8 @@ #ifndef EMP_VECTOR_UTILS_H #define EMP_VECTOR_UTILS_H +#include +#include #include #include #include @@ -72,13 +74,19 @@ namespace emp { return true; } - /// Return whether a value exists in a vector.s + /// Return whether a value exists in a vector template bool Has(const emp::vector & v, const T & val) { return FindValue(v, val) >= 0; } - /// Print the contects of a vector. + /// Return number of times a value occurs in a vector + template + int Count(const emp::vector & vec, const T & val) { + return std::count (vec.begin(), vec.end(), val); + } + + /// Print the contents of a vector. template void Print(const emp::vector & v, std::ostream & os=std::cout, const std::string & spacer=" ") { for (size_t id = 0; id < v.size(); id++) { @@ -220,6 +228,25 @@ namespace emp { return out_vv; } + + /// Returns a vector containing the numbers from @param N1 to @param N2 + // from https://stackoverflow.com/questions/13152252/is-there-a-compact-equivalent-to-python-range-in-c-stl + template + emp::vector NRange(T N1, T N2) { + emp::vector numbers(N2-N1); + std::iota(numbers.begin(), numbers.end(), N1); + return numbers; + } + + /// Return a new vector containing the same elements as @param v + /// with any duplicate elements removed. + /// Not guarunteed to preserve order + template + emp::vector RemoveDuplicates(const emp::vector & v) { + std::set temp_set(v.begin(), v.end()); + emp::vector new_vec(temp_set.begin(), temp_set.end()); + return new_vec; + } /// Tree manipulation in vectors. constexpr size_t tree_left(size_t id) { return id*2+1; } diff --git a/source/web/Input.h b/source/web/Input.h index 993d918386..e3b5f1c5f0 100644 --- a/source/web/Input.h +++ b/source/web/Input.h @@ -53,7 +53,7 @@ namespace web { std::string step = ""; std::string curr_val =""; - + bool show_value = false; bool autofocus; std::function callback; @@ -74,7 +74,15 @@ namespace web { std::string GetTypeName() const override { return "InputInfo"; } virtual void GetHTML(std::stringstream & HTML) override { + + // CSS from https://stackoverflow.com/questions/46695616/align-range-slider-and-label + HTML.str(""); // Clear the current text. + // HTML << "
"; // Needs to be part of form for label + output to work + if (label != "") { // Add label, if one exists + HTML << ""; + } HTML << "" << label << ""; // Close and label the Input. + HTML << ">" << ""; // Close the Input. + if (show_value) { + HTML << ""; // Add output to show value of slider + } + // HTML << "
"; + } + + virtual void TriggerJS() override { + + if (show_value) { + // Inspired by https://codepen.io/chriscoyier/pen/imdrE + EM_ASM_ARGS({ + + function modifyOffset() { + var el; + var newPlace; + var offset; + var siblings; + var k; + var width = this.offsetWidth; + var newPoint = (this.value - this.getAttribute("min")) / (this.getAttribute("max") - this.getAttribute("min")); + offset = -1; + if (newPoint < 0) { newPlace = 0; } + else if (newPoint > 1) { newPlace = width; } + else { newPlace = width * newPoint + offset; offset -= newPoint;} + siblings = this.parentNode.childNodes; + for (var i = 0; i < siblings.length; i++) { + sibling = siblings[i]; + if (sibling.id == this.id) { k = true; } + if ((k == true) && (sibling.nodeName == "OUTPUT")) { + outputTag = sibling; + } + } + outputTag.innerHTML = this.value; + } + + function modifyInputs() { + + var input_el = document.getElementById(UTF8ToString($0)); + input_el.addEventListener("input", modifyOffset); + // the following taken from http://stackoverflow.com/questions/2856513/trigger-onchange-event-manually + if ("fireEvent" in input_el) { + input_el.fireEvent("oninput"); + } else { + var evt = document.createEvent("HTMLEvents"); + evt.initEvent("input", false, true); + input_el.dispatchEvent(evt); + } + } + + modifyInputs(); + }, id.c_str()); + } } void DoChange(std::string new_val) { @@ -177,13 +237,14 @@ namespace web { /// @param in_label The label that should appear on the Input. /// @param in_id The HTML ID to use for this Input (leave blank for auto-generated) Input(const std::function & in_cb, const std::string & in_type, - const std::string & in_label, const std::string & in_id="") + const std::string & in_label, const std::string & in_id="", bool show_value=false) : WidgetFacet(in_id) { info = new InputInfo(in_id); Info()->label = in_label; Info()->type = in_type; + Info()->show_value = show_value; Info()->autofocus = false; Info()->curr_val = ""; diff --git a/source/web/Style.h b/source/web/Style.h index a2886afa7c..3acd9e7143 100644 --- a/source/web/Style.h +++ b/source/web/Style.h @@ -19,6 +19,7 @@ #include "../tools/string_utils.h" #include +#include #include namespace emp { @@ -29,6 +30,7 @@ namespace web { class Style { private: std::map settings; ///< CSS Setting values being tracked. + std::set classes; ///< CSS classes public: Style() { ; } @@ -39,6 +41,14 @@ namespace web { /// Return a count of the number of settings that have been set. size_t GetSize() const { return settings.size(); } + /// Return a count of the number of classes that have been added. + size_t GetNClasses() const { return classes.size(); } + + Style & AddClass(const std::string & in_clss) { + classes.insert(in_clss); + return *this; + } + Style & DoSet(const std::string & in_set, const std::string & in_val) { settings[in_set] = in_val; @@ -72,18 +82,29 @@ namespace web { return settings; } - /// Remove all setting values. - void Clear() { settings.clear(); } + const std::set & GetClasses() const { + return classes; + } + + + /// Remove all setting values and all classes. + void Clear() { settings.clear(); classes.clear(); } /// Remove a specific setting value. void Remove(const std::string & setting) { settings.erase(setting); } + /// Remove a specific class + void RemoveClass(const std::string & clss) { + classes.erase(clss); + } + /// Apply ALL of the style settings to a specified widget. void Apply(const std::string & widget_id) { + // Stop immediately if nothing to set. - if (settings.size() == 0) return; + if (settings.size() == 0 && classes.size() == 0) return; // Find the current object only once. #ifdef __EMSCRIPTEN__ @@ -106,6 +127,20 @@ namespace web { << "' to '" << css_pair.second << "'."; #endif } + + for (std::string clss : classes) { + +#ifdef EMSCRIPTEN + EM_ASM_ARGS({ + var name = UTF8ToString($0); + emp_i.cur_obj.addClass( name); + }, clss.c_str()); +#else + std::cout << "Adding class to '" << widget_id << "': '" << clss; +#endif + + } + } /// Apply only a SPECIFIC style setting from the setting library. @@ -141,6 +176,31 @@ namespace web { #endif } + static void ApplyClass(const std::string & widget_id, const std::string & clss){ + +#ifdef EMSCRIPTEN + EM_ASM_ARGS({ + var id = UTF8ToString($0); + var name = UTF8ToString($1); + $( '#' + id ).addClass( name); + }, widget_id.c_str(), clss.c_str()); +#else + std::cout << "Adding class to '" << widget_id << "': '" << clss; +#endif + } + + static void ApplyRemoveClass(const std::string & widget_id, const std::string & clss){ +#ifdef EMSCRIPTEN + EM_ASM_ARGS({ + var id = UTF8ToString($0); + var name = UTF8ToString($1); + $( '#' + id ).removeClass( name); + }, widget_id.c_str(), clss.c_str()); +#else + std::cout << "Adding class to '" << widget_id << "': '" << clss; +#endif + } + /// Have any settings be set? operator bool() const { return (bool) settings.size(); } }; diff --git a/source/web/Widget.h b/source/web/Widget.h index 8c31103356..25fee04e2d 100644 --- a/source/web/Widget.h +++ b/source/web/Widget.h @@ -573,6 +573,12 @@ namespace web { info->extras.style.DoSet(setting, value); if (IsActive()) Style::Apply(info->id, setting, value); } + + virtual void DoCSS(const std::string & clss) { + info->extras.style.AddClass(clss); + if (IsActive()) Style::ApplyClass(info->id, clss); + } + /// Attr-related options may be overridden in derived classes that have multiple attributes. /// By default DoAttr will track the new information and apply it (if active) to the widget. virtual void DoAttr(const std::string & setting, const std::string & value) { @@ -626,10 +632,14 @@ namespace web { /// Allow multiple CSS settings to be provided as a single object. /// (still go through DoCSS given need for virtual re-routing.) return_t & SetCSS(const Style & in_style) { + emp_assert(info != nullptr); for (const auto & s : in_style.GetMap()) { DoCSS(s.first, s.second); } + for (const auto & s : in_style.GetClasses()) { + DoCSS(s); + } return (return_t &) *this; } diff --git a/source/web/color_map.h b/source/web/color_map.h index 846da23853..d65bf26385 100644 --- a/source/web/color_map.h +++ b/source/web/color_map.h @@ -30,9 +30,9 @@ namespace emp { /// Generate a string to describe a JS color out of HSL values. std::string ColorHSL(double h, double s, double l) { - emp_assert(h >= 0 && h <= 360); - emp_assert(s >= 0 && s <= 100); - emp_assert(l >= 0 && l <= 100); + emp_assert(h >= 0 && h <= 360, h); + emp_assert(s >= 0 && s <= 100, s); + emp_assert(l >= 0 && l <= 100, l); std::stringstream ss; ss << "hsl(" << h << ',' << s << "%," << l << "%)"; return ss.str(); diff --git a/source/web/d3/scales.h b/source/web/d3/scales.h index 70579ec361..47ef7510ff 100644 --- a/source/web/d3/scales.h +++ b/source/web/d3/scales.h @@ -79,6 +79,24 @@ namespace D3 { return EM_ASM_INT({return js.objects[$0]($1);},this->id, input); } + /// Calculate the ouput for [input], based on the scale's scaling function + std::string ApplyScaleString(double input) { + //TODO: make this work for other types + char * buffer = (char *) EM_ASM_INT({ + result = js.objects[$0]($1); + // console.log(result); + var buffer = Module._malloc(result.length+1); + Module.stringToUTF8(result, buffer, result.length*4+1); + return buffer; + + },this->id, input); + + std::string result = std::string(buffer); + free(buffer); + return result; + + } + //TODO:Getters diff --git a/tests/Evolve/NK.cc b/tests/Evolve/NK.cc index ec3f6e73a6..f1edd2704e 100644 --- a/tests/Evolve/NK.cc +++ b/tests/Evolve/NK.cc @@ -21,12 +21,13 @@ TEST_CASE("Test NK Landscapes", "[Evolve]") nk0.SetState(0, 0, 1.0); nk0.SetState(1, 0, 1.0); nk0.SetState(2, 0, 1.0); - nk0.SetState(3, 0, 1.0); + nk0.SetState(3, 0, 2.0); nk0.SetState(4, 0, 1.0); emp::vector states(5, 0); - REQUIRE(nk0.GetFitness(states) == 5.0); + REQUIRE(nk0.GetFitness(states) == 6.0); emp::BitVector bv(5); - REQUIRE(nk0.GetFitness(bv) == 5.0); + REQUIRE(nk0.GetFitness(bv) == 6.0); + REQUIRE(nk0.GetSiteFitness(3, bv) == 2); nk0.RandomizeStates(rnd, 100); REQUIRE(nk0.GetFitness(0, 0) != 1.0); diff --git a/tests/Evolve/Resource.cc b/tests/Evolve/Resource.cc index f47812be54..01c0ffd432 100644 --- a/tests/Evolve/Resource.cc +++ b/tests/Evolve/Resource.cc @@ -61,11 +61,11 @@ TEST_CASE("Test resources", "[evo]") pop.SetFitFun([](BitOrg &org){ return 10; }); - emp::vector > fit_funs; + emp::vector > fit_funs; - fit_funs.push_back([](const BitOrg &org){ return org.CountOnes()/N; }); - fit_funs.push_back([](const BitOrg &org){ return org[0]; }); - fit_funs.push_back([](const BitOrg &org){ return 1 - org[0]; }); + fit_funs.push_back([](BitOrg &org){ return org.CountOnes()/N; }); + fit_funs.push_back([](BitOrg &org){ return org[0]; }); + fit_funs.push_back([](BitOrg &org){ return 1 - org[0]; }); emp::ResourceSelect(pop, fit_funs, resources, 5, POP_SIZE); diff --git a/tests/Evolve/Systematics.cc b/tests/Evolve/Systematics.cc index ec0f21eeb3..3d79f9e6ae 100644 --- a/tests/Evolve/Systematics.cc +++ b/tests/Evolve/Systematics.cc @@ -21,7 +21,7 @@ TEST_CASE("Test Systematics", "[Evolve]") { -{ + // Taxon emp::Taxon tx(0, "a"); REQUIRE(tx.GetID() == 0); @@ -39,7 +39,7 @@ TEST_CASE("Test Systematics", "[Evolve]") emp::Ptr< emp::Taxon > parentPtr(&tx); emp::Taxon tx_1(1, "b", parentPtr); REQUIRE(tx_1.GetParent() == parentPtr); - tx_1.AddOffspring(); + tx_1.AddTotalOffspring(); REQUIRE(tx_1.GetTotalOffspring() == 1); REQUIRE(tx.GetTotalOffspring() == 1); @@ -54,13 +54,13 @@ TEST_CASE("Test Systematics", "[Evolve]") REQUIRE(sys1.GetNumTaxa() == 0); sys1.SetTrackSynchronous(true); - sys1.AddOrg(15.0, 0, 0, false); + sys1.AddOrg(15.0, {0,0}, 0); REQUIRE(sys1.GetNumActive() == 1); REQUIRE(sys1.GetTaxonAt(0)->GetInfo() == "small"); - sys1.AddOrg(56.0, 1, 0, true); + sys1.AddOrg(56.0, {1,1}, 0); REQUIRE(sys1.GetNumActive() == 2); REQUIRE(sys1.GetNextTaxonAt(1)->GetInfo() == "large"); - sys1.RemoveNextOrg(1); + sys1.RemoveOrg({1,1}); REQUIRE(sys1.GetNumActive() == 1); // Base setters and getters @@ -127,24 +127,23 @@ TEST_CASE("Test Systematics", "[Evolve]") REQUIRE(emp::CountMuts(ptr2, "short") == 16); REQUIRE(emp::CountMutSteps(ptr1, "short") == 1); REQUIRE(emp::CountMutSteps(ptr2, "short") == 2); -} -{ - emp::Systematics sys([](int & i){return i;}, true, true, true); - std::cout << "\nAddOrg 25 (id1, no parent)\n"; - auto id1 = sys.AddOrg(25, nullptr, 0); - std::cout << "\nAddOrg -10 (id2; parent id1)\n"; - auto id2 = sys.AddOrg(-10, id1, 6); - std::cout << "\nAddOrg 26 (id3; parent id1)\n"; - auto id3 = sys.AddOrg(26, id1, 10); - std::cout << "\nAddOrg 27 (id4; parent id2)\n"; - auto id4 = sys.AddOrg(27, id2, 25); - std::cout << "\nAddOrg 28 (id5; parent id2)\n"; - auto id5 = sys.AddOrg(28, id2, 32); - std::cout << "\nAddOrg 29 (id6; parent id5)\n"; - auto id6 = sys.AddOrg(29, id5, 39); - std::cout << "\nAddOrg 30 (id7; parent id1)\n"; - auto id7 = sys.AddOrg(30, id1, 6); + emp::Systematics sys([](const int & i){return i;}, true, true, true, false); + + std::cout << "\nAddOrg 25 (id1, no parent)\n"; + auto id1 = sys.AddOrg(25, nullptr, 0); + std::cout << "\nAddOrg -10 (id2; parent id1)\n"; + auto id2 = sys.AddOrg(-10, id1, 6); + std::cout << "\nAddOrg 26 (id3; parent id1)\n"; + auto id3 = sys.AddOrg(26, id1, 10); + std::cout << "\nAddOrg 27 (id4; parent id2)\n"; + auto id4 = sys.AddOrg(27, id2, 25); + std::cout << "\nAddOrg 28 (id5; parent id2)\n"; + auto id5 = sys.AddOrg(28, id2, 32); + std::cout << "\nAddOrg 29 (id6; parent id5)\n"; + auto id6 = sys.AddOrg(29, id5, 39); + std::cout << "\nAddOrg 30 (id7; parent id1)\n"; + auto id7 = sys.AddOrg(30, id1, 6); std::cout << "\nRemoveOrg (id2)\n"; @@ -242,15 +241,289 @@ TEST_CASE("Test Systematics", "[Evolve]") std::cout << "id3 = " << id3 << std::endl; std::cout << "id4 = " << id4 << std::endl; - std::cout << "\nLineage:\n"; - sys.PrintLineage(id4); + std::stringstream result; + + sys.PrintLineage(id4, result); sys.PrintStatus(); + + REQUIRE(result.str() == "Lineage:\n27\n-10\n25\n"); + + CHECK(sys.GetStoreActive() == 1); + CHECK(sys.GetStoreAncestors() == 1); + CHECK(sys.GetStoreOutside() == 1); + CHECK(sys.GetArchive() == 1); + CHECK(sys.GetTrackSynchronous() == 0); + CHECK(sys.GetNextID() == 19); + CHECK(sys.GetNumActive() == 11); + CHECK(sys.GetNumAncestors() == 7); + CHECK(sys.GetNumOutside() == 1); + + auto ancestors = sys.GetAncestors(); + emp::vector>> ancestor_vec(ancestors.begin(), ancestors.end()); + emp::Sort(ancestor_vec, [](emp::Ptr> & a, emp::Ptr> & b){ + return a->GetID() < b->GetID(); + }); + + CHECK(ancestor_vec[0]->GetID() == 1); + CHECK(ancestor_vec[0]->GetNumOrgs() == 0); + CHECK(ancestor_vec[0]->GetNumOff() == 3); + CHECK(ancestor_vec[0]->GetParent() == nullptr); + + CHECK(ancestor_vec[1]->GetID() == 2); + CHECK(ancestor_vec[1]->GetNumOrgs() == 0); + CHECK(ancestor_vec[1]->GetNumOff() == 2); + CHECK(ancestor_vec[1]->GetParent()->GetID() == 1); + + CHECK(ancestor_vec[2]->GetID() == 7); + CHECK(ancestor_vec[2]->GetNumOrgs() == 0); + CHECK(ancestor_vec[2]->GetNumOff() == 1); + CHECK(ancestor_vec[2]->GetParent()->GetID() == 1); + + CHECK(ancestor_vec[3]->GetID() == 8); + CHECK(ancestor_vec[3]->GetNumOrgs() == 0); + CHECK(ancestor_vec[3]->GetNumOff() == 1); + CHECK(ancestor_vec[3]->GetParent()->GetID() == 7); + + CHECK(ancestor_vec[4]->GetID() == 9); + CHECK(ancestor_vec[4]->GetNumOrgs() == 0); + CHECK(ancestor_vec[4]->GetNumOff() == 1); + CHECK(ancestor_vec[4]->GetParent()->GetID() == 8); + + CHECK(ancestor_vec[5]->GetID() == 13); + CHECK(ancestor_vec[5]->GetNumOrgs() == 0); + CHECK(ancestor_vec[5]->GetNumOff() == 1); + CHECK(ancestor_vec[5]->GetParent()->GetID() == 12); + + CHECK(ancestor_vec[6]->GetID() == 14); + CHECK(ancestor_vec[6]->GetNumOrgs() == 0); + CHECK(ancestor_vec[6]->GetNumOff() == 1); + CHECK(ancestor_vec[6]->GetParent()->GetID() == 13); + + auto outside_taxon = *(sys.GetOutside().begin()); + CHECK(outside_taxon->GetID() == 10); + CHECK(outside_taxon->GetNumOrgs() == 0); + CHECK(outside_taxon->GetNumOff() == 0); + CHECK(outside_taxon->GetParent()->GetID() == 8); + + auto active = sys.GetActive(); + emp::vector>> active_vec(active.begin(), active.end()); + emp::Sort(active_vec, [](emp::Ptr> & a, emp::Ptr> & b){ + return a->GetID() < b->GetID(); + }); + + CHECK(active_vec[0]->GetID() == 3); + CHECK(active_vec[0]->GetNumOrgs() == 1); + CHECK(active_vec[0]->GetNumOff() == 0); + CHECK(active_vec[0]->GetParent()->GetID() == 1); + + CHECK(active_vec[1]->GetID() == 4); + CHECK(active_vec[1]->GetNumOrgs() == 1); + CHECK(active_vec[1]->GetNumOff() == 0); + CHECK(active_vec[1]->GetParent()->GetID() == 2); + + CHECK(active_vec[2]->GetID() == 5); + CHECK(active_vec[2]->GetNumOrgs() == 1); + CHECK(active_vec[2]->GetNumOff() == 1); + CHECK(active_vec[2]->GetParent()->GetID() == 2); + + CHECK(active_vec[3]->GetID() == 6); + CHECK(active_vec[3]->GetNumOrgs() == 1); + CHECK(active_vec[3]->GetNumOff() == 0); + CHECK(active_vec[3]->GetParent()->GetID() == 5); + + CHECK(active_vec[4]->GetID() == 11); + CHECK(active_vec[4]->GetNumOrgs() == 1); + CHECK(active_vec[4]->GetNumOff() == 3); + CHECK(active_vec[4]->GetParent()->GetID() == 9); + + CHECK(active_vec[5]->GetID() == 12); + CHECK(active_vec[5]->GetNumOrgs() == 1); + CHECK(active_vec[5]->GetNumOff() == 1); + CHECK(active_vec[5]->GetParent()->GetID() == 11); + + CHECK(active_vec[6]->GetID() == 15); + CHECK(active_vec[6]->GetNumOrgs() == 1); + CHECK(active_vec[6]->GetNumOff() == 0); + CHECK(active_vec[6]->GetParent()->GetID() == 14); + + CHECK(active_vec[7]->GetID() == 16); + CHECK(active_vec[7]->GetNumOrgs() == 1); + CHECK(active_vec[7]->GetNumOff() == 0); + CHECK(active_vec[7]->GetParent()->GetID() == 11); + + CHECK(active_vec[8]->GetID() == 17); + CHECK(active_vec[8]->GetNumOrgs() == 1); + CHECK(active_vec[8]->GetNumOff() == 2); + CHECK(active_vec[8]->GetParent()->GetID() == 11); + + CHECK(active_vec[9]->GetID() == 18); + CHECK(active_vec[9]->GetNumOrgs() == 1); + CHECK(active_vec[9]->GetNumOff() == 0); + CHECK(active_vec[9]->GetParent()->GetID() == 17); + + CHECK(active_vec[10]->GetID() == 19); + CHECK(active_vec[10]->GetNumOrgs() == 1); + CHECK(active_vec[10]->GetNumOff() == 0); + CHECK(active_vec[10]->GetParent()->GetID() == 17); + } + +TEST_CASE("Test not tracking ancestors", "[Evolve]") +{ + emp::Systematics sys([](const int & i){return i;}, true, false, false, false); + + std::cout << "\nAddOrg 25 (id1, no parent)\n"; + auto id1 = sys.AddOrg(25, nullptr, 0); + std::cout << "\nAddOrg -10 (id2; parent id1)\n"; + auto id2 = sys.AddOrg(-10, id1, 6); + std::cout << "\nAddOrg 26 (id3; parent id1)\n"; + auto id3 = sys.AddOrg(26, id1, 10); + std::cout << "\nAddOrg 27 (id4; parent id2)\n"; + auto id4 = sys.AddOrg(27, id2, 25); + std::cout << "\nAddOrg 28 (id5; parent id2)\n"; + auto id5 = sys.AddOrg(28, id2, 32); + std::cout << "\nAddOrg 29 (id6; parent id5)\n"; + auto id6 = sys.AddOrg(29, id5, 39); + std::cout << "\nAddOrg 30 (id7; parent id1)\n"; + auto id7 = sys.AddOrg(30, id1, 6); + + + std::cout << "\nRemoveOrg (id2)\n"; + sys.RemoveOrg(id1); + sys.RemoveOrg(id2); + + double mpd = sys.GetMeanPairwiseDistance(); + + std::cout << "\nAddOrg 31 (id8; parent id7)\n"; + auto id8 = sys.AddOrg(31, id7, 11); + std::cout << "\nAddOrg 32 (id9; parent id8)\n"; + auto id9 = sys.AddOrg(32, id8, 19); + + std::cout << "\nAddOrg 33 (id10; parent id8)\n"; + auto id10 = sys.AddOrg(33, id8, 19); + + sys.RemoveOrg(id7); + sys.RemoveOrg(id8); + + sys.RemoveOrg(id10); + + + std::cout << "\nAddOrg 34 (id11; parent id9)\n"; + auto id11 = sys.AddOrg(34, id9, 22); + std::cout << "\nAddOrg 35 (id12; parent id10)\n"; + auto id12 = sys.AddOrg(35, id11, 23); + + sys.RemoveOrg(id9); + + std::cout << "\nAddOrg 36 (id13; parent id12)\n"; + auto id13 = sys.AddOrg(36, id12, 27); + std::cout << "\nAddOrg 37 (id14; parent id13)\n"; + auto id14 = sys.AddOrg(37, id13, 30); + + sys.RemoveOrg(id13); + + std::cout << "\nAddOrg 38 (id15; parent id14)\n"; + auto id15 = sys.AddOrg(38, id14, 33); + + sys.RemoveOrg(id14); + + std::cout << "\nAddOrg 39 (id16; parent id11)\n"; + auto id16 = sys.AddOrg(39, id11, 35); + std::cout << "\nAddOrg 40 (id17; parent id11)\n"; + auto id17 = sys.AddOrg(40, id11, 35); + + std::cout << "\nAddOrg 41 (id18; parent id17)\n"; + auto id18 = sys.AddOrg(41, id17, 36); + + std::cout << "\nAddOrg 42 (id19; parent id17)\n"; + auto id19 = sys.AddOrg(42, id17, 37); + REQUIRE(id17->GetTotalOffspring() > 0); + + std::cout << "id3 = " << id3 << std::endl; + std::cout << "id4 = " << id4 << std::endl; + + std::stringstream result; + + sys.PrintLineage(id4, result); + sys.PrintStatus(); + CHECK(result.str() == "Lineage:\n27\n"); + + CHECK(sys.GetStoreActive() == 1); + CHECK(sys.GetStoreAncestors() == 0); + CHECK(sys.GetStoreOutside() == 0); + CHECK(sys.GetArchive() == 0); + CHECK(sys.GetTrackSynchronous() == 0); + CHECK(sys.GetNextID() == 19); + CHECK(sys.GetNumActive() == 11); + CHECK(sys.GetNumAncestors() == 0); + CHECK(sys.GetNumOutside() == 0); + + auto active = sys.GetActive(); + emp::vector>> active_vec(active.begin(), active.end()); + emp::Sort(active_vec, [](emp::Ptr> & a, emp::Ptr> & b){ + return a->GetID() < b->GetID(); + }); + + CHECK(active_vec[0]->GetID() == 3); + CHECK(active_vec[0]->GetNumOrgs() == 1); + CHECK(active_vec[0]->GetNumOff() == 0); + CHECK(active_vec[0]->GetParent() == nullptr); + + CHECK(active_vec[1]->GetID() == 4); + CHECK(active_vec[1]->GetNumOrgs() == 1); + CHECK(active_vec[1]->GetNumOff() == 0); + CHECK(active_vec[1]->GetParent() == nullptr); + + CHECK(active_vec[2]->GetID() == 5); + CHECK(active_vec[2]->GetNumOrgs() == 1); + CHECK(active_vec[2]->GetNumOff() == 1); + CHECK(active_vec[2]->GetParent() == nullptr); + + CHECK(active_vec[3]->GetID() == 6); + CHECK(active_vec[3]->GetNumOrgs() == 1); + CHECK(active_vec[3]->GetNumOff() == 0); + CHECK(active_vec[3]->GetParent()->GetID() == 5); + + CHECK(active_vec[4]->GetID() == 11); + CHECK(active_vec[4]->GetNumOrgs() == 1); + CHECK(active_vec[4]->GetNumOff() == 3); + CHECK(active_vec[4]->GetParent() == nullptr); + + CHECK(active_vec[5]->GetID() == 12); + CHECK(active_vec[5]->GetNumOrgs() == 1); + CHECK(active_vec[5]->GetNumOff() == 1); + CHECK(active_vec[5]->GetParent()->GetID() == 11); + + CHECK(active_vec[6]->GetID() == 15); + CHECK(active_vec[6]->GetNumOrgs() == 1); + CHECK(active_vec[6]->GetNumOff() == 0); + CHECK(active_vec[6]->GetParent() == nullptr); + + CHECK(active_vec[7]->GetID() == 16); + CHECK(active_vec[7]->GetNumOrgs() == 1); + CHECK(active_vec[7]->GetNumOff() == 0); + CHECK(active_vec[7]->GetParent()->GetID() == 11); + + CHECK(active_vec[8]->GetID() == 17); + CHECK(active_vec[8]->GetNumOrgs() == 1); + CHECK(active_vec[8]->GetNumOff() == 2); + CHECK(active_vec[8]->GetParent()->GetID() == 11); + + CHECK(active_vec[9]->GetID() == 18); + CHECK(active_vec[9]->GetNumOrgs() == 1); + CHECK(active_vec[9]->GetNumOff() == 0); + CHECK(active_vec[9]->GetParent()->GetID() == 17); + + CHECK(active_vec[10]->GetID() == 19); + CHECK(active_vec[10]->GetNumOrgs() == 1); + CHECK(active_vec[10]->GetNumOff() == 0); + CHECK(active_vec[10]->GetParent()->GetID() == 17); + } TEST_CASE("Pointer to systematics", "[evo]") { emp::Ptr> sys; - sys.New([](int & i){return i;}, true, true, true); + sys.New([](const int & i){return i;}, true, true, true); sys.Delete(); } @@ -258,7 +531,7 @@ TEST_CASE("Test Data Struct", "[evo]") { emp::Ptr >> sys; - sys.New([](int & i){return i;}, true, true, true); + sys.New([](const int & i){return i;}, true, true, true, false); auto id1 = sys->AddOrg(1, nullptr); id1->GetData().fitness.Add(2); id1->GetData().phenotype = 6; @@ -267,7 +540,7 @@ TEST_CASE("Test Data Struct", "[evo]") id2->GetData().mut_counts["substitution"] = 2; id2->GetData().fitness.Add(1); id2->GetData().phenotype = 6; - REQUIRE(id2->GetData().mut_counts["substitution"] == 2); + CHECK(id2->GetData().mut_counts["substitution"] == 2); auto id3 = sys->AddOrg(3, id1); id3->GetData().mut_counts["substitution"] = 5; @@ -285,20 +558,20 @@ TEST_CASE("Test Data Struct", "[evo]") id5->GetData().phenotype = 6; - REQUIRE(CountMuts(id4) == 3); - REQUIRE(CountDeleteriousSteps(id4) == 1); - REQUIRE(CountPhenotypeChanges(id4) == 1); - REQUIRE(CountUniquePhenotypes(id4) == 2); + CHECK(CountMuts(id4) == 3); + CHECK(CountDeleteriousSteps(id4) == 1); + CHECK(CountPhenotypeChanges(id4) == 1); + CHECK(CountUniquePhenotypes(id4) == 2); - REQUIRE(CountMuts(id3) == 5); - REQUIRE(CountDeleteriousSteps(id3) == 1); - REQUIRE(CountPhenotypeChanges(id3) == 0); - REQUIRE(CountUniquePhenotypes(id3) == 1); + CHECK(CountMuts(id3) == 5); + CHECK(CountDeleteriousSteps(id3) == 1); + CHECK(CountPhenotypeChanges(id3) == 0); + CHECK(CountUniquePhenotypes(id3) == 1); - REQUIRE(CountMuts(id5) == 4); - REQUIRE(CountDeleteriousSteps(id5) == 2); - REQUIRE(CountPhenotypeChanges(id5) == 2); - REQUIRE(CountUniquePhenotypes(id5) == 2); + CHECK(CountMuts(id5) == 4); + CHECK(CountDeleteriousSteps(id5) == 2); + CHECK(CountPhenotypeChanges(id5) == 2); + CHECK(CountUniquePhenotypes(id5) == 2); sys.Delete(); @@ -319,7 +592,7 @@ TEST_CASE("World systematics integration", "[evo]") { emp::World> world; emp::Ptr sys; - sys.New([](emp::vector & v){return v;}, true, true, true); + sys.New([](const emp::vector & v){return v;}, true, true, true); world.AddSystematics(sys); world.SetMutFun([](emp::vector & org, emp::Random & r){return 0;}); @@ -392,18 +665,19 @@ TEST_CASE("Run world", "[evo]") { world.SetPopStruct_Mixed(true); - std::function gene_fun = - [](emp::AvidaGP & org) { + std::function gene_fun = + [](const emp::AvidaGP & org) { return org.GetGenome(); }; - std::function(emp::AvidaGP &)> phen_fun = - [](emp::AvidaGP & org) { + std::function(const emp::AvidaGP &)> phen_fun = + [](const emp::AvidaGP & org) { emp::vector phen; + emp::AvidaGP org2 = org; for (int i = 0; i < 16; i++) { - org.ResetHardware(); - org.Process(20); - phen.push_back(org.GetOutput(i)); + org2.ResetHardware(); + org2.Process(20); + phen.push_back(org2.GetOutput(i)); } return phen; }; @@ -458,6 +732,8 @@ TEST_CASE("Run world", "[evo]") { return num_muts; }); + world.SetAutoMutate(); + // Setup the fitness function. std::function fit_fun = [](emp::AvidaGP & org) { @@ -531,3 +807,272 @@ TEST_CASE("Run world", "[evo]") { // } // std::cout << std::endl; } + + + +TEST_CASE("Test GetCanopy", "[evo]") +{ + emp::Systematics sys([](const int & i){return i;}, true, true, true, false); + + auto id1 = sys.AddOrg(1, nullptr, 0); + auto id2 = sys.AddOrg(2, id1, 2); + auto id3 = sys.AddOrg(3, id1, 3); + auto id4 = sys.AddOrg(4, id2, 3); + + sys.RemoveOrg(id1, 3); + sys.RemoveOrg(id2, 5); + + auto can_set = sys.GetCanopyExtantRoots(4); + + // Both 3 and 4 were alive at time point 4 so they are the canopy roots + CHECK(can_set.size() == 2); + CHECK(Has(can_set, id3)); + CHECK(Has(can_set, id4)); + + can_set = sys.GetCanopyExtantRoots(2); + + // Both 3 and 4 were not alive at time point 2, so the canopy roots + // will be 1 and 2. + CHECK(can_set.size() == 2); + CHECK(Has(can_set, id1)); + CHECK(Has(can_set, id2)); + + sys.RemoveOrg(id3, 7); + + can_set = sys.GetCanopyExtantRoots(2); + + // Only 4 is alive, but it wasn't alive at time point 2. 2 is the + // only canopy root because even though 1 is alive, because 4's + // lineage diverged from 1 when 2 was born. + CHECK(can_set.size() == 1); + CHECK(Has(can_set, id2)); + + auto id5 = sys.AddOrg(5, id4, 8); + sys.RemoveOrg(id4, 9); + auto id6 = sys.AddOrg(6, id5, 10); + sys.RemoveOrg(id5, 11); + + can_set = sys.GetCanopyExtantRoots(7); + // Should only be 4 + CHECK(can_set.size() == 1); + CHECK(Has(can_set, id4)); + + can_set = sys.GetCanopyExtantRoots(9); + // Should only be 5 + CHECK(can_set.size() == 1); + CHECK(Has(can_set, id5)); + + + auto id7 = sys.AddOrg(7, id6, 12); + auto id8 = sys.AddOrg(8, id7, 13); + auto id9 = sys.AddOrg(9, id8, 14); + auto id10 = sys.AddOrg(10, id9, 15); + + sys.RemoveOrg(id6, 20); + sys.RemoveOrg(id7, 20); + sys.RemoveOrg(id8, 20); + sys.RemoveOrg(id9, 20); + + can_set = sys.GetCanopyExtantRoots(22); + // Should only be 10 + CHECK(can_set.size() == 1); + CHECK(Has(can_set, id10)); + + can_set = sys.GetCanopyExtantRoots(14); + // Should only be 9, even though others were alive + CHECK(can_set.size() == 1); + CHECK(Has(can_set, id9)); + + can_set = sys.GetCanopyExtantRoots(13); + // Should only be 8, because 9 wasn't born yet + CHECK(can_set.size() == 1); + CHECK(Has(can_set, id8)); + + can_set = sys.GetCanopyExtantRoots(11); + CHECK(can_set.size() == 1); + CHECK(Has(can_set, id6)); + + can_set = sys.GetCanopyExtantRoots(12); + CHECK(can_set.size() == 1); + CHECK(Has(can_set, id7)); + + can_set = sys.GetCanopyExtantRoots(9); + CHECK(can_set.size() == 1); + CHECK(Has(can_set, id5)); + + + // auto id5 = sys.AddOrg(28, id2, 32); + // std::cout << "\nAddOrg 29 (id6; parent id5)\n"; + // auto id6 = sys.AddOrg(29, id5, 39); + // std::cout << "\nAddOrg 30 (id7; parent id1)\n"; + // auto id7 = sys.AddOrg(30, id1, 6); + +} + +// Tests from Shao 1990 tree balance paper +TEST_CASE("Tree balance", "[evo]") { + emp::Systematics tree1([](const int & i){return i;}, true, true, false, false); + + auto tree1org1 = tree1.AddOrg(1, nullptr); + auto tree1org2 = tree1.AddOrg(2, tree1org1); + auto tree1org3 = tree1.AddOrg(3, tree1org2); + auto tree1org4 = tree1.AddOrg(4, tree1org3); + auto tree1org5 = tree1.AddOrg(5, tree1org3); + auto tree1org6 = tree1.AddOrg(6, tree1org2); + auto tree1org7 = tree1.AddOrg(7, tree1org6); + auto tree1org8 = tree1.AddOrg(8, tree1org6); + auto tree1org9 = tree1.AddOrg(9, tree1org1); + auto tree1org10 = tree1.AddOrg(10, tree1org9); + auto tree1org11 = tree1.AddOrg(11, tree1org9); + tree1.RemoveOrg(tree1org1); + tree1.RemoveOrg(tree1org2); + tree1.RemoveOrg(tree1org3); + tree1.RemoveOrg(tree1org6); + tree1.RemoveOrg(tree1org9); + + CHECK(tree1.SackinIndex() == 16); + + emp::Systematics tree2([](const int & i){return i;}, true, true, false, false); + + auto tree2org1 = tree2.AddOrg(1, nullptr); + auto tree2org2 = tree2.AddOrg(2, tree2org1); + auto tree2org3 = tree2.AddOrg(3, tree2org2); + auto tree2org4 = tree2.AddOrg(4, tree2org3); + auto tree2org5 = tree2.AddOrg(5, tree2org3); + auto tree2org6 = tree2.AddOrg(6, tree2org2); + auto tree2org7 = tree2.AddOrg(7, tree2org1); + auto tree2org8 = tree2.AddOrg(8, tree2org7); + auto tree2org9 = tree2.AddOrg(9, tree2org7); + auto tree2org10 = tree2.AddOrg(10, tree2org9); + auto tree2org11 = tree2.AddOrg(11, tree2org9); + + tree2.RemoveOrg(tree2org1); + tree2.RemoveOrg(tree2org2); + tree2.RemoveOrg(tree2org3); + tree2.RemoveOrg(tree2org7); + tree2.RemoveOrg(tree2org9); + + CHECK(tree2.SackinIndex() == 16); + + emp::Systematics tree3([](const int & i){return i;}, true, true, false, false); + + auto tree3org1 = tree3.AddOrg(1, nullptr); + auto tree3org2 = tree3.AddOrg(2, tree3org1); + auto tree3org3 = tree3.AddOrg(3, tree3org2); + auto tree3org4 = tree3.AddOrg(4, tree3org2); + auto tree3org5 = tree3.AddOrg(5, tree3org4); + auto tree3org6 = tree3.AddOrg(6, tree3org4); + auto tree3org7 = tree3.AddOrg(7, tree3org6); + auto tree3org8 = tree3.AddOrg(8, tree3org6); + auto tree3org9 = tree3.AddOrg(9, tree3org1); + auto tree3org10 = tree3.AddOrg(10, tree3org9); + auto tree3org11 = tree3.AddOrg(11, tree3org9); + + tree3.RemoveOrg(tree3org1); + tree3.RemoveOrg(tree3org2); + tree3.RemoveOrg(tree3org4); + tree3.RemoveOrg(tree3org6); + tree3.RemoveOrg(tree3org9); + + CHECK(tree3.SackinIndex() == 17); + + emp::Systematics tree29([](const int & i){return i;}, true, true, false, false); + + auto tree29org1 = tree29.AddOrg(1, nullptr); + auto tree29org2 = tree29.AddOrg(2, tree29org1); + auto tree29org3 = tree29.AddOrg(3, tree29org1); + auto tree29org4 = tree29.AddOrg(4, tree29org3); + auto tree29org5 = tree29.AddOrg(5, tree29org3); + auto tree29org6 = tree29.AddOrg(6, tree29org3); + auto tree29org7 = tree29.AddOrg(7, tree29org3); + auto tree29org8 = tree29.AddOrg(8, tree29org3); + + tree29.RemoveOrg(tree29org1); + tree29.RemoveOrg(tree29org3); + + CHECK(tree29.SackinIndex() == 11); + + emp::Systematics tree30([](const int & i){return i;}, true, true, false, false); + + auto tree30org1 = tree30.AddOrg(1, nullptr); + auto tree30org2 = tree30.AddOrg(2, tree30org1); + auto tree30org3 = tree30.AddOrg(3, tree30org1); + auto tree30org4 = tree30.AddOrg(4, tree30org1); + auto tree30org5 = tree30.AddOrg(5, tree30org4); + auto tree30org6 = tree30.AddOrg(6, tree30org4); + auto tree30org7 = tree30.AddOrg(7, tree30org4); + auto tree30org8 = tree30.AddOrg(8, tree30org4); + + tree30.RemoveOrg(tree30org1); + tree30.RemoveOrg(tree30org4); + + CHECK(tree30.SackinIndex() == 10); + + emp::Systematics tree31([](const int & i){return i;}, true, true, false, false); + + auto tree31org1 = tree31.AddOrg(1, nullptr); + auto tree31org2 = tree31.AddOrg(2, tree31org1); + auto tree31org3 = tree31.AddOrg(3, tree31org1); + auto tree31org4 = tree31.AddOrg(4, tree31org1); + auto tree31org5 = tree31.AddOrg(5, tree31org1); + auto tree31org6 = tree31.AddOrg(6, tree31org5); + auto tree31org7 = tree31.AddOrg(7, tree31org5); + auto tree31org8 = tree31.AddOrg(8, tree31org5); + + tree31.RemoveOrg(tree31org1); + tree31.RemoveOrg(tree31org5); + + CHECK(tree31.SackinIndex() == 9); + + emp::Systematics tree32([](const int & i){return i;}, true, true, false, false); + + auto tree32org1 = tree32.AddOrg(1, nullptr); + auto tree32org2 = tree32.AddOrg(2, tree32org1); + auto tree32org3 = tree32.AddOrg(3, tree32org1); + auto tree32org4 = tree32.AddOrg(4, tree32org1); + auto tree32org5 = tree32.AddOrg(5, tree32org1); + auto tree32org6 = tree32.AddOrg(6, tree32org1); + auto tree32org7 = tree32.AddOrg(7, tree32org6); + auto tree32org8 = tree32.AddOrg(8, tree32org6); + + tree32.RemoveOrg(tree32org1); + tree32.RemoveOrg(tree32org6); + + CHECK(tree32.SackinIndex() == 8); + + emp::Systematics tree33([](const int & i){return i;}, true, true, false, false); + + auto tree33org1 = tree33.AddOrg(1, nullptr); + auto tree33org2 = tree33.AddOrg(2, tree33org1); + auto tree33org3 = tree33.AddOrg(3, tree33org1); + auto tree33org4 = tree33.AddOrg(4, tree33org1); + auto tree33org5 = tree33.AddOrg(5, tree33org1); + auto tree33org6 = tree33.AddOrg(6, tree33org1); + auto tree33org7 = tree33.AddOrg(7, tree33org1); + + tree33.RemoveOrg(tree33org1); + CHECK(tree33.SackinIndex() == 6); + + // From CollessLike metric paper + emp::Systematics treecl([](const int & i){return i;}, true, true, false, false); + auto treeclorg1 = treecl.AddOrg(1, nullptr); + auto treeclorg2 = treecl.AddOrg(2, treeclorg1); + auto treeclorg3 = treecl.AddOrg(3, treeclorg1); + auto treeclorg4 = treecl.AddOrg(4, treeclorg2); + auto treeclorg5 = treecl.AddOrg(5, treeclorg2); + auto treeclorg6 = treecl.AddOrg(6, treeclorg2); + auto treeclorg7 = treecl.AddOrg(7, treeclorg2); + auto treeclorg8 = treecl.AddOrg(8, treeclorg2); + auto treeclorg9 = treecl.AddOrg(9, treeclorg3); + auto treeclorg10 = treecl.AddOrg(10, treeclorg3); + auto treeclorg11 = treecl.AddOrg(11, treeclorg10); + auto treeclorg12 = treecl.AddOrg(12, treeclorg10); + + treecl.RemoveOrg(treeclorg1); + treecl.RemoveOrg(treeclorg2); + treecl.RemoveOrg(treeclorg3); + treecl.RemoveOrg(treeclorg10); + + CHECK(treecl.SackinIndex() == 18); + CHECK(treecl.CollessLikeIndex() == Approx(1.746074)); +} diff --git a/tests/config/config_web_interface.cc b/tests/config/config_web_interface.cc new file mode 100644 index 0000000000..7731f64a8f --- /dev/null +++ b/tests/config/config_web_interface.cc @@ -0,0 +1,11 @@ +#define CATCH_CONFIG_MAIN +#include "third-party/Catch/single_include/catch.hpp" + +#include "config/config_web_interface.h" +#include "assets/config_setup.h" + +TEST_CASE("Test config_web_interface", "[config]") +{ + MyConfig config; + emp::ConfigWebUI ui(config); +} diff --git a/tests/data/DataFile.cc b/tests/data/DataFile.cc index c9f55edb59..d04fabf3ca 100644 --- a/tests/data/DataFile.cc +++ b/tests/data/DataFile.cc @@ -45,50 +45,53 @@ int test_fun() { return val += 3; } -TEST_CASE("Test DataFile", "[data]") { - int test_int = 5; - - emp::DataFile dfile("new_test_file.dat"); - - REQUIRE(dfile.GetFilename() == "new_test_file.dat"); - - emp::DataMonitor data_fracs; - emp::DataMonitor data_squares; - emp::DataMonitor data_cubes; - - dfile.AddCurrent(data_fracs); - dfile.AddCurrent(data_squares); - dfile.AddCurrent(data_cubes); - dfile.AddMean(data_cubes); - dfile.AddTotal(data_cubes); - dfile.AddMin(data_cubes); - dfile.AddMax(data_cubes); - dfile.AddFun(test_fun); - dfile.AddVar(test_int); - - double frac = 0.0; - for (size_t i = 0; i < 10; i++) { - test_int += i; - data_fracs.Add(frac += 0.01); - data_squares.Add((int)(i*i)); - data_cubes.Add(i*i*i); - dfile.Update(); - - // std::cout << i << std::endl; - } - - dfile.SetupLine("[[",":", "]]\n"); - for (size_t i = 10; i < 20; i++) { - data_fracs.Add(frac += 0.01); - data_squares.Add((int)(i*i)); - data_cubes.Add(i*i*i); - dfile.Update(); - - // std::cout << i << std::endl; - } - - REQUIRE(compareFiles("new_test_file.dat", "test_file.dat")); -} +// TEST_CASE("Test DataFile", "[data]") { +// int test_int = 5; + +// emp::DataFile dfile("new_test_file.dat"); + +// REQUIRE(dfile.GetFilename() == "new_test_file.dat"); + +// emp::DataMonitor data_fracs; +// emp::DataMonitor data_squares; +// emp::DataMonitor data_cubes; + +// dfile.AddCurrent(data_fracs); +// dfile.AddCurrent(data_squares); +// dfile.AddCurrent(data_cubes); +// dfile.AddMean(data_cubes); +// dfile.AddTotal(data_cubes); +// dfile.AddMin(data_cubes); +// dfile.AddMax(data_cubes); +// dfile.AddStandardDeviation(data_cubes); +// dfile.AddSkew(data_cubes); +// dfile.AddKurtosis(data_cubes); +// dfile.AddFun(test_fun); +// dfile.AddVar(test_int); + +// double frac = 0.0; +// for (size_t i = 0; i < 10; i++) { +// test_int += i; +// data_fracs.Add(frac += 0.01); +// data_squares.Add((int)(i*i)); +// data_cubes.Add(i*i*i); +// dfile.Update(); + +// // std::cout << i << std::endl; +// } + +// dfile.SetupLine("[[",":", "]]\n"); +// for (size_t i = 10; i < 20; i++) { +// data_fracs.Add(frac += 0.01); +// data_squares.Add((int)(i*i)); +// data_cubes.Add(i*i*i); +// dfile.Update(); + +// // std::cout << i << std::endl; +// } + +// REQUIRE(compareFiles("new_test_file.dat", "test_file.dat")); +// } TEST_CASE("Test Container DataFile", "[data]") { @@ -137,73 +140,73 @@ TEST_CASE("Test Container DataFile", "[data]") { REQUIRE(compareFiles("new_test_make_container_file.dat", "test_make_container_file.dat")); } -TEST_CASE("Test timing", "[data]") { - int test_int = 5; +// TEST_CASE("Test timing", "[data]") { +// int test_int = 5; - emp::DataFile dfile("new_test_timing_file.dat"); +// emp::DataFile dfile("new_test_timing_file.dat"); - emp::DataMonitor data_fracs; - emp::DataMonitor data_squares; - emp::DataMonitor data_cubes; +// emp::DataMonitor data_fracs; +// emp::DataMonitor data_squares; +// emp::DataMonitor data_cubes; - dfile.AddVar(test_int); - dfile.AddCurrent(data_fracs); - dfile.AddCurrent(data_squares); - dfile.AddCurrent(data_cubes); - dfile.AddMean(data_cubes); - dfile.AddTotal(data_cubes); - dfile.AddMin(data_cubes); - dfile.AddMax(data_cubes); - dfile.AddFun(test_fun); +// dfile.AddVar(test_int); +// dfile.AddCurrent(data_fracs); +// dfile.AddCurrent(data_squares); +// dfile.AddCurrent(data_cubes); +// dfile.AddMean(data_cubes); +// dfile.AddTotal(data_cubes); +// dfile.AddMin(data_cubes); +// dfile.AddMax(data_cubes); +// dfile.AddFun(test_fun); - double frac = 0.0; +// double frac = 0.0; - dfile.SetTimingRepeat(2); +// dfile.SetTimingRepeat(2); - for (size_t i = 0; i < 10; i++) { - test_int = i; - data_fracs.Add(frac += 0.01); - data_squares.Add((int)(i*i)); - data_cubes.Add(i*i*i); - dfile.Update(i); +// for (size_t i = 0; i < 10; i++) { +// test_int = i; +// data_fracs.Add(frac += 0.01); +// data_squares.Add((int)(i*i)); +// data_cubes.Add(i*i*i); +// dfile.Update(i); - // std::cout << i << std::endl; - } +// // std::cout << i << std::endl; +// } - dfile.SetTimingOnce(5); +// dfile.SetTimingOnce(5); - for (size_t i = 0; i < 10; i++) { - test_int = i; - data_fracs.Add(frac += 0.01); - data_squares.Add((int)(i*i)); - data_cubes.Add(i*i*i); - dfile.Update(i); - // std::cout << i << std::endl; - } +// for (size_t i = 0; i < 10; i++) { +// test_int = i; +// data_fracs.Add(frac += 0.01); +// data_squares.Add((int)(i*i)); +// data_cubes.Add(i*i*i); +// dfile.Update(i); +// // std::cout << i << std::endl; +// } - dfile.SetTimingRange(2, 3, 9); +// dfile.SetTimingRange(2, 3, 9); - for (size_t i = 0; i < 10; i++) { - test_int = i; - data_fracs.Add(frac += 0.01); - data_squares.Add((int)(i*i)); - data_cubes.Add(i*i*i); - dfile.Update(i); +// for (size_t i = 0; i < 10; i++) { +// test_int = i; +// data_fracs.Add(frac += 0.01); +// data_squares.Add((int)(i*i)); +// data_cubes.Add(i*i*i); +// dfile.Update(i); - // std::cout << i << std::endl; - } +// // std::cout << i << std::endl; +// } - dfile.SetTiming([](size_t ud){return (bool)floor(sqrt((double)ud) == ceil(sqrt((double)ud)));}); +// dfile.SetTiming([](size_t ud){return (bool)floor(sqrt((double)ud) == ceil(sqrt((double)ud)));}); - for (size_t i = 0; i < 10; i++) { - test_int = i; - data_fracs.Add(frac += 0.01); - data_squares.Add((int)(i*i)); - data_cubes.Add(i*i*i); - dfile.Update(i); +// for (size_t i = 0; i < 10; i++) { +// test_int = i; +// data_fracs.Add(frac += 0.01); +// data_squares.Add((int)(i*i)); +// data_cubes.Add(i*i*i); +// dfile.Update(i); - // std::cout << i << std::endl; - } +// // std::cout << i << std::endl; +// } - REQUIRE(compareFiles("new_test_timing_file.dat", "test_timing_file.dat")); -} \ No newline at end of file +// REQUIRE(compareFiles("new_test_timing_file.dat", "test_timing_file.dat")); +// } \ No newline at end of file diff --git a/tests/data/DataNode.cc b/tests/data/DataNode.cc index 7360df0ef8..c754476648 100644 --- a/tests/data/DataNode.cc +++ b/tests/data/DataNode.cc @@ -165,10 +165,19 @@ TEST_CASE("Test DataRange", "[data]") { REQUIRE(data.GetPercentile(25) == 0); REQUIRE(data.GetPercentile(100) == 1600); - // std::cout << std::endl; - // data.PrintDebug(); + std::stringstream result; - // std::cout << std::endl; + data.PrintDebug(result); + REQUIRE(result.str() == "Main DataNode.\nDataNodeModule for data::Pull. (level 8)\nDataNodeModule for data::Range. (level 4)\nDataNodeModule for data::Log. (level 2)\nDataNodeModule for data::Current. (level 0)\nBASE DataNodeModule.\n"); + result.str(""); + + data.PrintLog(result); + REQUIRE(result.str() == "100, 200, 300, 400, 500, -800, -800, 1600, 0, 0\n"); + result.str(""); + + data.Add(5); + data.PrintCurrent(result); + REQUIRE(result.str() == "5"); } @@ -270,16 +279,23 @@ TEST_CASE("Test DataStats", "[data]") { TEST_CASE("Test histogram", "[data]") { emp::DataNode data; data.SetupBins(1,21,10); - data.Add(1,2,1,19); + data.Add(1,2,1,19, 0, -1, 49); REQUIRE(data.GetHistMin() == 1); REQUIRE(data.GetHistWidth(5) == 2); + REQUIRE(data.GetHistMax() == 21); REQUIRE(data.GetBinMins() == emp::vector({1,3,5,7,9,11,13,15,17,19})); - + + REQUIRE(data.GetOverflow() == 1); + REQUIRE(data.GetUnderflow() == 2); REQUIRE(data.GetHistCount(9) == 1); REQUIRE(data.GetHistCounts() == emp::vector({3,0,0,0,0,0,0,0,0,1})); data.Reset(); REQUIRE(data.GetHistCounts() == emp::vector({0,0,0,0,0,0,0,0,0,0})); + + std::stringstream result; + data.PrintDebug(result); + REQUIRE(result.str() == "Main DataNode.\nDataNodeModule for data::Pull. (level 8)\nDataNodeModule for data::Histogram. (level 5)\nDataNodeModule for data::Range. (level 4)\nDataNodeModule for data::Log. (level 2)\nDataNodeModule for data::Current. (level 0)\nBASE DataNodeModule.\n"); } diff --git a/tests/data/test_file.dat b/tests/data/test_file.dat index 6de3388b87..6dc26e85df 100644 --- a/tests/data/test_file.dat +++ b/tests/data/test_file.dat @@ -1,20 +1,20 @@ -0.01,0,0,0,0,0,0,13,5 -0.02,1,1,0.5,1,0,1,16,6 -0.03,4,8,3,9,0,8,19,8 -0.04,9,27,9,36,0,27,22,11 -0.05,16,64,20,100,0,64,25,15 -0.06,25,125,37.5,225,0,125,28,20 -0.07,36,216,63,441,0,216,31,26 -0.08,49,343,98,784,0,343,34,33 -0.09,64,512,144,1296,0,512,37,41 -0.1,81,729,202.5,2025,0,729,40,50 -[[0.11:100:1000:275:3025:0:1000:43:50]] -[[0.12:121:1331:363:4356:0:1331:46:50]] -[[0.13:144:1728:468:6084:0:1728:49:50]] -[[0.14:169:2197:591.5:8281:0:2197:52:50]] -[[0.15:196:2744:735:11025:0:2744:55:50]] -[[0.16:225:3375:900:14400:0:3375:58:50]] -[[0.17:256:4096:1088:18496:0:4096:61:50]] -[[0.18:289:4913:1300.5:23409:0:4913:64:50]] -[[0.19:324:5832:1539:29241:0:5832:67:50]] -[[0.2:361:6859:1805:36100:0:6859:70:50]] +0.01,0,0,0,0,0,0,0,-nan,-nan,13,5 +0.02,1,1,0.5,1,0,1,0.5,0,-2,16,6 +0.03,4,8,3,9,0,8,3.55903,0.665469,-1.5,19,8 +0.04,9,27,9,36,0,27,10.8397,0.90094,-0.906129,22,11 +0.05,16,64,20,100,0,64,24.0416,0.992224,-0.568543,25,15 +0.06,25,125,37.5,225,0,125,44.8655,1.0341,-0.382533,28,20 +0.07,36,216,63,441,0,216,75.0124,1.05548,-0.274217,31,26 +0.08,49,343,98,784,0,343,116.183,1.06714,-0.207518,34,33 +0.09,64,512,144,1296,0,512,170.078,1.07375,-0.1645,37,41 +0.1,81,729,202.5,2025,0,729,238.399,1.07755,-0.135707,40,50 +[[0.11:100:1000:275:3025:0:1000:322.847:1.07969:-0.115862:43:50]] +[[0.12:121:1331:363:4356:0:1331:425.121:1.08084:-0.101862:46:50]] +[[0.13:144:1728:468:6084:0:1728:546.924:1.08137:-0.0918054:49:50]] +[[0.14:169:2197:591.5:8281:0:2197:689.956:1.0815:-0.0844808:52:50]] +[[0.15:196:2744:735:11025:0:2744:855.917:1.08138:-0.0790924:55:50]] +[[0.16:225:3375:900:14400:0:3375:1046.51:1.0811:-0.0751027:58:50]] +[[0.17:256:4096:1088:18496:0:4096:1263.43:1.08071:-0.0721403:61:50]] +[[0.18:289:4913:1300.5:23409:0:4913:1508.39:1.08027:-0.0699431:64:50]] +[[0.19:324:5832:1539:29241:0:5832:1783.08:1.07978:-0.0683227:67:50]] +[[0.2:361:6859:1805:36100:0:6859:2089.2:1.07928:-0.0671417:70:50]] diff --git a/tests/test_OEE.cc b/tests/test_OEE.cc new file mode 100644 index 0000000000..8355c639d5 --- /dev/null +++ b/tests/test_OEE.cc @@ -0,0 +1,146 @@ +#define CATCH_CONFIG_MAIN + +#ifndef EMP_TRACK_MEM +#define EMP_TRACK_MEM +#endif + +#include "../third-party/Catch/single_include/catch.hpp" + +#include "Evolve/OEE.h" +#include "Evolve/World.h" +#include "Evolve/World_output.h" + +TEST_CASE("OEE", "[evo]") { + emp::Random random; + emp::World world(random, "OEEWorld"); + + emp::Ptr > sys_ptr; + sys_ptr.New([](int org){return org;}, true, true, false); + + emp::OEETracker> oee(sys_ptr, [](int org){return org;}, [](int org){return org;}, true); + oee.SetResolution(1); + oee.SetGenerationInterval(1); + + sys_ptr->AddOrg(1, 0, 0, false); + sys_ptr->AddOrg(2, 1, 0, false); + sys_ptr->AddOrg(3, 2, 0, false); + sys_ptr->PrintStatus(); + oee.Update(0); + + // Coalescence interval hasn't passed yet + CHECK(oee.CoalescenceFilter(0).size() == 0); + CHECK(oee.GetDataNode("change")->GetCurrent() == 0); + CHECK(oee.GetDataNode("novelty")->GetCurrent() == 0); + CHECK(oee.GetDataNode("diversity")->GetCurrent() == 0); + CHECK(oee.GetDataNode("complexity")->GetCurrent() == 0); + + sys_ptr->SetNextParent(0); + sys_ptr->RemoveOrgAfterRepro(2, 1); + sys_ptr->AddOrg(4, 2, 1, false); + sys_ptr->PrintStatus(); + oee.Update(1); + + // 1 and 2 should make it through filter + CHECK(oee.CoalescenceFilter(1).size() == 2); + CHECK(oee.GetDataNode("change")->GetCurrent() == 2); + CHECK(oee.GetDataNode("novelty")->GetCurrent() == 2); + CHECK(oee.GetDataNode("diversity")->GetCurrent() == 1); + CHECK(oee.GetDataNode("complexity")->GetCurrent() == 2); + + // If we change nothing, 4 will now pass filter + oee.Update(2); + CHECK(oee.CoalescenceFilter(2).size() == 3); + CHECK(oee.GetDataNode("change")->GetCurrent() == 1); + CHECK(oee.GetDataNode("novelty")->GetCurrent() == 1); + CHECK(oee.GetDataNode("diversity")->GetCurrent() == Approx(1.5853)); + CHECK(oee.GetDataNode("complexity")->GetCurrent() == 4); + + // If we change nothing again, change and novelty should drop to 0 + oee.Update(3); + CHECK(oee.CoalescenceFilter(3).size() == 3); + CHECK(oee.GetDataNode("change")->GetCurrent() == 0); + CHECK(oee.GetDataNode("novelty")->GetCurrent() == 0); + CHECK(oee.GetDataNode("diversity")->GetCurrent() == Approx(1.5853)); + CHECK(oee.GetDataNode("complexity")->GetCurrent() == 4); + + sys_ptr->SetNextParent(0); + sys_ptr->RemoveOrgAfterRepro(0,4); + sys_ptr->AddOrg(1, 0, 4, false); + sys_ptr->PrintStatus(); + + // Replacing 1 with a copy of itself should change nothing + oee.Update(4); + CHECK(oee.CoalescenceFilter(4).size() == 3); + CHECK(oee.GetDataNode("change")->GetCurrent() == 0); + CHECK(oee.GetDataNode("novelty")->GetCurrent() == 0); + CHECK(oee.GetDataNode("diversity")->GetCurrent() == Approx(1.5853)); + CHECK(oee.GetDataNode("complexity")->GetCurrent() == 4); + + sys_ptr->SetNextParent(0); + sys_ptr->RemoveOrgAfterRepro(0, 5); + sys_ptr->AddOrg(10, 0, 5, false); + sys_ptr->PrintStatus(); + + // Replacing 1 with a new descendant should change nothing at first + // because 1 still has descendants and 10 hasn't survived filter time + oee.Update(5); + CHECK(oee.CoalescenceFilter(5).size() == 3); + CHECK(oee.GetDataNode("change")->GetCurrent() == 0); + CHECK(oee.GetDataNode("novelty")->GetCurrent() == 0); + CHECK(oee.GetDataNode("diversity")->GetCurrent() == Approx(1.5853)); + CHECK(oee.GetDataNode("complexity")->GetCurrent() == 4); + + // 10 survives the filter and replaces 1 because 1 is no longer in the + // set being filtered + oee.Update(6); + CHECK(oee.CoalescenceFilter(6).size() == 3); + CHECK(oee.GetDataNode("change")->GetCurrent() == 1); + CHECK(oee.GetDataNode("novelty")->GetCurrent() == 1); + CHECK(oee.GetDataNode("diversity")->GetCurrent() == Approx(1.5853)); + CHECK(oee.GetDataNode("complexity")->GetCurrent() == 10); + + sys_ptr->SetNextParent(0); + sys_ptr->RemoveOrgAfterRepro(1, 7); + sys_ptr->AddOrg(2, 1, 7, false); + sys_ptr->PrintStatus(); + + // Adding an independent origin of 2 should increase change but not novelty + // (the time after this). For now, we're replacing 2, leaving it with + // no descendants, so it should go away immediately + oee.Update(7); + CHECK(oee.CoalescenceFilter(7).size() == 2); + CHECK(oee.GetDataNode("change")->GetCurrent() == 0); + CHECK(oee.GetDataNode("novelty")->GetCurrent() == 0); + CHECK(oee.GetDataNode("diversity")->GetCurrent() == Approx(1)); + CHECK(oee.GetDataNode("complexity")->GetCurrent() == 10); + + // Now we see the bump in change + oee.Update(8); + CHECK(oee.CoalescenceFilter(8).size() == 3); + CHECK(oee.GetDataNode("change")->GetCurrent() == 1); + CHECK(oee.GetDataNode("novelty")->GetCurrent() == 0); + CHECK(oee.GetDataNode("diversity")->GetCurrent() == Approx(1.5853)); + CHECK(oee.GetDataNode("complexity")->GetCurrent() == 10); + + sys_ptr->SetNextParent(0); + sys_ptr->AddOrg(10, 3, 9, false); + sys_ptr->PrintStatus(); + + // No effect this time + oee.Update(8); + CHECK(oee.GetDataNode("change")->GetCurrent() == 0); + CHECK(oee.GetDataNode("novelty")->GetCurrent() == 0); + CHECK(oee.GetDataNode("diversity")->GetCurrent() == Approx(1.5853)); + CHECK(oee.GetDataNode("complexity")->GetCurrent() == 10); + + // Now we should see diversity change + oee.Update(9); + CHECK(oee.GetDataNode("change")->GetCurrent() == 0); + CHECK(oee.GetDataNode("novelty")->GetCurrent() == 0); + CHECK(oee.GetDataNode("diversity")->GetCurrent() == Approx(1.5)); + CHECK(oee.GetDataNode("complexity")->GetCurrent() == 10); + + sys_ptr.Delete(); + + +} \ No newline at end of file diff --git a/tests/tools/Graph.cc b/tests/tools/Graph.cc index bd460fafdf..a2cf62bb7b 100644 --- a/tests/tools/Graph.cc +++ b/tests/tools/Graph.cc @@ -25,6 +25,22 @@ TEST_CASE("Test Graph", "[tools]") REQUIRE(!graph.HasEdge(2,4)); REQUIRE((graph.GetEdgeCount() == 1)); + // Labels + graph.SetLabel(1, "node 1"); + REQUIRE(graph.GetLabel(1) == "node 1"); + + // Degree + REQUIRE(graph.GetInDegree(1) == 1); + REQUIRE(graph.GetInDegree(0) == 0); + REQUIRE(graph.GetDegree(0) == 1); + + // Getters + emp::Graph::Node n = graph.GetNode(1); + REQUIRE(n.GetLabel() == "node 1"); + emp::vector nodes = graph.GetNodes(); + REQUIRE(nodes.size() == 10); + REQUIRE(nodes[1].GetLabel() == "node 1"); + // Assignment emp::Graph g2 = graph; REQUIRE((g2.GetEdgeCount() == 1)); @@ -40,7 +56,13 @@ TEST_CASE("Test Graph", "[tools]") graph.AddEdge(0,3); graph.AddEdge(0,6); REQUIRE((graph.GetDegree(0) == 3)); - + emp::BitVector bit_v(10); + REQUIRE(graph.GetMaskedDegree(0, bit_v) == 0); + bit_v.Set(3); + REQUIRE(graph.GetMaskedDegree(0, bit_v) == 1); + bit_v.Set(6); + REQUIRE(graph.GetMaskedDegree(0, bit_v) == 2); + // GetEdgeSet emp::BitVector bv = graph.GetEdgeSet(0); REQUIRE(!bv[0]); @@ -106,6 +128,9 @@ TEST_CASE("Test Graph", "[tools]") REQUIRE((wgraph.GetWeight(0,1) == 3.2)); REQUIRE((wgraph.GetDegree(0) == 1)); + emp::vector > weights = wgraph.GetWeights(); + REQUIRE(weights[0][1] == 3.2); + // AddEdgePair wgraph.AddEdgePair(3, 2, 1.5); REQUIRE(wgraph.HasEdgePair(3,2)); diff --git a/tests/tools/Random.cc b/tests/tools/Random.cc index be488adb76..aa23229b7c 100644 --- a/tests/tools/Random.cc +++ b/tests/tools/Random.cc @@ -34,32 +34,32 @@ TEST_CASE("Test Random", "[tools]") REQUIRE(rnd.GetSeed() == 1); rnd.ResetSeed(5); REQUIRE(rnd.GetSeed() == 5); - + // Get Double double r_d = rnd.GetDouble(emp::Range(0.0,5.0)); REQUIRE(r_d >= 0.0); REQUIRE(r_d < 5.0); - + // Get UInt size_t r_ui = rnd.GetUInt(emp::Range(0,5)); REQUIRE(r_ui < 5); - + // Get Int int r_i = rnd.GetInt(emp::Range(-5,5)); REQUIRE(r_i >= -5); REQUIRE(r_i < 5); - + // Get UInt64 uint64_t ui64 = rnd.GetUInt64(100); REQUIRE(ui64 < 100); ui64 = rnd.GetUInt64(100000000000); REQUIRE(ui64 < 100000000000); - + // Values are consistent when random seeded with 5 double rndNormal = rnd.GetRandNormal(5.0, 0.1); REQUIRE( abs(rndNormal - 5.0) < 0.5 ); - + REQUIRE(rnd.GetRandPoisson(1.0, 0.9) == 1.0); size_t b1_result = rnd.GetRandBinomial(3000, 0.1); @@ -69,9 +69,13 @@ TEST_CASE("Test Random", "[tools]") size_t b2_result = rnd.GetRandBinomial(100, 0.3); REQUIRE(b2_result > 15); REQUIRE(b2_result < 50); - + emp::RandomStdAdaptor randomStd(rnd); REQUIRE(randomStd(4) == 3); + + REQUIRE(rnd.GetRandGeometric(1) == 1); + REQUIRE(rnd.GetRandGeometric(0) == std::numeric_limits::infinity()); + // REQUIRE(rnd.GetRandGeometric(.25) == 8); } TEST_CASE("Another Test random", "[tools]") diff --git a/tests/tools/map_utils.cc b/tests/tools/map_utils.cc index d4d5e0b54b..cbf32ad076 100644 --- a/tests/tools/map_utils.cc +++ b/tests/tools/map_utils.cc @@ -3,6 +3,7 @@ #include "third-party/Catch/single_include/catch2/catch.hpp" #include "tools/map_utils.h" +#include "tools/vector_utils.h" #include #include @@ -36,4 +37,13 @@ TEST_CASE("Test map_utils", "[tools]") REQUIRE( emp::Find(test_123, "0", "nothing") == "nothing" ); REQUIRE( emp::Find(test_123, "1", "nothing") == "1" ); REQUIRE( emp::FindRef(test_123, "1", "nothing") == "1" ); + + // Test Keys + emp::vector key_v2 = emp::Keys(test_map); + REQUIRE( emp::Has(key_v2, 0) ); + REQUIRE( emp::Has(key_v2, 4) ); + REQUIRE( emp::Has(key_v2, 8) ); + REQUIRE( emp::Has(key_v2, 14) ); + REQUIRE( emp::Has(key_v2, 20) ); + REQUIRE( key_v2.size() == 5 ); } diff --git a/tests/tools/math.cc b/tests/tools/math.cc index 6a03d37a26..58f2fb9f8e 100644 --- a/tests/tools/math.cc +++ b/tests/tools/math.cc @@ -20,18 +20,18 @@ TEST_CASE("Test Math", "[tools]") REQUIRE(emp::Mod(5.5, 3.3) == 2.2); REQUIRE(emp::MinRef(0,4,-1,6,52) == -1); REQUIRE(emp::MaxRef(0,4,-1,6,52) == 52); - + REQUIRE(emp::Log10(100.0) == 2); REQUIRE(emp::Ln(emp::E) == 1); REQUIRE( emp::Abs(emp::Ln(emp::Exp(5)) - 5) < 0.01); - + REQUIRE(emp::IntLog2(10) == 3); REQUIRE(emp::CountOnes(15) == 4); REQUIRE(emp::CountOnes(255) == 8); - + unsigned long long large = 0x8000000000000000; REQUIRE(emp::MaskHigh(1) == large); - + REQUIRE(emp::Min(7,3,100,-50,62) == -50); REQUIRE(emp::Max(7,3,100,-50,62) == 100); } @@ -285,4 +285,8 @@ TEST_CASE("Another Test math", "[tools]") REQUIRE(emp::ToRange(12345678, 10, 20) == 20); REQUIRE(emp::ToRange(12345678, 10, 20.1) == 20.1); REQUIRE(emp::ToRange(12345678.0, 10.7, 20.1) == 20.1); + + REQUIRE(emp::Factorial(5) == 120); + REQUIRE(emp::Factorial(3) == 6); + } diff --git a/tests/tools/set_utils.cc b/tests/tools/set_utils.cc index 314142bd12..5d868e9222 100644 --- a/tests/tools/set_utils.cc +++ b/tests/tools/set_utils.cc @@ -41,6 +41,9 @@ TEST_CASE("Test set utils", "[tools]") { comp_set.insert(2); REQUIRE(emp::difference(s1, v1) == comp_set); comp_set.clear(); + comp_set.insert(3); + REQUIRE(emp::difference(v1, s1) == comp_set); + comp_set.clear(); comp_set.insert(1); REQUIRE(emp::intersection(s1, v1) == comp_set); REQUIRE(emp::intersection(v1, s1) == comp_set); @@ -84,4 +87,22 @@ TEST_CASE("Test set utils", "[tools]") { REQUIRE(emp::symmetric_difference(v1, s1) == comp_set); REQUIRE(emp::symmetric_difference(s1, v1) == comp_set); + std::multiset m_set; + m_set.insert(4); + m_set.insert(3); + REQUIRE(emp::Has(m_set, 3)); + REQUIRE(!emp::Has(m_set, 5)); + + std::unordered_multiset um_set; + um_set.insert(6); + um_set.insert(3); + REQUIRE(emp::Has(um_set, 3)); + REQUIRE(!emp::Has(um_set, 7)); + + std::unordered_set u_set; + u_set.insert(9); + u_set.insert(7); + REQUIRE(emp::Has(u_set, 7)); + REQUIRE(!emp::Has(u_set, 4)); + } diff --git a/tests/tools/string_utils.cc b/tests/tools/string_utils.cc index 522983f34b..457431a109 100644 --- a/tests/tools/string_utils.cc +++ b/tests/tools/string_utils.cc @@ -129,7 +129,25 @@ TEST_CASE("Test string_utils", "[tools]") REQUIRE(int_numbers[1] == 2); REQUIRE(int_numbers[2] == 3); - // TODO: try this with more arguments + REQUIRE(emp::is_digits("391830581734")); + REQUIRE(!emp::is_digits("3h91830581734")); + REQUIRE(emp::is_alphanumeric("39adg18af3tj05ykty81734")); + REQUIRE(!emp::is_alphanumeric("39adg18af?3tj05ykty81734")); + REQUIRE(emp::is_literal_char("'f'")); + REQUIRE(emp::is_literal_char("' '")); + REQUIRE(!emp::is_literal_char("f")); + REQUIRE(emp::is_literal_char("'\n'")); + REQUIRE(!emp::is_literal_char("'\\'")); + REQUIRE(emp::from_literal_char("'f'") == 'f'); + REQUIRE(emp::from_literal_char("'\n'") == '\n'); + REQUIRE(emp::is_literal_string("\"He llo!\"")); + REQUIRE(!emp::is_literal_string("\"He\"llo!\"")); + REQUIRE(emp::is_literal_string("\"Hel\n\t\r\\\'lo!\"")); + REQUIRE(emp::is_literal_string("\"Hel\n \t \r \'lo!\"")); + REQUIRE(emp::from_literal_string("\"Hello!\"") == "Hello!"); + REQUIRE(emp::from_literal_string("\"Hel\n \t \r \'lo!\"") == "Hel\n \t \r \'lo!"); + + // TODO: try this with more arguments int one; emp::from_string("1", one); REQUIRE(one == 1); @@ -340,6 +358,10 @@ TEST_CASE("Another Test string_utils", "[tools]") REQUIRE(cat_full == "ABC123"); emp::array test_arr({{ 4, 2, 5 }}); REQUIRE(emp::to_string(test_arr) == "[ 4 2 5 ]"); + REQUIRE(emp::count(emp::to_string(test_arr), ' ') == 4); + REQUIRE(emp::join(emp::vector({17,18,19}), ",") == "17,18,19"); + REQUIRE(emp::join(emp::vector({}), ",") == ""); + REQUIRE(emp::join(emp::vector({17}), ",") == "17"); // tests adapted from https://stackoverflow.com/questions/5288396/c-ostream-out-manipulation/5289170#5289170 std::string els[] = { "aap", "noot", "mies" }; @@ -394,6 +416,7 @@ TEST_CASE("Another Test string_utils", "[tools]") string_v.push_back("four"); REQUIRE( emp::to_english_list(string_v) == "one, two, three, and four" ); + REQUIRE( emp::to_quoted_list(string_v) == "'one', 'two', 'three', and 'four'"); emp::string_vec_t quoted_strings = emp::quote_strings(string_v); @@ -410,4 +433,5 @@ TEST_CASE("Another Test string_utils", "[tools]") REQUIRE( quoted_strings[0] == "([{}])" ); REQUIRE( quoted_strings[2] == "([{}])" ); + REQUIRE( emp::to_titlecase("Harry Potter and the pRisoner of azkaban") == "Harry Potter And The Prisoner Of Azkaban"); } diff --git a/tests/tools/vector_utils.cc b/tests/tools/vector_utils.cc index 8237a75d07..e6ae31ba6e 100644 --- a/tests/tools/vector_utils.cc +++ b/tests/tools/vector_utils.cc @@ -67,10 +67,10 @@ TEST_CASE("Test vector_utils", "[tools]") REQUIRE(v_d2.at(1) == 10.0); REQUIRE(v_d2.at(2) == 5.0); - // Heapify should change nothing - emp::vector v_d3 = emp::Slice(v_d2, 0, 2); - emp::Heapify(v_d3); - REQUIRE(v_d2.at(0) == 20.0); + // Heapify should change nothing + emp::vector v_d3 = emp::Slice(v_d2, 0, 2); + emp::Heapify(v_d3); + REQUIRE(v_d2.at(0) == 20.0); REQUIRE(v_d2.at(1) == 10.0); // HeapExtract @@ -82,17 +82,68 @@ TEST_CASE("Test vector_utils", "[tools]") emp::HeapInsert(v_d2, 35.0); REQUIRE(v_d2.at(0) == 35.0); - - - - - - + emp::vector range_vec = emp::NRange(4, 7); + REQUIRE(range_vec[0] == 4); + REQUIRE(range_vec[1] == 5); + REQUIRE(range_vec[2] == 6); + REQUIRE(range_vec.size() == 3); + + // RemoveDuplicates + range_vec.push_back(4); + REQUIRE(range_vec.size() == 4); + range_vec = emp::RemoveDuplicates(range_vec); + REQUIRE(range_vec.size() == 3); + + // Flatten + emp::vector> nested_v = {{2,1,6}, {4,5,3}}; + emp::vector flattened_v = emp::Flatten(nested_v); + REQUIRE(flattened_v[0] == 2); + REQUIRE(flattened_v[1] == 1); + REQUIRE(flattened_v[2] == 6); + REQUIRE(flattened_v[3] == 4); + REQUIRE(flattened_v[4] == 5); + REQUIRE(flattened_v[5] == 3); + + // FindMin and FindMax + REQUIRE(emp::FindMax(flattened_v) == 6); + REQUIRE(emp::FindMin(flattened_v) == 1); + + // Concat + nested_v = emp::Concat(nested_v, range_vec); + REQUIRE(nested_v[0][0] == 2); + REQUIRE(nested_v[0][1] == 1); + REQUIRE(nested_v[0][2] == 6); + REQUIRE(nested_v[1][0] == 4); + REQUIRE(nested_v[1][1] == 5); + REQUIRE(nested_v[1][2] == 3); + REQUIRE(nested_v[2][0] == 4); + REQUIRE(nested_v[2][1] == 5); + REQUIRE(nested_v[2][2] == 6); + + // FindEval + std::function is_even = [](int i){return ((i % 2) == 0);}; + REQUIRE(emp::FindEval(flattened_v, is_even, 1) == 2); + + // Scale + emp::Scale(range_vec, 2); + REQUIRE(range_vec[0] == 8); + REQUIRE(range_vec[1] == 10); + REQUIRE(range_vec[2] == 12); + + // Heapify on a larger vector + emp::Heapify(flattened_v); + REQUIRE(flattened_v.at(0) == 6); + REQUIRE(flattened_v.at(0) > flattened_v.at(1)); + REQUIRE(flattened_v.at(0) > flattened_v.at(2)); + REQUIRE(flattened_v.at(1) > flattened_v.at(3)); + REQUIRE(flattened_v.at(1) > flattened_v.at(4)); + REQUIRE(flattened_v.at(2) > flattened_v.at(5)); } TEST_CASE("Another Test vector utils", "[tools]") { emp::vector v1({6,2,5,1,3}); + emp::vector v2({7,6,7,1,7}); emp::Sort(v1); REQUIRE(v1 == emp::vector({1,2,3,5,6})); REQUIRE(emp::FindValue(v1, 3) == 2); @@ -101,6 +152,8 @@ TEST_CASE("Another Test vector utils", "[tools]") { REQUIRE(!emp::Has(v1, 4)); REQUIRE(emp::Product(v1) == 180); REQUIRE(emp::Slice(v1,1,3) == emp::vector({2,3})); + REQUIRE(emp::Count(v1, 2) == 1); + REQUIRE(emp::Count(v2, 7) == 3); // Test handling vector-of-vectors. using vv_int_t = emp::vector< emp::vector< int > >;