You can subscribe to this list here.
2011 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
(111) |
Jul
(77) |
Aug
(25) |
Sep
(12) |
Oct
(12) |
Nov
(24) |
Dec
(12) |
---|---|---|---|---|---|---|---|---|---|---|---|---|
2012 |
Jan
(7) |
Feb
(38) |
Mar
(152) |
Apr
(53) |
May
(4) |
Jun
(7) |
Jul
(13) |
Aug
(13) |
Sep
(51) |
Oct
(20) |
Nov
(4) |
Dec
(16) |
2013 |
Jan
(2) |
Feb
(19) |
Mar
(2) |
Apr
(12) |
May
(11) |
Jun
|
Jul
(1) |
Aug
|
Sep
|
Oct
|
Nov
|
Dec
(3) |
2014 |
Jan
(2) |
Feb
(46) |
Mar
(7) |
Apr
(8) |
May
(16) |
Jun
(25) |
Jul
(33) |
Aug
(4) |
Sep
(7) |
Oct
(23) |
Nov
(2) |
Dec
(1) |
2015 |
Jan
(2) |
Feb
(5) |
Mar
|
Apr
(17) |
May
(2) |
Jun
|
Jul
(1) |
Aug
(11) |
Sep
(1) |
Oct
(3) |
Nov
(1) |
Dec
|
2016 |
Jan
(11) |
Feb
(8) |
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
|
From: <bpe...@us...> - 2015-04-21 21:37:47
|
Revision: 4965 http://sourceforge.net/p/simupop/code/4965 Author: bpeng2000 Date: 2015-04-21 21:37:40 +0000 (Tue, 21 Apr 2015) Log Message: ----------- Fix a crash caused by incorrect fix of circular reference in 1.1.3 and 1.1.4 Modified Paths: -------------- trunk/src/utility.cpp trunk/src/utility.h trunk/test/test_03_operator.py Modified: trunk/src/utility.cpp =================================================================== --- trunk/src/utility.cpp 2015-04-14 04:33:16 UTC (rev 4964) +++ trunk/src/utility.cpp 2015-04-21 21:37:40 UTC (rev 4965) @@ -662,11 +662,12 @@ m_elems.push_back(defaultItems[i]); } + const vectorstr & stringList::elems(const GenoStruTrait * trait) const { - if (trait == NULL || ! allAvail()) + if (trait == NULL || !allAvail()) return m_elems; - + if (trait) { if (trait->genoStruIdx() == m_trait) return m_elems; @@ -678,6 +679,7 @@ return m_elems; } + intMatrix::intMatrix(PyObject * obj) : m_elems() { if (obj == NULL) @@ -835,7 +837,7 @@ } -pyFunc::pyFunc(PyObject * func) : m_func(func), m_numArgs(0) +pyFunc::pyFunc(PyObject * func) : m_func(func), m_numArgs(0), m_circular(false) { if (!m_func.isValid()) return; @@ -947,10 +949,9 @@ // check the super class (BaseOperator) because of the SWIG // interface if (PyObject_HasAttrString(self, "apply") && - PyObject_HasAttrString(self, "describe")) { - //std::cerr << "CURCULAR " << std::endl; - Py_DECREF(self); - } + PyObject_HasAttrString(self, "describe")) + //std::cerr << "CIRCULAR " << std::endl; + m_circular = true; Py_DECREF(self); } @@ -988,14 +989,55 @@ } Py_DECREF(co_varnames); // accepting arbitrary number of parameters? - PyObject * co_flag = PyObject_GetAttrString(code, "co_flags"); - DBG_ASSERT(co_flag, SystemError, "Invalid attribute co_flags for a function object"); - m_flags = static_cast<unsigned char>(PyInt_AsLong(co_flag)); - Py_DECREF(co_flag); + /* + PyObject * co_flag = PyObject_GetAttrString(code, "co_flags"); + DBG_ASSERT(co_flag, SystemError, "Invalid attribute co_flags for a function object"); + int flags = static_cast<unsigned char>(PyInt_AsLong(co_flag)); + Py_DECREF(co_flag); + */ Py_DECREF(code); } +pyFunc::pyFunc(const pyFunc & rhs) : m_func(rhs.m_func), m_name(rhs.m_name), m_numArgs(rhs.m_numArgs), + m_args(rhs.m_args), m_circular(rhs.m_circular) +{ + // prevent reference decrease of circular reference caused by removing of + // temporary rhs objects + rhs.m_circular = false; +} + + +pyFunc & pyFunc::operator=(const pyFunc & rhs) +{ + m_func = rhs.m_func; + m_name = rhs.m_name; + m_numArgs = rhs.m_numArgs; + m_args = rhs.m_args; + m_circular = rhs.m_circular; + // prevent reference decrease of circular reference caused by removing of + // temporary rhs objects + rhs.m_circular = false; +} + + +pyFunc::~pyFunc() +{ + // in the case of circular reference, we will have to + // manually reduce reference of self, then self.func in order + // for both of them to be removed from memory + if (m_circular) { + PyObject * self = PyObject_GetAttrString(m_func.object(), SELF_ATTR); + // remove circular reference of self + Py_XDECREF(self); + // remove circular reference of self.func (not sure if this is needed) + Py_XDECREF(m_func.object()); + // remove reference caused by PyObject_GetAttrString. + Py_XDECREF(self); + } +} + + void pyGenerator::set(PyObject * gen) { Py_XDECREF(m_iterator); Modified: trunk/src/utility.h =================================================================== --- trunk/src/utility.h 2015-04-14 04:33:16 UTC (rev 4964) +++ trunk/src/utility.h 2015-04-21 21:37:40 UTC (rev 4965) @@ -319,6 +319,12 @@ public: pyFunc(PyObject * func); + pyFunc(const pyFunc & rhs); + + pyFunc & operator=(const pyFunc & rhs); + + ~pyFunc(); + /// return number of arguments this function accepts. /// This function does not count tuple parameters. size_t numArgs() const @@ -433,7 +439,10 @@ vectorstr m_args; - int m_flags; + // if this function is a member function of a class + // that release of this function requires destroying + // self first. + mutable bool m_circular; }; /// CPPONLY Modified: trunk/test/test_03_operator.py =================================================================== --- trunk/test/test_03_operator.py 2015-04-14 04:33:16 UTC (rev 4964) +++ trunk/test/test_03_operator.py 2015-04-21 21:37:40 UTC (rev 4965) @@ -49,6 +49,16 @@ def __call__(self, x1, x2): return sum([x1, x2]) +class OutputCollector: + def __init__(self): + self.messages = [] + + def func1(self, msg): + self.messages.append(msg) + + def func2(self, msg): + self.messages.append(msg) + class TestOperator(unittest.TestCase): def setUp(self): @@ -217,6 +227,21 @@ PyOutput("func2", output=func2), ], gen=10) + + def testOutputMemberFunc(self): + '''Testing output to a function''' + + simu = Simulator(Population(), rep=5) + o = OutputCollector() + # each replicate + simu.evolve( + matingScheme=CloneMating(), + postOps = [ + PyOutput("func1", output=o.func1), + PyOutput("func2", output=o.func2), + ], gen=10) + self.assertEqual(o.messages, ['func1', 'func2'] * 50) + def testOutputFileHandler(self): '''Testing output to a file handler''' simu = Simulator(Population(), rep=5) @@ -240,7 +265,7 @@ ], gen=10) with open('test1.txt') as test1txt: self.assertEqual(test1txt.read(), 'func2'*50) - # + # runtime error only raise in python 3 simu = Simulator(Population(), rep=5) with open('test1.txt', 'wb') as out: # each replicate This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bpe...@us...> - 2015-04-14 04:33:18
|
Revision: 4964 http://sourceforge.net/p/simupop/code/4964 Author: bpeng2000 Date: 2015-04-14 04:33:16 +0000 (Tue, 14 Apr 2015) Log Message: ----------- Update online doc Modified Paths: -------------- trunk/src/simuPOP_doc.i Modified: trunk/src/simuPOP_doc.i =================================================================== --- trunk/src/simuPOP_doc.i 2015-04-14 04:22:46 UTC (rev 4963) +++ trunk/src/simuPOP_doc.i 2015-04-14 04:33:16 UTC (rev 4964) @@ -318,9 +318,9 @@ %ignore simuPOP::BaseOperator::isActive(ssize_t rep, ssize_t gen) const; -%ignore simuPOP::BaseOperator::infoSize() const; +%ignore simuPOP::BaseOperator::infoSize(const GenoStruTrait *trait=NULL) const; -%ignore simuPOP::BaseOperator::infoField(size_t idx) const; +%ignore simuPOP::BaseOperator::infoField(size_t idx, const GenoStruTrait *trait=NULL) const; %ignore simuPOP::BaseOperator::infoFields() const; @@ -1143,7 +1143,8 @@ Usage: - CombinedParentsChooser(fatherChooser, motherChooser) + CombinedParentsChooser(fatherChooser, motherChooser, + allowSelfing=True) Details: @@ -1155,7 +1156,10 @@ for motherChooser. Although these two parent choosers are supposed to return a father and a mother respectively, the sex of returned parents are not checked so it is possible to return parents with - the same sex using this parents chooser. + the same sex using this parents chooser. This choose by default + allows the selection of the same parents as father and mother + (self-fertilization), unless a parameter allowSelfing is used to + disable it. "; @@ -1559,7 +1563,7 @@ Dumper(genotype=True, structure=True, ancGens=UNSPECIFIED, width=1, max=100, loci=[], output=\">\", begin=0, end=-1, step=1, - at=[], reps=ALL_AVAIL, subPops=ALL_AVAIL, infoFields=[]) + at=[], reps=ALL_AVAIL, subPops=ALL_AVAIL, infoFields=ALL_AVAIL) Details: @@ -1569,10 +1573,13 @@ this operator will only output the first 100 (parameter max) individuals of the present generation (parameter ancGens). All loci will be outputed unless parameter loci are used to specify a - subset of loci. If a list of (virtual) subpopulations are - specified, this operator will only output individuals in these - outputs. Please refer to class BaseOperator for a detailed - explanation for common parameters such as output and stage. + subset of loci. This operator by default output values of all + information fields unless parameter infoFields is used to specify + a subset of info fields to display. If a list of (virtual) + subpopulations are specified, this operator will only output + individuals in these outputs. Please refer to class BaseOperator + for a detailed explanation for common parameters such as output + and stage. "; @@ -2922,7 +2929,7 @@ Usage: - indCompare(idx) + indCompare(idx, reverse=False) "; @@ -3304,7 +3311,7 @@ %ignore simuPOP::Individual::swap(Individual &ind, bool swapContent=true); -%ignore simuPOP::Individual::display(ostream &out, int width=1, const vectoru &loci=vectoru()); +%ignore simuPOP::Individual::display(ostream &out, int width=1, const vectoru &loci=vectoru(), const vectoru &infoIdx=vectoru()); %feature("docstring") simuPOP::IndividualIterator " @@ -6487,13 +6494,13 @@ Usage: - x.sortIndividuals(infoFields) + x.sortIndividuals(infoFields, reverse=False) Details: Sort individuals according to values at specified information fields (infoFields). Individuals will be sorted at an increasing - order. + order unless reverse is set to true. "; @@ -10916,7 +10923,7 @@ %ignore simuPOP::stringList::pushback(const string &str); -%ignore simuPOP::stringList::elems() const; +%ignore simuPOP::stringList::elems(const GenoStruTrait *trait=NULL) const; %feature("docstring") simuPOP::stringMatrix " This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bpe...@us...> - 2015-04-14 04:22:49
|
Revision: 4963 http://sourceforge.net/p/simupop/code/4963 Author: bpeng2000 Date: 2015-04-14 04:22:46 +0000 (Tue, 14 Apr 2015) Log Message: ----------- Add paramter reverse to function Population.sortIndividuals Modified Paths: -------------- trunk/src/individual.h trunk/src/population.cpp trunk/src/population.h trunk/test/test_02_population.py Modified: trunk/src/individual.h =================================================================== --- trunk/src/individual.h 2015-04-03 21:50:47 UTC (rev 4962) +++ trunk/src/individual.h 2015-04-14 04:22:46 UTC (rev 4963) @@ -779,15 +779,15 @@ /** CPPONLY - * A class used to compare two individuals by an information field. + * A class used to compare two individuals by one or more information fields. */ class indCompare { public: // accept the index to an information field - indCompare(size_t idx) : m_fields(1, idx) {} + indCompare(size_t idx, bool reverse=false) : m_fields(1, idx), m_reverse(reverse) {} - indCompare(const vectoru & idx) : m_fields(idx) {} + indCompare(const vectoru & idx, bool reverse=false) : m_fields(idx), m_reverse(reverse) {} bool operator()(const Individual & lhs, const Individual & rhs) const { @@ -799,7 +799,7 @@ continue; else // if < or >, we have a comparison - return v1 < v2; + return m_reverse ? v1 > v2 : v1 < v2; } // they are actually equal, but we are inplementing operator < return false; @@ -808,6 +808,8 @@ private: vectoru m_fields; + + bool m_reverse; }; Modified: trunk/src/population.cpp =================================================================== --- trunk/src/population.cpp 2015-04-03 21:50:47 UTC (rev 4962) +++ trunk/src/population.cpp 2015-04-14 04:22:46 UTC (rev 4963) @@ -1037,7 +1037,7 @@ } -void Population::sortIndividuals(const stringList & infoList) +void Population::sortIndividuals(const stringList & infoList, bool reverse) { const vectorstr & infoFields = infoList.elems(); @@ -1047,7 +1047,7 @@ for (size_t i = 0; i < infoFields.size(); ++i) fields[i] = infoIdx(infoFields[i]); for (size_t sp = 0; sp < numSubPop(); ++sp) - parallelSort(rawIndBegin(sp), rawIndEnd(sp), indCompare(fields)); + parallelSort(rawIndBegin(sp), rawIndEnd(sp), indCompare(fields, reverse=reverse)); setIndOrdered(false); } Modified: trunk/src/population.h =================================================================== --- trunk/src/population.h 2015-04-03 21:50:47 UTC (rev 4962) +++ trunk/src/population.h 2015-04-14 04:22:46 UTC (rev 4963) @@ -1110,9 +1110,9 @@ /** Sort individuals according to values at specified information * fields (\e infoFields). Individuals will be sorted at an increasing - * order. + * order unless \e reverse is set to \c true. */ - void sortIndividuals(const stringList & infoFields); + void sortIndividuals(const stringList & infoFields, bool reverse=false); /** Rearrange individuals to their new subpopulations according to their * integer values at information field \e field (value returned by Modified: trunk/test/test_02_population.py =================================================================== --- trunk/test/test_02_population.py 2015-04-03 21:50:47 UTC (rev 4962) +++ trunk/test/test_02_population.py 2015-04-14 04:22:46 UTC (rev 4963) @@ -1170,11 +1170,20 @@ 'Testing Population::sortIndividuals(infoFields)' pop = self.getPop(size=[1000, 2000], infoFields=['a', 'b']) initInfo(pop, lambda: random.randint(1, 5), infoFields=['a', 'b']) - pop.sortIndividuals('a') + print(pop.indInfo('a')) for sp in range(2): for i in range(1, pop.subPopSize(sp)): self.assertTrue(pop.individual(i-1, sp).a <= pop.individual(i, sp).a) self.assertTrue(pop.individual(999).a > pop.individual(0, 1).a) + # sorting in reverse order + initInfo(pop, lambda: random.randint(1, 5), infoFields=['a', 'b']) + pop.sortIndividuals('a', reverse=True) + for sp in range(2): + for i in range(1, pop.subPopSize(sp)): + self.assertTrue(pop.individual(i-1, sp).a >= pop.individual(i, sp).a) + self.assertTrue(pop.individual(999).a < pop.individual(0, 1).a) + + def testAddInfoFields(self): 'Testing Population::addInfoFields(fields, init=0)' This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bpe...@us...> - 2015-04-03 21:50:55
|
Revision: 4962 http://sourceforge.net/p/simupop/code/4962 Author: bpeng2000 Date: 2015-04-03 21:50:47 +0000 (Fri, 03 Apr 2015) Log Message: ----------- More testing, fix ScatterPlotter, and add rpy2 support to trajectory simulation plots Modified Paths: -------------- trunk/doc/userGuide.py trunk/src/plotter.py trunk/src/utils.py Modified: trunk/doc/userGuide.py =================================================================== --- trunk/doc/userGuide.py 2015-04-03 21:21:00 UTC (rev 4961) +++ trunk/doc/userGuide.py 2015-04-03 21:50:47 UTC (rev 4962) @@ -4728,10 +4728,8 @@ #end_file #begin_file log/forwardTrajectory.py -#begin_ignore import simuOpt simuOpt.setOptions(quiet=True, plotter='matplotlib') -#end_ignore import simuPOP as sim #begin_ignore sim.setRNG(seed=12345) @@ -4771,7 +4769,7 @@ #begin_file log/backTrajectory.py #begin_ignore import simuOpt -simuOpt.setOptions(quiet=True, plotter='matplotlib') +simuOpt.setOptions(quiet=True, plotter='rpy2') #end_ignore import simuPOP as sim #begin_ignore @@ -4791,11 +4789,11 @@ traj = simulateBackwardTrajectory(N=Nt, fitness=fitness, nLoci=2, endGen=1000, endFreq=[0.1, 0.2]) # matplotlib syntax -traj.plot('log/backTrajectory.png', set_ylim_top=0.3, set_ylim_bottom=0, - plot_c_loc=['r', 'b'], set_title_label='Simulated Trajectory (backward-time)') +#traj.plot('log/backTrajectory.png', set_ylim_top=0.3, set_ylim_bottom=0, +# plot_c_loc=['r', 'b'], set_title_label='Simulated Trajectory (backward-time)') # rpy syntax -# traj.plot('log/backTrajectory.png', plot_ylim=[0, 0.3], plot_xlim=[0, 1000], -# col_loc=['red', 'blue'], plot_main='Simulated Trajectory (backward-time)') +traj.plot('log/backTrajectory.png', plot_ylim=[0, 0.3], plot_xlim=[0, 1000], + col_loc=['red', 'blue'], plot_main='Simulated Trajectory (backward-time)') print('Trajectory simulated with length %s ' % len(traj.traj)) pop = sim.Population(size=Nt(0), loci=[1]*2) @@ -5125,7 +5123,7 @@ #begin_file log/varPlotByDim.py #begin_ignore import simuOpt -simuOpt.setOptions(quiet=True, plotter='matplotlib') +simuOpt.setOptions(quiet=True, plotter='rpy2') #end_ignore import simuPOP as sim #begin_ignore @@ -5159,35 +5157,33 @@ matingScheme=sim.RandomMating(), postOps=[ sim.Stat(alleleFreq=range(4)), - # rpy syntax + # rpy and rpy2 syntax + VarPlotter('[alleleFreq[x][0] for x in range(4)]', byDim=True, + update=10, saveAs='log/rpy_byDim.png', + legend=['Replicate %d' % x for x in range(3)], + ylab='Allele frequency', + ylim=[0, 1], + main_dim=['Genetic drift, freq=%.1f' % ((x+1)*0.10) for x in range(4)], + col_rep=['red', 'blue', 'black'], + lty_rep=[1, 2, 3], + # the default png dimension is 800x600 + dev_print_width=600, dev_print_height=500, + # do not draw axes in r.plot, leaving the job to drawFrame + plot_axes=False, + plotHook = rpy_drawFrame, + ), + # matplot lib syntax #VarPlotter('[alleleFreq[x][0] for x in range(4)]', byDim=True, - # update=10, saveAs='log/rpy_byDim.png', + # update=10, saveAs='log/varplot_byDim.png', # legend=['Replicate %d' % x for x in range(3)], - # ylab='Allele frequency', - # ylim=[0, 1], - # main_dim=['Genetic drift, freq=%.1f' % ((x+1)*0.10) for x in range(4)], - # col_rep=['red', 'blue', 'black'], - # lty_rep=[1, 2, 3], - # # the default png dimension is 800x600 - # dev_print_width=600, dev_print_height=500, - # # do not draw axes in r.plot, leaving the job to drawFrame - # plot_axes=False, - # # plot frame, grid etc after each r.plot call - # plotHook = rpy_drawFrame, + # set_ylabel_ylabel='Allele frequency', + # set_ylim_bottom=0, set_ylim_top=1, + # set_title_label_dim=['Genetic drift, freq=%.1f' % ((x+1)*0.10) for x in range(4)], + # plot_c_rep=['red', 'blue', 'black'], + # plot_linestyle_rep=['-', '-.', ':'], + # figure_figsize=(10,8), + # plotHook = mat_drawFrame, #), - # matplot lib syntax - VarPlotter('[alleleFreq[x][0] for x in range(4)]', byDim=True, - update=10, saveAs='log/varplot_byDim.png', - legend=['Replicate %d' % x for x in range(3)], - set_ylabel_ylabel='Allele frequency', - set_ylim_bottom=0, set_ylim_top=1, - set_title_label_dim=['Genetic drift, freq=%.1f' % ((x+1)*0.10) for x in range(4)], - plot_c_rep=['red', 'blue', 'black'], - plot_linestyle_rep=['-', '-.', ':'], - figure_figsize=(10,8), - # plot frame, grid etc after each r.plot call - plotHook = mat_drawFrame, - ), ], gen=100 ) @@ -5196,7 +5192,7 @@ #begin_file log/ScatterPlotter.py #begin_ignore import simuOpt -simuOpt.setOptions(quiet=True, plotter='matplotlib') +simuOpt.setOptions(quiet=True, plotter='rpy2') #end_ignore import simuPOP as sim #begin_ignore @@ -5229,25 +5225,25 @@ sim.PyTagger(passInfo)]), postOps=[ # rpy syntax + ScatterPlotter(['x', 'y'], + saveAs = 'log/ScatterPlotter.png', + subPops = [(0, 0), (0, 1), (0, 2), (0, 3), (0, 4)], + ylim = [0, 1.2], + main = "!'Ancestry distribution of individuals at generation %d' % gen", + legend = ['anc < 0.2', '0.2 <= anc < 0.4', '0.4 <= anc < 0.6', + '0.6 <= anc < 0.8', '0.8 <= anc'], + plot_axes = False, + par_mar = [0, 0, 2, 0], + ), + # matplotlib syntax #ScatterPlotter(['x', 'y'], # saveAs = 'log/ScatterPlotter.png', # subPops = [(0, 0), (0, 1), (0, 2), (0, 3), (0, 4)], - # ylim = [0, 1.2], - # main = "!'Ancestry distribution of individuals at generation %d' % gen", + # set_ylim_bottom = 0, set_ylim_top=1.2, + # set_title_label = "!'Ancestry distribution of individuals at generation %d' % gen", # legend = ['anc < 0.2', '0.2 <= anc < 0.4', '0.4 <= anc < 0.6', # '0.6 <= anc < 0.8', '0.8 <= anc'], - # plot_axes = False, - # par_mar = [0, 0, 2, 0], #), - # matplotlib syntax - ScatterPlotter(['x', 'y'], - saveAs = 'log/ScatterPlotter.png', - subPops = [(0, 0), (0, 1), (0, 2), (0, 3), (0, 4)], - set_ylim_bottom = 0, set_ylim_top=1.2, - set_title_label = "!'Ancestry distribution of individuals at generation %d' % gen", - legend = ['anc < 0.2', '0.2 <= anc < 0.4', '0.4 <= anc < 0.6', - '0.6 <= anc < 0.8', '0.8 <= anc'], - ), ], gen = 5, @@ -5257,7 +5253,7 @@ #begin_file log/HistPlotter.py #begin_ignore import simuOpt -simuOpt.setOptions(quiet=True, plotter='rpy') +simuOpt.setOptions(quiet=True, plotter='rpy2') #end_ignore import simuPOP as sim #begin_ignore Modified: trunk/src/plotter.py =================================================================== --- trunk/src/plotter.py 2015-04-03 21:21:00 UTC (rev 4961) +++ trunk/src/plotter.py 2015-04-03 21:50:47 UTC (rev 4962) @@ -584,7 +584,7 @@ self.gen = [] self.data = [] # when apply is called, self._rpy_plot is called. - PyOperator.__init__(self, func=self._rpy_plot if plotter == 'rpy' else self._mat_plot, + PyOperator.__init__(self, func=self._rpy_plot if 'rpy' in plotter else self._mat_plot, begin=begin, end=end, step=step, at=at, reps=reps, subPops=ALL_AVAIL, infoFields=[]) @@ -1017,7 +1017,7 @@ self.preHook = preHook self.postHook = postHook self.subPops = subPops - if plotter == 'rpy': + if 'rpy' in plotter: self.args = DerivedArgs( defaultFuncs = ['plot', 'points'], allFuncs = ['par', 'plot', 'points', 'dev_print', 'legend'], @@ -1045,7 +1045,7 @@ ) if len(self.subPops) > 1: - if plotter == 'rpy': + if 'rpy' in plotter: self.args.addDefault( pch_sp = range(1, len(self.subPops) + 1), col_sp = r.rainbow(len(self.subPops))) Modified: trunk/src/utils.py =================================================================== --- trunk/src/utils.py 2015-04-03 21:21:00 UTC (rev 4961) +++ trunk/src/utils.py 2015-04-03 21:50:47 UTC (rev 4962) @@ -856,19 +856,25 @@ if plotter is None: try: import rpy - self._rpy_plot(filename, **kwargs) + self._rpy_plot(filename, rpy2=False, **kwargs) except ImportError: try: - import matplotlib - self._mat_plot(filename, **kwargs) + import rpy2 + self._rpy_plot(filename, rpy2=True, **kwargs) except ImportError: - sys.stderr.write('Failed to draw figure: rpy or matplotlib is not available.') + try: + import matplotlib + self._mat_plot(filename, **kwargs) + except ImportError: + sys.stderr.write('Failed to draw figure: rpy or matplotlib is not available.') elif plotter == 'rpy': - self._rpy_plot(filename, **kwargs) + self._rpy_plot(filename, rpy2=False, **kwargs) + elif plotter == 'rpy2': + self._rpy_plot(filename, rpy2=True, **kwargs) else: self._mat_plot(filename, **kwargs) - def _rpy_plot(self, filename, **kwargs): + def _rpy_plot(self, filename, rpy2, **kwargs): import plotter # args = plotter.DerivedArgs( @@ -926,8 +932,11 @@ **args.getArgs('lines', None, sp=sp, loc=loc)) # plotter.saveFigure(**args.getArgs('dev_print', None, file=filename)) - plotter.r.dev_off() - + if rpy2: + plotter.r['dev.off']() + else: + plotter.r.dev_off() + def _mat_plot(self, filename, **kwargs): import matplotlib.pylab as plt import plotter This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bpe...@us...> - 2015-04-03 21:21:08
|
Revision: 4961 http://sourceforge.net/p/simupop/code/4961 Author: bpeng2000 Date: 2015-04-03 21:21:00 +0000 (Fri, 03 Apr 2015) Log Message: ----------- Add rpy2 as one of the supported plotters Modified Paths: -------------- trunk/simuOpt.py trunk/src/plotter.py Modified: trunk/simuOpt.py =================================================================== --- trunk/simuOpt.py 2015-02-03 01:56:06 UTC (rev 4960) +++ trunk/simuOpt.py 2015-04-03 21:21:00 UTC (rev 4961) @@ -196,8 +196,8 @@ commandline option ``--gui``. plotter - The plotting library used to draw figures. simuPOP will by default - use ``rpy`` (NOT ``rpy2``), and ``matplotlib`` if ``rpy`` is not + The plotting library used to draw figures. simuPOP will try to + use ``rpy``, ``rpy2``, and ``matplotlib`` if the module is available. You can set this parameter to one of them to force the use a particular library. @@ -273,7 +273,7 @@ simuOptions['NumThreads'] = numThreads elif numThreads is not None: raise TypeError('An integer number is expected for parameter numThreads.') - if plotter in ['rpy', 'matplotlib']: + if plotter in ['rpy', 'rpy2', 'matplotlib']: simuOptions['Plotter'] = plotter elif plotter is not None: raise TypeError('Only rpy and matplotlib are allowed for parameter plotter.') Modified: trunk/src/plotter.py =================================================================== --- trunk/src/plotter.py 2015-02-03 01:56:06 UTC (rev 4960) +++ trunk/src/plotter.py 2015-04-03 21:21:00 UTC (rev 4961) @@ -27,11 +27,10 @@ ''' This module defines several utility functions and Python operators that make -use of the Python rpy module (http://rpy.sourceforge.net) to plot expressions -and information fields of evolving populations using a popular statistical -analysis language R (http://www.r-project.org) or matplotlib. Note that -rpy2, the successor of rpy, is currently not supported. You can specify the -plotting library using simuOpt.setOptions(plotter='matplotlib'). +use of the Python rpy or rpy2 module (http://rpy.sourceforge.net) to plot +expressions and information fields of evolving populations using a popular +statistical analysis language R (http://www.r-project.org) or matplotlib. +You can specify the plotting library using simuOpt.setOptions(plotter='matplotlib'). Each operator calls a sequence of R or matplotlib functions to draw and save figures. A special parameter passing mechanism is used so that you can specify @@ -40,7 +39,8 @@ ``lty_rep=[1,2]`` to pass ``lty=1`` and ``lty=2`` to specify different line types for different replicates if ``rpy`` is used. The help message of each class will describe which and in what sequence these R or matplotlib functions are -called to help you figure out which parameters are allowed. +called to help you figure out which parameters are allowed. The syntax +will be different if different plotting module is used. ''' __all__ = [ @@ -53,8 +53,9 @@ 'HistPlotter', 'QQPlotter', 'BoxPlotter', - # export essential piece of rpy so that other modules could use them + # export essential piece of rpy or rpy2 so that other modules could use them 'r', + # exported with rpy. Not meaningful if rpy2 is used. 'with_mode', 'NO_CONVERSION' ] @@ -64,41 +65,81 @@ from simuOpt import simuOptions +# used by rpy2 +def my_py2ri(o): + '''Convert tuple to R object. This is defined because the default + py2ri function does not support tuple type.''' + if isinstance(o, (tuple, list)): + if all([isinstance(x, int) for x in o]): + res = ro.vectors.IntVector(o) + elif all([isinstance(x, (int, float)) for x in o]): + res = ro.vectors.FloatVector(o) + elif all([isinstance(x, str) for x in o]): + res = ro.vectors.StrVector(o) + elif isinstance(o, tuple): + raise RuntimeError('Failed to convert tuple {}'.format(o)) + else: + # default_py2ri can handle list of + # heterogeneous type, I guess + res = ro.default_py2ri(o) + else: + res = ro.default_py2ri(o) + return res + if simuOptions['Plotter'] is None: try: import rpy_options rpy_options.set_options(VERBOSE = False) from rpy import r, with_mode, NO_CONVERSION - use_rpy = True + plotter = 'rpy' except ImportError, e: try: - import matplotlib.pylab as plt - use_rpy = False - except: - print('Neither rpy nor matplotlib is available.') - raise e + import rpy2.robjects as ro + from rpy2.robjects import r, conversion + conversion.py2ri = my_py2ri + with_mode = None + NO_CONVERSION = None + plotter = 'rpy2' + except ImportError, e: + try: + import matplotlib.pylab as plt + plotter = 'matplotlib' + except: + print('Neither rpy nor matplotlib is available.') + raise e elif simuOptions['Plotter'] == 'rpy': import rpy_options rpy_options.set_options(VERBOSE = False) from rpy import r, with_mode, NO_CONVERSION - use_rpy = True + plotter = 'rpy' +elif simuOptions['Plotter'] == 'rpy2': + import rpy2.robjects as ro + from rpy2.robjects import r, conversion + conversion.py2ri = my_py2ri + with_mode = None + NO_CONVERSION = None + plotter = 'rpy2' else: import matplotlib.pylab as plt - use_rpy = False + plotter = 'matplotlib' # if under windows, fix a bug with rpy which uses blocking i/o so # R figure will not be refreshed in time. See # https://Stat.ethz.ch/pipermail/r-devel/2006-January/036049.html # for details. -if use_rpy and os.name == 'nt': - r.options(windowsBuffered=False) +if os.name == 'nt': # In addition to options(windowsBuffered=False), I find that I also need to # call windows.options(buffered=False) to make functions such as hist work. # # This function is only available for R 2.9.0 (rev 48333) - if int(r.R_Version()['svn rev']) >= 48333: - r.windows_options(buffered=False) + if plotter == 'rpy': + r.options(windowsBuffered=False) + if int(r.R_Version()['svn rev']) >= 48333: + r.windows_options(buffered=False) + elif plotter == 'rpy2': + if int(r('R.Version()$"svn rev"')[0]) >= 48333: + r("windows.options(buffered=False)") from simuPOP import PyOperator, ALL_AVAIL @@ -107,7 +148,7 @@ '''Create a new graphics window and return its device number in R. This function essentially calls ``getOption('device')()`` in R. ''' - if use_rpy: + if plotter == 'rpy': # open a new window try: # 46754 is the revision number for R 2.8.0 @@ -124,6 +165,24 @@ if device == 0: raise RuntimeError('Can not open new device') return device + elif plotter == 'rpy2': + # open a new window + try: + # 46754 is the revision number for R 2.8.0 + if int(r('R.Version()$"svn rev"')[0]) < 46754: + # For R < 2.8.0, getOption('device') returns a string (such as 'X11') + r(r.getOption('device') + '()') + else: + # For R >= 2.8.0, getOption('device') returns a function + r('getOption("device")()') + except: + raise RuntimeError("Failed to get R version to start a graphical device"); + # get device number + device = r['dev.cur']() + if device == 0: + raise RuntimeError('Can not open new device') + return device + else: return plt.figure() @@ -139,7 +198,7 @@ ''' if file is None: return - if use_rpy: + if 'rpy' in plotter: filename, ext = os.path.splitext(file) dirname = os.path.dirname(file) if dirname != '' and not os.path.isdir(dirname): @@ -175,7 +234,10 @@ print 'Can not determine which device to use to save file %s. A postscript driver is used.' % name device = r.postscript params.update(kwargs) - r.dev_print(file=file, device=device, **params) + if plotter == 'rpy': + r.dev_print(file=file, device=device, **params) + else: + r['dev.print'](file=file, device=device, **params) else: plt.savefig(file) @@ -482,7 +544,7 @@ self.preHook = preHook self.postHook = postHook self.plotHook = plotHook - if use_rpy: + if 'rpy' in plotter: self.args = DerivedArgs( defaultFuncs = ['plot', 'lines'], allFuncs = ['par', 'plot', 'lines', 'dev_print', 'legend'], @@ -522,15 +584,18 @@ self.gen = [] self.data = [] # when apply is called, self._rpy_plot is called. - PyOperator.__init__(self, func=self._rpy_plot if use_rpy else self._mat_plot, + PyOperator.__init__(self, func=self._rpy_plot if plotter == 'rpy' else self._mat_plot, begin=begin, end=end, step=step, at=at, reps=reps, subPops=ALL_AVAIL, infoFields=[]) def __del__(self): # Close the device if needed. - if use_rpy: + if plotter == 'rpy': if not self.leaveOpen and hasattr(self, 'device'): r.dev_off(self.device) + elif plotter == 'rpy2': + if not self.leaveOpen and hasattr(self, 'device'): + r['dev.off'](self.device) else: if not self.leaveOpen: plt.close() @@ -604,7 +669,10 @@ if not hasattr(self, 'device'): self.device = newDevice() else: # if there are multiple devices, set it back - r.dev_set(self.device) + if plotter == 'rpy': + r.dev_set(self.device) + else: + r['dev.set'](self.device) # call the preHook function if given if self.preHook is not None: self.preHook(r) @@ -949,7 +1017,7 @@ self.preHook = preHook self.postHook = postHook self.subPops = subPops - if use_rpy: + if plotter == 'rpy': self.args = DerivedArgs( defaultFuncs = ['plot', 'points'], allFuncs = ['par', 'plot', 'points', 'dev_print', 'legend'], @@ -977,7 +1045,7 @@ ) if len(self.subPops) > 1: - if use_rpy: + if plotter == 'rpy': self.args.addDefault( pch_sp = range(1, len(self.subPops) + 1), col_sp = r.rainbow(len(self.subPops))) @@ -986,14 +1054,17 @@ self.args.addDefault( scatter_c_sp=[cm(i*1.0/len(self.subPops)) for i in range(len(self.subPops))]) # when apply is called, self._rpy_plot is called. - PyOperator.__init__(self, func=self._rpy_plot if use_rpy else self._mat_plot, + PyOperator.__init__(self, func=self._rpy_plot if 'rpy' in plotter else self._mat_plot, begin=begin, end=end, step=step, at=at, reps=reps) def __del__(self): # Close the device if needed. if not self.leaveOpen and hasattr(self, 'device'): - r.dev_off(self.device) + if plotter == 'rpy': + r.dev_off(self.device) + elif plotter == 'rpy2': + r['dev.off'](self.device) def _rpy_plot(self, pop): "Evaluate expression in pop and save result. Plot all data if needed" @@ -1003,7 +1074,10 @@ if not hasattr(self, 'device'): self.device = newDevice() else: # if there are multiple devices, set it back - r.dev_set(self.device) + if plotter == 'rpy': + r.dev_set(self.device) + else: + r['dev.set'](self.device) # call the preHook function if given if self.preHook is not None: self.preHook(r) @@ -1195,9 +1269,12 @@ self.postHook = postHook self.plotHook = plotHook self.subPops = subPops - if use_rpy: + if 'rpy' in plotter: if self.func is not None: - self.rfunc = r(self.func) + if plotter == 'rpy': + self.rfunc = r(self.func) + else: + self.rfunc = r[self.func] self.args = DerivedArgs( defaultFuncs = [] if self.func is None else [self.func], allFuncs = ['par', 'dev_print', 'legend'] + ([] if self.func is None else [self.func]), @@ -1213,9 +1290,12 @@ def __del__(self): # Close the device if needed. - if use_rpy: + if plotter == 'rpy': if not self.leaveOpen and hasattr(self, 'device'): r.dev_off(self.device) + elif plotter == 'rpy2': + if not self.leaveOpen and hasattr(self, 'device'): + r['dev.off'](self.device) else: if not self.leaveOpen: plt.close() @@ -1228,7 +1308,10 @@ if not hasattr(self, 'device'): self.device = newDevice() else: # if there are multiple devices, set it back - r.dev_set(self.device) + if plotter == 'rpy': + r.dev_set(self.device) + else: + r['dev.set'](self.device) # call the preHook function if given if self.preHook is not None: self.preHook(r) @@ -1388,7 +1471,7 @@ self.infoFields = infoFields if len(self.infoFields) == 0: raise RuntimeError('At least one information field should be given') - if not use_rpy: + if 'rpy' not in plotter: raise RuntimeError('BoxPlotter function does not support matplotlib plotter') self.saveAs = saveAs self.leaveOpen = leaveOpen @@ -1416,8 +1499,12 @@ def __del__(self): # Close the device if needed. if not self.leaveOpen and hasattr(self, 'device'): - r.dev_off(self.device) + if plotter == 'rpy': + r.dev_off(self.device) + elif plotter == 'rpy2': + r['dev.off'](self.device) + def _rpy_plot(self, pop): "Evaluate expression in pop and save result. Plot all data if needed" gen = pop.dvars().gen @@ -1426,7 +1513,10 @@ if not hasattr(self, 'device'): self.device = newDevice() else: # if there are multiple devices, set it back - r.dev_set(self.device) + if plotter == 'rpy': + r.dev_set(self.device) + else: + r['dev.set'](self.device) # call the preHook function if given if self.preHook is not None: self.preHook(r) @@ -1470,9 +1560,15 @@ data.extend(spData) owner.extend([pop.subPopName(sp)]*len(spData)) # - r.boxplot(r('data ~ owner'), data=r.data_frame(data=data, owner=owner), - **self.args.getArgs('boxplot', pop, fld=fldIdx, - main='Field %s at gen %d' % (fld, gen))) + if plotter == 'rpy': + r.boxplot(r('data ~ owner'), data=r.data_frame(data=data, owner=owner), + **self.args.getArgs('boxplot', pop, fld=fldIdx, + main='Field %s at gen %d' % (fld, gen))) + else: + r.boxplot(ro.Formula('data ~ owner'), + data=ro.vectors.DataFrame({'data': data, 'owner': owner}), + **self.args.getArgs('boxplot', pop, fld=fldIdx, + main='Field %s at gen %d' % (fld, gen))) if self.plotHook is not None: self.plotHook(r=r, field=fld) elif not self.byField and self.bySubPop: @@ -1485,9 +1581,15 @@ data.extend(fldData) owner.extend([fld]*len(fldData)) # - r.boxplot(r('data ~ owner'), data=r.data_frame(data=data, owner=owner), - **self.args.getArgs('boxplot', pop, sp=spIdx, - main='Subpop %s at gen %d' % (pop.subPopName(sp), gen))) + if plotter == 'rpy': + r.boxplot(r('data ~ owner'), data=r.data_frame(data=data, owner=owner), + **self.args.getArgs('boxplot', pop, sp=spIdx, + main='Subpop %s at gen %d' % (pop.subPopName(sp), gen))) + else: + r.boxplot(ro.Formula('data ~ owner'), + data=ro.vectors.DataFrame({'data': data, 'owner': owner}), + **self.args.getArgs('boxplot', pop, sp=spIdx, + main='Subpop %s at gen %d' % (pop.subPopName(sp), gen))) if self.plotHook is not None: self.plotHook(r=r, subPop=sp) else: @@ -1508,8 +1610,13 @@ else: owner.extend(['%s, %s' % (fld, pop.subPopName(sp))]*len(spData)) # - r.boxplot(r("data ~ owner"), data=r.data_frame(data=data, owner=owner), - **self.args.getArgs('boxplot', pop)) + if plotter == 'rpy': + r.boxplot(r("data ~ owner"), data=r.data_frame(data=data, owner=owner), + **self.args.getArgs('boxplot', pop)) + else: + r.boxplot(ro.Formula('data ~ owner'), + data=ro.vectors.DataFrame({'data': data, 'owner': owner}), + **self.args.getArgs('boxplot', pop)) if self.plotHook is not None: self.plotHook(r=r) # call the postHook function if given This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: Bo P. <be...@gm...> - 2015-02-03 15:45:33
|
Hi, Roman, Implementation is easy but I need to consider side effects such as compatibility (some users process dump output) and duplicated information. I have implemented a feature that makes use of existing but ignored parameter infoFields that allows you to specify the information fields to display. This will hopefully make the output easier to read, especially when your application has many info fields. The change is in the subversion version of simuPOP so you will have to compile from source to use it. Cheers, Bo On Tue, Feb 3, 2015 at 4:50 AM, Roman Lustrik <rom...@bi...> wrote: > Hello Bo, > > if you can't implement the solution "on-the-fly", don't worry about it. If > it's more work than worth, then I'll be able to live with setting > `structure = True`. > > Cheers, > Roman > > ---- > In god we trust, all others bring data. > > ------------------------------ > *From: *"Bo Peng" <be...@gm...> > *To: *"Roman Lustrik" <rom...@bi...> > *Cc: *sim...@li... > *Sent: *Monday, February 2, 2015 10:40:24 PM > *Subject: *Re: [simuPOP-devel] feature request for displaying column names > > > Hi, Roman, > > The info field information is already available if you do not use option > 'structure=False'. I do not know if it is a good idea to repeat that > information in the genotype sections. However, the Dumper() operator (dump > function) currently ignore parameter infoFields. Will it help if the > function lists only the information fields you specify, and in the order > they are specified? > > Cheers, > Bo > > > On Fri, Jan 30, 2015 at 5:19 AM, Roman Lustrik <rom...@bi... > > wrote: > >> Dear developer(s), >> >> I have been using simuPOP for the last month or so and have grown not >> only to like it, but love it. While developing increasingly complex >> simulation, I have learned that having a bit more informed dump() result >> would be helpful. Would you consider implementing an option to turn on/off >> printing of infoField names? For example, I have a number of columns and it >> would be a bit easier if infoField names were displayed on the spot. >> >> SubPopulation 0 (), 23 Individuals: >> 0: FU 699 | 911 | 0 0 0 0 1 3 5 >> 1: FU 175 | 238 | 0 0 0 0 1 3 5 >> 2: MU 699 | 238 | 0 0 0 0 1 3 5 >> 3: MU 364 | 912 | 0 0 0 0 1 1 0 >> 4: FU 323 | 836 | 0 0 0 0 1 1 0 >> >> If infoFields were displayed, I imagine the output would look something >> like this: >> >> SubPopulation 0 (), 23 Individuals: >> >> age ywc fitness hc iscub >> father_idx mother_idx >> 0: FU 699 | 911 | 0 0 0 0 1 3 5 >> 1: FU 175 | 238 | 0 0 0 0 1 3 5 >> 2: MU 699 | 238 | 0 0 0 0 1 3 5 >> 3: MU 364 | 912 | 0 0 0 0 1 1 0 >> 4: FU 323 | 836 | 0 0 0 0 1 1 0 >> >> >> Cheers, >> >> Roman >> >> >> >> ---- >> In god we trust, all others bring data. >> >> >> >> ------------------------------------------------------------------------------ >> Dive into the World of Parallel Programming. The Go Parallel Website, >> sponsored by Intel and developed in partnership with Slashdot Media, is >> your >> hub for all things parallel software development, from weekly thought >> leadership blogs to news, videos, case studies, tutorials and more. Take a >> look and join the conversation now. http://goparallel.sourceforge.net/ >> _______________________________________________ >> simuPOP-devel mailing list >> sim...@li... >> https://lists.sourceforge.net/lists/listinfo/simupop-devel >> >> > > |
From: Roman L. <rom...@bi...> - 2015-02-03 10:47:29
|
Hello Bo, if you can't implement the solution "on-the-fly", don't worry about it. If it's more work than worth, then I'll be able to live with setting `structure = True`. Cheers, Roman ---- In god we trust, all others bring data. ----- Original Message ----- From: "Bo Peng" <be...@gm...> To: "Roman Lustrik" <rom...@bi...> Cc: sim...@li... Sent: Monday, February 2, 2015 10:40:24 PM Subject: Re: [simuPOP-devel] feature request for displaying column names Hi, Roman, The info field information is already available if you do not use option 'structure=False'. I do not know if it is a good idea to repeat that information in the genotype sections. However, the Dumper() operator (dump function) currently ignore parameter infoFields. Will it help if the function lists only the information fields you specify, and in the order they are specified? Cheers, Bo On Fri, Jan 30, 2015 at 5:19 AM, Roman Lustrik < rom...@bi... > wrote: Dear developer(s), I have been using simuPOP for the last month or so and have grown not only to like it, but love it. While developing increasingly complex simulation, I have learned that having a bit more informed dump() result would be helpful. Would you consider implementing an option to turn on/off printing of infoField names? For example, I have a number of columns and it would be a bit easier if infoField names were displayed on the spot. SubPopulation 0 (), 23 Individuals: 0: FU 699 | 911 | 0 0 0 0 1 3 5 1: FU 175 | 238 | 0 0 0 0 1 3 5 2: MU 699 | 238 | 0 0 0 0 1 3 5 3: MU 364 | 912 | 0 0 0 0 1 1 0 4: FU 323 | 836 | 0 0 0 0 1 1 0 If infoFields were displayed, I imagine the output would look something like this: SubPopulation 0 (), 23 Individuals: age ywc fitness hc iscub father_idx mother_idx 0: FU 699 | 911 | 0 0 0 0 1 3 5 1: FU 175 | 238 | 0 0 0 0 1 3 5 2: MU 699 | 238 | 0 0 0 0 1 3 5 3: MU 364 | 912 | 0 0 0 0 1 1 0 4: FU 323 | 836 | 0 0 0 0 1 1 0 Cheers, Roman ---- In god we trust, all others bring data. ------------------------------------------------------------------------------ Dive into the World of Parallel Programming. The Go Parallel Website, sponsored by Intel and developed in partnership with Slashdot Media, is your hub for all things parallel software development, from weekly thought leadership blogs to news, videos, case studies, tutorials and more. Take a look and join the conversation now. http://goparallel.sourceforge.net/ _______________________________________________ simuPOP-devel mailing list sim...@li... https://lists.sourceforge.net/lists/listinfo/simupop-devel |
From: <bpe...@us...> - 2015-02-03 01:56:09
|
Revision: 4960 http://sourceforge.net/p/simupop/code/4960 Author: bpeng2000 Date: 2015-02-03 01:56:06 +0000 (Tue, 03 Feb 2015) Log Message: ----------- 1. allow stringList.elems to take population as a parameter so that ALL_AVAIL is honored. 2. BaseOperator.infoSize() now takes a trait as a parameter so that it can get ALL_AVAIL fields. 3. Adjust Dumper(infoFields) so that it accepts ALL_AVAIL as its default parameter. 4. Make use of parameter infoFields of Dumper() to output just specified info fields. Modified Paths: -------------- trunk/src/individual.cpp trunk/src/individual.h trunk/src/operator.h trunk/src/outputer.cpp trunk/src/outputer.h trunk/src/utility.cpp trunk/src/utility.h trunk/test/test_03_operator.py Modified: trunk/src/individual.cpp =================================================================== --- trunk/src/individual.cpp 2015-02-03 01:31:54 UTC (rev 4959) +++ trunk/src/individual.cpp 2015-02-03 01:56:06 UTC (rev 4960) @@ -1006,7 +1006,7 @@ } -void Individual::display(ostream & out, int width, const vectoru & loci) +void Individual::display(ostream & out, int width, const vectoru & loci, const vectoru & infoIdx) { out << (sex() == MALE ? 'M' : 'F') << (affected() ? 'A' : 'U') << " "; int pEnd = ploidy(); @@ -1028,8 +1028,9 @@ } if (infoSize() != 0) { out << "| "; - for (vectorf::const_iterator info = infoBegin(); info != infoEnd(); ++info) - out << " " << *info; + for (vectoru::const_iterator info = infoIdx.begin(); + info != infoIdx.end(); ++info) + out << " " << *(infoBegin() + *info); } } Modified: trunk/src/individual.h =================================================================== --- trunk/src/individual.h 2015-02-03 01:31:54 UTC (rev 4959) +++ trunk/src/individual.h 2015-02-03 01:56:06 UTC (rev 4960) @@ -716,7 +716,8 @@ //@{ /// CPPONLY - void display(ostream & out, int width = 1, const vectoru & loci = vectoru()); + void display(ostream & out, int width = 1, const vectoru & loci = vectoru(), + const vectoru & infoIdx = vectoru()); //@} Modified: trunk/src/operator.h =================================================================== --- trunk/src/operator.h 2015-02-03 01:31:54 UTC (rev 4959) +++ trunk/src/operator.h 2015-02-03 01:56:06 UTC (rev 4960) @@ -200,17 +200,16 @@ /// get the length of information fields for this operator /// CPPONLY - size_t infoSize() const + size_t infoSize(const GenoStruTrait * trait = NULL) const { - return m_infoFields.elems().size(); + return m_infoFields.elems(trait).size(); } - - + /// get the information field specified by user (or by default) /// CPPONLY - string infoField(size_t idx) const + string infoField(size_t idx, const GenoStruTrait * trait = NULL) const { - DBG_ASSERT(idx < m_infoFields.elems().size(), IndexError, + DBG_ASSERT(idx < m_infoFields.elems(trait).size(), IndexError, (boost::format("Given info index %1% is out of range of 0 ~ %2%") % idx % m_infoFields.elems().size()).str()); return m_infoFields.elems()[idx]; } Modified: trunk/src/outputer.cpp =================================================================== --- trunk/src/outputer.cpp 2015-02-03 01:31:54 UTC (rev 4959) +++ trunk/src/outputer.cpp 2015-02-03 01:56:06 UTC (rev 4960) @@ -105,6 +105,11 @@ subPopList::const_iterator sp = subPops.begin(); subPopList::const_iterator spEnd = subPops.end(); + vectoru infoIdx; + if (infoSize(&pop) > 0) + for (size_t i = 0; i < infoSize(); ++i) + infoIdx.push_back(pop.infoIdx(infoField(i))); + for ( ; sp != spEnd; ++sp) { size_t spSize = pop.subPopSize(*sp); out << "SubPopulation " << *sp << " (" << pop.subPopName(*sp) << "), " @@ -114,7 +119,7 @@ IndIterator ind = const_cast<Population &>(pop).indIterator(sp->subPop()); for ( ; ind.valid(); ++ind, ++count) { out << setw(4) << (&*ind - &*pop.rawIndBegin()) << ": "; - ind->display(out, m_width, m_loci); + ind->display(out, m_width, m_loci, infoIdx); out << endl; if (m_max > 0 && count + 1 >= m_max && count < pop.popSize()) break; Modified: trunk/src/outputer.h =================================================================== --- trunk/src/outputer.h 2015-02-03 01:31:54 UTC (rev 4959) +++ trunk/src/outputer.h 2015-02-03 01:56:06 UTC (rev 4960) @@ -100,7 +100,9 @@ * this operator will only output the first 100 (parameter \e max) * individuals of the present generation (parameter \e ancGens). All loci * will be outputed unless parameter \e loci are used to specify a subset - * of loci. If a list of (virtual) subpopulations are specified, this + * of loci. This operator by default output values of all information fields + * unless parameter \e infoFields is used to specify a subset of info fields + * to display. If a list of (virtual) subpopulations are specified, this * operator will only output individuals in these outputs. Please refer to * class \c BaseOperator for a detailed explanation for common parameters * such as \e output and \e stage. @@ -109,7 +111,7 @@ int width = 1, UINT max = 100, const uintList & loci = vectoru(), const stringFunc & output = ">", int begin = 0, int end = -1, int step = 1, const intList & at = vectori(), const intList & reps = intList(), const subPopList & subPops = subPopList(), - const stringList & infoFields = vectorstr()) : + const stringList & infoFields = stringList()) : BaseOperator(output, begin, end, step, at, reps, subPops, infoFields), m_showGenotype(genotype), m_showStructure(structure), m_ancGens(ancGens), m_width(width), m_loci(loci.elems()), m_max(max) Modified: trunk/src/utility.cpp =================================================================== --- trunk/src/utility.cpp 2015-02-03 01:31:54 UTC (rev 4959) +++ trunk/src/utility.cpp 2015-02-03 01:56:06 UTC (rev 4960) @@ -587,7 +587,7 @@ } -stringList::stringList(PyObject * obj) : m_elems(), m_allAvail(false) +stringList::stringList(PyObject * obj) : m_elems(), m_allAvail(false), m_trait(MaxTraitIndex) { if (obj == NULL || obj == Py_None) m_allAvail = true; @@ -662,6 +662,21 @@ m_elems.push_back(defaultItems[i]); } +const vectorstr & stringList::elems(const GenoStruTrait * trait) const +{ + if (trait == NULL || ! allAvail()) + return m_elems; + + if (trait) { + if (trait->genoStruIdx() == m_trait) + return m_elems; + m_elems.clear(); + const vectorstr fields = trait->infoFields(); + m_elems.insert(m_elems.end(), fields.begin(), fields.end()); + m_trait = trait->genoStruIdx(); + } + return m_elems; +} intMatrix::intMatrix(PyObject * obj) : m_elems() { Modified: trunk/src/utility.h =================================================================== --- trunk/src/utility.h 2015-02-03 01:31:54 UTC (rev 4959) +++ trunk/src/utility.h 2015-02-03 01:56:06 UTC (rev 4960) @@ -706,13 +706,13 @@ stringList(PyObject * str = NULL); /// CPPONLY - stringList(const string & str) : m_elems(1, str), m_allAvail(false) + stringList(const string & str) : m_elems(1, str), m_allAvail(false), m_trait(MaxTraitIndex) { } /// CPPONLY - stringList(const string & str1, const string & str2) : m_elems(), m_allAvail(false) + stringList(const string & str1, const string & str2) : m_elems(), m_allAvail(false), m_trait(MaxTraitIndex) { m_elems.push_back(str1); m_elems.push_back(str2); @@ -720,7 +720,7 @@ /// CPPONLY - stringList(const vectorstr & str) : m_elems(str), m_allAvail(false) + stringList(const vectorstr & str) : m_elems(str), m_allAvail(false), m_trait(MaxTraitIndex) { } @@ -758,18 +758,18 @@ /// CPPONLY - const vectorstr & elems() const - { - return m_elems; - } + /// Return lists. Population from information fields of a population + /// if allAvail is true. + const vectorstr & elems(const GenoStruTrait * trait = NULL) const; private: void addString(PyObject * str); protected: - vectorstr m_elems; + mutable vectorstr m_elems; bool m_allAvail; + mutable TraitIndexType m_trait; }; Modified: trunk/test/test_03_operator.py =================================================================== --- trunk/test/test_03_operator.py 2015-02-03 01:31:54 UTC (rev 4959) +++ trunk/test/test_03_operator.py 2015-02-03 01:56:06 UTC (rev 4960) @@ -361,6 +361,15 @@ ) + def testDumper(self): + '''Testing operator Dumper''' + pop = Population(100, loci=10, infoFields=('a', 'b')) + dump(pop, genotype=False, output='a.dump') + dump(pop, structure = False, output='a.dump') + dump(pop, width=2, output='a.dump') + dump(pop, infoFields='a', output='a.dump') + dump(pop, infoFields=('b', 'a', 'b'), output='a.dump') + def testcloseOutput(self): '''Testing global function closeOutput''' pop = Population(100, loci=[2]) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bpe...@us...> - 2015-02-03 01:31:57
|
Revision: 4959 http://sourceforge.net/p/simupop/code/4959 Author: bpeng2000 Date: 2015-02-03 01:31:54 +0000 (Tue, 03 Feb 2015) Log Message: ----------- Update comments about importing 0 alleles with negative adjustment values for function import Modified Paths: -------------- trunk/src/utils.py Modified: trunk/src/utils.py =================================================================== --- trunk/src/utils.py 2015-01-15 14:47:50 UTC (rev 4958) +++ trunk/src/utils.py 2015-02-03 01:31:54 UTC (rev 4959) @@ -3183,15 +3183,20 @@ adjust Adjust alleles by specified value (default to 0 for no adjustment). This parameter is mostly used to convert alleles 1 and 2 in a GenePop file to - alleles 0 and 1 (with adjust=-1) in simuPOP. + alleles 0 and 1 (with adjust=-1) in simuPOP. Negative allele (e.g. missing + value 0) will be imported as regular allele with module-dependent values + (e.g. -1 imported as 255 for standard module). + FSTAT (http://www2.unil.ch/popgen/softwares/fstat.htm). This format accepts the following parameters: adjust Adjust alleles by specified value (default to 0 for no adjustment). This parameter is mostly used to convert alleles 1 and 2 in a GenePop file to - alleles 0 and 1 (with adjust=-1) in simuPOP. + alleles 0 and 1 (with adjust=-1) in simuPOP. Negative allele (e.g. missing + value 0) will be imported as regular allele with module-dependent values + (e.g. -1 imported as 255 for standard module). Phylip (Joseph Felsenstein's Phylip format). This function ignores sequence names and import sequences in a haploid (default) or diploid population (if This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: Bo P. <be...@gm...> - 2015-02-02 21:40:52
|
Hi, Roman, The info field information is already available if you do not use option 'structure=False'. I do not know if it is a good idea to repeat that information in the genotype sections. However, the Dumper() operator (dump function) currently ignore parameter infoFields. Will it help if the function lists only the information fields you specify, and in the order they are specified? Cheers, Bo On Fri, Jan 30, 2015 at 5:19 AM, Roman Lustrik <rom...@bi...> wrote: > Dear developer(s), > > I have been using simuPOP for the last month or so and have grown not only > to like it, but love it. While developing increasingly complex simulation, > I have learned that having a bit more informed dump() result would be > helpful. Would you consider implementing an option to turn on/off printing > of infoField names? For example, I have a number of columns and it would be > a bit easier if infoField names were displayed on the spot. > > SubPopulation 0 (), 23 Individuals: > 0: FU 699 | 911 | 0 0 0 0 1 3 5 > 1: FU 175 | 238 | 0 0 0 0 1 3 5 > 2: MU 699 | 238 | 0 0 0 0 1 3 5 > 3: MU 364 | 912 | 0 0 0 0 1 1 0 > 4: FU 323 | 836 | 0 0 0 0 1 1 0 > > If infoFields were displayed, I imagine the output would look something > like this: > > SubPopulation 0 (), 23 Individuals: > > age ywc fitness hc iscub > father_idx mother_idx > 0: FU 699 | 911 | 0 0 0 0 1 3 5 > 1: FU 175 | 238 | 0 0 0 0 1 3 5 > 2: MU 699 | 238 | 0 0 0 0 1 3 5 > 3: MU 364 | 912 | 0 0 0 0 1 1 0 > 4: FU 323 | 836 | 0 0 0 0 1 1 0 > > > Cheers, > > Roman > > > > ---- > In god we trust, all others bring data. > > > > ------------------------------------------------------------------------------ > Dive into the World of Parallel Programming. The Go Parallel Website, > sponsored by Intel and developed in partnership with Slashdot Media, is > your > hub for all things parallel software development, from weekly thought > leadership blogs to news, videos, case studies, tutorials and more. Take a > look and join the conversation now. http://goparallel.sourceforge.net/ > _______________________________________________ > simuPOP-devel mailing list > sim...@li... > https://lists.sourceforge.net/lists/listinfo/simupop-devel > > |
From: Roman L. <rom...@bi...> - 2015-01-30 11:16:05
|
Dear developer(s), I have been using simuPOP for the last month or so and have grown not only to like it, but love it. While developing increasingly complex simulation, I have learned that having a bit more informed dump() result would be helpful. Would you consider implementing an option to turn on/off printing of infoField names? For example, I have a number of columns and it would be a bit easier if infoField names were displayed on the spot. SubPopulation 0 (), 23 Individuals: 0: FU 699 | 911 | 0 0 0 0 1 3 5 1: FU 175 | 238 | 0 0 0 0 1 3 5 2: MU 699 | 238 | 0 0 0 0 1 3 5 3: MU 364 | 912 | 0 0 0 0 1 1 0 4: FU 323 | 836 | 0 0 0 0 1 1 0 If infoFields were displayed, I imagine the output would look something like this: SubPopulation 0 (), 23 Individuals: age ywc fitness hc iscub father_idx mother_idx 0: FU 699 | 911 | 0 0 0 0 1 3 5 1: FU 175 | 238 | 0 0 0 0 1 3 5 2: MU 699 | 238 | 0 0 0 0 1 3 5 3: MU 364 | 912 | 0 0 0 0 1 1 0 4: FU 323 | 836 | 0 0 0 0 1 1 0 Cheers, Roman ---- In god we trust, all others bring data. |
From: <bpe...@us...> - 2015-01-15 14:47:58
|
Revision: 4958 http://sourceforge.net/p/simupop/code/4958 Author: bpeng2000 Date: 2015-01-15 14:47:50 +0000 (Thu, 15 Jan 2015) Log Message: ----------- Update two examples on users guide Modified Paths: -------------- trunk/doc/userGuide.py Modified: trunk/doc/userGuide.py =================================================================== --- trunk/doc/userGuide.py 2014-12-08 20:21:32 UTC (rev 4957) +++ trunk/doc/userGuide.py 2015-01-15 14:47:50 UTC (rev 4958) @@ -6485,9 +6485,10 @@ yield males[rank][randint(0, len(males[rank]) - 1)], \ females[rank][randint(0, len(females[rank]) - 1)] -def setRank(off, dad): +def setRank(rank): 'The rank of offspring can increase or drop to zero randomly' - off.rank = (dad.rank + randint(-1, 1)) % 3 + # only use rank of the father + return (rank[0] + randint(-1, 1)) % 3 pop = sim.Population(size=[1000, 2000], loci=1, infoFields='rank') pop.evolve( @@ -6497,7 +6498,10 @@ ], matingScheme=sim.HomoMating( sim.PyParentsChooser(randomChooser), - sim.OffspringGenerator(ops=sim.MendelianGenoTransmitter()) + sim.OffspringGenerator(ops=[ + sim.MendelianGenoTransmitter(), + sim.PyTagger(setRank), + ]) ), gen = 5 ) @@ -6531,9 +6535,10 @@ # find a female in the same rank yield m, females[rank][randint(0, len(females[rank]) - 1)] -def setRank(off, dad): +def setRank(rank): 'The rank of offspring can increase or drop to zero randomly' - off.rank = (dad.rank + randint(-1, 1)) % 3 + # only use rank of the father + return (rank[0] + randint(-1, 1)) % 3 pop = sim.Population(size=[1000, 2000], loci=1, infoFields='rank') pop.evolve( @@ -6543,7 +6548,10 @@ ], matingScheme=sim.HomoMating( sim.PyParentsChooser(randomChooser), - sim.OffspringGenerator(ops=sim.MendelianGenoTransmitter()) + sim.OffspringGenerator(ops=[ + sim.MendelianGenoTransmitter(), + sim.PyTagger(setRank), + ]) ), gen = 5 ) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bpe...@us...> - 2014-12-08 20:21:40
|
Revision: 4957 http://sourceforge.net/p/simupop/code/4957 Author: bpeng2000 Date: 2014-12-08 20:21:32 +0000 (Mon, 08 Dec 2014) Log Message: ----------- Add parameter allowSelfing to CombinedParentsChooser. Add HermaphroditicMating. Modified Paths: -------------- trunk/doc/refManual.lyx trunk/src/__init__.py trunk/src/mating.cpp trunk/src/mating.h trunk/test/test_05_matings.py Modified: trunk/doc/refManual.lyx =================================================================== --- trunk/doc/refManual.lyx 2014-11-12 02:06:43 UTC (rev 4956) +++ trunk/doc/refManual.lyx 2014-12-08 20:21:32 UTC (rev 4957) @@ -679,6 +679,13 @@ \backslash +simuPOPHermaphroditicMatingRef +\end_layout + +\begin_layout Plain Layout + + +\backslash simuPOPControlledRandomMatingRef \end_layout Modified: trunk/src/__init__.py =================================================================== --- trunk/src/__init__.py 2014-11-12 02:06:43 UTC (rev 4956) +++ trunk/src/__init__.py 2014-12-08 20:21:32 UTC (rev 4957) @@ -143,6 +143,7 @@ 'RandomSelection', 'MonogamousMating', 'SelfMating', + 'HermaphroditicMating', 'CloneMating', 'HaplodiploidMating', #'consanguineousMating', @@ -852,6 +853,33 @@ subPops = subPops, weight = weight) +class HermaphroditicMating(HomoMating): + '''A hermaphroditic mating scheme that chooses two parents randomly + from the population regardless of sex. The parents could be chosen + with or without replacement (parameter *replacement*). Selfing (if + the same parents are chosen) is allowed unless *allowSelfing* is + set to *False* ''' + def __init__(self, replacement=True, allowSelfing=True, numOffspring = 1, + sexMode = RANDOM_SEX, ops = MendelianGenoTransmitter(), + subPopSize = [], subPops = ALL_AVAIL, weight = 0, + selectionField = 'fitness'): + '''Creates a hermaphroditic mating scheme where individuals can + serve as father or mother, or both (self-fertilization). Please + refer to class ``CombinedParentsChooser`` for parameter *allowSelfing``, + to ``RandomParentChooser`` for parameter *replacement* and + *selectionField*, to class ``OffspringGenerator`` for parameters *ops*, + *sexMode* and *numOffspring*, and to class ``HomoMating`` for parameters + *subPopSize*, *subPops* and *weight*. ''' + HomoMating.__init__(self, + chooser = CombinedParentsChooser( + RandomParentChooser(replacement, selectionField), + RandomParentChooser(replacement, selectionField), + allowSelfing=allowSelfing), + generator = OffspringGenerator(ops, numOffspring, sexMode), + subPopSize = subPopSize, + subPops = subPops, + weight = weight) + ## ## def consanguineousMating(infoFields = [], func = None, param = None, ## replacement = False, numOffspring = 1., sexMode = RANDOM_SEX, Modified: trunk/src/mating.cpp =================================================================== --- trunk/src/mating.cpp 2014-11-12 02:06:43 UTC (rev 4956) +++ trunk/src/mating.cpp 2014-12-08 20:21:32 UTC (rev 4957) @@ -1079,9 +1079,10 @@ */ CombinedParentsChooser::CombinedParentsChooser(const ParentChooser & fatherChooser, - const ParentChooser & motherChooser) + const ParentChooser & motherChooser, bool allowSelfing) : m_fatherChooser(fatherChooser.clone()), - m_motherChooser(motherChooser.clone()) + m_motherChooser(motherChooser.clone()), + m_allowSelfing(allowSelfing) { } @@ -1095,12 +1096,20 @@ ParentChooser::IndividualPair CombinedParentsChooser::chooseParents() { - ParentChooser::IndividualPair p1 = m_fatherChooser->chooseParents(); - ParentChooser::IndividualPair p2 = m_motherChooser->chooseParents(); - Individual * dad = p1.first != NULL ? p1.first : p1.second; - Individual * mom = p2.second != NULL ? p2.second : p2.first; + size_t attempts = 0; + while (attempts < 100) { + ParentChooser::IndividualPair p1 = m_fatherChooser->chooseParents(); + ParentChooser::IndividualPair p2 = m_motherChooser->chooseParents(); + Individual * dad = p1.first != NULL ? p1.first : p1.second; + Individual * mom = p2.second != NULL ? p2.second : p2.first; - return ParentChooser::IndividualPair(dad, mom); + if (dad != mom || m_allowSelfing) + return ParentChooser::IndividualPair(dad, mom); + if (attempts ++ == 100) + throw RuntimeError("Failed to select distinct parents using CombinedParentsChooser."); + } + // just make the compiler happy + return ParentChooser::IndividualPair((Individual *)(0), (Individual *)(0)); } Modified: trunk/src/mating.h =================================================================== --- trunk/src/mating.h 2014-11-12 02:06:43 UTC (rev 4956) +++ trunk/src/mating.h 2014-12-08 20:21:32 UTC (rev 4957) @@ -1267,15 +1267,19 @@ * for \e motherChooser. Although these two parent choosers are supposed * to return a father and a mother respectively, the sex of returned * parents are not checked so it is possible to return parents with the - * same sex using this parents chooser. + * same sex using this parents chooser. This choose by default allows + * the selection of the same parents as father and mother + * (self-fertilization), unless a parameter \e allowSelfing is used to + * disable it. */ CombinedParentsChooser(const ParentChooser & fatherChooser, - const ParentChooser & motherChooser); + const ParentChooser & motherChooser, bool allowSelfing=true); /// CPPONLY CombinedParentsChooser(const CombinedParentsChooser & rhs) : ParentChooser("fitness"), m_fatherChooser(rhs.m_fatherChooser->clone()), - m_motherChooser(rhs.m_motherChooser->clone()) + m_motherChooser(rhs.m_motherChooser->clone()), + m_allowSelfing(rhs.m_allowSelfing) { m_initialized = false; } @@ -1322,6 +1326,7 @@ private: ParentChooser * m_fatherChooser; ParentChooser * m_motherChooser; + bool m_allowSelfing; }; /** This parent chooser accepts a Python generator function that repeatedly Modified: trunk/test/test_05_matings.py =================================================================== --- trunk/test/test_05_matings.py 2014-11-12 02:06:43 UTC (rev 4956) +++ trunk/test/test_05_matings.py 2014-12-08 20:21:32 UTC (rev 4957) @@ -734,6 +734,28 @@ self.assertLess(m.a, 10) self.assertGreaterEqual(m.a, 0) + def testHermaphroditicMating(self): + 'Test HermaphroditicMating' + pop = Population(100, loci=10) + pop.evolve( + matingScheme=HermaphroditicMating(), + gen=10) + # + def theSame(dad, mom): + if dad.ind_id == mom.ind_id: + raise RuntimeError() + return True + # + pop = Population(100, loci=10, infoFields='ind_id') + pop.evolve( + preOps=IdTagger(), + matingScheme=HermaphroditicMating(allowSelfing=False, + ops=[MendelianGenoTransmitter(), + IdTagger(), + PyOperator(func=theSame), + ]), + gen=100) + if __name__ == '__main__': unittest.main() sys.exit(0) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bpe...@us...> - 2014-11-12 02:06:46
|
Revision: 4956 http://sourceforge.net/p/simupop/code/4956 Author: bpeng2000 Date: 2014-11-12 02:06:43 +0000 (Wed, 12 Nov 2014) Log Message: ----------- Add conda release command Modified Paths: -------------- trunk/development/conda/meta.yaml trunk/development/release.py trunk/setup.py Modified: trunk/development/conda/meta.yaml =================================================================== --- trunk/development/conda/meta.yaml 2014-11-11 03:05:43 UTC (rev 4955) +++ trunk/development/conda/meta.yaml 2014-11-12 02:06:43 UTC (rev 4956) @@ -1,6 +1,6 @@ package: name: simupop - version: !!str 1.0.0 + version: !!str 1.1.4 source: fn: simuPOP-1.1.4.tar.gz @@ -38,7 +38,7 @@ test: # Python imports imports: - - simupop + - simuPOP - src # commands: Modified: trunk/development/release.py =================================================================== --- trunk/development/release.py 2014-11-11 03:05:43 UTC (rev 4955) +++ trunk/development/release.py 2014-11-12 02:06:43 UTC (rev 4956) @@ -332,6 +332,10 @@ .format(dmg, ver, dest, pyver)) +def condaRelease(release): + '''Build conda binary''' + run_command('conda build conda') + def tagRelease(release): ''' if there are changes, commit it ''' if subprocess.check_output('svn diff', shell=True) != '': @@ -386,3 +390,5 @@ createMacImage(ver, pyver='2.7') if 'tag' in args.actions: tagRelease(ver) + if 'conda' in args.actions: + condaRelease() Modified: trunk/setup.py =================================================================== --- trunk/setup.py 2014-11-11 03:05:43 UTC (rev 4955) +++ trunk/setup.py 2014-11-12 02:06:43 UTC (rev 4956) @@ -835,7 +835,7 @@ url = "http://simupop.sourceforge.net", description = "Forward-time population genetics simulation environment", long_description = DESCRIPTION, - download_url = 'http://simupop.sourceforge.net/Main/Download', + download_url = 'http://sourceforge.net/projects/simupop/files/simupop/', classifiers = [ 'Development Status :: 5 - Production/Stable', 'Environment :: Console', This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bpe...@us...> - 2014-11-11 03:05:46
|
Revision: 4955 http://sourceforge.net/p/simupop/code/4955 Author: bpeng2000 Date: 2014-11-11 03:05:43 +0000 (Tue, 11 Nov 2014) Log Message: ----------- Add conda distributing scripts Added Paths: ----------- trunk/development/conda/ trunk/development/conda/bld.bat trunk/development/conda/build.sh trunk/development/conda/meta.yaml Added: trunk/development/conda/bld.bat =================================================================== --- trunk/development/conda/bld.bat (rev 0) +++ trunk/development/conda/bld.bat 2014-11-11 03:05:43 UTC (rev 4955) @@ -0,0 +1,8 @@ +"%PYTHON%" setup.py install +if errorlevel 1 exit 1 + +:: Add more build steps here, if they are necessary. + +:: See +:: http://docs.continuum.io/conda/build.html +:: for a list of environment variables that are set during the build process. Added: trunk/development/conda/build.sh =================================================================== --- trunk/development/conda/build.sh (rev 0) +++ trunk/development/conda/build.sh 2014-11-11 03:05:43 UTC (rev 4955) @@ -0,0 +1,9 @@ +#!/bin/bash + +$PYTHON setup.py install + +# Add more build steps here, if they are necessary. + +# See +# http://docs.continuum.io/conda/build.html +# for a list of environment variables that are set during the build process. Added: trunk/development/conda/meta.yaml =================================================================== --- trunk/development/conda/meta.yaml (rev 0) +++ trunk/development/conda/meta.yaml 2014-11-11 03:05:43 UTC (rev 4955) @@ -0,0 +1,63 @@ +package: + name: simupop + version: !!str 1.0.0 + +source: + fn: simuPOP-1.1.4.tar.gz + url: https://pypi.python.org/packages/source/s/simuPOP/simuPOP-1.1.4.tar.gz#md5=f157cd3902747d72c59fd8460a78e3fc + md5: +# patches: + # List any patch files here + # - fix.patch + +# build: + # preserve_egg_dir: True + # entry_points: + # Put any entry points (scripts to be generated automatically) here. The + # syntax is module:function. For example + # + # - simupop = simupop:main + # + # Would create an entry point called simupop that calls simupop.main() + + + # If this is a new build for the same version, increment the build + # number. If you do not include this key, it defaults to 0. + # number: 1 + +requirements: + build: + - python + - setuptools + - distribute + + run: + - python + - distribute + +test: + # Python imports + imports: + - simupop + - src + + # commands: + # You can put test commands to be run here. Use this to test that the + # entry points work. + + + # You can also put a file called run_test.py in the recipe that will be run + # at test time. + + # requires: + # Put any additional test requirements here. For example + # - nose + +about: + home: http://simupop.sourceforge.net + license: GNU General Public License (GPL) + summary: 'Forward-time population genetics simulation environment' + +# See +# http://docs.continuum.io/conda/build.html for +# more information about meta.yaml This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bpe...@us...> - 2014-10-17 15:40:10
|
Revision: 4954 http://sourceforge.net/p/simupop/code/4954 Author: bpeng2000 Date: 2014-10-17 15:40:08 +0000 (Fri, 17 Oct 2014) Log Message: ----------- Small fix for PRECISION Modified Paths: -------------- trunk/src/utility.h Modified: trunk/src/utility.h =================================================================== --- trunk/src/utility.h 2014-10-16 21:21:49 UTC (rev 4953) +++ trunk/src/utility.h 2014-10-17 15:40:08 UTC (rev 4954) @@ -90,7 +90,7 @@ // this is used to compare loci positions when loci are provided by // (chr, pos) pair. -#define PRECISION(f) (double(int(f)) == f ? f : (int((f) * 10000. + 0.5) / 10000.)) +#define PRECISION(f) (double(size_t(f)) == f ? f : (size_t((f) * 10000. + 0.5) / 10000.)) namespace simuPOP { This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bpe...@us...> - 2014-10-16 21:21:53
|
Revision: 4953 http://sourceforge.net/p/simupop/code/4953 Author: bpeng2000 Date: 2014-10-16 21:21:49 +0000 (Thu, 16 Oct 2014) Log Message: ----------- Import the initialization time for long genomic region Modified Paths: -------------- trunk/src/transmitter.cpp Modified: trunk/src/transmitter.cpp =================================================================== --- trunk/src/transmitter.cpp 2014-10-15 14:57:44 UTC (rev 4952) +++ trunk/src/transmitter.cpp 2014-10-16 21:21:49 UTC (rev 4953) @@ -674,15 +674,15 @@ // get loci distance * m_rates and then recombinant points for (size_t loc = chBegin; loc < chEnd - 1; ++loc) { // if this locus will be recombined. - vectoru::const_iterator pos = find(loci.begin(), loci.end(), loc); - if (pos != loci.end()) { + size_t pos = m_loci.indexOf(loc); + if (pos != NOT_FOUND) { double r = 0; if (useLociDist) r = m_intensity > 0 ? ((ind.locusPos(loc + 1) - ind.locusPos(loc)) * m_intensity) : r; else if (m_rates.size() == 1 && !useLociDist) r = max(m_rates[0], 0.); else - r = m_rates[pos - loci.begin()]; + r = m_rates[pos]; m_recBeforeLoci.push_back(loc + 1); vecP.push_back(min(0.5, r)); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bpe...@us...> - 2014-10-15 14:57:47
|
Revision: 4952 http://sourceforge.net/p/simupop/code/4952 Author: bpeng2000 Date: 2014-10-15 14:57:44 +0000 (Wed, 15 Oct 2014) Log Message: ----------- Update reference manual Modified Paths: -------------- trunk/doc/refManual.lyx trunk/doc/refManual.pdf trunk/doc/userGuide.pdf Modified: trunk/doc/refManual.lyx =================================================================== --- trunk/doc/refManual.lyx 2014-10-15 14:46:55 UTC (rev 4951) +++ trunk/doc/refManual.lyx 2014-10-15 14:57:44 UTC (rev 4952) @@ -2049,13 +2049,6 @@ \backslash -simuPOPdemographyBaseDemographicModelRef -\end_layout - -\begin_layout Plain Layout - - -\backslash simuPOPdemographyEventBasedModelRef \end_layout Modified: trunk/doc/refManual.pdf =================================================================== (Binary files differ) Modified: trunk/doc/userGuide.pdf =================================================================== (Binary files differ) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bpe...@us...> - 2014-10-15 14:47:05
|
Revision: 4951 http://sourceforge.net/p/simupop/code/4951 Author: bpeng2000 Date: 2014-10-15 14:46:55 +0000 (Wed, 15 Oct 2014) Log Message: ----------- Prepare a 1.1.4 release Modified Paths: -------------- trunk/doc/refManual.lyx trunk/doc/userGuide.lyx trunk/simuPOP_version.py trunk/src/demography.py trunk/test/test_19_demography.py Modified: trunk/doc/refManual.lyx =================================================================== --- trunk/doc/refManual.lyx 2014-10-15 05:19:11 UTC (rev 4950) +++ trunk/doc/refManual.lyx 2014-10-15 14:46:55 UTC (rev 4951) @@ -5,7 +5,7 @@ \textclass manual \begin_preamble \renewcommand{\py@ptsize}{12pt} -\setreleaseinfo{Release 1.1.3 (\mbox{Rev: 4913})} +\setreleaseinfo{Release 1.1.4 (\mbox{Rev: 4950})} % file revision $Rev: 3372$ \authoraddress{ {\bf Department of Epidemiology, U.T. M.D. Anderson Cancer Center}\\ Modified: trunk/doc/userGuide.lyx =================================================================== --- trunk/doc/userGuide.lyx 2014-10-15 05:19:11 UTC (rev 4950) +++ trunk/doc/userGuide.lyx 2014-10-15 14:46:55 UTC (rev 4951) @@ -6,7 +6,7 @@ \begin_preamble \renewcommand{\py@ptsize}{12pt} -\setreleaseinfo{Release 1.1.3 (\mbox{Rev: 4913})} +\setreleaseinfo{Release 1.1.4 (\mbox{Rev: 4950})} % file revision $Rev: 3372$ \authoraddress{ {\bf Department of Epidemiology, U.T. M.D. Anderson Cancer Center}\\ @@ -26113,27 +26113,27 @@ \end_layout \begin_layout Standard -Finally, an AdmixEvent mix two or more subpopulations by certain proportions, - and either create a new subpopulation or replace an existing subpopulation. +Finally, an +\family typewriter +AdmixtureEvent +\family default + mix two or more subpopulations by certain proportions, and either create + a new subpopulation or replace an existing subpopulation. In particular, \end_layout \begin_layout Standard - -\family typewriter \begin_inset listings inline false status open \begin_layout Plain Layout -AdmixEvent(subPops=['MX', 'EU'], at=-10, sizes=[0.4, 0.6], names='MXL') +AdmixtureEvent(subPops=['MX', 'EU'], at=-10, sizes=[0.4, 0.6], name='MXL') \end_layout \end_inset - -\family default creates a new admixed population called MXL with 40% of individuals from the MX population, and the rest from the EU population. The admixture process happens once and follows an Hybrid Isolation model. @@ -26147,7 +26147,7 @@ \begin_layout Plain Layout -AdmixEvent(subPops=['MX', 'EU'], begin=-10, sizes=[0.8, 0.2], toSubPop='MX') +AdmixtureEvent(subPops=['MX', 'EU'], begin=-10, sizes=[0.8, 0.2], toSubPop='MX') \end_layout \end_inset @@ -26156,9 +26156,41 @@ 10 generations. Because 20% of the admixed population will be replaced by individuals from the EU population, this models a continuous gene flow model of admixture. + If you would like to control the exact size of the admixed population, + you can specify the number of individuals as integer numbers instead of + proportions: \end_layout \begin_layout Standard +\begin_inset listings +inline false +status open + +\begin_layout Plain Layout + +AdmixtureEvent(subPops=['MX', 'EU'], begin=-10, sizes=[int(1400*0.8), int(1400*0.2 +)], toSubPop='MX') +\end_layout + +\end_inset + +Note that the type of elements in parameter +\family typewriter +sizes +\family default + is important, +\family typewriter +1. + +\family default + stands for all subpopulation and +\family typewriter +1 +\family default + stands for one individual from it. +\end_layout + +\begin_layout Standard Example \begin_inset CommandInset ref LatexCommand ref Modified: trunk/simuPOP_version.py =================================================================== --- trunk/simuPOP_version.py 2014-10-15 05:19:11 UTC (rev 4950) +++ trunk/simuPOP_version.py 2014-10-15 14:46:55 UTC (rev 4951) @@ -1,2 +1,2 @@ -SIMUPOP_VER="1.1.4svn" -SIMUPOP_REV="4913" +SIMUPOP_VER="1.1.4" +SIMUPOP_REV="4950" Modified: trunk/src/demography.py =================================================================== --- trunk/src/demography.py 2014-10-15 05:19:11 UTC (rev 4950) +++ trunk/src/demography.py 2014-10-15 14:46:55 UTC (rev 4951) @@ -454,29 +454,30 @@ plt.savefig(filename) plt.close() - def _checkSize(self, pop): + def _checkSize(self, pop, param): gen = pop.dvars().gen - if gen in self.intended_size: - sz = self.intended_size[gen] + if gen - param in self.intended_size: + sz = self.intended_size[gen - param] if isinstance(sz, int): sz = (sz,) else: sz = tuple(sz) if sz != pop.subPopSizes(): - raise ValueError('Mismatch population size at generation %s: observed=%s, intended=%s' % \ - (gen, pop.subPopSizes(), sz)) + raise ValueError('Mismatch population size at generation %s (with starting gen %s): observed=%s, intended=%s' % \ + (gen - param, param, pop.subPopSizes(), sz)) return True - def _assertSize(self, sizes, initSize=[]): + def _assertSize(self, sizes, initSize=[], startGen=0): '''This function is provided for testing purpose. ''' self.intended_size = sizes pop = Population(size=initSize if initSize else self.init_size, infoFields=self.info_fields) + pop.dvars().gen = startGen pop.evolve( matingScheme=RandomSelection(subPopSize=self), - postOps=PyOperator(self._checkSize), - finalOps=PyOperator(self._checkSize), + postOps=PyOperator(self._checkSize, param=startGen), + finalOps=PyOperator(self._checkSize, param=startGen), gen=self.num_gens ) Modified: trunk/test/test_19_demography.py =================================================================== --- trunk/test/test_19_demography.py 2014-10-15 05:19:11 UTC (rev 4950) +++ trunk/test/test_19_demography.py 2014-10-15 14:46:55 UTC (rev 4951) @@ -23,49 +23,50 @@ sys.argv=new_argv from simuPOP import * from time import sleep +import random from simuPOP.demography import * class TestDemography(unittest.TestCase): def testInstantChangeModel(self): - InstantChangeModel(T=100, N0=200)._assertSize( + InstantChangeModel(T=100, N0=200)._assertSize(startGen=random.randint(0,100), sizes= { 0: 200, 99: 200, }) - InstantChangeModel(T=100, N0=[200, 300])._assertSize( + InstantChangeModel(T=100, N0=[200, 300])._assertSize(startGen=random.randint(0,100), sizes= { 0: [200, 300], 99: [200, 300] }) - InstantChangeModel(T=100, N0=200, G=10, NG=30)._assertSize( + InstantChangeModel(T=100, N0=200, G=10, NG=30)._assertSize(startGen=random.randint(0,100), sizes= { 9: 200, 10: 30, 99: 30 }) - InstantChangeModel(T=100, N0=200, G=10, NG=30)._assertSize( + InstantChangeModel(T=100, N0=200, G=10, NG=30)._assertSize(startGen=random.randint(0,100), sizes= { 0: 200, 9: 200, 10: 30, 99: 30 }, initSize=500) - InstantChangeModel(T=100, N0=200, G=10, NG=[30, 40])._assertSize( + InstantChangeModel(T=100, N0=200, G=10, NG=[30, 40])._assertSize(startGen=random.randint(0,100), sizes= { 0: 200, 9: 200, 10: (30, 40), 99: (30, 40), }, initSize=500) - InstantChangeModel(T=100, N0=[200, 300], G=10, NG=[30, 40])._assertSize( + InstantChangeModel(T=100, N0=[200, 300], G=10, NG=[30, 40])._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300), 9: (200, 300), 10: (30, 40), 99: (30, 40), }, initSize=500) - InstantChangeModel(T=100, N0=[200, 300], G=10, NG=[[30, 40]])._assertSize( + InstantChangeModel(T=100, N0=[200, 300], G=10, NG=[[30, 40]])._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300), 9: (200, 300), @@ -74,7 +75,7 @@ }, initSize=500) self.assertRaises(ValueError, InstantChangeModel, T=100, N0=[200, 300], G=[10], NG=[30, 40]) - InstantChangeModel(T=100, N0=[200, 300], G=(10, 20), NG=[30, 40])._assertSize( + InstantChangeModel(T=100, N0=[200, 300], G=(10, 20), NG=[30, 40])._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300), 9: (200, 300), @@ -82,7 +83,7 @@ 20: 40, 99: 40, }, initSize=500) - InstantChangeModel(T=100, N0=[200, 300], G=(10, 20), NG=[30, [40, 50]])._assertSize( + InstantChangeModel(T=100, N0=[200, 300], G=(10, 20), NG=[30, [40, 50]])._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300), 9: (200, 300), @@ -93,7 +94,7 @@ self.assertRaises(ValueError, InstantChangeModel, T=100, N0=[200, 300], G=(10, 20), NG=[30, [40, 50], 40], ) # testing the use of proportions - InstantChangeModel(T=100, N0=[200, 300], G=(10, 20), NG=[0.5, [0.4, 0.5]])._assertSize( + InstantChangeModel(T=100, N0=[200, 300], G=(10, 20), NG=[0.5, [0.4, 0.5]])._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300), 9: (200, 300), @@ -101,7 +102,7 @@ 20: (100, 125), 99: (100, 125) }) - InstantChangeModel(T=100, G=(10, 20), NG=[0.5, [0.4, 0.5]])._assertSize( + InstantChangeModel(T=100, G=(10, 20), NG=[0.5, [0.4, 0.5]])._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300), 9: (200, 300), @@ -109,7 +110,7 @@ 20: (100, 125), 99: (100, 125) }, initSize=[200, 300]) - InstantChangeModel(T=100, G=(10, 20), NG=[0.5, [0.4, None]])._assertSize( + InstantChangeModel(T=100, G=(10, 20), NG=[0.5, [0.4, None]])._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300), 9: (200, 300), @@ -117,7 +118,7 @@ 20: (100, 250), 99: (100, 250) }, initSize=[200, 300]) - InstantChangeModel(T=100, G=(10, 20), NG=[0.5, [None, 50]])._assertSize( + InstantChangeModel(T=100, G=(10, 20), NG=[0.5, [None, 50]])._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300), 9: (200, 300), @@ -126,39 +127,39 @@ 99: (250, 50) }, initSize=[200, 300]) # test the removeEmptySubPops parameter - InstantChangeModel(T=0, N0=[100, 0])._assertSize( + InstantChangeModel(T=0, N0=[100, 0])._assertSize(startGen=random.randint(0,100), sizes= { 0: (100, 0), }, initSize=[200, 300]) - InstantChangeModel(T=0, N0=[0, 100], removeEmptySubPops=True)._assertSize( + InstantChangeModel(T=0, N0=[0, 100], removeEmptySubPops=True)._assertSize(startGen=random.randint(0,100), sizes= { 0: (100,), }, initSize=[200, 300]) - InstantChangeModel(T=0, removeEmptySubPops=True)._assertSize( + InstantChangeModel(T=0, removeEmptySubPops=True)._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300, 400, 500), }, initSize=[200, 300, 400, 0, 0, 500]) def testAdmixtureModel(self): - AdmixtureModel(T=10, model=('HI', 1, 2, 0.5))._assertSize( + AdmixtureModel(T=10, model=('HI', 1, 2, 0.5))._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300, 400, 600), 4: (200, 300, 400, 600), 5: (200, 300, 400, 600), }, initSize=[200, 300, 400]) - AdmixtureModel(T=10, model=('HI', 0, 2, 0.5))._assertSize( + AdmixtureModel(T=10, model=('HI', 0, 2, 0.5))._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300, 400, 400), 4: (200, 300, 400, 400), 5: (200, 300, 400, 400), }, initSize=[200, 300, 400]) - AdmixtureModel(T=10, model=('CGF', 0, 2, 0.8))._assertSize( + AdmixtureModel(T=10, model=('CGF', 0, 2, 0.8))._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300, 400), 4: (200, 300, 400), 5: (200, 300, 400), }, initSize=[200, 300, 400]) - AdmixtureModel(T=10, model=('CGF', 1, 2, 0.8))._assertSize( + AdmixtureModel(T=10, model=('CGF', 1, 2, 0.8))._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300, 400), 4: (200, 300, 400), @@ -169,46 +170,46 @@ self.assertRaises(LinearGrowthModel) self.assertRaises(LinearGrowthModel, T=100, N0=200) # - LinearGrowthModel(T=100, N0=200, NT=400)._assertSize( + LinearGrowthModel(T=100, N0=200, NT=400)._assertSize(startGen=random.randint(0,100), sizes= { 0: 202, 1: 204, 99: 400, }, initSize=500) - LinearGrowthModel(T=100, N0=200, r=0.01)._assertSize( + LinearGrowthModel(T=100, N0=200, r=0.01)._assertSize(startGen=random.randint(0,100), sizes= { 0: 202, 2: 206, 99: 400, }, initSize=500) - LinearGrowthModel(N0=200, r=0.01, NT=400)._assertSize( + LinearGrowthModel(N0=200, r=0.01, NT=400)._assertSize(startGen=random.randint(0,100), sizes= { 0: 202, 2: 206, 99: 400, }, initSize=500) # - LinearGrowthModel(T=100, N0=[200, 200], NT=[400, 800])._assertSize( + LinearGrowthModel(T=100, N0=[200, 200], NT=[400, 800])._assertSize(startGen=random.randint(0,100), sizes= { 0: [202, 206], 1: [204, 212], 99: [400, 800], }, initSize=500) - LinearGrowthModel(T=100, N0=[200, 200], r=0.01)._assertSize( + LinearGrowthModel(T=100, N0=[200, 200], r=0.01)._assertSize(startGen=random.randint(0,100), sizes= { 0: [202, 202], 2: [206, 206], 99: [400, 400], }, initSize=500) # unknown NT - LinearGrowthModel(T=100, N0=[200, 200], r=[0.01, 0.03])._assertSize( + LinearGrowthModel(T=100, N0=[200, 200], r=[0.01, 0.03])._assertSize(startGen=random.randint(0,100), sizes= { 0: [202, 206], 2: [206, 218], 99: [400, 800], }, initSize=500) # unknown T - LinearGrowthModel(N0=[200, 200], r=0.01, NT=[400, 800])._assertSize( + LinearGrowthModel(N0=[200, 200], r=0.01, NT=[400, 800])._assertSize(startGen=random.randint(0,100), sizes= { 0: [202, 202], 2: [206, 206], @@ -218,7 +219,7 @@ 300: [400, 800], }, initSize=500) # N0 size ... - LinearGrowthModel(N0=[1., 1.], r=0.01, NT=[400, 800])._assertSize( + LinearGrowthModel(N0=[1., 1.], r=0.01, NT=[400, 800])._assertSize(startGen=random.randint(0,100), sizes= { 0: [202, 202], 2: [206, 206], @@ -227,7 +228,7 @@ 298: [400, 798], 300: [400, 800], }, initSize=200) - LinearGrowthModel(N0=[(None, 'name'), 1.], r=0.01, NT=[400, 800])._assertSize( + LinearGrowthModel(N0=[(None, 'name'), 1.], r=0.01, NT=[400, 800])._assertSize(startGen=random.randint(0,100), sizes= { 0: [202, 202], 2: [206, 206], @@ -240,7 +241,7 @@ self.assertRaises(LinearGrowthModel, N0=[1., 1.], r=0.01, NT=[2000]) self.assertRaises(LinearGrowthModel, N0=[1., 1.], r=[0.01], NT=[2000, 3000]) # carrying capacity - LinearGrowthModel(T=100, N0=[200, 200], r=0.01, NT=[300, 350])._assertSize( + LinearGrowthModel(T=100, N0=[200, 200], r=0.01, NT=[300, 350])._assertSize(startGen=random.randint(0,100), sizes= { 0: [202, 202], 2: [206, 206], @@ -252,59 +253,59 @@ self.assertRaises(ExponentialGrowthModel) self.assertRaises(ExponentialGrowthModel, T=100, N0=200) # - ExponentialGrowthModel(T=100, N0=200, NT=400)._assertSize( + ExponentialGrowthModel(T=100, N0=200, NT=400)._assertSize(startGen=random.randint(0,100), sizes= { 0: 201, 1: 203, 99: 400, }, initSize=500) - ExponentialGrowthModel(T=100, N0=200, r=0.01)._assertSize( + ExponentialGrowthModel(T=100, N0=200, r=0.01)._assertSize(startGen=random.randint(0,100), sizes= { 0: 202, 2: 206, 99: 544, }, initSize=500) - ExponentialGrowthModel(N0=200, r=0.01, NT=400)._assertSize( + ExponentialGrowthModel(N0=200, r=0.01, NT=400)._assertSize(startGen=random.randint(0,100), sizes= { 0: 202, 2: 206, 99: 544, }, initSize=500) # - ExponentialGrowthModel(T=100, N0=[200, 200], NT=[400, 800])._assertSize( + ExponentialGrowthModel(T=100, N0=[200, 200], NT=[400, 800])._assertSize(startGen=random.randint(0,100), sizes= { 0: [201, 203], 1: [203, 206], 99: [400, 800], }, initSize=500) - ExponentialGrowthModel(T=100, N0=[200, 200], r=0.01)._assertSize( + ExponentialGrowthModel(T=100, N0=[200, 200], r=0.01)._assertSize(startGen=random.randint(0,100), sizes= { 0: [202, 202], 2: [206, 206], 99: [544, 544], }, initSize=500) # unknown NT - ExponentialGrowthModel(T=100, N0=[200, 200], r=[0.01, 0.03])._assertSize( + ExponentialGrowthModel(T=100, N0=[200, 200], r=[0.01, 0.03])._assertSize(startGen=random.randint(0,100), sizes= { 0: [202, 206], 2: [206, 219], 99: [544, 4017], }, initSize=500) # unknown T - ExponentialGrowthModel(N0=[200, 200], r=0.01, NT=[400, 800])._assertSize( + ExponentialGrowthModel(N0=[200, 200], r=0.01, NT=[400, 800])._assertSize(startGen=random.randint(0,100), sizes= { 2: [206, 206], 99: [400, 544], 100: [400, 549], }, initSize=500) # N0 size ... - ExponentialGrowthModel(N0=[1., 1.], r=0.01, NT=[400, 800])._assertSize( + ExponentialGrowthModel(N0=[1., 1.], r=0.01, NT=[400, 800])._assertSize(startGen=random.randint(0,100), sizes= { 2: [206, 206], 99: [400, 544], 100: [400, 549], }, initSize=200) - ExponentialGrowthModel(N0=[(None, 'name'), 1.], r=0.01, NT=[400, 800])._assertSize( + ExponentialGrowthModel(N0=[(None, 'name'), 1.], r=0.01, NT=[400, 800])._assertSize(startGen=random.randint(0,100), sizes= { 2: [206, 206], 99: [400, 544], @@ -314,7 +315,7 @@ self.assertRaises(ExponentialGrowthModel, N0=[1., 1.], r=0.01, NT=[2000]) self.assertRaises(ExponentialGrowthModel, N0=[1., 1.], r=[0.01], NT=[2000, 3000]) # carrying capacity - ExponentialGrowthModel(T=100, N0=[200, 200], r=0.01, NT=[300, 350])._assertSize( + ExponentialGrowthModel(T=100, N0=[200, 200], r=0.01, NT=[300, 350])._assertSize(startGen=random.randint(0,100), sizes= { 0: [202, 202], 2: [206, 206], @@ -326,7 +327,7 @@ MultiStageModel([ InstantChangeModel(T=1000, N0=1000, G=[500, 600], NG=[100, 1000]), ExponentialGrowthModel(T=100, NT=10000) - ])._assertSize( + ])._assertSize(startGen=random.randint(0,100), sizes= { 100: 1000, 500: 100, @@ -340,15 +341,15 @@ LinearGrowthModel(T=100, N0=1000, r=0.01), ExponentialGrowthModel(T=100, N0=[0.4, 0.6], r=0.001), ExponentialGrowthModel(r=0.01, NT=[2000, 4000]) - ])._assertSize( + ])._assertSize(startGen=random.randint(0,100), sizes= { 0: 1010, 99: 2000, - 100: [800, 1201], - 102: [802, 1203], + 100: [801, 1201], + 102: [802, 1204], 199: [884, 1326], 250: [1472, 2208], - 300: [2000, 3640], + 300: [2000, 3641], 310: [2000, 4000], }) MultiStageModel([ @@ -356,15 +357,15 @@ ExponentialGrowthModel(T=100, N0=[0.4, 0.6], r=0.001), ExponentialGrowthModel(r=0.01, NT=[2000, 4000]), InstantChangeModel(N0=[1000, 0], T=100) - ])._assertSize( + ])._assertSize(startGen=random.randint(0,100), sizes= { 0: 1010, 99: 2000, - 100: [800, 1201], - 102: [802, 1203], + 100: [801, 1201], + 102: [802, 1204], 199: [884, 1326], 250: [1472, 2208], - 300: [2000, 3640], + 300: [2000, 3641], 310: [2000, 4000], 350: [1000, 0] }) @@ -374,15 +375,15 @@ ExponentialGrowthModel(r=0.01, NT=[2000, 4000]), AdmixtureModel(model=['HI', 0, 1, 0.3], T=1), InstantChangeModel(N0=[0, 0, None], T=100) - ])._assertSize( + ])._assertSize(startGen=random.randint(0,100), sizes= { 0: 1010, 99: 2000, - 100: [800, 1201], - 102: [802, 1203], + 100: [801, 1201], + 102: [802, 1204], 199: [884, 1326], 250: [1472, 2208], - 300: [2000, 3640], + 300: [2000, 3641], 310: [2000, 4000], 311: [2000, 4000, 5714], 350: [0, 0, 5714] @@ -393,15 +394,15 @@ ExponentialGrowthModel(r=0.01, NT=[2000, 4000]), AdmixtureModel(model=['HI', 0, 1, 0.3], T=1), InstantChangeModel(N0=[0, 0, None], removeEmptySubPops=True, T=1) - ])._assertSize( + ])._assertSize(startGen=random.randint(0,100), sizes= { 0: 1010, 99: 2000, - 100: [800, 1201], - 102: [802, 1203], + 100: [801, 1201], + 102: [802, 1204], 199: [884, 1326], 250: [1472, 2208], - 300: [2000, 3640], + 300: [2000, 3641], 310: [2000, 4000], 311: [2000, 4000, 5714], 312: [5714] @@ -410,22 +411,22 @@ LinearGrowthModel(T=100, N0=1000, r=0.01), InstantChangeModel(N0=[1000, 1000], T=0), ExponentialGrowthModel(T=100, N0=[0.4, 0.6], r=0.001), - ])._assertSize( + ])._assertSize(startGen=random.randint(0,100), sizes= { 0: 1010, 99: 2000, - 100: [400, 600], + 100: [400, 601], }) MultiStageModel([ LinearGrowthModel(T=100, N0=1000, r=0.01), InstantChangeModel(N0=[1000, 1000], T=0), ExponentialGrowthModel(T=100, N0=[0.4, 0.6], r=0.001), InstantChangeModel(N0=[1000, 1000], T=1), - ])._assertSize( + ])._assertSize(startGen=random.randint(0,100), sizes= { 0: 1010, 99: 2000, - 100: [400, 600], + 100: [400, 601], 199: [442, 663], 200: [1000, 1000], }) @@ -529,7 +530,7 @@ def testEventBasedModel(self): 'Test event based implementation of demographic models' - EventBasedModel(T=100, N0=200)._assertSize( + EventBasedModel(T=100, N0=200)._assertSize(startGen=random.randint(0,100), sizes= { 0: 200, 99: 200, @@ -539,7 +540,7 @@ ResizeEvent(sizes=300, subPops=1, names='AF', at=2), MergeEvent(subPops=[1,2], at=3), ] - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 200, 200), 2: (200, 300, 200), @@ -551,7 +552,7 @@ ResizeEvent(sizes=400, subPops=2, at=2), ResizeEvent(sizes=350, subPops=1, at=2), ] - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 200, 200), 2: (200, 350, 400), @@ -563,7 +564,7 @@ MergeEvent(subPops=[1,0], at=2), ResizeEvent(sizes=350, subPops=1, at=2), ] - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 200, 200), 2: (500, 350), @@ -575,7 +576,7 @@ def testResizeEvent(self): EventBasedModel(T=100, N0=(200, 200, 200), events=ResizeEvent(sizes=300, subPops=1, names='AF', at=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 200, 200), 2: (200, 300, 200), @@ -583,7 +584,7 @@ }, initSize=500) EventBasedModel(T=100, N0=(200, 200, 200), events=ResizeEvent(sizes=1.5, subPops=1, names='AF', at=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 200, 200), 2: (200, 300, 200), @@ -591,7 +592,7 @@ }, initSize=500) EventBasedModel(T=100, N0=(200, 200, 200), events=ResizeEvent(sizes=(1.5, 1.5, 1.8), names=('AF', 'A', 'B'), at=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 200, 200), 2: (300, 300, 360), @@ -599,7 +600,7 @@ }, initSize=500) EventBasedModel(T=100, N0=(200, 200, 200), events=ResizeEvent(sizes=(300, 1.5), subPops=[0, 1], at=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 200, 200), 2: (300, 300, 200), @@ -607,7 +608,7 @@ }, initSize=500) EventBasedModel(T=100, N0=(200, 200, 200), events=ResizeEvent(sizes=(300, 0.9), subPops=[0, 1], begin=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 200, 200), 2: (300, 180, 200), @@ -615,7 +616,7 @@ }, initSize=500) EventBasedModel(T=100, N0=(200, 200, 200), events=ResizeEvent(sizes=(300, 0), subPops=[0, 1], at=2, removeEmptySubPops=True) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 200, 200), 2: (300, 200), @@ -626,7 +627,7 @@ def testMergeEvent(self): EventBasedModel(T=100, N0=(200, 200, 200), events=MergeEvent(name='AF', at=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 200, 200), 2: (600,), @@ -634,7 +635,7 @@ }, initSize=500) EventBasedModel(T=100, N0=((100, 'AF'), (200, 'EU'), (300, 'AS')), events=MergeEvent(name='AF', subPops=('AF', 'AS'), at=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (100, 200, 300), 2: (400, 200), @@ -642,7 +643,7 @@ }, initSize=500) EventBasedModel(T=100, N0=((100, 'AF'), (200, 'EU'), (300, 'AS')), events=MergeEvent(name='EU', subPops=('EU', 'AS'), at=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (100, 200, 300), 2: (100, 500), @@ -652,7 +653,7 @@ def testExpansionEvent(self): EventBasedModel(T=100, N0=200, events=ExpansionEvent(rates=0.01) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: 202, 2: 206, @@ -660,7 +661,7 @@ }, initSize=500) EventBasedModel(T=100, N0=(200, 'AF'), events=ExpansionEvent(rates=0.01, subPops='AF') - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: 202, 2: 206, @@ -668,7 +669,7 @@ }, initSize=500) EventBasedModel(T=100, N0=[200, 200], events=ExpansionEvent(rates=0.01) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: [202, 202], 2: [206, 206], @@ -676,7 +677,7 @@ }, initSize=500) EventBasedModel(T=100, N0=[(200, 'A'), (200, 'B')], events=ExpansionEvent(rates=0.01, subPops=('A', 'B')) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: [202, 202], 2: [206, 206], @@ -685,7 +686,7 @@ # unknown NT EventBasedModel(T=100, N0=[200, 200], events=ExpansionEvent(rates=[0.01, 0.03]) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: [202, 206], 2: [206, 219], @@ -693,7 +694,7 @@ }, initSize=500) EventBasedModel(T=100, N0=[(200, 'AF'), (200, 'EU'), (200, 'AS')], events=ExpansionEvent(rates=[0.01, 0.02], subPops=['EU', 'AS']) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: [200, 202, 204], 2: [200, 206, 212], @@ -701,7 +702,7 @@ }, initSize=500) EventBasedModel(T=100, N0=200, events=ExpansionEvent(slopes=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: 202, 2: 206, @@ -709,7 +710,7 @@ }, initSize=500) EventBasedModel(T=100, N0=[200, 200], events=ExpansionEvent(slopes=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: [202, 202], 2: [206, 206], @@ -717,7 +718,7 @@ }, initSize=500) EventBasedModel(T=100, N0=[200, 200], events=ExpansionEvent(slopes=[2, 6]) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: [202, 206], 2: [206, 218], @@ -728,7 +729,7 @@ def testSplitEvent(self): EventBasedModel(T=100, N0=200, events=SplitEvent(sizes=(0.4, 0.6), at=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: [200], 2: [80, 120], @@ -736,7 +737,7 @@ }, initSize=500) EventBasedModel(T=100, N0=200, events=SplitEvent(sizes=(0.4, 250), at=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: [200], 2: [80, 250], @@ -747,7 +748,7 @@ )._assertSize, { }, initSize=500) EventBasedModel(T=100, N0=((100, 'A'), (200, 'B')), events=SplitEvent(sizes=(0.4, 250), at=2, subPops='B') - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: [100, 200], 2: [100, 80, 250], @@ -757,7 +758,7 @@ def testAdmixtureEvent(self): EventBasedModel(T=10, events=AdmixtureEvent(subPops=[1, 2], sizes=[0.5, 0.5], at=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300, 400), 2: (200, 300, 400, 600), @@ -766,7 +767,7 @@ }, initSize=[200, 300, 400]) EventBasedModel(T=10, events=AdmixtureEvent(subPops=[0, 2], sizes=[0.5, 0.5], at=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300, 400), 2: (200, 300, 400, 400), @@ -775,7 +776,7 @@ }, initSize=[200, 300, 400]) EventBasedModel(T=10, events=AdmixtureEvent(subPops=[0, 2], sizes=[0.5, 0.1], at=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300, 400), 2: (200, 300, 400, 400), @@ -784,7 +785,7 @@ }, initSize=[200, 300, 400]) EventBasedModel(T=10, events=AdmixtureEvent(subPops=[0, 1, 2], sizes=[0.2, 0.3, 0.5], at=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300, 400), # take 160, 240, 400 @@ -795,7 +796,7 @@ # EventBasedModel(T=10, events=AdmixtureEvent(subPops=[0, 1, 2], sizes=[20, 30, 50], at=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300, 400), # take 160, 240, 400 @@ -805,7 +806,7 @@ }, initSize=[200, 300, 400]) EventBasedModel(T=10, events=AdmixtureEvent(subPops=[0, 1, 2], sizes=[20, 30, 50], begin=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300, 400), 2: (200, 300, 400, 100), @@ -815,7 +816,7 @@ # with toSubPop ... EventBasedModel(T=10, events=AdmixtureEvent(subPops=[1, 2], sizes=[0.5, 0.5], at=2, toSubPop=1) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300, 400), 2: (200, 300, 400), @@ -824,7 +825,7 @@ }, initSize=[200, 300, 400]) EventBasedModel(T=10, events=AdmixtureEvent(subPops=[0, 1, 2], sizes=[20, 30, 50], toSubPop=1, at=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300, 400), 2: (200, 100, 400), @@ -833,7 +834,7 @@ }, initSize=[200, 300, 400]) EventBasedModel(T=10, events=AdmixtureEvent(subPops=[0, 1, 2], sizes=[20, 30, 50], toSubPop=1, begin=2) - )._assertSize( + )._assertSize(startGen=random.randint(0,100), sizes= { 0: (200, 300, 400), 2: (200, 100, 400), This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bpe...@us...> - 2014-10-15 05:19:16
|
Revision: 4950 http://sourceforge.net/p/simupop/code/4950 Author: bpeng2000 Date: 2014-10-15 05:19:11 +0000 (Wed, 15 Oct 2014) Log Message: ----------- Fine-tune exponential growth models to make sure two implementations give identical results Modified Paths: -------------- trunk/src/demography.py trunk/test/test_19_demography.py Modified: trunk/src/demography.py =================================================================== --- trunk/src/demography.py 2014-10-15 03:30:30 UTC (rev 4949) +++ trunk/src/demography.py 2014-10-15 05:19:11 UTC (rev 4950) @@ -630,7 +630,7 @@ # elif self.NT is None: # get final population size from T and r - self.NT = [int(round(x*((1.+self.r)**self.num_gens))) for x in self.init_size] + self.NT = [x*math.exp(self.r*self.num_gens) for x in self.init_size] elif None in self.NT or \ any([isinstance(x, float) for x in self.NT]): raise ValueError('Relative ending population size is not allowed' @@ -641,7 +641,7 @@ raise ValueError('Please specify growth rate for each subpopulation ' 'if multiple growth rates are specified.') if self.NT is None: - self.NT = [int(round(x*math.exp(y*self.num_gens))) for x,y in zip(self.init_size, self.r)] + self.NT = [x*math.exp(y*self.num_gens) for x,y in zip(self.init_size, self.r)] elif None in self.NT or \ any([isinstance(x, float) for x in self.NT]): raise ValueError('Relative ending population size is not allowed' @@ -665,7 +665,7 @@ for (n0, nt) in zip(self.init_size, self.NT)]) else: # with r ... - return self._save_size(pop.dvars().gen, [min(nt, int(n0 * math.exp(r * (self._gen + 1)))) + return self._save_size(pop.dvars().gen, [min(int(round(nt)), int(round(n0 * math.exp(r * (self._gen + 1))))) for (n0, nt, r) in zip(self.init_size, self.NT, self.r)]) @@ -1301,6 +1301,8 @@ of subpopulation index or names. ``capacity`` can be empty (no limit on carrying capacity), or one or more numbers for each of the subpopulations. ''' + self._N0 = None + self._T0 = None self.rates = rates self.slopes = slopes if self.rates and self.slopes: @@ -1332,17 +1334,48 @@ sz = pop.dvars()._expected_size else: sz = list(pop.subPopSizes()) + # + if self._N0 is None or len(self._N0) != len(subPops): + self._N0 = [sz[x] for x in subPops] + self._T0 = pop.dvars().gen for idx, sp in enumerate(subPops): if self.rates: if type(self.rates) in [list, tuple]: - sz[sp] = int(round(sz[sp] * (1 + self.rates[idx]))) + # if current size match, use it to do next + if sz[sp] == int(round(self._N0[idx] * math.exp((pop.dvars().gen - self._T0) * self.rates[idx]))): + sz[sp] = int(round(self._N0[idx] * math.exp((pop.dvars().gen - self._T0 + 1) * self.rates[idx]))) + # otherwise reset ... + else: + self._N0 = [sz[x] for x in subPops] + self._T0 = pop.dvars().gen + sz[sp] = int(round(sz[sp] * math.exp(self.rates[idx]))) else: - sz[sp] = int(round(sz[sp] * (1 + self.rates))) + if sz[sp] == int(round(self._N0[idx] * math.exp((pop.dvars().gen - self._T0) * self.rates))): + sz[sp] = int(round(self._N0[idx] * math.exp((pop.dvars().gen - self._T0 + 1) * self.rates))) + # otherwise reset ... + else: + self._N0 = [sz[x] for x in subPops] + self._T0 = pop.dvars().gen + sz[sp] = int(round(sz[sp] * math.exp(self.rates))) else: if type(self.slopes) in [list, tuple]: - sz[sp] = int(round(sz[sp] + self.slopes[idx])) + # if current size match, use it to do next + if sz[sp] == int(round(self._N0[idx] + (pop.dvars().gen - self._T0) * self.slopes[idx])): + sz[sp] = int(round(self._N0[idx] + (pop.dvars().gen - self._T0 + 1) * self.slopes[idx])) + # otherwise reset ... + else: + self._N0 = [sz[x] for x in subPops] + self._T0 = pop.dvars().gen + sz[sp] = int(round(sz[sp] + self.slopes[idx])) else: - sz[sp] = int(round(sz[sp] + self.slopes)) + # if current size match, use it to do next + if sz[sp] == int(round(self._N0[idx] + (pop.dvars().gen - self._T0) * self.slopes)): + sz[sp] = int(round(self._N0[idx] + (pop.dvars().gen - self._T0 + 1) * self.slopes)) + # otherwise reset ... + else: + self._N0 = [sz[x] for x in subPops] + self._T0 = pop.dvars().gen + sz[sp] = int(round(sz[sp] + self.slopes)) if self.capacity: for idx, sp in enumerate(subPops): if isinstance(self.capacity, (list, tuple)): Modified: trunk/test/test_19_demography.py =================================================================== --- trunk/test/test_19_demography.py 2014-10-15 03:30:30 UTC (rev 4949) +++ trunk/test/test_19_demography.py 2014-10-15 05:19:11 UTC (rev 4950) @@ -262,13 +262,13 @@ { 0: 202, 2: 206, - 99: 541, + 99: 544, }, initSize=500) ExponentialGrowthModel(N0=200, r=0.01, NT=400)._assertSize( { 0: 202, 2: 206, - 99: 541, + 99: 544, }, initSize=500) # ExponentialGrowthModel(T=100, N0=[200, 200], NT=[400, 800])._assertSize( @@ -281,33 +281,33 @@ { 0: [202, 202], 2: [206, 206], - 99: [541, 541], + 99: [544, 544], }, initSize=500) # unknown NT ExponentialGrowthModel(T=100, N0=[200, 200], r=[0.01, 0.03])._assertSize( { 0: [202, 206], - 2: [206, 218], - 99: [543, 4017], + 2: [206, 219], + 99: [544, 4017], }, initSize=500) # unknown T ExponentialGrowthModel(N0=[200, 200], r=0.01, NT=[400, 800])._assertSize( { 2: [206, 206], - 99: [400, 543], + 99: [400, 544], 100: [400, 549], }, initSize=500) # N0 size ... ExponentialGrowthModel(N0=[1., 1.], r=0.01, NT=[400, 800])._assertSize( { 2: [206, 206], - 99: [400, 543], + 99: [400, 544], 100: [400, 549], }, initSize=200) ExponentialGrowthModel(N0=[(None, 'name'), 1.], r=0.01, NT=[400, 800])._assertSize( { 2: [206, 206], - 99: [400, 543], + 99: [400, 544], 100: [400, 549], }, initSize=200) self.assertRaises(ExponentialGrowthModel, N0=[1., 1.], r=0.01, NT=[2, 3.]) @@ -656,7 +656,7 @@ { 0: 202, 2: 206, - 99: 532, + 99: 544, }, initSize=500) EventBasedModel(T=100, N0=(200, 'AF'), events=ExpansionEvent(rates=0.01, subPops='AF') @@ -664,7 +664,7 @@ { 0: 202, 2: 206, - 99: 532, + 99: 544, }, initSize=500) EventBasedModel(T=100, N0=[200, 200], events=ExpansionEvent(rates=0.01) @@ -672,7 +672,7 @@ { 0: [202, 202], 2: [206, 206], - 99: [532, 532], + 99: [544, 544], }, initSize=500) EventBasedModel(T=100, N0=[(200, 'A'), (200, 'B')], events=ExpansionEvent(rates=0.01, subPops=('A', 'B')) @@ -680,7 +680,7 @@ { 0: [202, 202], 2: [206, 206], - 99: [532, 532], + 99: [544, 544], }, initSize=500) # unknown NT EventBasedModel(T=100, N0=[200, 200], @@ -688,8 +688,8 @@ )._assertSize( { 0: [202, 206], - 2: [206, 218], - 99: [532, 3858], + 2: [206, 219], + 99: [544, 4017], }, initSize=500) EventBasedModel(T=100, N0=[(200, 'AF'), (200, 'EU'), (200, 'AS')], events=ExpansionEvent(rates=[0.01, 0.02], subPops=['EU', 'AS']) @@ -697,7 +697,7 @@ { 0: [200, 202, 204], 2: [200, 206, 212], - 99: [200, 532, 1445], + 99: [200, 544, 1478], }, initSize=500) EventBasedModel(T=100, N0=200, events=ExpansionEvent(slopes=2) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bpe...@us...> - 2014-10-15 03:30:40
|
Revision: 4949 http://sourceforge.net/p/simupop/code/4949 Author: bpeng2000 Date: 2014-10-15 03:30:30 +0000 (Wed, 15 Oct 2014) Log Message: ----------- Update demographic model module Modified Paths: -------------- trunk/src/demography.py Modified: trunk/src/demography.py =================================================================== --- trunk/src/demography.py 2014-10-15 01:46:39 UTC (rev 4948) +++ trunk/src/demography.py 2014-10-15 03:30:30 UTC (rev 4949) @@ -1412,7 +1412,7 @@ raise ValueError('Subpopulation index {} out of range'.format(self.toSubPop)) toSubPop = self.toSubPop else: - if self.toSubPop != pop.subPopNames(): + if self.toSubPop not in pop.subPopNames(): raise ValueError('No subpopulation with name {} can be located'.format(self.toSubPop)) toSubPop = pop.subPopNames().index(self.toSubPop) # @@ -1428,7 +1428,7 @@ num_migrants = [] for sp in range(pop.numSubPop()): if sp in subPops: - num_migrants.append(min(pop.subPopSize(sp), self.numbers[subPops.index(sp)])) + num_migrants.append(self.numbers[subPops.index(sp)]) else: # not involved num_migrants.append(0) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bpe...@us...> - 2014-10-15 01:46:44
|
Revision: 4948 http://sourceforge.net/p/simupop/code/4948 Author: bpeng2000 Date: 2014-10-15 01:46:39 +0000 (Wed, 15 Oct 2014) Log Message: ----------- Add more tests to demographic models Modified Paths: -------------- trunk/src/demography.py trunk/test/test_19_demography.py Modified: trunk/src/demography.py =================================================================== --- trunk/src/demography.py 2014-10-14 21:37:37 UTC (rev 4947) +++ trunk/src/demography.py 2014-10-15 01:46:39 UTC (rev 4948) @@ -1303,6 +1303,8 @@ ''' self.rates = rates self.slopes = slopes + if self.rates and self.slopes: + raise ValueError('Please specify one and only one of parameter rates and slopes.') self.capacity = capacity DemographicEvent.__init__(self, ops, output, begin, end, step, at, reps, subPops, infoFields) Modified: trunk/test/test_19_demography.py =================================================================== --- trunk/test/test_19_demography.py 2014-10-14 21:37:37 UTC (rev 4947) +++ trunk/test/test_19_demography.py 2014-10-15 01:46:39 UTC (rev 4948) @@ -534,7 +534,44 @@ 0: 200, 99: 200, }) + EventBasedModel(T=100, N0=(200, 200, 200), + events=[ + ResizeEvent(sizes=300, subPops=1, names='AF', at=2), + MergeEvent(subPops=[1,2], at=3), + ] + )._assertSize( + { + 0: (200, 200, 200), + 2: (200, 300, 200), + 99: (200, 500), + }, initSize=500) + EventBasedModel(T=100, N0=(200, 200, 200), + events=[ + ResizeEvent(sizes=300, subPops=1, at=2), + ResizeEvent(sizes=400, subPops=2, at=2), + ResizeEvent(sizes=350, subPops=1, at=2), + ] + )._assertSize( + { + 0: (200, 200, 200), + 2: (200, 350, 400), + 99: (200, 350, 400), + }, initSize=500) + EventBasedModel(T=100, N0=(200, 200, 200), + events=[ + ResizeEvent(sizes=300, subPops=1, at=2), + MergeEvent(subPops=[1,0], at=2), + ResizeEvent(sizes=350, subPops=1, at=2), + ] + )._assertSize( + { + 0: (200, 200, 200), + 2: (500, 350), + 99: (500, 350), + }, initSize=500) + + def testResizeEvent(self): EventBasedModel(T=100, N0=(200, 200, 200), events=ResizeEvent(sizes=300, subPops=1, names='AF', at=2) @@ -686,7 +723,36 @@ 2: [206, 218], 99: [400, 800], }, initSize=500) + self.assertRaises(ValueError, ExpansionEvent, rates=[0.2, 0.4], slopes=[2, 6]) + def testSplitEvent(self): + EventBasedModel(T=100, N0=200, + events=SplitEvent(sizes=(0.4, 0.6), at=2) + )._assertSize( + { + 0: [200], + 2: [80, 120], + 99: [80, 120], + }, initSize=500) + EventBasedModel(T=100, N0=200, + events=SplitEvent(sizes=(0.4, 250), at=2) + )._assertSize( + { + 0: [200], + 2: [80, 250], + 99: [80, 250], + }, initSize=500) + self.assertRaises(ValueError, EventBasedModel(T=100, N0=200, + events=SplitEvent(sizes=(0.4, 250), begin=2) + )._assertSize, { }, initSize=500) + EventBasedModel(T=100, N0=((100, 'A'), (200, 'B')), + events=SplitEvent(sizes=(0.4, 250), at=2, subPops='B') + )._assertSize( + { + 0: [100, 200], + 2: [100, 80, 250], + 99: [100, 80, 250], + }, initSize=500) def testAdmixtureEvent(self): EventBasedModel(T=10, This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bpe...@us...> - 2014-10-14 21:37:40
|
Revision: 4947 http://sourceforge.net/p/simupop/code/4947 Author: bpeng2000 Date: 2014-10-14 21:37:37 +0000 (Tue, 14 Oct 2014) Log Message: ----------- More fix Modified Paths: -------------- trunk/src/demography.py Modified: trunk/src/demography.py =================================================================== --- trunk/src/demography.py 2014-10-14 21:34:25 UTC (rev 4946) +++ trunk/src/demography.py 2014-10-14 21:37:37 UTC (rev 4947) @@ -1761,6 +1761,8 @@ raise ValueError('Length of evolution T0=%d should be more than T_AF=%d' % (T0, T_AF)) # scale = float(scale) + if isinstance(outcome, str): + outcome = [outcome] EventBasedModel.__init__(self, T = int(T0/scale), N0 = (int(N_A/scale), 'Ancestral'), This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bpe...@us...> - 2014-10-14 21:34:30
|
Revision: 4946 http://sourceforge.net/p/simupop/code/4946 Author: bpeng2000 Date: 2014-10-14 21:34:25 +0000 (Tue, 14 Oct 2014) Log Message: ----------- More tests and more bug fixes of the event based implementation of models Modified Paths: -------------- trunk/doc/userGuide.py trunk/src/demography.py trunk/test/test_19_demography.py Modified: trunk/doc/userGuide.py =================================================================== --- trunk/doc/userGuide.py 2014-10-14 16:55:00 UTC (rev 4945) +++ trunk/doc/userGuide.py 2014-10-14 21:34:25 UTC (rev 4946) @@ -2263,55 +2263,6 @@ #end_file -#begin_file log/ReuseBurnIn.py -#begin_ignore -import simuOpt -simuOpt.setOptions(quiet=True) -#end_ignore -import os -import simuPOP as sim -#begin_ignore -if os.path.isfile('burnin.pop'): - os.remove('burnin.pop') - -sim.setRNG(seed=1234) -#end_ignore -from simuPOP.demography import LinearGrowthModel - -for rep in range(3): - pop = sim.Population(1000, loci=1) - evolved = pop.evolve( - initOps=sim.InitSex(), - preOps=[ - sim.RevertIf(os.path.isfile('burnin.pop'), 'burnin.pop', at=0), - sim.IfElse(not os.path.isfile('burnin.pop'), - [sim.PyOutput('Save burnin at gen=10\n'), - sim.SavePopulation('burnin.pop')], at=10), - ], - matingScheme=sim.RandomMating(), - gen=20 - ) - print('Evolved {} generations'.format(evolved)) - -for rep in range(3): - pop = sim.Population(1000, loci=1) - evolved = pop.evolve( - initOps=sim.InitSex(), - preOps=[ - sim.RevertIf(os.path.isfile('burnin.pop'), 'burnin.pop', at=0), - sim.IfElse(not os.path.isfile('burnin.pop'), - [sim.PyOutput('Save burnin at gen=10\n'), - sim.SavePopulation('burnin.pop')], at=10), - ], - matingScheme=sim.RandomMating( - subPopSize=LinearGrowthModel(T=20, r=0.1)), - ) - print('Evolved {} generations'.format(evolved)) - - - -#end_file - #begin_file log/DiscardIf.py #begin_ignore import simuOpt Modified: trunk/src/demography.py =================================================================== --- trunk/src/demography.py 2014-10-14 16:55:00 UTC (rev 4945) +++ trunk/src/demography.py 2014-10-14 21:34:25 UTC (rev 4946) @@ -489,7 +489,7 @@ raise ValueError('Generation number %d out of bound (0<= t < %d is expected' % (x, T)) else: - return int(math.exp(((x+1-T0)*math.log(NT) + (T-x-1)*math.log(N0))/(T-T0))) + return int(round(math.exp(((x+1-T0)*math.log(NT) + (T-x-1)*math.log(N0))/(T-T0)))) def _linearIntepolate(self, N0, NT, T, x, T0=0): '''x=0, ... T-1 @@ -500,7 +500,7 @@ raise ValueError('Generation number %d out of bound (0<= t < %d is expected)' % (x, T)) else: - return int(((x+1-T0)*NT + (T-x-1)*N0)/(T-T0)) + return int(round(((x+1-T0)*NT + (T-x-1)*N0)/(T-T0))) def _setup(self, pop): return True @@ -623,14 +623,14 @@ # what is the number of generations if self.r == 0: raise ValueError('Cannot reach destination size with r=0') - T = [int((math.log(y) - math.log(x)) / self.r) + 1 for (x,y) in zip(self.init_size, self.NT)] + T = [int(round((math.log(y) - math.log(x)) / self.r)) + 1 for (x,y) in zip(self.init_size, self.NT)] if max(T) < 0: raise ValueError('Cannot reach destination size in this configuraiton.') self.num_gens = max(T + [1]) # elif self.NT is None: # get final population size from T and r - self.NT = [int(x*((1.+self.r)**self.num_gens)) for x in self.init_size] + self.NT = [int(round(x*((1.+self.r)**self.num_gens))) for x in self.init_size] elif None in self.NT or \ any([isinstance(x, float) for x in self.NT]): raise ValueError('Relative ending population size is not allowed' @@ -641,7 +641,7 @@ raise ValueError('Please specify growth rate for each subpopulation ' 'if multiple growth rates are specified.') if self.NT is None: - self.NT = [int(x*math.exp(y*self.num_gens)) for x,y in zip(self.init_size, self.r)] + self.NT = [int(round(x*math.exp(y*self.num_gens))) for x,y in zip(self.init_size, self.r)] elif None in self.NT or \ any([isinstance(x, float) for x in self.NT]): raise ValueError('Relative ending population size is not allowed' @@ -1069,7 +1069,7 @@ raise ValueError('Cannot apply to a population without variables _gen or _num_gens') # gen = pop.dvars()._gen - end = pop.dvars()._num_gens + end = pop.dvars()._num_gens - 1 rep = pop.dvars().rep # if self.reps != ALL_AVAIL and (\ @@ -1111,7 +1111,7 @@ return False def _identifySubPops(self, pop): - if self.subPops == ALL_AVAIL: + if self.subPops is ALL_AVAIL: return range(pop.numSubPop()) # ret = [] @@ -1171,18 +1171,24 @@ if isinstance(self.sizes, (tuple, list)) and len(subPops) != len(self.sizes): raise ValueError('If a list of sizes are specified, please specify them for all subpopulations.') if self.names and len(self.names) != len(subPops): - raise ValueError('If a list of names are specified, please assign to all specified subpopulations.') + raise ValueError('Legnth of names (%d) and subPops mismatch (%d)' % (len(self.names), len(subPops))) + # + if '_expected_size' in pop.vars(): + # if there was an intended size, realized it. + pop.resize(pop.vars().pop('_expected_size')) # new size ... sz = list(pop.subPopSizes()) for idx, sp in enumerate(subPops): if isinstance(self.sizes, int): sz[sp] = self.sizes - elif isinstance(sizes, float): + elif isinstance(self.sizes, float): sz[sp] = int(round(sz[sp] * self.sizes)) - elif isinstance(sp, (list, tuple)) and isinstance(self.sizes[idx], int): + elif isinstance(self.sizes, (list, tuple)) and isinstance(self.sizes[idx], int): sz[sp] = self.sizes[idx] - elif isinstance(sp, (list, tuple)) and isinstance(self.sizes[idx], float): + elif isinstance(self.sizes, (list, tuple)) and isinstance(self.sizes[idx], float): sz[sp] = int(round(self.sizes[idx] * sz[sp])) + else: + raise ValueError('Unrecognized sizes parameter.') # pop.resize(sz, propagate=True) if self.removeEmptySubPops: @@ -1201,12 +1207,27 @@ '''A demographic event that merges subpopulations into a single subpopulation. The merged subpopulation will have the name of the first merged subpopulation unless a separate ``name`` is supported.''' - DemographicEvent.__init__(self, - ops=[MergeSubPops(subPops=subPops, name=name)] + ops, - output=output, begin=begin, end=end, step=step, at=at, - reps=reps, subPops=subPops, infoFields=infoFields) + self.name = name + DemographicEvent.__init__(self, ops, output, begin, end, step, at, reps, + subPops, infoFields) + def apply(self, pop): + if not DemographicEvent.apply(self, pop): + return False + # + if not self._applicable(pop): + return True + # + if '_expected_size' in pop.vars(): + # if there was an intended size, realized it. + pop.resize(pop.vars().pop('_expected_size')) + # identify applicable subpopulations + subPops = self._identifySubPops(pop) + # + pop.mergeSubPops(subPops=subPops, name=self.name) + return True + class SplitEvent(DemographicEvent): '''A demographic event that splits a specified population into two or more subpopulations.''' @@ -1236,6 +1257,10 @@ if len(subPops) != 1: raise ValueError('Please specify one and only one subpopulation for event CopyEvent.') subPop = subPops[0] + # + if '_expected_size' in pop.vars(): + # if there was an intended size, realized it. + pop.resize(pop.vars().pop('_expected_size')) # new size ... sz = [] for x in self.sizes: @@ -1301,24 +1326,27 @@ raise ValueError('If specified, please specify carrying capacity for all ' 'subpopulations or each of the {} subpopulations'.format(len(subPops))) # - sz = list(pop.subPopSizes()) + if '_expected_size' in pop.vars(): + sz = pop.dvars()._expected_size + else: + sz = list(pop.subPopSizes()) for idx, sp in enumerate(subPops): if self.rates: if type(self.rates) in [list, tuple]: - sz[idx] = int(round(sz[idx] * (1 + self.rates[idx]))) + sz[sp] = int(round(sz[sp] * (1 + self.rates[idx]))) else: - sz[idx] = int(round(sz[idx] * (1 + self.rates))) + sz[sp] = int(round(sz[sp] * (1 + self.rates))) else: if type(self.slopes) in [list, tuple]: - sz[idx] = int(round(sz[idx] + self.slopes[idx])) + sz[sp] = int(round(sz[sp] + self.slopes[idx])) else: - sz[idx] = int(round(sz[idx] + self.slopes)) + sz[sp] = int(round(sz[sp] + self.slopes)) if self.capacity: for idx, sp in enumerate(subPops): if isinstance(self.capacity, (list, tuple)): - sz[idx] = min(sz[idx], self.capacity[idx]) + sz[sp] = min(sz[sp], self.capacity[idx]) else: - sz[idx] = min(sz[idx], self.capacity) + sz[sp] = min(sz[sp], self.capacity) pop.dvars()._expected_size = sz return True @@ -1369,6 +1397,10 @@ # identify applicable subpopulations subPops = self._identifySubPops(pop) # + if '_expected_size' in pop.vars(): + # if there was an intended size, realized it. + pop.resize(pop.vars().pop('_expected_size')) + # if self.toSubPop is None: toSubPop = self.toSubPop else: @@ -1465,7 +1497,7 @@ return True -class OutOfAfricaModel(EventBasedModel): +class _OutOfAfricaModel_event(EventBasedModel): '''A dempgraphic model for the CHB, CEU, and YRI populations, as defined in Gutenkunst 2009, Plos Genetics. The model is depicted in Figure 2, and the default parameters are listed in Table 1 of this paper. ''' @@ -1529,7 +1561,7 @@ sizes=int(N_AF/scale), names='AF'), # split B from AF SplitEvent(at=-int(T_B/scale), subPops='AF', - sizes=(1.0, int(N_AF/scale)), names=('AF', 'B')), + sizes=(1.0, int(N_B/scale)), names=('AF', 'B')), # migration between AF and B DemographicEvent( ops=Migrator(rate=[ @@ -1564,44 +1596,113 @@ ] ) - # # for python 2.x and 3.x compatibility - # scale = float(scale) - # MultiStageModel.__init__(self, [ - # InstantChangeModel( - # T=int((T0-T_B)/scale), - # N0=(int(N_A/scale), 'Ancestral'), - # # change population size twice, one at T_AF, one at T_B - # G=[int((T0-T_AF)/scale)], - # NG=[(int(N_AF/scale), 'AF')] - # ), - # # - # # at T_B, split to population B from subpopulation 1 - # InstantChangeModel( - # T=int((T_B - T_EU_AS)/scale), - # # change population size twice, one at T_AF, one at T_B - # N0=[None, (int(N_B/scale), 'B')], - # ops=Migrator(rate=[ - # [m_AF_B, 0], - # [0, m_AF_B]]) - # ), - # ExponentialGrowthModel( - # T=int(T_EU_AS/scale), - # N0=[None, - # # split B into EU and AS at the beginning of this - # # exponential growth stage - # [(int(N_EU0/scale), 'EU'), (int(N_AS0/scale), 'AS')]], - # r=[0, r_EU*scale, r_AS*scale], - # infoFields='migrate_to', - # ops=Migrator(rate=[ - # [0, m_AF_EU, m_AF_AS], - # [m_AF_EU, 0, m_EU_AS], - # [m_AF_AS, m_EU_AS, 0] - # ]) - # ), - # ] + finalStage, ops=ops, infoFields=infoFields - # ) -class SettlementOfNewWorldModel(EventBasedModel): +class _OutOfAfricaModel_outcome(MultiStageModel): + '''A dempgraphic model for the CHB, CEU, and YRI populations, as defined in + Gutenkunst 2009, Plos Genetics. The model is depicted in Figure 2, and the + default parameters are listed in Table 1 of this paper. ''' + def __init__(self, + T0, + N_A=7300, + N_AF=12300, + N_B=2100, + N_EU0=1000, + r_EU=0.004, + N_AS0=510, + r_AS=0.0055, + m_AF_B=0.00025, + m_AF_EU=0.00003, + m_AF_AS=0.000019, + m_EU_AS=0.000096, + T_AF=220000//25, + T_B=140000//25, + T_EU_AS=21200//25, + ops=[], + infoFields=[], + outcome=['AF', 'EU', 'AS'], + scale=1 + ): + '''Counting **backward in time**, this model evolves a population for ``T0`` + generations (required parameter). The ancient population ``A`` started at + size ``N_A`` and expanded at ``T_AF`` generations from now, to pop ``AF`` + with size ``N_AF``. Pop ``B`` split from pop ``AF`` at ``T_B`` generations + from now, with size ``N_B``; Pop ``AF`` remains as ``N_AF`` individuals. + Pop ``EU`` and ``AS`` split from pop ``B`` at ``T_EU_AS`` generations + from now; with size ``N_EU0`` individuals and ``N_ASO`` individuals, + respectively. Pop ``EU`` grew exponentially with rate ``r_EU``; Pop + ``AS`` grew exponentially with rate ``r_AS``. The ``YRI``, ``CEU`` and + ``CHB`` samples are drawn from ``AF``, ``EU`` and ``AS`` populations + respectively. Additional operators could be added to ``ops``. Information + fields required by these operators should be passed to ``infoFields``. If + a scaling factor ``scale`` is specified, all population sizes and + generation numbers will be divided by a factor of ``scale``. This demographic + model by default returns all populations (``AF``, ``EU``, ``AS``) but + you can choose to keep only selected subpopulations using parameter + ``outcome`` (e.g. ``outcome=['EU', 'AS']``). + + This model merges all subpopulations if it is applied to an initial + population with multiple subpopulation. + ''' + # + if T0 < T_AF: + raise ValueError('Length of evolution T0=%d should be more than T_AF=%d' % (T0, T_AF)) + # + if isinstance(outcome, str): + outcome = [outcome] + # + final_subpops = [None, None, None] + for (idx, name) in enumerate(['AF', 'EU', 'AS']): + if name not in outcome: + final_subpops[idx] = 0 + # + if 0 in final_subpops: + finalStage = [ + InstantChangeModel(T=1, + N0=final_subpops, + removeEmptySubPops=True) + ] + else: + finalStage = [] + # for python 2.x and 3.x compatibility + scale = float(scale) + MultiStageModel.__init__(self, [ + InstantChangeModel( + T=int((T0-T_B)/scale), + N0=(int(N_A/scale), 'Ancestral'), + # change population size twice, one at T_AF, one at T_B + G=[int((T0-T_AF)/scale)], + NG=[(int(N_AF/scale), 'AF')] + ), + # + # at T_B, split to population B from subpopulation 1 + InstantChangeModel( + T=int((T_B - T_EU_AS)/scale), + # change population size twice, one at T_AF, one at T_B + N0=[None, (int(N_B/scale), 'B')], + ops=Migrator(rate=[ + [m_AF_B, 0], + [0, m_AF_B]]) + ), + ExponentialGrowthModel( + T=int(T_EU_AS/scale), + N0=[None, + # split B into EU and AS at the beginning of this + # exponential growth stage + [(int(N_EU0/scale), 'EU'), (int(N_AS0/scale), 'AS')]], + r=[0, r_EU*scale, r_AS*scale], + infoFields='migrate_to', + ops=Migrator(rate=[ + [0, m_AF_EU, m_AF_AS], + [m_AF_EU, 0, m_EU_AS], + [m_AF_AS, m_EU_AS, 0] + ]) + ), + ] + finalStage, ops=ops, infoFields=infoFields + ) + +OutOfAfricaModel = _OutOfAfricaModel_outcome + +class _SettlementOfNewWorldModel_event(EventBasedModel): '''A dempgraphic model for settlement of the new world of Americans, as defined in Gutenkunst 2009, Plos Genetics. The model is depicted in Figure 3, and the default parameters are listed in Table 2 of this paper. ''' @@ -1670,7 +1771,7 @@ sizes=int(N_AF/scale), names='AF'), # split B from AF SplitEvent(at=-int(T_B/scale), subPops='AF', - sizes=(1.0, int(N_AF/scale)), names=('AF', 'B')), + sizes=(1.0, int(N_B/scale)), names=('AF', 'B')), # migration between AF and B DemographicEvent( ops=Migrator(rate=[ @@ -1723,102 +1824,157 @@ ] ) - ## - ## The following is the old outcome-based implementation. The new version - ## uses event-based model, which is easier to implement and change. - ## - ## # try to figure out how to mix two populations - ## N_EU=int(N_EU0*math.exp(r_EU*T_EU_AS)) - ## N_MX=int(N_MX0*math.exp(r_MX*T_MX)) - ## # - ## # - ## # for python 2.x and 3.x compatibility - ## if isinstance(outcome, str): - ## outcome = [outcome] - ## if 'MXL' in outcome: - ## # with admixture - ## final_subpops = [None, None, None, None, None] - ## for (idx, name) in enumerate(['AF', 'EU', 'AS', 'MX', 'MXL']): - ## if name not in outcome: - ## final_subpops[idx] = 0 - ## # - ## admixtureStage = [ - ## AdmixtureModel(T=1, - ## N0=[None, None, None, None], - ## # mixing European and Mexican population - ## model=['HI', 1, 3, 1-f_MX, 'MXL']), - ## InstantChangeModel(T=1, - ## N0=final_subpops, - ## removeEmptySubPops=True) - ## ] - ## else: - ## final_subpops = [None, None, None, None] - ## for (idx, name) in enumerate(['AF', 'EU', 'AS', 'MX']): - ## if name not in outcome: - ## final_subpops[idx] = 0 - ## admixtureStage = [ - ## InstantChangeModel(T=1, - ## N0=final_subpops, - ## removeEmptySubPops=True) - ## ] - ## # - ## scale = float(scale) - ## MultiStageModel.__init__(self, [ - ## InstantChangeModel( - ## T=int((T0-T_B)/scale), - ## N0=(int(N_A/scale), 'Ancestral'), - ## # change population size twice, one at T_AF, one at T_B - ## G=[int((T0-T_AF)/scale)], - ## NG=[(int(N_AF/scale), 'AF')] - ## ), - ## # - ## # at T_B, split to population B from subpopulation 1 - ## InstantChangeModel( - ## T=int((T_B - T_EU_AS)/scale), - ## # change population size twice, one at T_AF, one at T_B - ## N0=[None, (int(N_B/scale), 'B')], - ## ops=Migrator(rate=[ - ## [m_AF_B, 0], - ## [0, m_AF_B]]) - ## ), - ## ExponentialGrowthModel( - ## T=int((T_EU_AS - T_MX)/scale), - ## N0=[None, - ## # split B into EU and AS at the beginning of this - ## # exponential growth stage - ## [(int(N_EU0/scale), 'EU'), (int(N_AS0/scale), 'AS')]], - ## r=[0, r_EU*scale, r_AS*scale], - ## infoFields='migrate_to', - ## ops=Migrator(rate=[ - ## [0, m_AF_EU, m_AF_AS], - ## [m_AF_EU, 0, m_EU_AS], - ## [m_AF_AS, m_EU_AS, 0] - ## ]) - ## ), - ## ExponentialGrowthModel(T=int(T_MX/scale), - ## N0=[None, - ## # initial population size has to be calculated - ## # because we only know the final population size of - ## # EU and AS - ## None, - ## # split MX from AS - ## [(None, 'AS'), (int(N_MX0//scale), 'MX')]], - ## r=[0, r_EU*scale, r_AS*scale, r_MX*scale], - ## infoFields='migrate_to', - ## ops=Migrator(rate=[ - ## [0, m_AF_EU, m_AF_AS], - ## [m_AF_EU, 0, m_EU_AS], - ## [m_AF_AS, m_EU_AS, 0] - ## ], - ## # the last MX population does not involve in - ## # migration - ## subPops=[0, 1, 2], - ## toSubPops=[0, 1, 2]) - ## ) - ## ] + admixtureStage, - ## ops=ops, infoFields=infoFields - ## ) +class _SettlementOfNewWorldModel_outcome(MultiStageModel): + '''A dempgraphic model for settlement of the new world of Americans, as defined + in Gutenkunst 2009, Plos Genetics. The model is depicted in Figure 3, and the + default parameters are listed in Table 2 of this paper. ''' + def __init__(self, + T0, + N_A=7300, + N_AF=12300, + N_B=2100, + N_EU0=1500, + r_EU=0.0023, + N_AS0=590, + r_AS=0.0037, + N_MX0=800, + r_MX=0.005, + m_AF_B=0.00025, + m_AF_EU=0.00003, + m_AF_AS=0.000019, + m_EU_AS=0.0000135, + T_AF=220000//25, + T_B=140000//25, + T_EU_AS=26400//25, + T_MX=21600//25, + f_MX=0.48, + ops=[], + infoFields=[], + outcome='MXL', + scale=1 + ): + '''Counting **backward in time**, this model evolves a population for ``T0`` + generations. The ancient population ``A`` started at size ``N_A`` and + expanded at ``T_AF`` generations from now, to pop ``AF`` with size ``N_AF``. + Pop ``B`` split from pop ``AF`` at ``T_B`` generations from now, with + size ``N_B``; Pop ``AF`` remains as ``N_AF`` individuals. Pop ``EU`` and + ``AS`` split from pop ``B`` at ``T_EU_AS`` generations from now; with + size ``N_EU0`` individuals and ``N_ASO`` individuals, respectively. Pop + ``EU`` grew exponentially with final population size ``N_EU``; Pop + ``AS`` grew exponentially with final populaiton size ``N_AS``. Pop ``MX`` + split from pop ``AS`` at ``T_MX`` generations from now with size ``N_MX0``, + grew exponentially to final size ``N_MX``. Migrations are allowed between + populations with migration rates ``m_AF_B``, ``m_EU_AS``, ``m_AF_EU``, + and ``m_AF_AS``. At the end of the evolution, the ``AF`` and ``CHB`` + populations are removed, and the ``EU`` and ``MX`` populations are merged + with ``f_MX`` proportion for ``MX``. The Mexican American<F19> sample could + be sampled from the last single population. Additional operators could + be added to ``ops``. Information fields required by these operators + should be passed to ``infoFields``. If a scaling factor ``scale`` + is specified, all population sizes and generation numbers will be divided by + a factor of ``scale``. This demographic model by default only returns the + mixed Mexican America model (``outputcom='MXL'``) but you can specify any + combination of ``AF``, ``EU``, ``AS``, ``MX`` and ``MXL``. + This model merges all subpopulations if it is applied to an initial population + with multiple subpopulation. + ''' + if T0 < T_AF: + raise ValueError('Length of evolution T0=%d should be more than T_AF=%d' % (T0, T_AF)) + # try to figure out how to mix two populations + N_EU=int(N_EU0*math.exp(r_EU*T_EU_AS)) + N_MX=int(N_MX0*math.exp(r_MX*T_MX)) + # + # for python 2.x and 3.x compatibility + if isinstance(outcome, str): + outcome = [outcome] + if 'MXL' in outcome: + # with admixture + final_subpops = [None, None, None, None, None] + for (idx, name) in enumerate(['AF', 'EU', 'AS', 'MX', 'MXL']): + if name not in outcome: + final_subpops[idx] = 0 + # + admixtureStage = [ + AdmixtureModel(T=1, + N0=[None, None, None, None], + # mixing European and Mexican population + model=['HI', 1, 3, 1-f_MX, 'MXL']), + InstantChangeModel(T=1, + N0=final_subpops, + removeEmptySubPops=True) + ] + else: + final_subpops = [None, None, None, None] + for (idx, name) in enumerate(['AF', 'EU', 'AS', 'MX']): + if name not in outcome: + final_subpops[idx] = 0 + admixtureStage = [ + InstantChangeModel(T=1, + N0=final_subpops, + removeEmptySubPops=True) + ] + # + scale = float(scale) + MultiStageModel.__init__(self, [ + InstantChangeModel( + T=int((T0-T_B)/scale), + N0=(int(N_A/scale), 'Ancestral'), + # change population size twice, one at T_AF, one at T_B + G=[int((T0-T_AF)/scale)], + NG=[(int(N_AF/scale), 'AF')] + ), + # + # at T_B, split to population B from subpopulation 1 + InstantChangeModel( + T=int((T_B - T_EU_AS)/scale), + # change population size twice, one at T_AF, one at T_B + N0=[None, (int(N_B/scale), 'B')], + ops=Migrator(rate=[ + [m_AF_B, 0], + [0, m_AF_B]]) + ), + ExponentialGrowthModel( + T=int((T_EU_AS - T_MX)/scale), + N0=[None, + # split B into EU and AS at the beginning of this + # exponential growth stage + [(int(N_EU0/scale), 'EU'), (int(N_AS0/scale), 'AS')]], + r=[0, r_EU*scale, r_AS*scale], + infoFields='migrate_to', + ops=Migrator(rate=[ + [0, m_AF_EU, m_AF_AS], + [m_AF_EU, 0, m_EU_AS], + [m_AF_AS, m_EU_AS, 0] + ]) + ), + ExponentialGrowthModel(T=int(T_MX/scale), + N0=[None, + # initial population size has to be calculated + # because we only know the final population size of + # EU and AS + None, + # split MX from AS + [(None, 'AS'), (int(N_MX0//scale), 'MX')]], + r=[0, r_EU*scale, r_AS*scale, r_MX*scale], + infoFields='migrate_to', + ops=Migrator(rate=[ + [0, m_AF_EU, m_AF_AS], + [m_AF_EU, 0, m_EU_AS], + [m_AF_AS, m_EU_AS, 0] + ], + # the last MX population does not involve in + # migration + subPops=[0, 1, 2], + toSubPops=[0, 1, 2]) + ) + ] + admixtureStage, + ops=ops, infoFields=infoFields + ) + + +SettlementOfNewWorldModel = _SettlementOfNewWorldModel_outcome + class CosiModel(MultiStageModel): '''A dempgraphic model for Africa, Asia and Europe, as described in Schaffner et al, Genome Research, 2005, and implemented in the coalescent Modified: trunk/test/test_19_demography.py =================================================================== --- trunk/test/test_19_demography.py 2014-10-14 16:55:00 UTC (rev 4945) +++ trunk/test/test_19_demography.py 2014-10-14 21:34:25 UTC (rev 4946) @@ -255,33 +255,33 @@ ExponentialGrowthModel(T=100, N0=200, NT=400)._assertSize( { 0: 201, - 1: 202, + 1: 203, 99: 400, }, initSize=500) ExponentialGrowthModel(T=100, N0=200, r=0.01)._assertSize( { 0: 202, 2: 206, - 99: 540, + 99: 541, }, initSize=500) ExponentialGrowthModel(N0=200, r=0.01, NT=400)._assertSize( { 0: 202, 2: 206, - 99: 540, + 99: 541, }, initSize=500) # ExponentialGrowthModel(T=100, N0=[200, 200], NT=[400, 800])._assertSize( { - 0: [201, 202], - 1: [202, 205], + 0: [201, 203], + 1: [203, 206], 99: [400, 800], }, initSize=500) ExponentialGrowthModel(T=100, N0=[200, 200], r=0.01)._assertSize( { 0: [202, 202], 2: [206, 206], - 99: [540, 540], + 99: [541, 541], }, initSize=500) # unknown NT ExponentialGrowthModel(T=100, N0=[200, 200], r=[0.01, 0.03])._assertSize( @@ -535,6 +535,83 @@ 99: 200, }) + def testResizeEvent(self): + EventBasedModel(T=100, N0=(200, 200, 200), + events=ResizeEvent(sizes=300, subPops=1, names='AF', at=2) + )._assertSize( + { + 0: (200, 200, 200), + 2: (200, 300, 200), + 99: (200, 300, 200), + }, initSize=500) + EventBasedModel(T=100, N0=(200, 200, 200), + events=ResizeEvent(sizes=1.5, subPops=1, names='AF', at=2) + )._assertSize( + { + 0: (200, 200, 200), + 2: (200, 300, 200), + 99: (200, 300, 200), + }, initSize=500) + EventBasedModel(T=100, N0=(200, 200, 200), + events=ResizeEvent(sizes=(1.5, 1.5, 1.8), names=('AF', 'A', 'B'), at=2) + )._assertSize( + { + 0: (200, 200, 200), + 2: (300, 300, 360), + 99: (300, 300, 360), + }, initSize=500) + EventBasedModel(T=100, N0=(200, 200, 200), + events=ResizeEvent(sizes=(300, 1.5), subPops=[0, 1], at=2) + )._assertSize( + { + 0: (200, 200, 200), + 2: (300, 300, 200), + 99: (300, 300, 200), + }, initSize=500) + EventBasedModel(T=100, N0=(200, 200, 200), + events=ResizeEvent(sizes=(300, 0.9), subPops=[0, 1], begin=2) + )._assertSize( + { + 0: (200, 200, 200), + 2: (300, 180, 200), + 99: (300, 5, 200), + }, initSize=500) + EventBasedModel(T=100, N0=(200, 200, 200), + events=ResizeEvent(sizes=(300, 0), subPops=[0, 1], at=2, removeEmptySubPops=True) + )._assertSize( + { + 0: (200, 200, 200), + 2: (300, 200), + 99: (300, 200), + }, initSize=500) + + + def testMergeEvent(self): + EventBasedModel(T=100, N0=(200, 200, 200), + events=MergeEvent(name='AF', at=2) + )._assertSize( + { + 0: (200, 200, 200), + 2: (600,), + 99: (600,), + }, initSize=500) + EventBasedModel(T=100, N0=((100, 'AF'), (200, 'EU'), (300, 'AS')), + events=MergeEvent(name='AF', subPops=('AF', 'AS'), at=2) + )._assertSize( + { + 0: (100, 200, 300), + 2: (400, 200), + 99: (400, 200), + }, initSize=500) + EventBasedModel(T=100, N0=((100, 'AF'), (200, 'EU'), (300, 'AS')), + events=MergeEvent(name='EU', subPops=('EU', 'AS'), at=2) + )._assertSize( + { + 0: (100, 200, 300), + 2: (100, 500), + 99: (100, 500), + }, initSize=500) + def testExpansionEvent(self): EventBasedModel(T=100, N0=200, events=ExpansionEvent(rates=0.01) @@ -542,16 +619,32 @@ { 0: 202, 2: 206, - 99: 466, + 99: 532, }, initSize=500) + EventBasedModel(T=100, N0=(200, 'AF'), + events=ExpansionEvent(rates=0.01, subPops='AF') + )._assertSize( + { + 0: 202, + 2: 206, + 99: 532, + }, initSize=500) EventBasedModel(T=100, N0=[200, 200], events=ExpansionEvent(rates=0.01) )._assertSize( { 0: [202, 202], 2: [206, 206], - 99: [466, 466], + 99: [532, 532], }, initSize=500) + EventBasedModel(T=100, N0=[(200, 'A'), (200, 'B')], + events=ExpansionEvent(rates=0.01, subPops=('A', 'B')) + )._assertSize( + { + 0: [202, 202], + 2: [206, 206], + 99: [532, 532], + }, initSize=500) # unknown NT EventBasedModel(T=100, N0=[200, 200], events=ExpansionEvent(rates=[0.01, 0.03]) @@ -559,10 +652,16 @@ { 0: [202, 206], 2: [206, 218], - 99: [466, 3574], + 99: [532, 3858], }, initSize=500) - - def testExpansionEvent(self): + EventBasedModel(T=100, N0=[(200, 'AF'), (200, 'EU'), (200, 'AS')], + events=ExpansionEvent(rates=[0.01, 0.02], subPops=['EU', 'AS']) + )._assertSize( + { + 0: [200, 202, 204], + 2: [200, 206, 212], + 99: [200, 532, 1445], + }, initSize=500) EventBasedModel(T=100, N0=200, events=ExpansionEvent(slopes=2) )._assertSize( @@ -579,7 +678,6 @@ 2: [206, 206], 99: [400, 400], }, initSize=500) - # unknown NT EventBasedModel(T=100, N0=[200, 200], events=ExpansionEvent(slopes=[2, 6]) )._assertSize( @@ -589,6 +687,7 @@ 99: [400, 800], }, initSize=500) + def testAdmixtureEvent(self): EventBasedModel(T=10, events=AdmixtureEvent(subPops=[1, 2], sizes=[0.5, 0.5], at=2) @@ -676,19 +775,6 @@ 5: (200, 100, 400), }, initSize=[200, 300, 400]) - - def testAsEvent(self): - EventBasedModel(T=10, - events=MergeEvent(subPops=[0, 2], at=2) - )._assertSize( - { - 0: (200, 300, 400), - 2: (600, 300), - 4: (600, 300), - 5: (600, 300), - }, initSize=[200, 300, 400]) - - def testRevertAndDemo(self): 'Test the use of demographic model with RevertIf operator' pop = Population(100, loci=1) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <bpe...@us...> - 2014-10-14 16:55:05
|
Revision: 4945 http://sourceforge.net/p/simupop/code/4945 Author: bpeng2000 Date: 2014-10-14 16:55:00 +0000 (Tue, 14 Oct 2014) Log Message: ----------- Use cached mechanism to allow retriving past population sizes of demographic models Modified Paths: -------------- trunk/doc/userGuide.lyx trunk/src/demography.py trunk/test/test_19_demography.py Modified: trunk/doc/userGuide.lyx =================================================================== --- trunk/doc/userGuide.lyx 2014-10-14 12:45:31 UTC (rev 4944) +++ trunk/doc/userGuide.lyx 2014-10-14 16:55:00 UTC (rev 4945) @@ -17493,8 +17493,8 @@ \end_layout \begin_layout Subsection -Conditionally revert or forward an evolutionary process to a saved state - (operator +Conditionally revert an evolutionary process to a saved state (operator + \family typewriter RevertIf \family default @@ -17513,13 +17513,17 @@ \family default file) and will revert the current evolving population to the saved population if the condition is met. - There are at least two secenarios for which this operator can be very useful. + \end_layout \begin_layout Standard -One of the biggest problem with introducing a disease allele to an evolving - population (using an operator PointMutate) is that the disease allele will - very likely get lost because of genetic drift. +For example, one of the biggest problem with introducing a disease allele + to an evolving population (using an operator +\family typewriter +PointMutator +\family default +) is that the disease allele will very likely get lost because of genetic + drift. It is possible to simulate the allele frequency trajectory backward in time and follow the trajectory during the forward-time simulation phase ( @@ -17548,7 +17552,7 @@ \begin_layout Standard A natural way to simulate the introduction of disease alleles is therefore to terminate and restart the simulation whenever the disease allele gets - lost. + lost (or fixed). This could be done by splitting the evolutionary process into two stages. The disease allele is introduced at the second stage and the simulation will be terminated as soon as the introduced allele is lost (using a @@ -17612,59 +17616,15 @@ \end_layout \begin_layout Standard -The second application is to skip certain part of the evolutionary process - if it has been executed before. - For example, many forward-time simulation scenarios have a long burn-in - process. - It you simulate several replicates of the evolutionary process, you can - save the burned in initial population of the first replicate, and reuse - it for subsequent simulations. - In Example -\begin_inset CommandInset ref -LatexCommand ref -reference "ReuseBurnin" - -\end_inset - -, we save the population at generation 10 for the first replicate. - The next few replicates would load the saved file and jump directly to - generation 10. - Note that although the +Although operator \family typewriter RevertIf \family default - is applied at generation 0, it cannot be put in parameter -\family typewriter -initOps -\family default - because it otherwise evolve another 20 generations from generation 10. - However, combining the -\family typewriter -RevertIf -\family default - operator with complex demographic models can have unexpected results. - As shown in the second part of the example, because the demographic model - starts from the generation when the model is first called (intialized), - the evolutionary model will evolve from generation 10 to 30. - The problem can be avoided by calling -\family typewriter -RevertIf -\family default - at the second generation where the demographic model is already initialized - at generation 0. + can also be used to fast-forward an evolutionary process (skip to a previously + saved state directly), it is tricky to make it work with complex demographic + models and therefore not recommended. \end_layout -\begin_layout Standard -\begin_inset CommandInset include -LatexCommand lstinputlisting -filename "log/ReuseBurnIn.log" -lstparams "caption={Jump to the middle of an evolutionary process by uing a previously saved population.},label={ReuseBurnin}" - -\end_inset - - -\end_layout - \begin_layout Subsection Conditional during mating operator (operator \family typewriter @@ -25403,6 +25363,15 @@ properly. \end_layout +\begin_layout Itemize +saves population sizes of evolved generations, which makes it possible to + revert an evolutionary process to an previous state using operator +\family typewriter +RevertIf +\family default +. +\end_layout + \begin_layout Standard A demographic model can be defined in two ways. The first approach is to specify the size of subpopulations at each generation, Modified: trunk/src/demography.py =================================================================== --- trunk/src/demography.py 2014-10-14 12:45:31 UTC (rev 4944) +++ trunk/src/demography.py 2014-10-14 16:55:00 UTC (rev 4945) @@ -132,7 +132,7 @@ rates[-1][x[0] * n + x[1]] = r * 1.0 / len(neighbors) return rates -class BaseDemographicModel: +class DemographicModel: '''This class is the base class for all demographic models and provides common interface and utility functions for derived classes. A demographic model is essentially a callable Python object that encapsulates @@ -187,6 +187,9 @@ self.ops = list(ops) else: self.ops = [ops] + # + # history + self.size_cache = {} def _reset(self): if hasattr(self, '_start_gen'): @@ -270,8 +273,6 @@ # if size is None or size is [], return unchanged if not size: return - if '__locked__' in pop.vars() and pop.dvars().__locked__: - raise RuntimeError('Change population size of a locked population is not allowed.') named_size = self._convertToNamedSize(size) if pop.numSubPop() > 1: # merge subpopualtions @@ -504,6 +505,30 @@ def _setup(self, pop): return True + def _save_size(self, gen, sz): + if self.size_cache: + prev = [x for x in self.size_cache.keys() if x < gen] + if prev and self.size_cache[max(prev)] == sz: + return sz + self.size_cache[gen] = sz + return sz + + def _cached_size(self, gen): + if gen in self.size_cache: + return self.size_cache[gen] + if not self.size_cache: + raise RuntimeError('Failed to determine size for generation {}' + .format(gen)) + # if we look further, keep contant size + if gen > max(self.size_cache.keys()): + return self.size_cache[max(self.size_cache.keys())] + prev = [x for x in self.size_cache.keys() if x < gen] + if prev: + return self.size_cache[max(prev)] + else: + raise RuntimeError('Failed to determine size for generation {}' + .format(gen)) + def __call__(self, pop): # When calling the demographic function, there are two quite separate scenarios # @@ -512,51 +537,30 @@ # be two cases, the first one for initialization, the other one for instant # population change (in the case of InstantChangeModel). # - # 2. When the demographic function is called randomly, that is to say we might - # have a population already at the destination size, and we just want to start - # from that generation. In this case, all _fitTo_Size calls from the initialization - # part should be disallowed. The other call from InstantChangeModel can be - # called but the population is assumed to be already in place, so it should - # not have any impact. + # 2. When the demographic function is called randomly, only a previously evolved + # generation is allowed. This is because many demographic models depend on + # size of previous generations to determine new population size, which makes + # jumping forward impossible. # - # It is not quite clear how to tell these two cases apart. For a single-model or - # multi-stage model case, the first time the demographic model is called, - # _start_gen should be set and fitSize will be called. - # - # Then, for any subsequent calls, a single-model will not call initizliation again. - # For prograss-like calling, events will be triggered to change population size, - # for random access, we only assume that population already has the right size - # and do not call fitSize. - # # In the case of multi-model, random access will not call fitSize for any # intermediate steps. # # the first time this demographic function is called - self._randomAccess = False + self._use_cached = False if not hasattr(self, '_start_gen'): self._reset() self._start_gen = pop.dvars().gen self._last_gen = pop.dvars().gen # resize populations if necessary - if '__locked__' in pop.vars() and pop.dvars().__locked__: - # pop is assumed to be in destination population size - # the initial population size has to be determined from user input - # size, which is not possible if dynamic population size was - # specified. - self.init_size = self._extractSize(self._raw_init_size) - for sz in self.init_size: - if sz is None or isinstance(sz, (float, list, tuple)): - raise RuntimeError('Random access to an uninitialized demographic model with ' - 'dynamic population size is not allowed.') - else: - self._fitToSize(pop, self._raw_init_size) - # by this time, we should know the initial population size - self.init_size = pop.subPopSizes() + self._fitToSize(pop, self._raw_init_size) + # by this time, we should know the initial population size + self.init_size = pop.subPopSizes() # then we can set up model if the model depends on initial # population size self._setup(pop) elif pop.dvars().gen != self._last_gen + 1: - self._randomAccess = True + self._use_cached = True + return self._cached_size(pop.dvars().gen) # self._gen = pop.dvars().gen - self._start_gen self._last_gen = pop.dvars().gen @@ -568,11 +572,11 @@ self._reset() return [] if '_expected_size' in pop.vars(): - return pop.vars().pop('_expected_size') + return self._save_size(pop.dvars().gen, pop.vars().pop('_expected_size')) else: - return pop.subPopSizes() + return self._save_size(pop.dvars().gen, pop.subPopSizes()) -class ExponentialGrowthModel(BaseDemographicModel): +class ExponentialGrowthModel(DemographicModel): '''A model for exponential population growth with carry capacity''' def __init__(self, T=None, N0=[], NT=None, r=None, ops=[], infoFields=[]): '''An exponential population growth model that evolves a population from size @@ -586,7 +590,7 @@ population will keep contant after it reaches size ``NT``. Optionally, one or more operators (e.g. a migrator) ``ops`` can be applied to population. ''' - BaseDemographicModel.__init__(self, T, N0, ops, infoFields) + DemographicModel.__init__(self, T, N0, ops, infoFields) # if [x is None or x == [] for x in [T, NT, r]].count(True) > 1: raise ValueError('Please specify at least two parameters of T, NT and r') @@ -647,25 +651,25 @@ 'is expected') def __call__(self, pop): - if not BaseDemographicModel.__call__(self, pop): - return [] + if not DemographicModel.__call__(self, pop): + return [] # - # this model does not need differntiation between _randomAccess or not - # because it does not call fitSize to change population size. - # pop passed to _setup() is not used. + if self._use_cached: + return self._cached_size(pop.dvars().gen) # if self._gen == self.num_gens: return [] elif self.r is None: - return [self._expIntepolate(n0, nt, self.num_gens, self._gen) - for (n0, nt) in zip(self.init_size, self.NT)] + return self._save_size(pop.dvars().gen, + [self._expIntepolate(n0, nt, self.num_gens, self._gen) + for (n0, nt) in zip(self.init_size, self.NT)]) else: # with r ... - return [min(nt, int(n0 * math.exp(r * (self._gen + 1)))) - for (n0, nt, r) in zip(self.init_size, self.NT, self.r)] + return self._save_size(pop.dvars().gen, [min(nt, int(n0 * math.exp(r * (self._gen + 1)))) + for (n0, nt, r) in zip(self.init_size, self.NT, self.r)]) -class LinearGrowthModel(BaseDemographicModel): +class LinearGrowthModel(DemographicModel): '''A model for linear population growth with carry capacity.''' def __init__(self, T=None, N0=[], NT=None, r=None, ops=[], infoFields=[]): '''An linear population growth model that evolves a population from size @@ -679,7 +683,7 @@ population will keep contant after it reaches size ``NT``. Optionally, one or more operators (e.g. a migrator) ``ops`` can be applied to population. ''' - BaseDemographicModel.__init__(self, T, N0, ops, infoFields) + DemographicModel.__init__(self, T, N0, ops, infoFields) # if [x is None or x == [] for x in [T, NT, r]].count(True) > 1: raise ValueError('Please specify at least two parameters of T, NT and r') @@ -736,25 +740,25 @@ 'is expected') def __call__(self, pop): - if not BaseDemographicModel.__call__(self, pop): + if not DemographicModel.__call__(self, pop): return [] # - # this model does not need differntiation between _randomAccess or not - # because it does not call fitSize to change population size. + if self._use_cached: + return self._cached_size(pop.dvars().gen) # if self._gen == self.num_gens: return [] elif self.r is None: # no r use intepolation - return [self._linearIntepolate(n0, nt, self.num_gens, self._gen) - for (n0, nt) in zip(self.init_size, self.NT)] + return self._save_size(pop.dvars().gen, [self._linearIntepolate(n0, nt, self.num_gens, self._gen) + for (n0, nt) in zip(self.init_size, self.NT)]) else: # with r ... - return [min(nt, int(n0 + n0 * (self._gen + 1.) * r)) - for (n0, nt, r) in zip(self.init_size, self.NT, self.r)] + return self._save_size(pop.dvars().gen, [min(nt, int(n0 + n0 * (self._gen + 1.) * r)) + for (n0, nt, r) in zip(self.init_size, self.NT, self.r)]) -class InstantChangeModel(BaseDemographicModel): +class InstantChangeModel(DemographicModel): '''A model for instant population change (growth, resize, merge, split).''' def __init__(self, T=None, N0=[], G=[], NG=[], ops=[], infoFields=[], removeEmptySubPops=False): '''An instant population growth model that evolves a population @@ -769,7 +773,7 @@ operators should be passed to parameter ``infoFields``. If ``removeEmpty`` option is set to ``True``, empty subpopulation will be removed. This option can be used to remove subpopulations.''' - BaseDemographicModel.__init__(self, T, N0, ops, infoFields) + DemographicModel.__init__(self, T, N0, ops, infoFields) if isinstance(G, int): self.G = [G] if isinstance(NG, int) : @@ -794,37 +798,23 @@ def __call__(self, pop): # this one fixes N0... (self.init_size) - if not BaseDemographicModel.__call__(self, pop): + if not DemographicModel.__call__(self, pop): return [] + # + if self._use_cached: + return self._cached_size(pop.dvars().gen) # - if self._randomAccess or '__locked__' in pop.vars(): - # for random access, we just return calculate size - if not self.G or self._gen < self.G[0]: - sz = self.init_size - else: - # we do not expect pop.gen() to come in sequence - # so that we can change the population size at exactly - # self.G generations, therefore we have to check which - # internal self._gen falls into. - for i in range(len(self.G)): - if self._gen >= self.G[i] and (i == len(self.G) - 1 or self._gen < self.G[i+1]): - sz = self.NG[i] - break - if not isinstance(sz, (list, tuple)): - sz = [sz] - else: - # this is a sequential call, we do fitSize - if self._gen in self.G: - sz = self.NG[self.G.index(self._gen)] - self._fitToSize(pop, sz) - if self.removeEmptySubPops: - pop.removeSubPops([idx for idx,x in enumerate(pop.subPopSizes()) if x==0]) - sz = pop.subPopSizes() - return sz + if self._gen in self.G: + sz = self.NG[self.G.index(self._gen)] + self._fitToSize(pop, sz) + if self.removeEmptySubPops: + pop.removeSubPops([idx for idx,x in enumerate(pop.subPopSizes()) if x==0]) + sz = pop.subPopSizes() + return self._save_size(pop.dvars().gen, sz) -class AdmixtureModel(BaseDemographicModel): +class AdmixtureModel(DemographicModel): '''An population admixture model that mimicks the admixing of two subpopulations using either an HI (hybrid isolation) or CGF (continuous gene flow) model. In the HI model, an admixed population @@ -851,7 +841,7 @@ 1-mu percent of individuals in receipient population will be replaced by individuals in the doner population. ''' - BaseDemographicModel.__init__(self, T, N0, ops, infoFields) + DemographicModel.__init__(self, T, N0, ops, infoFields) if not model or model[0] not in ('HI', 'CGF') or len(model) < 4 or model[1] == model[2] or model[3] > 1: raise ValueError('model should be either ("HI", parent1, paren2, mu)' 'or ("CGF", receipient, doner, alpha)') @@ -878,11 +868,12 @@ return (int(N*mu + 0.5), int(N*(1-mu) + 0.5)) def __call__(self, pop): - if not BaseDemographicModel.__call__(self, pop): + if not DemographicModel.__call__(self, pop): return [] - if self._randomAccess: - raise ValueError('simuPOP does not yet support random access to ' - 'an admixture demographic model') + # + if self._use_cached: + return self._cached_size(pop.dvars().gen) + # if self.model[0] == 'HI': if self._gen != 0: return pop.subPopSizes() @@ -924,10 +915,10 @@ pop.splitSubPop(self.model[2], [request, pop.subPopSize(self.model[2]) - request]) # step3, merge the split indiiduals to subpopulation self.model[1] pop.mergeSubPops([self.model[1], self.model[2]]) - return pop.subPopSizes() + return self._save_size(pop.dvars().gen, pop.subPopSizes()) -class MultiStageModel(BaseDemographicModel): +class MultiStageModel(DemographicModel): '''A multi-stage demographic model that connects a number of demographic models. ''' def __init__(self, models, ops=[], infoFields=[]): @@ -946,7 +937,7 @@ total_gens = sum(gens) else: total_gens = -1 - BaseDemographicModel.__init__(self, numGens=total_gens, + DemographicModel.__init__(self, numGens=total_gens, initSize=models[0].init_size, ops=ops, infoFields=flds+infoFields) # self.models = models @@ -989,69 +980,40 @@ # so that the starting gen is after the current gen. We will # need to handle this case. # determines generation number internally as self.gen - if not BaseDemographicModel.__call__(self, pop): + if not DemographicModel.__call__(self, pop): return [] - if self._randomAccess: - # in this case, we assume that the population is already - # at the generation and is in need of size for next - pop.dvars().__locked__ = True - # self._start_gen works, so is self._gen - # but self._model_idx cannot be used, try to find it - g = 0 - for idx,model in enumerate(self.models): - if model.num_gens == 0: - continue - if model.num_gens < 0 or self._gen - g < model.num_gens: - self._model_idx = idx - # if we are jumping to an uninitialized model, we cannot - # initialize it with a later generation and has to forcefully - # initialize it with the correct starting generation. - if self._gen > g and not hasattr(model, '_last_gen'): - model._start_gen = g - model._last_gen = g - model.init_size = model._extractSize(model._raw_init_size) - for sz in model.init_size: - if sz is None or isinstance(sz, (float, list, tuple)): - raise RuntimeError('Random access to an uninitialized demographic model with ' - 'dynamic population size is not allowed.') - model._setup(pop) - model._randomAccess = True - sz = model.__call__(pop) - pop.vars().pop('__locked__') - return sz - g += model.num_gens - raise RuntimeError('Failed to jump to generation {}'.format(pop.dvars().gen)) - else: - # sequential access to the demographic model - # - # There are three cases - # 1. within the current model, a valid step is returned - # --> return - # 2. within the current model, a [] is returned by a - # terminator. - # --> proceed to the next available model, call, and return - # 3. at the end of the current model, - # --> proceed to the next available model, call, and return - # 4. at the beginning of a zero-step model - # --> call - # --> process to the next available model, call, and return - # - # in the middle - if self.models[self._model_idx].num_gens < 0 or \ - self.models[self._model_idx].num_gens > self._gen - self._model_start_gen: - sz = self.models[self._model_idx].__call__(pop) - if not sz: - sz = self._advance(pop) - elif self.models[self._model_idx].num_gens == 0: - sz = self.models[self._model_idx].__call__(pop) + # + if self._use_cached: + return self._cached_size(pop.dvars().gen) + # + # There are three cases + # 1. within the current model, a valid step is returned + # --> return + # 2. within the current model, a [] is returned by a + # terminator. + # --> proceed to the next available model, call, and return + # 3. at the end of the current model, + # --> proceed to the next available model, call, and return + # 4. at the beginning of a zero-step model + # --> call + # --> process to the next available model, call, and return + # + # in the middle + if self.models[self._model_idx].num_gens < 0 or \ + self.models[self._model_idx].num_gens > self._gen - self._model_start_gen: + sz = self.models[self._model_idx].__call__(pop) + if not sz: sz = self._advance(pop) - elif self.models[self._model_idx].num_gens <= self._gen - self._model_start_gen: - sz = self._advance(pop) - return sz + elif self.models[self._model_idx].num_gens == 0: + sz = self.models[self._model_idx].__call__(pop) + sz = self._advance(pop) + elif self.models[self._model_idx].num_gens <= self._gen - self._model_start_gen: + sz = self._advance(pop) + return self._save_size(pop.dvars().gen, sz) -class EventBasedModel(BaseDemographicModel): +class EventBasedModel(DemographicModel): '''An event based demographic model in which the demographic changes are triggered by demographic events such as population growth, split, join, and admixture. The population size will be kept constant if no event is applied @@ -1066,7 +1028,7 @@ ''' if isinstance(events, DemographicEvent): events = [events] - BaseDemographicModel.__init__(self, numGens=T, initSize=N0, + DemographicModel.__init__(self, numGens=T, initSize=N0, ops=ops + events, infoFields=infoFields) @@ -1740,8 +1702,9 @@ [0, m_AF_EU, m_AF_AS], [m_AF_EU, 0, m_EU_AS], [m_AF_AS, m_EU_AS, 0] - ]), - begin=-int(T_EU_AS/scale) + ], + subPops=['AF', 'EU', 'AS']), + begin=-int(T_EU_AS/scale), ), # admixture AdmixtureEvent( Modified: trunk/test/test_19_demography.py =================================================================== --- trunk/test/test_19_demography.py 2014-10-14 12:45:31 UTC (rev 4944) +++ trunk/test/test_19_demography.py 2014-10-14 16:55:00 UTC (rev 4945) @@ -467,16 +467,19 @@ # number. pop = Population(100, loci=1) m = LinearGrowthModel(N0=[0.5, 0.5], r=[0.1, 0.2], T=100) - for gen in [0, 10, 5, 7, 6, 8, 20]: + for gen in list(range(11)) + [5, 7, 6, 8]: pop.dvars().gen = gen self.assertEqual(m(pop), [50 + 5 * (gen+1), 50 + 10*(gen+1)]) # m = InstantChangeModel(N0=100, T=100, G=[10, 40], NG=[30, 50]) - for gen in [0, 5, 60, 10, 45, 15, 70]: + for gen in list(range(71)) + [0, 5, 60, 10, 45, 15, 70]: # the object has to be initialized with gen = 0 pop.dvars().gen = gen - self.assertEqual(list(m(pop)), - [{0:100, 5: 100, 10: 30, 15: 30, 45: 50, 60:50, 70:50}[gen]]) + if gen in [0, 5, 60, 10, 45, 15, 70]: + self.assertEqual(list(m(pop)), + [{0:100, 5: 100, 10: 30, 15: 30, 45: 50, 60:50, 70:50}[gen]]) + else: + m(pop) # # for multiple stage stuff m = MultiStageModel([ @@ -484,17 +487,20 @@ InstantChangeModel(N0=[1000, 1000], T=1), LinearGrowthModel(T=100, N0=[400, 600], r=0.01), ]) - for gen in [0, 10, 40, 20, 50, 110, 120, 10, 40, 50]: + for gen in list(range(121)) + [0, 10, 40, 20, 50, 110, 120, 10, 40, 50]: # the object has to be initialized with gen = 0 pop.dvars().gen = gen - self.assertEqual(m(pop), - {0:[1000 + 10], 50: [1000 + 10*51], - 10:[1000 + 10*11], - 20:[1000 + 10*21], - 40:[1000 + 10*41], - 110:[400 + 4*10, 600+6*10], - 120:[400 + 4*20, 600+6*20], - }[gen]) + if gen in [0, 10, 40, 20, 50, 110, 120, 10, 40, 50]: + self.assertEqual(list(m(pop)), + {0:[1000 + 10], 50: [1000 + 10*51], + 10:[1000 + 10*11], + 20:[1000 + 10*21], + 40:[1000 + 10*41], + 110:[400 + 4*10, 600+6*10], + 120:[400 + 4*20, 600+6*20], + }[gen]) + else: + m(pop) # # for multiple stage stuff m = MultiStageModel([ @@ -502,20 +508,23 @@ InstantChangeModel(N0=[1000, 1000], T=1), LinearGrowthModel(T=100, N0=[0.4, 0.6], r=0.01), ]) - for gen in [0, 10, 40, 20, 50, 100, 101, 110, 120, 10, 40, 50]: + for gen in list(range(121)) + [0, 10, 40, 20, 50, 100, 101, 110, 120, 10, 40, 50]: # the object has to be initialized with gen = 0 pop.dvars().gen = gen - self.assertEqual(m(pop), - {0:[1000 + 10], 50: [1000 + 10*51], - 10:[1000 + 10*11], - 20:[1000 + 10*21], - 40:[1000 + 10*41], - 99:[1000 + 10*100], - 100:[1000, 1000], - 101:[400 + 4*1, 600+6*1], - 110:[400 + 4*10, 600+6*10], - 120:[400 + 4*20, 600+6*20], - }[gen]) + if gen in [0, 10, 40, 20, 50, 100, 101, 110, 120, 10, 40, 50]: + self.assertEqual(list(m(pop)), + {0:[1000 + 10], 50: [1000 + 10*51], + 10:[1000 + 10*11], + 20:[1000 + 10*21], + 40:[1000 + 10*41], + 99:[1000 + 10*100], + 100:[1000, 1000], + 101:[400 + 4*1, 600+6*1], + 110:[400 + 4*10, 600+6*10], + 120:[400 + 4*20, 600+6*20], + }[gen]) + else: + m(pop) # def testEventBasedModel(self): This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |