Recomdination

src/events.pxi

@@ -7,6 +7,7 @@ DEF SUSCCHANGE = 4
 DEF MIGRATION = 5
 DEF MULTITYPE = 6
 
+#To delete this
 cdef class Event:
     cdef:
         Py_ssize_t type_, haplotype, population, newHaplotype, newPopulation
@@ -26,6 +27,7 @@ cdef class Events:
         Py_ssize_t size, ptr
 
         Py_ssize_t[::1] types, haplotypes, populations, newHaplotypes, newPopulations
+        # Py_ssize_t[:,::1] events
         double[::1] times
 
     def __init__(self):
@@ -35,7 +37,9 @@ cdef class Events:
     @cython.boundscheck(False)
     @cython.wraparound(False)
     cdef void AddEvent(self, double time_, Py_ssize_t type_, Py_ssize_t haplotype, Py_ssize_t population, Py_ssize_t newHaplotype, Py_ssize_t newPopulation):
-        self.times[ self.ptr ] = time_
+    # cdef void AddEvents(self, double time_, Py_ssize_t[::1] new_event):
+        self.times[self.ptr] = time_
+        # self.events[self.ptr] = new_event
         self.types[ self.ptr ] = type_
         self.haplotypes[ self.ptr ] = haplotype
         self.populations[ self.ptr ] = population
@@ -48,25 +52,29 @@ cdef class Events:
     cdef Event GetEvent(self, Py_ssize_t e_id):
         ev = Event( self.times[ e_id ], self.types[ e_id ], self.haplotypes[ e_id ], self.populations[ e_id ], self.newHaplotypes[ e_id ], self.newPopulations[ e_id ])
         return( ev )
+        # return self.times[e_id], self.events[e_id]
 
     cdef void CreateEvents(self, Py_ssize_t iterations):
         if self.ptr == 0:
             self.size += iterations
             self.times = np.zeros(self.size, dtype=float)
+            # self.events = np.zeros((self.size, 7), dtype=np.int64)
             self.types = np.zeros(self.size, dtype=np.int64)
             self.haplotypes = np.zeros(self.size, dtype=np.int64)
             self.populations = np.zeros(self.size, dtype=np.int64)
             self.newHaplotypes = np.zeros(self.size, dtype=np.int64)
             self.newPopulations = np.zeros(self.size, dtype=np.int64)
         elif iterations + self.ptr - self.size > 0:
                 self.times = np.concatenate((self.times, np.zeros(iterations + self.ptr - self.size, dtype=float)))
+                # self.events = np.concatenate((self.events, np.zeros((iterations + self.ptr - self.size, 7), dtype=np.int64)))
                 self.types = np.concatenate((self.types, np.zeros(iterations + self.ptr - self.size, dtype=np.int64)))
                 self.haplotypes = np.concatenate((self.haplotypes, np.zeros(iterations + self.ptr - self.size, dtype=np.int64)))
                 self.populations = np.concatenate((self.populations, np.zeros(iterations + self.ptr - self.size, dtype=np.int64)))
                 self.newHaplotypes = np.concatenate((self.newHaplotypes, np.zeros(iterations + self.ptr - self.size, dtype=np.int64)))
                 self.newPopulations = np.concatenate((self.newPopulations, np.zeros(iterations + self.ptr - self.size, dtype=np.int64)))
                 self.size = iterations + self.ptr
 
+#To delete this
 cdef class multiEvent:
     cdef:
         Py_ssize_t num, type_, haplotype, population, newHaplotype, newPopulation
@@ -107,6 +115,7 @@ cdef class multiEvents:
         Py_ssize_t size, ptr
 
         Py_ssize_t[::1] num, types, haplotypes, populations, newHaplotypes, newPopulations
+        # Py_ssize_t[:,::1] events
         double[::1] times
 
     def __init__(self):
@@ -116,8 +125,10 @@ cdef class multiEvents:
     @cython.boundscheck(False)
     @cython.wraparound(False)
     cdef void AddEvents(self, Py_ssize_t num, double time_, Py_ssize_t type_, Py_ssize_t haplotype, Py_ssize_t population, Py_ssize_t newHaplotype, Py_ssize_t newPopulation):
-        self.num[ self.ptr ] = num
+    # cdef void AddEvents(self, double time_, Py_ssize_t[::1] new_event):
         self.times[ self.ptr ] = time_
