Adjustment to the way innovation rate is implemented, to remove syste…

…matic bias compared to Ewens 2004 predictions for K_n

analytics/ctmixtures-export-data.py

@@ -85,6 +85,7 @@ def export_population_stats():
     pop_fields = data.mixture_model_stats.pop_columns_to_export()
 
     # adjust the fields for the new summary statistics
+    pop_fields.append('innovation_rate')
     pop_fields.append('slatkin_locus_max')
     pop_fields.append('slatkin_locus_min')
     pop_fields.append('entropy_locus_max')
@@ -108,6 +109,7 @@ def export_population_stats():
         row['model_class_label'] = sample['model_class_label']
         row['num_trait_configurations'] = sample['num_trait_configurations']
         row['configuration_slatkin'] = sample['configuration_slatkin']
+        row['innovation_rate'] = sample['innovation_rate']
 
         # slatkin exact
         slatkin_values = sample['slatkin_exact']

conf/allneutral-wellmixed.json

@@ -16,6 +16,6 @@
     "TRAIT_FACTORY_CLASS" : "ctmixtures.traits.LocusAlleleTraitFactory",
     "TIME_AVERAGING_CLASS" : "pytransmission.aggregation.MoranCumulativeTimeAverager",
     "DYNAMICS_CLASS" : "ctmixtures.dynamics.MoranDynamics",
-    "MODEL_CLASS_LABEL" : "pureneutral"
+    "MODEL_CLASS_LABEL" : "allneutral"
 }
 

setup.py

@@ -15,7 +15,7 @@
       scripts = [
           'admin/ctmixtures-planner.py',
           'admin/ctmixtures-priorsampler-runbuilder.py',
-          'admin/ctmixtures-export-data.py',
+          'analytics/ctmixtures-export-data.py',
           'simulations/sim-ctmixture-notimeaveraging.py',
           'simulations/sim-ctmixture-timeaveraging.py'
       ],

simulations/sim-ctmixture-notimeaveraging.py

@@ -74,10 +74,9 @@ def setup():
     simconfig.num_features = int(args.numloci)
     simconfig.num_traits = int(args.maxinittraits)
     simconfig.popsize = int(args.popsize)
-    mut = pg.moran_mutation_rate_from_theta(float(args.popsize), float(args.innovationrate))
-    simconfig.innovation_rate = float(args.numloci) * mut
+    simconfig.innovation_rate = pg.moran_mutation_rate_from_theta(float(args.popsize), float(args.innovationrate)) * simconfig.num_features
     simconfig.configured_innovation_rate = float(args.innovationrate)
-    log.debug("configured theta = %s, using numloci %s * per-locus mutation rate %s = all-loci innovation rate: %s", args.innovationrate, args.numloci, mut, simconfig.innovation_rate)
+    log.debug("configured theta = %s, using all-loci innovation rate: %s", args.innovationrate,simconfig.innovation_rate)
 
 
     simconfig.maxtime = int(args.simulationendtime)

simulations/sim-ctmixture-timeaveraging.py

@@ -73,10 +73,9 @@ def setup():
     simconfig.num_features = int(args.numloci)
     simconfig.num_traits = int(args.maxinittraits)
     simconfig.popsize = int(args.popsize)
-    mut = pg.moran_mutation_rate_from_theta(float(args.popsize), float(args.innovationrate))
-    simconfig.innovation_rate = float(args.numloci) * mut
+    simconfig.innovation_rate = pg.moran_mutation_rate_from_theta(float(args.popsize), float(args.innovationrate)) * simconfig.num_features
     simconfig.configured_innovation_rate = float(args.innovationrate)
-    log.debug("configured theta = %s, using numloci %s * per-locus mutation rate %s = all-loci innovation rate: %s", args.innovationrate, args.numloci, mut, simconfig.innovation_rate)
+    log.debug("configured theta = %s, using all-loci innovation rate: %s", args.innovationrate,simconfig.innovation_rate)
 
 
     simconfig.maxtime = int(args.simulationendtime)