+        # self.events[self.ptr] = new_event
+        self.num[ self.ptr ] = num
         self.types[ self.ptr ] = type_
         self.haplotypes[ self.ptr ] = haplotype
         self.populations[ self.ptr ] = population
@@ -130,20 +141,24 @@ cdef class multiEvents:
     cdef multiEvent GetEvent(self, Py_ssize_t e_id):
         ev = multiEvent( self.num[ e_id ], self.times[ e_id ], self.types[ e_id ], self.haplotypes[ e_id ], self.populations[ e_id ], self.newHaplotypes[ e_id ], self.newPopulations[ e_id ])
         return( ev )
+        # return self.times[e_id], self.events[e_id]
 
     cdef void CreateEvents(self, Py_ssize_t iterations):
         if self.ptr == 0:
             self.size = iterations
-            self.num = np.zeros(self.size, dtype=np.int64)
             self.times = np.zeros(self.size, dtype=float)
+            # self.events = np.zeros((self.size, 8), dtype=np.int64)
+            self.num = np.zeros(self.size, dtype=np.int64)
             self.types = np.zeros(self.size, dtype=np.int64)
             self.haplotypes = np.zeros(self.size, dtype=np.int64)
             self.populations = np.zeros(self.size, dtype=np.int64)
             self.newHaplotypes = np.zeros(self.size, dtype=np.int64)
             self.newPopulations = np.zeros(self.size, dtype=np.int64)
         else:
-            self.num = np.concatenate((self.num, np.zeros(iterations + self.ptr - self.size, dtype=np.int64)))
+        # elif iterations + self.ptr - self.size > 0:
             self.times = np.concatenate((self.times, np.zeros(iterations + self.ptr - self.size, dtype=float)))
+            # self.events = np.concatenate((self.events, np.zeros((iterations + self.ptr - self.size, 8), dtype=np.int64)))
+            self.num = np.concatenate((self.num, np.zeros(iterations + self.ptr - self.size, dtype=np.int64)))
             self.types = np.concatenate((self.types, np.zeros(iterations + self.ptr - self.size, dtype=np.int64)))
             self.haplotypes = np.concatenate((self.haplotypes, np.zeros(iterations + self.ptr - self.size, dtype=np.int64)))
             self.populations = np.concatenate((self.populations, np.zeros(iterations + self.ptr - self.size, dtype=np.int64)))

src/models.pxi

@@ -1,3 +1,5 @@
+from libcpp.map cimport map
+
 cdef class Mutations:
     cdef:
         vector[Py_ssize_t] nodeId, AS, DS, site
@@ -68,29 +70,30 @@ cdef class Lockdowns:
 
 cdef class Recombination:
     cdef:
-        vector[Py_ssize_t] positions
-        # vector[Py_ssize_t] positions, hi1s, hi2s
-        # vector[double] times
+        map[double, Py_ssize_t] positions, time
+        map[Py_ssize_t, Py_ssize_t] branch
+        vector[Py_ssize_t] haplotypes
 
-        vector[Py_ssize_t] idevents, his, hi2s, nhis, posRecombs
+        Py_ssize_t counter
 
     def __init__(self):
-        pass
+        self.counter = 0
 
     @cython.boundscheck(False)
     @cython.wraparound(False)
-    cdef void AddRecombination_forward(self, Py_ssize_t Id, Py_ssize_t hi, Py_ssize_t hi2, Py_ssize_t posRecomb, Py_ssize_t nhi):
-        self.idevents.push_back(Id)
-        self.his.push_back(hi)
-        self.hi2s.push_back(hi2)
-        self.nhis.push_back(nhi)
-        self.posRecombs.push_back(posRecomb)
+    cdef void AddRecombination(self, double time, Py_ssize_t position, Py_ssize_t haplotype):
+        self.positions[time] = position
+        self.haplotypes.push_back(haplotype)
+
+    cdef Py_ssize_t get_next_recombination(self):
+        self.counter -= 1
+        return self.haplotypes[self.counter]
+
+    cdef void AddBranch(self, double time, Py_ssize_t child, Py_ssize_t parent):
+        self.time[time] = parent
+        self.branch[child] = parent
+
+    cdef void UpdateBranch(self):
+        self.branch.clear()
+        self.time.clear()
 
-    @cython.boundscheck(False)
-    @cython.wraparound(False)
-    cdef void AddRecombination(self, Py_ssize_t position):
-    # cdef void AddRecombination(self, Py_ssize_t position, Py_ssize_t hi1, Py_ssize_t hi2, double time):
-        self.positions.push_back(position)
-        # self.hi1s.push_back(hi1)
-        # self.hi2s.push_back(hi2)
-        # self.times.push_back(time)

testing/check_simulator.py

@@ -80,17 +80,17 @@ def check_data(num):
 check_data(name)
 
 #5 model
-# simulator = VGsim.Simulator(populations_number=2, seed=my_seed)
-
-# name = 5
-# simulator.set_population_size(2000000)
-# simulator.set_contact_density(1.3, population=0)
-# simulator.set_contact_density(0.8, population=1)
-# simulator.set_migration_probability(probability=0.01, source=0, target=1)
-# simulator.set_migration_probability(probability=0.005, source=1, target=0)
-# simulator.simulate(iterations)
-# simulator.export_chain_events('test_' + str(name))
-# check_data(name)
+simulator = VGsim.Simulator(populations_number=2, seed=my_seed)
+
+name = 5
+simulator.set_population_size(2000000)
+simulator.set_contact_density(1.3, population=0)
+simulator.set_contact_density(0.8, population=1)
+simulator.set_migration_probability(probability=0.01, source=0, target=1)
+simulator.set_migration_probability(probability=0.005, source=1, target=0)
+simulator.simulate(iterations)
+simulator.export_chain_events('test_' + str(name))
+check_data(name)
 
 #6 model
 simulator = VGsim.Simulator(populations_number=3, seed=my_seed)

testing/example.py

@@ -16,7 +16,7 @@
 #Set mutation rates
 mutation_rate=0.00003
 substitution_weights=[1,1,1,2]#ATCG
-simulator.set_mutation_rate(mutation_rate)
+simulator.set_mutation_rate(mutation_rate, haplotype='G*', mutation=1)
 simulator.set_mutation_probabilities(substitution_weights)
 simulator.set_mutation_rate(3*mutation_rate, haplotype='G*', mutation=1)
 
@@ -69,13 +69,13 @@
 	os.mkdir('example_output')
 os.chdir('example_output')
 
-#Generate a plot with epidemiological curves
-simulator.add_plot_infectious(population=0, haplotype='AA', step_num=100)
-simulator.add_plot_infectious(population=1, haplotype='AA', step_num=100)
-simulator.add_plot_infectious(population=2, haplotype='AA', step_num=100)
-simulator.add_title(name="AA haplotype")
-simulator.add_legend()
-simulator.plot('plot')
+# #Generate a plot with epidemiological curves
+# simulator.add_plot_infectious(population=0, haplotype='AA', step_num=100)
+# simulator.add_plot_infectious(population=1, haplotype='AA', step_num=100)
+# simulator.add_plot_infectious(population=2, haplotype='AA', step_num=100)
+# simulator.add_title(name="AA haplotype")
+# simulator.add_legend()
+# simulator.plot('plot')
 
 #Save the genealogy with mutations and migrations
 simulator.export_newick()

testing/getting_reference.py

@@ -78,7 +78,7 @@
 simulator = VGsim.Simulator(number_of_sites=2, seed=2020)
 
 name = 8
-simulator.set_mutation_rate(probabilities=[1, 0, 0, 1])
+simulator.set_mutation_probabilities([1, 0, 0, 1])
 simulator.simulate(100000)
 simulator.export_chain_events('reference_' + str(name))
 
@@ -89,7 +89,8 @@
 simulator.set_transmission_rate(5.0, haplotype=12)
 simulator.set_recovery_rate(1.5)
 simulator.set_sampling_rate(0.3)
-simulator.set_mutation_rate(0.01, probabilities=[1, 0, 0, 1])
+simulator.set_mutation_rate(0.01)
+simulator.set_mutation_probabilities([1, 0, 0, 1])
 simulator.set_migration_probability(0.01, source=0, target=1)
 simulator.set_migration_probability(0.005, source=2, target=1)
 simulator.set_sampling_multiplier(2.5, population=1)

tests/test_interface.py

@@ -187,11 +187,14 @@ def test_set_coinfection_parameters(recombination, error, text):
 
 #TRANSMISSION RATE
 # @pytest.mark.parametrize('rate, haplotype, condition', [()])
-@pytest.mark.parametrize('rate, haplotype, answer', [(0    , None, [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
-													 (0.002, 'A*', [0.002, 0.002, 0.002, 0.002, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]),
-													 (0.003, 'AT', [2, 0.003, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]),
-													 (0.004, 0   , [0.004, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]),
-													 (0.005, 15  , [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0.005])])
+@pytest.mark.parametrize('rate , haplotype    , answer', 
+						[(0    , None         , [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
+						 (0.002, 'A*'         , [0.002, 0.002, 0.002, 0.002, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]),
+						 (0.003, 'AT'         , [2, 0.003, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]),
+						 (0.004, 0            , [0.004, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]),
+						 (0.005, 15           , [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0.005]),
+						 (0.006, [0, 15, 'T*'], [0.006, 2, 2, 2, 0.006, 0.006, 0.006, 0.006, 2, 2, 2, 2, 2, 2, 2, 0.006]),
+						 (0.007, []           , [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2])])
 def test_set_transmisstion_rate(rate, haplotype, answer):
 	model = Simulator(number_of_sites=2)
 	model.set_transmission_rate(rate=rate, haplotype=haplotype)
@@ -204,7 +207,7 @@ def test_set_transmisstion_rate(rate, haplotype, answer):
 @pytest.mark.parametrize('rate, haplotype, error, text', [(None , None  , TypeError , 'Incorrect type of transmission rate. Type should be int or float.'),
 														  ('str', None  , TypeError , 'Incorrect type of transmission rate. Type should be int or float.'),
 														  (-1   , None  , ValueError, 'Incorrect value of transmission rate. Value should be more or equal 0.'),
-														  (0.01 , [1, 2], TypeError , 'Incorrect type of haplotype. Type should be int or str or None.'),
+														  # (0.01 , [1, 2], TypeError , 'Incorrect type of haplotype. Type should be int or str or None.'),
 														  (0.01 , -1    , IndexError, 'There are no such haplotype!'),
 														  (0.01 , 16    , IndexError, 'There are no such haplotype!'),
 														  (0.01 , 'str' , ValueError, r'Incorrect haplotype. Haplotype should contain only \"A\", \"T\", \"C\", \"G\", \"\*\" and length of haplotype should be equal number of mutations sites.'),
@@ -229,7 +232,7 @@ def test_set_recovery_rate(rate, haplotype, answer):
 @pytest.mark.parametrize('rate, haplotype, error, text', [(None , None  , TypeError , 'Incorrect type of recovery rate. Type should be int or float.'),
 														  ('str', None  , TypeError , 'Incorrect type of recovery rate. Type should be int or float.'),
 														  (-1   , None  , ValueError, 'Incorrect value of recovery rate. Value should be more or equal 0.'),
-														  (0.01 , [1, 2], TypeError , 'Incorrect type of haplotype. Type should be int or str or None.'),
+														  # (0.01 , [1, 2], TypeError , 'Incorrect type of haplotype. Type should be int or str or None.'),
 														  (0.01 , -1    , IndexError, 'There are no such haplotype!'),
 														  (0.01 , 16    , IndexError, 'There are no such haplotype!'),
 														  (0.01 , 'str' , ValueError, r'Incorrect haplotype. Haplotype should contain only \"A\", \"T\", \"C\", \"G\", \"\*\" and length of haplotype should be equal number of mutations sites.'),
@@ -253,7 +256,7 @@ def test_set_sampling_rate(rate, haplotype, answer):
 @pytest.mark.parametrize('rate, haplotype, error, text', [(None , None  , TypeError , 'Incorrect type of sampling rate. Type should be int or float.'),
 														  ('str', None  , TypeError , 'Incorrect type of sampling rate. Type should be int or float.'),
 														  (-1   , None  , ValueError, 'Incorrect value of sampling rate. Value should be more or equal 0.'),
-														  (0.01 , [1, 2], TypeError , 'Incorrect type of haplotype. Type should be int or str or None.'),
+														  # (0.01 , [1, 2], TypeError , 'Incorrect type of haplotype. Type should be int or str or None.'),
 														  (0.01 , -1    , IndexError, 'There are no such haplotype!'),
 														  (0.01 , 16    , IndexError, 'There are no such haplotype!'),
 														  (0.01 , 'str' , ValueError, r'Incorrect haplotype. Haplotype should contain only \"A\", \"T\", \"C\", \"G\", \"\*\" and length of haplotype should be equal number of mutations sites.'),
@@ -282,7 +285,7 @@ def test_set_mutation_rate(rate, haplotype, mutation, answer):
 @pytest.mark.parametrize('rate, haplotype, mutation, error, text', [(None , None  , None , TypeError , 'Incorrect type of mutation rate. Type should be int or float.'),
 																	('str', None  , None , TypeError , 'Incorrect type of mutation rate. Type should be int or float.'),
 																	(-1   , None  , None , ValueError, 'Incorrect value of mutation rate. Value should be more or equal 0.'),
-																	(2    , [1, 2], None , TypeError , 'Incorrect type of haplotype. Type should be int or str or None.'),
+																	# (2    , [1, 2], None , TypeError , 'Incorrect type of haplotype. Type should be int or str or None.'),
 																	(2    , -1    , None , IndexError, 'There are no such haplotype!'),
 																	(2    , 16    , None , IndexError, 'There are no such haplotype!'),
 																	(2    , 'str' , None , ValueError, r'Incorrect haplotype. Haplotype should contain only \"A\", \"T\", \"C\", \"G\", \"\*\" and length of haplotype should be equal number of mutations sites.'),
@@ -319,7 +322,7 @@ def test_set_mutation_probabilities(probabilities, haplotype, mutation, answer):
 																			 ([0, 0, 0, 5]    , None  , None , ValueError, 'Incorrect probabilities list. The sum of three elements without mutation allele should be more 0.'),
 																			 ([None, 0, 0, 5] , None  , None , TypeError , 'Incorrect type of mutation probabilities. Type should be int or float.'),
 																			 (['str', 0, 0, 5], None  , None , TypeError , 'Incorrect type of mutation probabilities. Type should be int or float.'),
-																			 ([2, 3, 4, 5]    , [1, 2], None , TypeError , 'Incorrect type of haplotype. Type should be int or str or None.'),
+																			 # ([2, 3, 4, 5]    , [1, 2], None , TypeError , 'Incorrect type of haplotype. Type should be int or str or None.'),
 																			 ([2, 3, 4, 5]    , -1    , None , IndexError, 'There are no such haplotype!'),
 																			 ([2, 3, 4, 5]    , 16    , None , IndexError, 'There are no such haplotype!'),
 																			 ([2, 3, 4, 5]    , 'str' , None , ValueError, r'Incorrect haplotype. Haplotype should contain only \"A\", \"T\", \"C\", \"G\", \"\*\" and length of haplotype should be equal number of mutations sites.'),
@@ -512,7 +515,7 @@ def test_set_susceptibility_type(susceptibility_type, haplotype, answer):
 																		 ('str', None  , TypeError , 'Incorrect type of susceptibility type. Type should be int.'),
 																		 (-1   , None  , IndexError, 'There are no such susceptibility type!'),
 																		 (4    , None  , IndexError, 'There are no such susceptibility type!'),
-																		 (2    , [1, 2], TypeError , 'Incorrect type of haplotype. Type should be int or str or None.'),
+																		 # (2    , [1, 2], TypeError , 'Incorrect type of haplotype. Type should be int or str or None.'),
 																		 (2    , -1    , IndexError, 'There are no such haplotype!'),
 																		 (2    , 16    , IndexError, 'There are no such haplotype!'),
 																		 (2    , 'str' , ValueError, r'Incorrect haplotype. Haplotype should contain only \"A\", \"T\", \"C\", \"G\", \"\*\" and length of haplotype should be equal number of mutations sites.'),
@@ -538,7 +541,7 @@ def test_set_susceptibility(rate, haplotype, susceptibility_type, answer):
 @pytest.mark.parametrize('rate, haplotype, susceptibility_type, error, text', [(None , None  , None , TypeError , 'Incorrect type of susceptibility rate. Type should be int or float.'),
 																			   ('str', None  , None , TypeError , 'Incorrect type of susceptibility rate. Type should be int or float.'),
 																			   (-1   , None  , None , ValueError, 'Incorrect value of susceptibility rate. Value should be more or equal 0.'),
-																			   (2    , [1, 2], None , TypeError , 'Incorrect type of haplotype. Type should be int or str or None.'),
+																			   # (2    , [1, 2], None , TypeError , 'Incorrect type of haplotype. Type should be int or str or None.'),
 																			   (2    , -1    , None , IndexError, 'There are no such haplotype!'),
 																			   (2    , 16    , None , IndexError, 'There are no such haplotype!'),
 																			   (2    , 'str' , None , ValueError, r'Incorrect haplotype. Haplotype should contain only \"A\", \"T\", \"C\", \"G\", \"\*\" and length of haplotype should be equal number of mutations sites.'),