Menu
▾
▴
simupop-devel
|
From: <si...@us...> - 2011-08-23 19:38:21
|
Revision: 4299
http://simupop.svn.sourceforge.net/simupop/?rev=4299&view=rev
Author: simupop
Date: 2011-08-23 19:38:15 +0000 (Tue, 23 Aug 2011)
Log Message:
-----------
Fix a documentation for recombinator
Modified Paths:
--------------
trunk/src/simuPOP_doc.i
trunk/src/transmitter.h
Modified: trunk/src/simuPOP_doc.i
===================================================================
--- trunk/src/simuPOP_doc.i 2011-08-23 17:44:47 UTC (rev 4298)
+++ trunk/src/simuPOP_doc.i 2011-08-23 19:38:15 UTC (rev 4299)
@@ -7704,7 +7704,7 @@
can be limited to a given set of loci, which can be a list of loci
indexes, names or ALL_AVAIL. Each locus in this list specifies a
recombination point between the locus and the locus immediately
- before it. Loci that are the first locus on each chromosome are
+ after it. Loci that are the last locus on each chromosome are
ignored. If a single recombination rate (parameter rates) is
specified, it will used for all loci (all loci or loci specified
by parameter loci), regardless of physical distances between
Modified: trunk/src/transmitter.h
===================================================================
--- trunk/src/transmitter.h 2011-08-23 17:44:47 UTC (rev 4298)
+++ trunk/src/transmitter.h 2011-08-23 19:38:15 UTC (rev 4299)
@@ -476,7 +476,7 @@
* be limited to a given set of \e loci, which can be a list of loci
* indexes, names or \c ALL_AVAIL. Each locus in this list specifies
* a recombination point between the locus and the locus immediately
- * \b before it. Loci that are the first locus on each chromosome are
+ * \b after it. Loci that are the last locus on each chromosome are
* ignored.
*
* If a single recombination rate (parameter \e rates) is specified, it
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2011-11-16 16:36:17
|
Revision: 4340
http://simupop.svn.sourceforge.net/simupop/?rev=4340&view=rev
Author: simupop
Date: 2011-11-16 16:36:08 +0000 (Wed, 16 Nov 2011)
Log Message:
-----------
Merge more files
Modified Paths:
--------------
trunk/src/__init__.py
trunk/src/main.cpp
trunk/src/qtrait.h
trunk/src/sandbox.py
Property Changed:
----------------
trunk/src/__init__.py
trunk/src/main.cpp
trunk/src/qtrait.h
trunk/src/sandbox.py
Modified: trunk/src/__init__.py
===================================================================
--- trunk/src/__init__.py 2011-11-16 16:23:48 UTC (rev 4339)
+++ trunk/src/__init__.py 2011-11-16 16:36:08 UTC (rev 4340)
@@ -286,6 +286,8 @@
from simuPOP.simuPOP_laop import *
elif simuOptions['AlleleType'] == 'binary':
from simuPOP.simuPOP_baop import *
+ elif simuOptions['AlleleType'] == 'mutant':
+ from simuPOP.simuPOP_muop import *
else:
from simuPOP.simuPOP_op import *
else:
@@ -295,6 +297,8 @@
from simuPOP.simuPOP_la import *
elif simuOptions['AlleleType'] == 'binary':
from simuPOP.simuPOP_ba import *
+ elif simuOptions['AlleleType'] == 'mutant':
+ from simuPOP.simuPOP_mu import *
else:
from simuPOP.simuPOP_std import *
Property changes on: trunk/src/__init__.py
___________________________________________________________________
Added: svn:mergeinfo
+ /branches/mutant/src/__init__.py:4262-4317
/branches/mutant_new/src/__init__.py:4318-4337
/branches/openMP/src/__init__.py:3999-4241
Modified: trunk/src/main.cpp
===================================================================
--- trunk/src/main.cpp 2011-11-16 16:23:48 UTC (rev 4339)
+++ trunk/src/main.cpp 2011-11-16 16:36:08 UTC (rev 4340)
@@ -87,6 +87,8 @@
#else
# ifdef BINARYALLELE
<< "binary"
+# elif MUTANTALLELE
+ << "mutant"
# else
<< "short"
# endif
Property changes on: trunk/src/main.cpp
___________________________________________________________________
Added: svn:mergeinfo
+ /branches/mutant/src/main.cpp:4262-4317
/branches/mutant_new/src/main.cpp:4318-4337
/branches/openMP/src/main.cpp:3999-4241
Modified: trunk/src/qtrait.h
===================================================================
--- trunk/src/qtrait.h 2011-11-16 16:23:48 UTC (rev 4339)
+++ trunk/src/qtrait.h 2011-11-16 16:36:08 UTC (rev 4340)
@@ -144,7 +144,7 @@
* field.
*/
PyQuanTrait(PyObject * func, const lociList & loci = vectoru(),
- const uintList ancGens = uintList(false), int begin = 0, int end = -1, int step = 1,
+ const uintList ancGens = uintList(NULL), 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()) :
BaseQuanTrait(ancGens, begin, end, step, at, reps, subPops, infoFields),
Property changes on: trunk/src/qtrait.h
___________________________________________________________________
Added: svn:mergeinfo
+ /branches/mutant/src/qtrait.h:4262-4317
/branches/mutant_new/src/qtrait.h:4318-4337
/branches/openMP/src/qtrait.h:3999-4241
Modified: trunk/src/sandbox.py
===================================================================
--- trunk/src/sandbox.py 2011-11-16 16:23:48 UTC (rev 4339)
+++ trunk/src/sandbox.py 2011-11-16 16:36:08 UTC (rev 4340)
@@ -57,6 +57,8 @@
from simuPOP_laop import *
elif simuOptions['AlleleType'] == 'binary':
from simuPOP_baop import *
+ elif simuOptions['AlleleType'] == 'mutant':
+ from simuPOP_muop import *
else:
from simuPOP_op import *
else:
@@ -66,6 +68,8 @@
from simuPOP_la import *
elif simuOptions['AlleleType'] == 'binary':
from simuPOP_ba import *
+ elif simuOptions['AlleleType'] == 'mutant':
+ from simuPOP_mu import *
else:
from simuPOP_std import *
Property changes on: trunk/src/sandbox.py
___________________________________________________________________
Added: svn:mergeinfo
+ /branches/mutant/src/sandbox.py:4262-4317
/branches/mutant_new/src/sandbox.py:4318-4337
/branches/openMP/src/sandbox.py:3999-4241
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2011-12-25 20:39:33
|
Revision: 4356
http://simupop.svn.sourceforge.net/simupop/?rev=4356&view=rev
Author: simupop
Date: 2011-12-25 20:39:27 +0000 (Sun, 25 Dec 2011)
Log Message:
-----------
Fix compiling errors
Modified Paths:
--------------
trunk/src/sandbox.cpp
trunk/src/stator.cpp
Modified: trunk/src/sandbox.cpp
===================================================================
--- trunk/src/sandbox.cpp 2011-12-25 20:26:08 UTC (rev 4355)
+++ trunk/src/sandbox.cpp 2011-12-25 20:39:27 UTC (rev 4356)
@@ -1190,7 +1190,7 @@
}
-void MutSpaceRecombinator::transmitGenotype0(Population & offPop, const Individual & parent,
+void MutSpaceRecombinator::transmitGenotype0(Population & pop, Population & offPop, const Individual & parent,
size_t offIndex, int ploidy) const
{
#ifdef BINARYALLELE
Modified: trunk/src/stator.cpp
===================================================================
--- trunk/src/stator.cpp 2011-12-25 20:26:08 UTC (rev 4355)
+++ trunk/src/stator.cpp 2011-12-25 20:39:27 UTC (rev 4356)
@@ -773,8 +773,8 @@
IndIterator ind = pop.indIterator(sp->subPop());
for (; ind.valid(); ++ind) {
GenoIterator it = ind->genoBegin();
+#ifdef MUTANTALLELE
GenoIterator it_end = ind->genoEnd();
-#ifdef MUTANTALLELE
compressed_vector<Allele>::index_array_type::iterator index_it = it.getIndexIterator();
compressed_vector<Allele>::index_array_type::iterator index_it_end = it_end.getIndexIterator();
compressed_vector<Allele>::value_array_type::iterator value_it = it.getValueIterator();
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ti...@us...> - 2012-02-16 22:08:40
|
Revision: 4373
http://simupop.svn.sourceforge.net/simupop/?rev=4373&view=rev
Author: tiwcpe8
Date: 2012-02-16 22:08:27 +0000 (Thu, 16 Feb 2012)
Log Message:
-----------
Fix error from copy funcition when runing test_14_transmitter.py
Modified Paths:
--------------
trunk/src/initializer.cpp
trunk/src/mutant_vector.h
Modified: trunk/src/initializer.cpp
===================================================================
--- trunk/src/initializer.cpp 2012-02-10 21:02:45 UTC (rev 4372)
+++ trunk/src/initializer.cpp 2012-02-16 22:08:27 UTC (rev 4373)
@@ -241,12 +241,10 @@
#else
IndIterator it = pop.indIterator(sp->subPop());
#endif
- std::cout << "Start init!!!!!!!\n";
for (size_t ii=0; it.valid(); ++it,++ii )
for (vectoru::iterator p = ploidy.begin(); p != ploidy.end(); ++p)
for (vectoru::const_iterator loc = loci.begin(); loc != loci.end(); ++loc, ++idx)
it->setAllele(ToAllele(m_genotype[idx % sz]), *loc, static_cast<int>(*p));
- std::cout << "End init!!!!!!!\n";
}
#ifdef _OPENMP
Modified: trunk/src/mutant_vector.h
===================================================================
--- trunk/src/mutant_vector.h 2012-02-10 21:02:45 UTC (rev 4372)
+++ trunk/src/mutant_vector.h 2012-02-16 22:08:27 UTC (rev 4373)
@@ -302,7 +302,7 @@
return *this;
if (m_com_index == -1)
m_com_index = 0;
- typename compressed_vector<T>::index_array_type::iterator lower =
+ typename compressed_vector<T>::index_array_type::const_iterator lower =
std::lower_bound(m_container->index_data().begin() + m_com_index, m_container->index_data().begin() + m_container->filled(), m_index);
m_com_index = lower - m_container->index_data().begin();
return *this;
@@ -317,7 +317,7 @@
return *this;
if (m_com_index == -1)
m_com_index = 0;
- typename compressed_vector<T>::index_array_type::iterator lower =
+ typename compressed_vector<T>::index_array_type::const_iterator lower =
std::lower_bound(m_container->index_data().begin() + m_com_index, m_container->index_data().begin() + m_container->filled(), m_index);
m_com_index = lower - m_container->index_data().begin();
return *this;
@@ -332,7 +332,7 @@
return *this;
if (m_com_index == -1)
m_com_index = 0;
- typename compressed_vector<T>::index_array_type::iterator lower =
+ typename compressed_vector<T>::index_array_type::const_iterator lower =
std::lower_bound(m_container->index_data().begin() + m_com_index, m_container->index_data().begin() + m_container->filled(), m_index);
m_com_index = lower - m_container->index_data().begin();
return *this;
@@ -347,7 +347,7 @@
return *this;
if (m_com_index == -1)
m_com_index = 0;
- typename compressed_vector<T>::index_array_type::iterator lower =
+ typename compressed_vector<T>::index_array_type::const_iterator lower =
std::lower_bound(m_container->index_data().begin() + m_com_index, m_container->index_data().begin() + m_container->filled(), m_index);
m_com_index = lower - m_container->index_data().begin();
return *this;
@@ -364,7 +364,10 @@
{
iterator result = *this;
result.m_index += iter.m_index;
- result.m_com_index = findPositionIndexData(result);
+ typename compressed_vector<T>::index_array_type::const_iterator lower =
+ std::lower_bound(m_container->index_data().begin(), m_container->index_data().begin() + m_container->filled(), result.m_index);
+ result.m_com_index = lower - m_container->index_data().begin();
+
return result;
}
@@ -372,12 +375,15 @@
{
iterator result = *this;
result.m_index += size;
- result.m_com_index = findPositionIndexData(result);
+ typename compressed_vector<T>::index_array_type::const_iterator lower =
+ std::lower_bound(m_container->index_data().begin(), m_container->index_data().begin() + m_container->filled(), result.m_index);
+ result.m_com_index = lower - m_container->index_data().begin();
return result;
}
size_t operator - (iterator & iter)
{
+
return (*this).m_index - iter.m_index;
}
@@ -386,7 +392,11 @@
{
iterator result = *this;
result.m_index -= iter.m_index;
- result.m_com_index = findPositionIndexData(result);
+ if (m_com_index == -1)
+ return result;
+ typename compressed_vector<T>::index_array_type::const_reverse_iterator upper =
+ std::upper_bound(m_container->index_data().rbegin() + (m_container->index_data().size() - m_com_index), m_container->index_data().rend(), result.m_index);
+ result.m_com_index = m_container->index_data().rend() - upper;
return result;
}
@@ -394,7 +404,11 @@
{
iterator result = *this;
result.m_index -= size;
- result.m_com_index = findPositionIndexData(result);
+ if (m_com_index == -1)
+ return result;
+ typename compressed_vector<T>::index_array_type::const_reverse_iterator upper =
+ std::upper_bound(m_container->index_data().rbegin() + (m_container->index_data().size() - m_com_index), m_container->index_data().rend(), result.m_index);
+ result.m_com_index = m_container->index_data().rend() - upper;
return result;
}
@@ -413,8 +427,8 @@
return 0;
}
- size_t findPositionIndexData (iterator & it) const {
- for (size_t i = 0; i < it.m_container->filled(); i++) {
+ size_t findPositionIndexData (iterator & it, size_t idxStart = 0) const {
+ for (size_t i = idxStart; i < it.m_container->filled(); i++) {
if (i < it.m_container->filled() - 1 && it.m_index > it.m_container->index_data()[i]) {
continue;
} else if (i == it.m_container->filled() - 1 && it.m_index > it.m_container->index_data()[i]) {
@@ -434,18 +448,19 @@
if (m_com_index == -1 )
return m_container->index_data().begin();
else if (m_com_index > (ssize_t)m_container->filled())
- return m_container->index_data().end();
+ return m_container->index_data().begin() + m_container->filled();
else
return m_container->index_data().begin() + m_com_index;
}
typename compressed_vector<T>::value_array_type::iterator getValueIterator () {
//return m_container->value_data().begin() + findPositionIndexData();
- //return m_container->value_data().begin() + (getIndexIterator() - m_container->index_data().begin());
+ //size_t com_index = getIndexIterator() - m_container->index_data().begin();
+ //return m_container->value_data().begin() + com_index;
if (m_com_index == -1 )
return m_container->value_data().begin();
else if (m_com_index > (ssize_t)m_container->filled())
- return m_container->value_data().end();
+ return m_container->value_data().begin() + m_container->filled();
else
return m_container->value_data().begin() + m_com_index;
}
@@ -560,12 +575,12 @@
size_t size = end - begin;
mutant_vectora::iterator it_end = it + size;
if (it_end.getIndex() <= it_end.getContainer()->size()) {
- compressed_vector<Allele>::value_array_type::iterator src_value_begin = begin.getValueIterator();
- compressed_vector<Allele>::value_array_type::iterator src_value_end = end.getValueIterator();
- compressed_vector<Allele>::index_array_type::iterator src_index_begin = begin.getIndexIterator();
+ compressed_vector<Allele>::index_array_type::const_iterator src_index_begin = begin.getIndexIterator();
+ compressed_vector<Allele>::value_array_type::const_iterator src_value_begin = begin.getValueIterator();
+ compressed_vector<Allele>::value_array_type::const_iterator src_value_end = end.getValueIterator();
compressed_vector<Allele>::value_array_type::iterator dest_value_begin = it.getValueIterator();
compressed_vector<Allele>::value_array_type::iterator dest_value_end = it_end.getValueIterator();
- compressed_vector<Allele>::index_array_type::iterator dest_index_begin;
+ compressed_vector<Allele>::index_array_type::iterator dest_index_begin = it.getIndexIterator();
size_t src_size = src_value_end - src_value_begin;
size_t dest_size = dest_value_end - dest_value_begin;
size_t src_idx_num_begin = begin.getIndex();
@@ -577,20 +592,20 @@
it.getContainer()->reserve(2 * filled_size + diff_size, true);
//After using reserve(), dest_value_begin no longer valid, get a new one
dest_value_begin = it.getValueIterator();
+ dest_index_begin = it.getIndexIterator();
+ src_index_begin = begin.getIndexIterator();
}
- dest_index_begin = it.getIndexIterator();
std::copy_backward(dest_index_begin, it.getContainer()->index_data().begin() + filled_size, it.getContainer()->index_data().begin() + filled_size + diff_size);
for (size_t i = 0; i < src_size; i++) {
size_t range = *(src_index_begin + i) - src_idx_num_begin;
*(dest_index_begin+i) = dest_idx_num_begin + range;
}
std::copy_backward(dest_value_begin, it.getContainer()->value_data().begin() + filled_size, it.getContainer()->value_data().begin() + filled_size + diff_size);
- std::copy(src_value_begin, src_value_end, dest_value_begin);
+ std::copy(begin.getValueIterator(), end.getValueIterator(), dest_value_begin);
it.getContainer()->set_filled(filled_size + diff_size);
} else if (src_size < dest_size) {
size_t diff_size = dest_size - src_size;
size_t filled_size = it.getContainer()->filled();
- dest_index_begin = it.getIndexIterator();
std::copy(dest_index_begin + diff_size, it.getContainer()->index_data().begin() + filled_size, dest_index_begin);
for (size_t i = 0; i < src_size; i++) {
size_t range = *(src_index_begin + i) - src_idx_num_begin;
@@ -601,7 +616,6 @@
it.getContainer()->set_filled(filled_size - diff_size);
} else {
- dest_index_begin = it.getIndexIterator();
for (size_t i = 0; i < src_size; i++) {
size_t range = *(src_index_begin + i) - src_idx_num_begin;
*(dest_index_begin + i) = dest_idx_num_begin + range;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2012-02-24 05:40:04
|
Revision: 4389
http://simupop.svn.sourceforge.net/simupop/?rev=4389&view=rev
Author: simupop
Date: 2012-02-24 05:39:58 +0000 (Fri, 24 Feb 2012)
Log Message:
-----------
Minor warning fixes for MS VC compiler.
Modified Paths:
--------------
trunk/src/genoStru.cpp
trunk/src/sandbox.cpp
trunk/src/stator.cpp
Modified: trunk/src/genoStru.cpp
===================================================================
--- trunk/src/genoStru.cpp 2012-02-24 04:19:33 UTC (rev 4388)
+++ trunk/src/genoStru.cpp 2012-02-24 05:39:58 UTC (rev 4389)
@@ -258,7 +258,7 @@
//
m_mitochondrial = -1;
size_t mt_len = 0;
- for (size_t i = 0; i < m_chromTypes.size(); ++i) {
+ for (int i = 0; i < static_cast<int>(m_chromTypes.size()); ++i) {
if (m_chromTypes[i] == MITOCHONDRIAL) {
DBG_ASSERT(m_mitochondrial == -1, ValueError,
"Only one mitochondrial chromosome can be specified");
Modified: trunk/src/sandbox.cpp
===================================================================
--- trunk/src/sandbox.cpp 2012-02-24 04:19:33 UTC (rev 4388)
+++ trunk/src/sandbox.cpp 2012-02-24 05:39:58 UTC (rev 4389)
@@ -484,7 +484,7 @@
mutants.insert(mutLoc);
}
*/
- ind.setAllele(1, mutLoc, p, ch);
+ ind.setAllele(1, mutLoc, int(p), int(ch));
} // while
} // each individual
} // each subpopulation
Modified: trunk/src/stator.cpp
===================================================================
--- trunk/src/stator.cpp 2012-02-24 04:19:33 UTC (rev 4388)
+++ trunk/src/stator.cpp 2012-02-24 05:39:58 UTC (rev 4389)
@@ -697,13 +697,13 @@
const vectoru & loci = m_loci.elems(&pop);
DBG_DO(DBG_STATOR, cerr << "Count number of segregating sites for " << loci.size() << " loci " << endl);
- std::set<ULONG> allSegSites;
+ std::set<size_t> allSegSites;
// for each subpopulation.
subPopList subPops = m_subPops.expandFrom(pop);
subPopList::const_iterator sp = subPops.begin();
subPopList::const_iterator spEnd = subPops.end();
for (; sp != spEnd; ++sp) {
- std::set<ULONG> segSites;
+ std::set<size_t> segSites;
pop.activateVirtualSubPop(*sp);
// go through all loci
@@ -2936,8 +2936,8 @@
for (size_t idx = 0; idx < loci.size(); ++idx) {
size_t chromType = pop.chromType(pop.chromLocusPair(loci[idx]).first);
if (idx == 0)
- use_observed_het = pop.ploidy() == 2 and (chromType == AUTOSOME || chromType == CUSTOMIZED);
- else if ((pop.ploidy() == 2 and (chromType == AUTOSOME || chromType == CUSTOMIZED)) != use_observed_het)
+ use_observed_het = pop.ploidy() == 2 && (chromType == AUTOSOME || chromType == CUSTOMIZED);
+ else if ((pop.ploidy() == 2 && (chromType == AUTOSOME || chromType == CUSTOMIZED)) != use_observed_het)
throw ValueError("Structure statistics can only be estimated from loci on chromosomes of the same type, because other wise the observed number of alleles will be different.");
}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2012-02-27 04:52:12
|
Revision: 4396
http://simupop.svn.sourceforge.net/simupop/?rev=4396&view=rev
Author: simupop
Date: 2012-02-27 04:52:06 +0000 (Mon, 27 Feb 2012)
Log Message:
-----------
Fix a compiling error caused by the addition of a warning message
Modified Paths:
--------------
trunk/src/mating.cpp
trunk/src/simuPOP_doc.i
Modified: trunk/src/mating.cpp
===================================================================
--- trunk/src/mating.cpp 2012-02-27 04:26:12 UTC (rev 4395)
+++ trunk/src/mating.cpp 2012-02-27 04:52:06 UTC (rev 4396)
@@ -1647,6 +1647,7 @@
DBG_WARNIF(spSize == 0, "WARNING: One of the parental (virtual) subpopulation is empty and will not "
"produce any offspring.");
w_pos[i] = static_cast<double>(spSize);
+ }
}
DBG_DO(DBG_DEVEL, cerr << "Positive mating scheme weights: " << w_pos << '\n'
<< "Negative mating scheme weights: " << w_neg << endl);
Modified: trunk/src/simuPOP_doc.i
===================================================================
--- trunk/src/simuPOP_doc.i 2012-02-27 04:26:12 UTC (rev 4395)
+++ trunk/src/simuPOP_doc.i 2012-02-27 04:52:06 UTC (rev 4396)
@@ -3781,7 +3781,7 @@
aabb for diploid two-locus cases,
* and in general 2**n for diploid and 3**n for haploid cases if
there are n loci. This operator does not support haplodiploid
- populations and sex chromosomes.
+ populations, sex and mitochondrial chromosomes.
";
@@ -9370,7 +9370,13 @@
Wright's fixation index but it is extended for multi-allele (more
than two alleles) and multi-loci cases. This statistics should be
used if you would like to obtain a true Fst value of a large
- Population.
+ Population. Nei's Gst uses only allele frequency information so it
+ is available for all population type (haploid, diploid etc). Weir
+ and Cockerham's Fst uses heterozygosity frequency so it is best
+ for autosome of diploid populations. For non-diploid population,
+ sex, and mitochondrial DNAs, simuPOP uses expected heterozygosity
+ (1 - sum p_i^2) when heterozygosity is needed. These statistics
+ output the following variables:
* F_st (default) The WC84 Fst statistic estimated for all
specified loci.
* F_is The WC84 Fis statistic estimated for all specified loci.
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2012-02-28 20:23:59
|
Revision: 4398
http://simupop.svn.sourceforge.net/simupop/?rev=4398&view=rev
Author: simupop
Date: 2012-02-28 20:23:52 +0000 (Tue, 28 Feb 2012)
Log Message:
-----------
Add a map type to return mutants
Modified Paths:
--------------
trunk/src/simuPOP_cfg.h
trunk/src/simuPOP_common.i
Modified: trunk/src/simuPOP_cfg.h
===================================================================
--- trunk/src/simuPOP_cfg.h 2012-02-27 22:45:20 UTC (rev 4397)
+++ trunk/src/simuPOP_cfg.h 2012-02-28 20:23:52 UTC (rev 4398)
@@ -367,6 +367,7 @@
typedef std::map<int, double> intDict;
typedef std::map<size_t, double> uintDict;
typedef std::map<vectori, double> tupleDict;
+typedef std::map<size_t, size_t> alleleDict;
#define ValidPyObject(obj) (obj != NULL && obj != Py_None)
#define InvalidPyObject(obj) (obj == NULL || obj == Py_None)
Modified: trunk/src/simuPOP_common.i
===================================================================
--- trunk/src/simuPOP_common.i 2012-02-27 22:45:20 UTC (rev 4397)
+++ trunk/src/simuPOP_common.i 2012-02-28 20:23:52 UTC (rev 4398)
@@ -171,6 +171,7 @@
%template() vector<double>; /* e.g. lociPos */
%template() vector<long>; /* e.g. vspID(vectori) */
%template() map<vector<long>, double>; /* e.g. MapSelector */
+ %template() map<size_t, size_t>; /* e.g. mutant list */
}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2012-02-28 20:38:40
|
Revision: 4399
http://simupop.svn.sourceforge.net/simupop/?rev=4399&view=rev
Author: simupop
Date: 2012-02-28 20:38:31 +0000 (Tue, 28 Feb 2012)
Log Message:
-----------
Cannot use a map because there can be several alleles for each location
Modified Paths:
--------------
trunk/src/simuPOP_cfg.h
trunk/src/simuPOP_common.i
Modified: trunk/src/simuPOP_cfg.h
===================================================================
--- trunk/src/simuPOP_cfg.h 2012-02-28 20:23:52 UTC (rev 4398)
+++ trunk/src/simuPOP_cfg.h 2012-02-28 20:38:31 UTC (rev 4399)
@@ -367,7 +367,7 @@
typedef std::map<int, double> intDict;
typedef std::map<size_t, double> uintDict;
typedef std::map<vectori, double> tupleDict;
-typedef std::map<size_t, size_t> alleleDict;
+typedef std::vector<std::pair<size_t, size_t> > mutantList;
#define ValidPyObject(obj) (obj != NULL && obj != Py_None)
#define InvalidPyObject(obj) (obj == NULL || obj == Py_None)
Modified: trunk/src/simuPOP_common.i
===================================================================
--- trunk/src/simuPOP_common.i 2012-02-28 20:23:52 UTC (rev 4398)
+++ trunk/src/simuPOP_common.i 2012-02-28 20:38:31 UTC (rev 4399)
@@ -171,7 +171,7 @@
%template() vector<double>; /* e.g. lociPos */
%template() vector<long>; /* e.g. vspID(vectori) */
%template() map<vector<long>, double>; /* e.g. MapSelector */
- %template() map<size_t, size_t>; /* e.g. mutant list */
+ %template() vector<pair<size_t, size_t> >; /* e.g. mutant list */
}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ti...@us...> - 2012-02-28 21:24:52
|
Revision: 4400
http://simupop.svn.sourceforge.net/simupop/?rev=4400&view=rev
Author: tiwcpe8
Date: 2012-02-28 21:24:46 +0000 (Tue, 28 Feb 2012)
Log Message:
-----------
Add mutants function in individual.cpp and .h
Modified Paths:
--------------
trunk/src/individual.cpp
trunk/src/individual.h
Modified: trunk/src/individual.cpp
===================================================================
--- trunk/src/individual.cpp 2012-02-28 20:38:31 UTC (rev 4399)
+++ trunk/src/individual.cpp 2012-02-28 21:24:46 UTC (rev 4400)
@@ -238,7 +238,82 @@
m_genoPtr + beginP * totNumLoci() + chromEnd(endCh - 1));
}
+#ifdef MUTANTALLELE
+mutantList Individual::mutants(const uintList & ply, const uintList & ch) {
+ mutantList mutantAllele;
+ DBG_WARNIF(true, "The returned object of function Individual.genotype() is a special "
+ "carray object that reflects the underlying genotype of an individual. "
+ "It will become invalid once the population changes. Please use "
+ "list(ind.genotype()) if you would like to keep a copy of genotypes");
+ size_t beginP = 0;
+ size_t endP = 0;
+ size_t beginCh = 0;
+ size_t endCh = 0;
+
+ if (ply.allAvail())
+ endP = ploidy();
+ else {
+ const vectoru & ploidys = ply.elems();
+ if (ploidys.empty()) {
+ return mutantAllele;
+ }
+ beginP = ploidys[0];
+ endP = ploidys[0];
+ CHECKRANGEPLOIDY(static_cast<size_t>(beginP));
+ for (size_t i = 1; i < ploidys.size(); ++i) {
+ CHECKRANGEPLOIDY(static_cast<size_t>(ploidys[i]));
+ if (beginP > ploidys[i])
+ beginP = ploidys[i];
+ if (endP < ploidys[i])
+ endP = ploidys[i];
+ }
+ ++endP;
+ }
+ if (ch.allAvail())
+ endCh = numChrom();
+ else {
+ const vectoru & chroms = ch.elems();
+ if (chroms.empty()) {
+ return mutantAllele;
+ }
+ beginCh = chroms[0];
+ endCh = chroms[0];
+ CHECKRANGECHROM(static_cast<size_t>(beginCh));
+ for (size_t i = 1; i < chroms.size(); ++i) {
+ CHECKRANGECHROM(static_cast<size_t>(chroms[i]));
+ if (beginCh > chroms[i])
+ beginCh = chroms[i];
+ if (endCh < chroms[i])
+ endCh = chroms[i];
+ }
+ ++endCh;
+ }
+
+ if (endP > beginP + 1) {
+ // has to be all chromosomes
+ DBG_FAILIF(beginCh != 0 || endCh != numChrom(), ValueError,
+ "If multiple ploidy are chosen, all chromosomes has to be chosen.");
+ compressed_vector<Allele>::index_array_type::iterator idx_beginP = (m_genoPtr + beginP * totNumLoci()).getIndexIterator();
+ compressed_vector<Allele>::index_array_type::iterator idx_endP = (m_genoPtr + endP * totNumLoci()).getIndexIterator();
+ compressed_vector<Allele>::value_array_type::iterator value_beginP = (m_genoPtr + beginP * totNumLoci()).getValueIterator();
+ mutantAllele.reserve(idx_endP - idx_beginP);
+ for(;idx_beginP != idx_endP; ++idx_beginP, ++value_beginP)
+ mutantAllele.push_back(std::pair<size_t, size_t>(*idx_beginP, *value_beginP));
+ return mutantAllele;
+ } else {
+ compressed_vector<Allele>::index_array_type::iterator idx_beginP = (m_genoPtr + beginP * totNumLoci() + chromBegin(beginCh)).getIndexIterator();
+ compressed_vector<Allele>::index_array_type::iterator idx_endP = (m_genoPtr + endP * totNumLoci() + chromEnd(endCh -1)).getIndexIterator();
+ compressed_vector<Allele>::value_array_type::iterator value_beginP = (m_genoPtr + beginP * totNumLoci() + chromBegin(beginCh)).getValueIterator();
+ mutantAllele.reserve(idx_endP - idx_beginP);
+ for(;idx_beginP != idx_endP; ++idx_beginP, ++value_beginP)
+ mutantAllele.push_back(std::pair<size_t, size_t>(*idx_beginP, *value_beginP));
+ //mutantAllele[*idx_beginP] = *value_beginP;//.insert(std::pair<size_t, size_t>(*idx_beginP, *value_beginP);
+ return mutantAllele;
+ }
+}
+#endif
+
// Fix me: This one has to optimize
PyObject * Individual::genoAtLoci(const lociList & lociList)
{
Modified: trunk/src/individual.h
===================================================================
--- trunk/src/individual.h 2012-02-28 20:38:31 UTC (rev 4399)
+++ trunk/src/individual.h 2012-02-28 21:24:46 UTC (rev 4400)
@@ -197,7 +197,9 @@
*/
PyObject * genotype(const uintList & ploidy = uintList(), const uintList & chroms = uintList());
-
+#ifdef MUTANTALLELE
+ mutantList mutants(const uintList & ploidy = uintList(), const uintList & chroms = uintList());
+#endif
/** CPPONLY
* Return a Python object with alleles at specified loci. This function
* is usually used to collect alleles to send to a user-provided function.
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ti...@us...> - 2012-03-01 21:17:56
|
Revision: 4404
http://simupop.svn.sourceforge.net/simupop/?rev=4404&view=rev
Author: tiwcpe8
Date: 2012-03-01 21:17:49 +0000 (Thu, 01 Mar 2012)
Log Message:
-----------
Fix bug in copy function and removeMarkedIndividuals
Modified Paths:
--------------
trunk/src/mutant_vector.h
trunk/src/population.cpp
Modified: trunk/src/mutant_vector.h
===================================================================
--- trunk/src/mutant_vector.h 2012-03-01 17:15:33 UTC (rev 4403)
+++ trunk/src/mutant_vector.h 2012-03-01 21:17:49 UTC (rev 4404)
@@ -576,13 +576,14 @@
mutant_vectora::iterator it_end = it + size;
if (it_end.getIndex() <= it_end.getContainer()->size()) {
compressed_vector<Allele>::index_array_type::const_iterator src_index_begin = begin.getIndexIterator();
+ compressed_vector<Allele>::index_array_type::const_iterator src_index_end = end.getIndexIterator();
compressed_vector<Allele>::value_array_type::const_iterator src_value_begin = begin.getValueIterator();
compressed_vector<Allele>::value_array_type::const_iterator src_value_end = end.getValueIterator();
+ compressed_vector<Allele>::index_array_type::iterator dest_index_begin = it.getIndexIterator();
+ compressed_vector<Allele>::index_array_type::iterator dest_index_end = it_end.getIndexIterator();
compressed_vector<Allele>::value_array_type::iterator dest_value_begin = it.getValueIterator();
- compressed_vector<Allele>::value_array_type::iterator dest_value_end = it_end.getValueIterator();
- compressed_vector<Allele>::index_array_type::iterator dest_index_begin = it.getIndexIterator();
- size_t src_size = src_value_end - src_value_begin;
- size_t dest_size = dest_value_end - dest_value_begin;
+ size_t src_size = src_index_end - src_index_begin;
+ size_t dest_size = dest_index_end - dest_index_begin;
size_t src_idx_num_begin = begin.getIndex();
size_t dest_idx_num_begin = it.getIndex();
if (src_size > dest_size) {
Modified: trunk/src/population.cpp
===================================================================
--- trunk/src/population.cpp 2012-03-01 17:15:33 UTC (rev 4403)
+++ trunk/src/population.cpp 2012-03-01 21:17:49 UTC (rev 4404)
@@ -995,7 +995,7 @@
if (oldInd != newInd) {
*newInd = *oldInd;
#ifdef MUTANTALLELE
- simuPOP::copy(oldPtr, oldPtr + step, newPtr);
+ simuPOP::copy(oldPtr + 0, oldPtr + step, newPtr);
#else
copy(oldPtr, oldPtr + step, newPtr);
#endif
@@ -1021,7 +1021,7 @@
// do not remove.
if (oldPtr != newPtr) {
#ifdef MUTANTALLELE
- simuPOP::copy(oldPtr, oldPtr + step * spSize, newPtr);
+ simuPOP::copy(oldPtr + 0, oldPtr + step * spSize, newPtr);
#else
copy(oldPtr, oldPtr + step * spSize, newPtr);
#endif
@@ -1076,7 +1076,7 @@
if (oldInd != newInd) {
*newInd = *oldInd;
#ifdef MUTANTALLELE
- simuPOP::copy(oldPtr, oldPtr + step, newPtr);
+ simuPOP::copy(oldPtr + 0, oldPtr + step, newPtr);
#else
copy(oldPtr, oldPtr + step, newPtr);
#endif
@@ -1763,7 +1763,7 @@
//
copy(oldInd, oldInd + spSize, newInd);
#ifdef MUTANTALLELE
- simuPOP::copy(oldPtr, oldPtr + step * spSize, newPtr);
+ simuPOP::copy(oldPtr + 0, oldPtr + step * spSize, newPtr);
#else
copy(oldPtr, oldPtr + step * spSize, newPtr);
#endif
@@ -1793,7 +1793,7 @@
++newSize;
*newInd = *oldInd;
#ifdef MUTANTALLELE
- simuPOP::copy(oldPtr, oldPtr + step, newPtr);
+ simuPOP::copy(oldPtr + 0, oldPtr + step, newPtr);
#else
copy(oldPtr, oldPtr + step, newPtr);
#endif
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2012-03-02 15:15:37
|
Revision: 4406
http://simupop.svn.sourceforge.net/simupop/?rev=4406&view=rev
Author: simupop
Date: 2012-03-02 15:15:25 +0000 (Fri, 02 Mar 2012)
Log Message:
-----------
Allow continuous use of RevertFixedSites in the sandbox module
Modified Paths:
--------------
trunk/src/mutator.h
trunk/src/sandbox.cpp
trunk/src/sandbox.h
trunk/src/sandbox.py
trunk/src/simuPOP_common.i
Modified: trunk/src/mutator.h
===================================================================
--- trunk/src/mutator.h 2012-03-02 05:52:54 UTC (rev 4405)
+++ trunk/src/mutator.h 2012-03-02 15:15:25 UTC (rev 4406)
@@ -645,17 +645,19 @@
};
-/** This operator checks all loci of a population and revert a mutant to
- * wildtype allele if it is fixed in the population. If a list of (virtual)
- * subpopulations are specified, alleles are reverted in each subpopulation,
- * regardless if the allele is fixed in other subpopulations.
+/** This operator checks all or specified loci of a population and revert all
+ * mutants at loci to wildtype alleles if they are fixed in the population.
+ * If a list of (virtual) subpopulations are specified, alleles are reverted
+ * if they are fixed in each subpopulation, regardless if the alleles are
+ * fixed in other subpopulations.
*/
class RevertFixedSites : public BaseOperator
{
public:
- /** Create an operator to revert alleles at fixed loci from value non-zero to zero.
- * If parameter \e subPops are specified, only individuals in these subpopulations
- * are considered.
+ /** Create an operator to set all alleles to zero at specified (parameter
+ * \e loci) or all loci if they are fixed (having no zero-allele) at these
+ * loci. If parameter \e subPops are specified, only individuals in these
+ * subpopulations are considered.
*/
RevertFixedSites(const lociList & loci = lociList(),
const stringFunc & output = "", int begin = 0, int end = -1, int step = 1,
@@ -688,7 +690,7 @@
string describe(bool format = true) const
{
(void)format; // avoid warning about unused parameter
- return "Revert fixed alleles to wildtype allele if it is fixed in the population.";
+ return "Revert fixed alleles to wildtype allele 0 if they are fixed in the population.";
}
private:
Modified: trunk/src/sandbox.cpp
===================================================================
--- trunk/src/sandbox.cpp 2012-03-02 05:52:54 UTC (rev 4405)
+++ trunk/src/sandbox.cpp 2012-03-02 15:15:25 UTC (rev 4406)
@@ -27,10 +27,11 @@
#include "sandbox.h"
namespace simuPOP {
+namespace sandbox {
#ifdef MUTANTALLELE
-bool SB_RevertFixedSites::apply(Population & pop) const
+bool RevertFixedSites::apply(Population & pop) const
{
if (pop.popSize() == 0 || pop.totNumLoci() == 0)
return true;
@@ -712,7 +713,7 @@
// ########################### THE OLD VERSION FOR LONG ALLELE TYPE ####
#else
-bool SB_RevertFixedSites::apply(Population & pop) const
+bool RevertFixedSites::apply(Population & pop) const
{
if (pop.popSize() == 0 || pop.totNumLoci() == 0)
return true;
@@ -1404,5 +1405,5 @@
}
#endif
-
}
+}
Modified: trunk/src/sandbox.h
===================================================================
--- trunk/src/sandbox.h 2012-03-02 05:52:54 UTC (rev 4405)
+++ trunk/src/sandbox.h 2012-03-02 15:15:25 UTC (rev 4406)
@@ -43,6 +43,7 @@
#include <set>
namespace simuPOP {
+namespace sandbox {
#ifdef MUTANTALLELE
@@ -51,13 +52,13 @@
* to wildtype allele if it is fixed in the population. If a valid output is
* specifieid, fixed alleles will be outputed with a leading generation number.
*/
-class SB_RevertFixedSites : public BaseOperator
+class RevertFixedSites : public BaseOperator
{
public:
/** Create an operator to revert alleles at fixed loci from value 1 to 0.
* Parameter \e subPops is ignored.
*/
- SB_RevertFixedSites(const stringFunc & output = "", int begin = 0, int end = -1, int step = 1,
+ RevertFixedSites(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())
@@ -67,7 +68,7 @@
/// destructor
- virtual ~SB_RevertFixedSites()
+ virtual ~RevertFixedSites()
{
}
@@ -75,7 +76,7 @@
/// HIDDEN Deep copy of a Migrator
virtual BaseOperator * clone() const
{
- return new SB_RevertFixedSites(*this);
+ return new RevertFixedSites(*this);
}
@@ -394,13 +395,13 @@
* to wildtype allele if it is fixed in the population. If a valid output is
* specifieid, fixed alleles will be outputed with a leading generation number.
*/
-class SB_RevertFixedSites : public BaseOperator
+class RevertFixedSites : public BaseOperator
{
public:
/** Create an operator to revert alleles at fixed loci from value 1 to 0.
* Parameter \e subPops is ignored.
*/
- SB_RevertFixedSites(const stringFunc & output = "", int begin = 0, int end = -1, int step = 1,
+ RevertFixedSites(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())
@@ -410,7 +411,7 @@
/// destructor
- virtual ~SB_RevertFixedSites()
+ virtual ~RevertFixedSites()
{
}
@@ -418,7 +419,7 @@
/// HIDDEN Deep copy of a Migrator
virtual BaseOperator * clone() const
{
- return new SB_RevertFixedSites(*this);
+ return new RevertFixedSites(*this);
}
@@ -716,7 +717,7 @@
};
#endif
-
}
+}
#endif
Modified: trunk/src/sandbox.py
===================================================================
--- trunk/src/sandbox.py 2012-03-02 05:52:54 UTC (rev 4405)
+++ trunk/src/sandbox.py 2012-03-02 15:15:25 UTC (rev 4406)
@@ -34,9 +34,11 @@
version is no longer available.
"""
-__all__ = [
- 'sb_revertFixedSites',
- 'SB_RevertFixedSites',
+__all__ = []
+
+__long__all__ = [
+ 'revertFixedSites',
+ 'RevertFixedSites',
'MutSpaceMutator',
'MutSpaceSelector',
'MutSpaceRecombinator',
@@ -48,31 +50,39 @@
#
# These classes are available in binary modules but are not exposed
-# in simuPOP.
+# in simuPOP. Existing functions are only available in long module.
#
if simuOptions['Optimized']:
if simuOptions['AlleleType'] == 'short':
- from simuPOP_op import *
+ pass
elif simuOptions['AlleleType'] == 'long':
- from simuPOP_laop import *
+ from simuPOP_laop import sb_RevertFixedSites as RevertFixedSites
+ from simuPOP_laop import sb_MutSpaceMutator as MutSpaceMutator
+ from simuPOP_laop import sb_MutSpaceSelector as MutSpaceSelector
+ from simuPOP_laop import sb_MutSpaceRecombinator as MutSpaceRecombinator
+ __all__.extend(__long__all__)
elif simuOptions['AlleleType'] == 'binary':
- from simuPOP_baop import *
+ pass
elif simuOptions['AlleleType'] == 'mutant':
- from simuPOP_muop import *
+ pass
else:
- from simuPOP_op import *
+ pass
else:
if simuOptions['AlleleType'] == 'short':
- from simuPOP_std import *
+ pass
elif simuOptions['AlleleType'] == 'long':
- from simuPOP_la import *
+ from simuPOP_la import sb_RevertFixedSites as RevertFixedSites
+ from simuPOP_la import sb_MutSpaceMutator as MutSpaceMutator
+ from simuPOP_la import sb_MutSpaceSelector as MutSpaceSelector
+ from simuPOP_la import sb_MutSpaceRecombinator as MutSpaceRecombinator
+ __all__.extend(__long__all__)
elif simuOptions['AlleleType'] == 'binary':
- from simuPOP_ba import *
+ pass
elif simuOptions['AlleleType'] == 'mutant':
- from simuPOP_mu import *
+ pass
else:
- from simuPOP_std import *
+ pass
-def sb_revertFixedSites(pop):
+def revertFixedSites(pop):
'''Apply operator ``RevertFixedSites`` to ``pop``'''
- SB_RevertFixedSites().apply(pop)
+ RevertFixedSites().apply(pop)
Modified: trunk/src/simuPOP_common.i
===================================================================
--- trunk/src/simuPOP_common.i 2012-03-02 05:52:54 UTC (rev 4405)
+++ trunk/src/simuPOP_common.i 2012-03-02 15:15:25 UTC (rev 4406)
@@ -309,5 +309,11 @@
%include "qtrait.h"
%include "penetrance.h"
%include "pedigree.h"
+
+%rename(sb_RevertFixedSites) simuPOP::sandbox::RevertFixedSites;
+%rename(sb_MutSpaceSelector) simuPOP::sandbox::MutSpaceSelector;
+%rename(sb_MutSpaceMutator) simuPOP::sandbox::MutSpaceMutator;
+%rename(sb_MutSpaceRecombinator) simuPOP::sandbox::MutSpaceRecombinator;
+
%include "sandbox.h"
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2012-03-02 15:44:27
|
Revision: 4407
http://simupop.svn.sourceforge.net/simupop/?rev=4407&view=rev
Author: simupop
Date: 2012-03-02 15:44:20 +0000 (Fri, 02 Mar 2012)
Log Message:
-----------
Remove MutSpaceSelector from sandbox for mutant modual. The operator that is available in the core can be used for all mudules
Modified Paths:
--------------
trunk/src/sandbox.h
trunk/src/selector.cpp
trunk/src/selector.h
Modified: trunk/src/sandbox.h
===================================================================
--- trunk/src/sandbox.h 2012-03-02 15:15:25 UTC (rev 4406)
+++ trunk/src/sandbox.h 2012-03-02 15:44:20 UTC (rev 4407)
@@ -47,52 +47,6 @@
#ifdef MUTANTALLELE
-
-/** This operator looks into a population in mutational space and revert a mutant
- * to wildtype allele if it is fixed in the population. If a valid output is
- * specifieid, fixed alleles will be outputed with a leading generation number.
- */
-class RevertFixedSites : public BaseOperator
-{
-public:
- /** Create an operator to revert alleles at fixed loci from value 1 to 0.
- * Parameter \e subPops is ignored.
- */
- RevertFixedSites(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())
- : BaseOperator(output, begin, end, step, at, reps, subPops, infoFields)
- {
- }
-
-
- /// destructor
- virtual ~RevertFixedSites()
- {
- }
-
-
- /// HIDDEN Deep copy of a Migrator
- virtual BaseOperator * clone() const
- {
- return new RevertFixedSites(*this);
- }
-
-
- /// HIDDEN apply the Migrator to populaiton \e pop.
- virtual bool apply(Population & pop) const;
-
- /// HIDDEN
- string describe(bool format = true) const
- {
- (void)format; // avoid warning about unused parameter
- return "Revert fixed alleles to wildtype allele if it is fixed in the population.";
- }
-
-
-};
-
/** This selector assumes that alleles are mutant locations in the mutational
* space and assign fitness values to them according to a random distribution.
* The overall individual fitness is determined by either an additive, an
Modified: trunk/src/selector.cpp
===================================================================
--- trunk/src/selector.cpp 2012-03-02 15:15:25 UTC (rev 4406)
+++ trunk/src/selector.cpp 2012-03-02 15:44:20 UTC (rev 4407)
@@ -268,6 +268,246 @@
}
+#ifdef MUTANTALLELE
+double RandomFitnessSelector::indFitness(Population & /* pop */, Individual * ind) const
+{
+ if (m_mode == MULTIPLICATIVE) {
+ return randomSelMulFitnessExt(ind->genoBegin(), ind->genoEnd());
+ } else if (m_mode == ADDITIVE) {
+ if (m_additive)
+ return randomSelAddFitness(ind->genoBegin(), ind->genoEnd());
+ else
+ return randomSelAddFitnessExt(ind->genoBegin(), ind->genoEnd());
+ } else if (m_mode == EXPONENTIAL) {
+ if (m_additive)
+ return randomSelExpFitness(ind->genoBegin(), ind->genoEnd());
+ else
+ return randomSelExpFitnessExt(ind->genoBegin(), ind->genoEnd());
+ }
+ return 0;
}
+bool RandomFitnessSelector::apply(Population & pop) const
+{
+ m_newMutants.clear();
+ if (!BaseSelector::apply(pop))
+ return false;
+ // output NEW mutant...
+ if (!m_newMutants.empty() && !noOutput()) {
+ ostream & out = getOstream(pop.dict());
+ vectoru::const_iterator it = m_newMutants.begin();
+ vectoru::const_iterator it_end = m_newMutants.end();
+ for (; it != it_end; ++it) {
+ SelCoef s = m_selFactory[*it];
+ out << *it << '\t' << s.first << '\t' << s.second << '\n';
+ }
+ closeOstream();
+ }
+ return true;
+}
+
+
+RandomFitnessSelector::SelCoef RandomFitnessSelector::getFitnessValue(size_t mutant) const
+{
+ size_t sz = m_selDist.size();
+ double s = 0;
+ double h = 0.5;
+
+ if (sz == 0) {
+ // call a function
+ const pyFunc & func = m_selDist.func();
+ PyObject * res;
+ if (func.numArgs() == 0)
+ res = func("()");
+ else {
+ DBG_FAILIF(func.arg(0) != "loc", ValueError,
+ "Only parameter loc is accepted for this user-defined function.");
+ res = func("(i)", mutant);
+ }
+ if (PyNumber_Check(res)) {
+ s = PyFloat_AsDouble(res);
+ } else if (PySequence_Check(res)) {
+ size_t sz = PySequence_Size(res);
+ DBG_FAILIF(sz == 0, RuntimeError, "Function return an empty list.");
+ PyObject * item = PySequence_GetItem(res, 0);
+ s = PyFloat_AsDouble(item);
+ Py_DECREF(item);
+ if (sz > 1) {
+ item = PySequence_GetItem(res, 1);
+ h = PyFloat_AsDouble(item);
+ Py_DECREF(item);
+ }
+ }
+ Py_DECREF(res);
+ } else {
+ int mode = static_cast<int>(m_selDist[0]);
+ if (mode == CONSTANT) {
+ // constant
+ s = m_selDist[1];
+ if (m_selDist.size() > 2)
+ h = m_selDist[2];
+ } else {
+ // a gamma distribution
+ s = getRNG().randGamma(m_selDist[1], m_selDist[2]);
+ if (m_selDist.size() > 3)
+ h = m_selDist[3];
+ }
+ }
+ m_selFactory[mutant] = SelCoef(s, h);
+ m_newMutants.push_back(mutant);
+ if (m_additive && h != 0.5)
+ m_additive = false;
+ return SelCoef(s, h);
+}
+
+
+double RandomFitnessSelector::randomSelAddFitness(GenoIterator it, GenoIterator it_end) const
+{
+ double s = 0;
+ compressed_vector<size_t>::index_array_type::iterator index_it = it.getIndexIterator();
+ compressed_vector<size_t>::index_array_type::iterator index_it_end = it_end.getIndexIterator();
+ for (; index_it != index_it_end; ++index_it) {
+ SelMap::iterator sit = m_selFactory.find(static_cast<unsigned int>(*index_it));
+ if (sit == m_selFactory.end())
+ s += getFitnessValue(*index_it).first / 2.;
+ else
+ s += sit->second.first / 2;
+ }
+ return 1 - s > 0 ? 1 - s : 0;
+}
+
+
+double RandomFitnessSelector::randomSelExpFitness(GenoIterator it, GenoIterator it_end) const
+{
+ double s = 0;
+
+ compressed_vector<size_t>::index_array_type::iterator index_it = it.getIndexIterator();
+ compressed_vector<size_t>::index_array_type::iterator index_it_end = it_end.getIndexIterator();
+ for (; index_it != index_it_end; ++index_it) {
+ if (*index_it == 0)
+ continue;
+ SelMap::iterator sit = m_selFactory.find(static_cast<unsigned int>(*index_it));
+ if (sit == m_selFactory.end())
+ s += getFitnessValue(*index_it).first / 2.;
+ else
+ s += sit->second.first / 2;
+ }
+ return exp(-s);
+}
+
+
+double RandomFitnessSelector::randomSelMulFitnessExt(GenoIterator it, GenoIterator it_end) const
+{
+ MutCounter cnt;
+ compressed_vector<size_t>::index_array_type::iterator index_it = it.getIndexIterator();
+ compressed_vector<size_t>::index_array_type::iterator index_it_end = it_end.getIndexIterator();
+ for (; index_it != index_it_end; ++index_it) {
+ if (*index_it == 0)
+ continue;
+ MutCounter::iterator mit = cnt.find(*index_it);
+ if (mit == cnt.end())
+ cnt[*index_it] = 1;
+ else
+ ++mit->second;
+ }
+
+ double s = 1;
+ MutCounter::iterator mit = cnt.begin();
+ MutCounter::iterator mit_end = cnt.end();
+ for (; mit != mit_end; ++mit) {
+ SelMap::iterator sit = m_selFactory.find(mit->first);
+ if (sit == m_selFactory.end()) {
+ SelCoef sf = getFitnessValue(mit->first);
+ if (mit->second == 1)
+ s *= 1 - sf.first * sf.second;
+ else
+ s *= 1 - sf.first;
+ } else {
+ if (mit->second == 1)
+ s *= 1 - sit->second.first * sit->second.second;
+ else
+ s *= 1 - sit->second.first;
+ }
+ }
+ return s;
+}
+
+
+double RandomFitnessSelector::randomSelAddFitnessExt(GenoIterator it, GenoIterator it_end) const
+{
+ MutCounter cnt;
+ compressed_vector<size_t>::index_array_type::iterator index_it = it.getIndexIterator();
+ compressed_vector<size_t>::index_array_type::iterator index_it_end = it_end.getIndexIterator();
+ for (; index_it != index_it_end; ++index_it) {
+ if (*index_it == 0)
+ continue;
+ MutCounter::iterator mit = cnt.find(*index_it);
+ if (mit == cnt.end())
+ cnt[*index_it] = 1;
+ else
+ ++mit->second;
+ }
+
+ double s = 0;
+ MutCounter::iterator mit = cnt.begin();
+ MutCounter::iterator mit_end = cnt.end();
+ for (; mit != mit_end; ++mit) {
+ SelMap::iterator sit = m_selFactory.find(mit->first);
+ if (sit == m_selFactory.end()) {
+ SelCoef sf = getFitnessValue(mit->first);
+ if (mit->second == 1)
+ s += sf.first * sf.second;
+ else
+ s += sf.first;
+ } else {
+ if (mit->second == 1)
+ s += sit->second.first * sit->second.second;
+ else
+ s += sit->second.first;
+ }
+ }
+ return 1 - s > 0 ? 1 - s : 0;
+}
+
+
+double RandomFitnessSelector::randomSelExpFitnessExt(GenoIterator it, GenoIterator it_end) const
+{
+ MutCounter cnt;
+ compressed_vector<size_t>::index_array_type::iterator index_it = it.getIndexIterator();
+ compressed_vector<size_t>::index_array_type::iterator index_it_end = it_end.getIndexIterator();
+ for (; index_it != index_it_end; ++index_it) {
+ if (*index_it == 0)
+ continue;
+ MutCounter::iterator mit = cnt.find(*index_it);
+ if (mit == cnt.end())
+ cnt[*index_it] = 1;
+ else
+ ++mit->second;
+ }
+
+ double s = 0;
+ MutCounter::iterator mit = cnt.begin();
+ MutCounter::iterator mit_end = cnt.end();
+ for (; mit != mit_end; ++mit) {
+ SelMap::iterator sit = m_selFactory.find(mit->first);
+ if (sit == m_selFactory.end()) {
+ SelCoef sf = getFitnessValue(mit->first);
+ if (mit->second == 1)
+ s += sf.first * sf.second;
+ else
+ s += sf.first;
+ } else {
+ if (mit->second == 1)
+ s += sit->second.first * sit->second.second;
+ else
+ s += sit->second.first;
+ }
+ }
+ return exp(-s);
+}
+
+
+}
+
+
Modified: trunk/src/selector.h
===================================================================
--- trunk/src/selector.h 2012-03-02 15:15:25 UTC (rev 4406)
+++ trunk/src/selector.h 2012-03-02 15:44:20 UTC (rev 4407)
@@ -463,5 +463,121 @@
};
+
+/** This selector assumes that alleles are mutant locations in the mutational
+ * space and assign fitness values to them according to a random distribution.
+ * The overall individual fitness is determined by either an additive, an
+ * multiplicative or an exponential model.
+ */
+class RandomFitnessSelector : public BaseSelector
+{
+public:
+ /** Create a selector that assigns individual fitness values according to
+ * random fitness effects. \e selDist can be
+ * \li <tt>(CONSTANT, s, h)</tt> where s will be used for all mutants. The
+ * fitness value for genotypes AA, Aa and aa will be (1, 1-hs, 1-s).
+ * If h is unspecified, a default value h=0.5 (additive model) will
+ * be used.
+ * \li <tt>(GAMMA_DISTRIBUTION, theta, k, h</tt> where s follows a gamma
+ * distribution with scale parameter theta and shape parameter k.
+ * Fitness values for genotypes AA, Aa and aa will be 1, 1-hs and 1-s.
+ * A default value h=0.5 will be used if h is unspecified.
+ * \li a Python function, which will be called when selection coefficient
+ * of a new mutant is needed. This function should return a single
+ * value s (with default value h=0.5) or a sequence of (h, s). Mutant
+ * location will be passed to this function if it accepts a parameter
+ * \c loc. This allows the definition of site-specific selection
+ * coefficients.
+ * Individual fitness will be combined in \c ADDITIVE,
+ * \c MULTIPLICATIVE or \c EXPONENTIAL mode. (See \c MlSelector for
+ * details).
+ * If an output is given, mutants and their fitness values will be written
+ * to the output, in the form of 'mutant s h'.
+ */
+ RandomFitnessSelector(const floatListFunc & selDist, int mode = EXPONENTIAL,
+ 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 = stringList("fitness")) :
+ BaseSelector(output, begin, end, step, at, reps, subPops, infoFields),
+ m_selDist(selDist), m_mode(mode), m_selFactory(), m_additive(true)
+ {
+ if (m_selDist.size() == 0) {
+ DBG_FAILIF(!m_selDist.func().isValid(), ValueError,
+ "Please specify either a distribution with parameter or a function.");
+ } else if (static_cast<int>(m_selDist[0]) == CONSTANT) {
+ DBG_FAILIF(m_selDist.size() < 2, ValueError, "At least one parameter is needed for constant selection coefficient.");
+ } else if (static_cast<int>(m_selDist[0]) == GAMMA_DISTRIBUTION) {
+ DBG_FAILIF(m_selDist.size() < 3, ValueError, "At least two parameters are needed for gamma distribution.");
+ }
+ }
+
+
+ virtual ~RandomFitnessSelector()
+ {
+ }
+
+
+ /// HIDDEN Deep copy of a map selector
+ virtual BaseOperator * clone() const
+ {
+ return new RandomFitnessSelector(*this);
+ }
+
+
+ /** CPPONLY
+ * calculate/return the fitness value, currently assuming diploid
+ */
+ virtual double indFitness(Population & pop, Individual * ind) const;
+
+ /// HIDDEN
+ string describe(bool format = true) const
+ {
+ (void)format; // avoid warning about unused parameter
+ return "<simuPOP.RandomFitnessSelector>" ;
+ }
+
+
+ /// CPPONLY
+ bool apply(Population & pop) const;
+
+ typedef std::pair<double, double> SelCoef;
+
+private:
+ SelCoef getFitnessValue(size_t mutant) const;
+
+
+ double randomSelAddFitness(GenoIterator it, GenoIterator it_end) const;
+
+ double randomSelExpFitness(GenoIterator it, GenoIterator it_end) const;
+
+ // extended models does not assume additivity (h != 0.5)
+ double randomSelMulFitnessExt(GenoIterator it, GenoIterator it_end) const;
+
+ double randomSelAddFitnessExt(GenoIterator it, GenoIterator it_end) const;
+
+ double randomSelExpFitnessExt(GenoIterator it, GenoIterator it_end) const;
+
+private:
+ ///
+ floatListFunc m_selDist;
+
+ int m_mode;
+ ///
+#if TR1_SUPPORT == 0
+ typedef std::map<unsigned int, SelCoef> SelMap;
+ typedef std::map<unsigned int, int> MutCounter;
+#else
+ // this is faster than std::map
+ typedef std::tr1::unordered_map<size_t, SelCoef> SelMap;
+ typedef std::tr1::unordered_map<size_t, size_t> MutCounter;
+#endif
+ mutable SelMap m_selFactory;
+ mutable vectoru m_newMutants;
+ // whether or not all markers are additive.
+ mutable bool m_additive;
+};
+
+
}
#endif
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2012-03-02 15:48:25
|
Revision: 4409
http://simupop.svn.sourceforge.net/simupop/?rev=4409&view=rev
Author: simupop
Date: 2012-03-02 15:48:18 +0000 (Fri, 02 Mar 2012)
Log Message:
-----------
Move things out of sandbox
Modified Paths:
--------------
trunk/src/mutator.h
trunk/src/sandbox.cpp
trunk/src/sandbox.h
Modified: trunk/src/mutator.h
===================================================================
--- trunk/src/mutator.h 2012-03-02 15:44:54 UTC (rev 4408)
+++ trunk/src/mutator.h 2012-03-02 15:48:18 UTC (rev 4409)
@@ -700,5 +700,105 @@
+/** This is an infite site mutation model in mutational space. The alleles
+ * in the population is assumed to be locations of mutants. A mutation
+ * rate is given that mutate alleles in 'regions'. If number of mutants
+ * for an individual exceed the number of loci, 10 loci will be added
+ * to everyone in the population.
+ */
+class MutSpaceMutator : public BaseOperator
+{
+public:
+ /** This operator accepts a list of ranges which is the 'real range' of
+ * each chromosome. Mutation happens with muation rate \e rate and mutants
+ * will be recorded to the population (instead of alleles). By default,
+ * this mutator assumes an finite-allele model where all mutations are
+ * allowed and if a mutant (allele 1) is mutated, it will be mutated to
+ * allele 0 (back mutation). Alternatively (\e model = 2), an
+ * infinite-sites mutation model can be used where mutations can happen
+ * only at a new locus. Mutations happen at a locus with existing mutants
+ * will be moved to a random locus without existing mutant. A warning
+ * message will be printed if there is no vacant locus available. If a
+ * valid \e output is given, mutants will be outputted in the format of
+ * "gen mutant ind type" where type is 0 for forward (0->1), 1 for
+ * backward (1->0), 2 for relocated mutations, and 3 for ignored mutation
+ * because no vacent locus is available. The second mode has the
+ * advantage that all mutants in the simulated population can be traced
+ * to a single mutation event. If the regions are reasonably wide and
+ * mutation rates are low, these two mutation models should yield
+ * similar results.
+ */
+ MutSpaceMutator(double rate,
+ // FIXME: we should not have ranges, because ranges are just
+ // chromosomes, so, ranges=[1, 63000], [5000, 50000] shouldbe
+ // loci=[63000, 450000], lociPos=range(63000) + range(5000, 500000)
+ //
+ // The problem is that you need to optimize the storage of loci positions
+ // by saving ranges, not every locations in genoStru.h/cpp.
+ //
+ const intMatrix & ranges,
+ int model = 1,
+ 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()) :
+ BaseOperator(output, begin, end, step, at, reps, subPops, infoFields),
+ m_rate(rate), m_ranges(ranges), m_model(model)
+ {
+ const matrixi & rngs = m_ranges.elems();
+
+ for (size_t i = 0; i < rngs.size(); ++i) {
+ DBG_FAILIF(rngs[i].size() != 2, ValueError, "Ranges should have two elements");
+ for (size_t j = i + 1; j < rngs.size(); ++j) {
+ DBG_FAILIF(rngs[j].size() != 2, ValueError, "Ranges should have two elements");
+ if (i == j)
+ continue;
+ if (rngs[i][0] >= rngs[j][0] && rngs[i][0] <= rngs[j][1])
+ throw ValueError("Overlapping ranges are currently not supported because of potential conflict of mutant locations on different chromosomes.");
+ if (rngs[i][1] >= rngs[j][0] && rngs[i][1] <= rngs[j][1])
+ throw ValueError("Overlapping ranges are currently not supported because of potential conflict of mutant locations on different chromosomes.");
+ }
+ }
+#ifdef BINARYALLELE
+ DBG_FAILIF(true, ValueError, "This operator does not work in binary allele type.");
+#endif
+ }
+
+
+ /// destructor.
+ ~MutSpaceMutator()
+ {
+ }
+
+
+ virtual bool apply(Population & pop) const;
+
+ /// HIDDEN Deep copy of a \c MutSpaceMutator
+ virtual BaseOperator * clone() const
+ {
+ return new MutSpaceMutator(*this);
+ }
+
+
+ /// HIDDEN
+ string describe(bool format = true) const
+ {
+ (void)format; // avoid warning about unused parameter
+ return "<simuPOP.MutSpaceMutator>";
+ }
+
+
+private:
+ size_t locateVacantLocus(Population & pop, size_t beg, size_t end, std::set<size_t> & mutants) const;
+
+private:
+ const double m_rate;
+
+ const intMatrix m_ranges;
+
+ const int m_model;
+};
+
+
}
#endif
Modified: trunk/src/sandbox.cpp
===================================================================
--- trunk/src/sandbox.cpp 2012-03-02 15:44:54 UTC (rev 4408)
+++ trunk/src/sandbox.cpp 2012-03-02 15:48:18 UTC (rev 4409)
@@ -30,380 +30,6 @@
namespace sandbox {
#ifdef MUTANTALLELE
-
-
-size_t MutSpaceMutator::locateVacantLocus(Population & /* pop */, size_t beg, size_t end, std::set<size_t> & mutants) const
-{
- // FIXME: IGNORE this for now
- // this get a random locations
- size_t loc = getRNG().randInt(static_cast<ULONG>(end - beg)) + beg;
-
- // see if it exists in existing mutants.
- std::set<size_t>::iterator it= std::find(mutants.begin(), mutants.end(), loc);
- if (it == mutants.end())
- return loc;
- //
- // FIXME:
- //
- // assume mutants have 100 105 106 107
- // and loc = 106,
- // goes back and 104, and return
- //
- // look forward and backward
- size_t loc1 = loc + 1;
- std::set<size_t>::iterator it1(it);
- ++it1;
- for (; it1 != mutants.end() && loc1 != end; ++it1, ++loc1) {
- if (*it1 != loc1)
- return loc1;
- }
- size_t loc2 = loc - 1;
- std::set<size_t>::reverse_iterator it2(it);
- --it2;
- for (; it2 != mutants.rend() && loc2 != beg; --it2, --loc2) {
- if (*it2 != loc2)
- return loc2;
- }
- // still cannot find
- return 0;
-}
-
-
-bool MutSpaceMutator::apply(Population & pop) const
-{
- // FIXME:
- //
- const matrixi & ranges = m_ranges.elems();
- vectoru width(ranges.size());
-
- width[0] = ranges[0][1] - ranges[0][0];
- for (size_t i = 1; i < width.size(); ++i)
- width[i] = ranges[i][1] - ranges[i][0] + width[i - 1];
- //
- // FIXME: width[i] is the accumulative number of loci on each chromosome
- //
- // for three ranges(chromsomes) with loci: 500, 1000, 2000
- // width[i] = 0, 500, 1500, 35000
- //
- // ploidy with = totalNumLoci()
- //size_t ploidyWidth = width.back();
- //
- //size_t indWidth = pop.ploidy() * ploidyWidth;
-
- size_t indWidth = pop.genoSize();
- size_t ploidyWidth = pop.totNumLoci();
-
- ostream * out = NULL;
- if (!noOutput())
- out = &getOstream(pop.dict());
-
- // build a set of existing mutants
- std::set<size_t> mutants;
- //bool saturated = false;
-
- subPopList subPops = applicableSubPops(pop);
- subPopList::const_iterator sp = subPops.begin();
- subPopList::const_iterator spEnd = subPops.end();
- for (; sp != spEnd; ++sp) {
- DBG_FAILIF(sp->isVirtual(), ValueError, "This operator does not support virtual subpopulation.");
- for (size_t indIndex = 0; indIndex < pop.subPopSize(sp->subPop()); ++indIndex) {
- size_t loc = 0;
- while (true) {
- // using a geometric distribution to determine mutants
- loc += getRNG().randGeometric(m_rate);
- if (loc > indWidth)
- break;
- Individual & ind = pop.individual(indIndex);
- size_t p = (loc - 1) / ploidyWidth;
- // chromosome and position on chromosome?
- size_t mutLoc = (loc - 1) - p * ploidyWidth;
- size_t ch = 0;
- for (size_t reg = 0; reg < width.size(); ++reg) {
- if (mutLoc < width[reg]) {
- ch = reg;
- break;
- }
- }
- mutLoc += ranges[ch][0];
- if (ch > 0)
- mutLoc -= width[ch - 1];
-
- // FIXME: as the first step, ignore model 2, so you do not have to figure out
- // the function locateVacantLocus...
- /*
- // mutLoci is the location of the mutant, for individual ..., chromsome ...
- if (m_model == 2) {
- // under an infinite-site model
- if (saturated) {
- if (out)
- (*out) << pop.gen() << '\t' << mutLoc << '\t' << indIndex
- << "\t3\n";
- continue;
- }
- bool ok = false;
- //
- // if the first time
- if (mutants.empty()) {
- // first try our luck...
- // FIXME:
- // Here we are checing all genotypes (mutant), if mutLoc exists.
- // for the new module, yo uneed to check pop.mutBegin()??? (iterate through all
- // loci with non-zero allele)
- ok = find(pop.genoBegin(false), pop.genoEnd(false), ToAllele(mutLoc)) == pop.genoEnd(false);
- if (!ok) {
- std::set<size_t> existing(pop.genoBegin(false), pop.genoEnd(false));
- mutants.swap(existing);
- mutants.erase(0);
- saturated = mutants.size() == ploidyWidth;
- if (saturated)
- cerr << "Failed to introduce new mutants at generation " << pop.gen() << " because all loci have existing mutants." << endl;
- }
- }
- if (!ok && mutants.find(mutLoc) != mutants.end()) {
-
- size_t newLoc = locateVacantLocus(pop, ranges[ch][0], ranges[ch][1], mutants);
- // nothing is found
- if (out)
- (*out) << pop.gen() << '\t' << mutLoc << '\t' << indIndex
- << (newLoc == 0 ? "\t3\n" : "\t2\n");
- if (newLoc != 0)
- mutLoc = newLoc;
- else {
- cerr << "Failed to introduce a new mutant at generation " << pop.gen() << " because all loci have existing mutants." << endl;
- // ignore this mutation, and subsequent mutations...
- saturated = true;
- continue;
- }
- // if there is no existing mutant, new mutant is allowed
- }
- mutants.insert(mutLoc);
- }
- */
- ind.setAllele(1, mutLoc, int(p), int(ch));
- } // while
- } // each individual
- } // each subpopulation
- if (out)
- closeOstream();
- return true;
-}
-
-
-void MutSpaceRecombinator::transmitGenotype0(Population &, Population & offPop, const Individual & parent,
- size_t offIndex, int ploidy) const
-{
- size_t nCh = parent.numChrom();
-
- // count duplicates...
- for (size_t ch = 0; ch < parent.numChrom(); ++ch) {
- MutCounter cnt;
- vectoru alleles;
- alleles.reserve(parent.numLoci(ch));
- if (nCh == 1) {
- // this is faster... for a most common case
- GenoIterator it = parent.genoBegin();
- GenoIterator it_end = parent.genoEnd();
- for (; it != it_end; ++it) {
- if (*it == 0u)
- continue;
- MutCounter::iterator mit = cnt.find(*it);
- if (mit == cnt.end())
- cnt[*it] = 1;
- else
- ++mit->second;
- }
- } else {
- GenoIterator it = parent.genoBegin(0, ch);
- GenoIterator it_end = parent.genoEnd(0, ch);
- for (; it != it_end; ++it) {
- if (*it == 0u)
- break;
- MutCounter::iterator mit = cnt.find(*it);
- if (mit == cnt.end())
- cnt[*it] = 1;
- else
- ++mit->second;
- }
- it = parent.genoBegin(1, ch);
- it_end = parent.genoEnd(1, ch);
- for (; it != it_end; ++it) {
- if (*it == 0u)
- break;
- MutCounter::iterator mit = cnt.find(*it);
- if (mit == cnt.end())
- cnt[*it] = 1;
- else
- ++mit->second;
- }
- }
- // no valid allele
- if (cnt.empty()) {
- GenoIterator it = offPop.individual(offIndex).genoBegin(ploidy, ch);
- GenoIterator it_end = offPop.individual(offIndex).genoEnd(ploidy, ch);
- std::fill(it, it_end, 0);
- continue;
- }
- // keep 1 count with probability 0.5, keep 2 count with probability 1
- MutCounter::iterator mit = cnt.begin();
- MutCounter::iterator mit_end = cnt.end();
- for (; mit != mit_end; ++mit) {
- if (mit->second == 2 || getRNG().randBit())
- alleles.push_back(mit->first);
- }
- // not enough size
- if (alleles.size() + 1 > offPop.numLoci(ch)) {
- DBG_DO(DBG_TRANSMITTER, cerr << "Extending size of chromosome " << ch <<
- " to " << alleles.size() + 2 << endl);
- size_t sz = alleles.size() - offPop.numLoci(ch) + 2;
- vectorf added(sz);
- for (size_t j = 0; j < sz; ++j)
- added[j] = static_cast<double>(offPop.numLoci(ch) + j + 1);
- vectoru addedChrom(sz, ch);
- offPop.addLoci(addedChrom, added);
- }
- //
- GenoIterator it = offPop.individual(offIndex).genoBegin(ploidy, ch);
- GenoIterator it_end = offPop.individual(offIndex).genoEnd(ploidy, ch);
- for (size_t i = 0; i < alleles.size(); ++i, ++it)
- *it = ToAllele(alleles[i]);
- // fill the rest with 0.
- std::fill(it, it_end, 0);
- }
-}
-
-
-void MutSpaceRecombinator::transmitGenotype1(Population &, Population & offPop, const Individual & parent,
- size_t offIndex, int ploidy) const
-{
- const matrixi & ranges = m_ranges.elems();
-
- for (size_t ch = 0; ch < parent.numChrom(); ++ch) {
- size_t width = ranges[ch][1] - ranges[ch][0];
- size_t beg = 0;
- size_t end = getRNG().randGeometric(m_rate);
- int p = getRNG().randBit() ? 0 : 1;
- // no recombination
- if (end >= width) {
- copyChromosome(parent, p, offPop.individual(offIndex), ploidy, ch);
- continue;
- }
- // we are in trouble... get some properties of alleles to reduce comparison
- vectoru alleles;
- size_t minAllele[2];
- size_t maxAllele[2];
- size_t cnt[2];
- cnt[0] = 0;
- cnt[1] = 0;
- minAllele[0] = ranges[ch][1];
- minAllele[1] = ranges[ch][1];
- maxAllele[0] = ranges[ch][0];
- maxAllele[1] = ranges[ch][0];
- GenoIterator it = parent.genoBegin(0, ch);
- GenoIterator it_end = parent.genoEnd(0, ch);
- for (; it != it_end; ++it) {
- if (*it == 0u)
- break;
- ++cnt[0];
- if (*it < minAllele[0])
- minAllele[0] = *it;
- if (*it > maxAllele[0])
- maxAllele[0] = *it;
- }
- it = parent.genoBegin(1, ch);
- it_end = parent.genoEnd(1, ch);
- for (; it != it_end; ++it) {
- if (*it == 0u)
- break;
- ++cnt[1];
- if (*it < minAllele[1])
- minAllele[1] = *it;
- if (*it > maxAllele[1])
- maxAllele[1] = *it;
- }
- minAllele[0] -= ranges[ch][0];
- minAllele[1] -= ranges[ch][0];
- maxAllele[0] -= ranges[ch][0];
- maxAllele[1] -= ranges[ch][0];
- do {
- // copy piece
- // this algorithm is NOT efficient, but we count the rareness of recombination. :-)
- if (cnt[p] > 0 && end >= minAllele[p] && beg <= maxAllele[p]) {
- it = parent.genoBegin(p, ch);
- it_end = parent.genoEnd(p, ch);
- for (; it != it_end; ++it) {
- if (*it == 0u)
- break;
- if (*it >= beg + ranges[ch][0] && *it < end + ranges[ch][0]) {
- alleles.push_back(*it);
- --cnt[p];
- }
- }
- }
- // change ploidy
- p = (p + 1) % 2;
- // next step
- beg = end;
- end += getRNG().randGeometric(m_rate);
- } while (end < width && (cnt[0] > 0 || cnt[1] > 0));
- // last piece
- if (cnt[0] > 0 || cnt[1] > 0) {
- it = parent.genoBegin(p, ch);
- it_end = parent.genoEnd(p, ch);
- for (; it != it_end; ++it) {
- if (*it >= beg + static_cast<size_t>(ranges[ch][0]) &&
- *it < static_cast<size_t>(ranges[ch][1]))
- alleles.push_back(*it);
- }
- }
- // set alleles
- // not enough size
- if (alleles.size() + 1 > offPop.numLoci(ch)) {
- DBG_DO(DBG_TRANSMITTER, cerr << "Extending size of chromosome " << ch <<
- " to " << alleles.size() + 2 << endl);
- size_t sz = alleles.size() - offPop.numLoci(ch) + 2;
- vectorf added(sz);
- for (size_t j = 0; j < sz; ++j)
- added[j] = static_cast<double>(offPop.numLoci(ch) + j + 1);
- vectoru addedChrom(sz, ch);
- offPop.addLoci(addedChrom, added);
- }
- //
- it = offPop.individual(offIndex).genoBegin(ploidy, ch);
- it_end = offPop.individual(offIndex).genoEnd(ploidy, ch);
- for (size_t i = 0; i < alleles.size(); ++i, ++it)
- *it = ToAllele(alleles[i]);
- // fill the rest with 0.
- std::fill(it, it_end, 0);
- }
-}
-
-
-bool MutSpaceRecombinator::applyDuringMating(Population & pop, Population & offPop,
- RawIndIterator offspring,
- Individual * dad, Individual * mom) const
-{
- // if offspring does not belong to subPops, do nothing, but does not fail.
- if (!applicableToAllOffspring() && !applicableToOffspring(offPop, offspring))
- return true;
-
- initializeIfNeeded(*offspring);
-
- // standard genotype transmitter
- if (m_rate == 0) {
- for (int ch = 0; static_cast<size_t>(ch) < pop.numChrom(); ++ch) {
- copyChromosome(*mom, getRNG().randBit(), *offspring, 0, ch);
- copyChromosome(*dad, getRNG().randBit(), *offspring, 1, ch);
- }
- } else if (m_rate == 0.5) {
- transmitGenotype0(pop, offPop, *mom, offspring - offPop.rawIndBegin(), 0);
- transmitGenotype0(pop, offPop, *dad, offspring - offPop.rawIndBegin(), 1);
- } else {
- transmitGenotype1(pop, offPop, *mom, offspring - offPop.rawIndBegin(), 0);
- transmitGenotype1(pop, offPop, *dad, offspring - offPop.rawIndBegin(), 1);
- }
- return true;
-}
-
// ########################### THE OLD VERSION FOR LONG ALLELE TYPE ####
#else
Modified: trunk/src/sandbox.h
===================================================================
--- trunk/src/sandbox.h 2012-03-02 15:44:54 UTC (rev 4408)
+++ trunk/src/sandbox.h 2012-03-02 15:48:18 UTC (rev 4409)
@@ -47,301 +47,6 @@
#ifdef MUTANTALLELE
-/** This selector assumes that alleles are mutant locations in the mutational
- * space and assign fitness values to them according to a random distribution.
- * The overall individual fitness is determined by either an additive, an
- * multiplicative or an exponential model.
- */
-class MutSpaceSelector : public BaseSelector
-{
-public:
- /** Create a selector that assigns individual fitness values according to
- * random fitness effects. \e selDist can be
- * \li <tt>(CONSTANT, s, h)</tt> where s will be used for all mutants. The
- * fitness value for genotypes AA, Aa and aa will be (1, 1-hs, 1-s).
- * If h is unspecified, a default value h=0.5 (additive model) will
- * be used.
- * \li <tt>(GAMMA_DISTRIBUTION, theta, k, h</tt> where s follows a gamma
- * distribution with scale parameter theta and shape parameter k.
- * Fitness values for genotypes AA, Aa and aa will be 1, 1-hs and 1-s.
- * A default value h=0.5 will be used if h is unspecified.
- * \li a Python function, which will be called when selection coefficient
- * of a new mutant is needed. This function should return a single
- * value s (with default value h=0.5) or a sequence of (h, s). Mutant
- * location will be passed to this function if it accepts a parameter
- * \c loc. This allows the definition of site-specific selection
- * coefficients.
- * Individual fitness will be combined in \c ADDITIVE,
- * \c MULTIPLICATIVE or \c EXPONENTIAL mode. (See \c MlSelector for
- * details).
- * If an output is given, mutants and their fitness values will be written
- * to the output, in the form of 'mutant s h'.
- */
- MutSpaceSelector(const floatListFunc & selDist, int mode = EXPONENTIAL,
- 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 = stringList("fitness")) :
- BaseSelector(output, begin, end, step, at, reps, subPops, infoFields),
- m_selDist(selDist), m_mode(mode), m_selFactory(), m_additive(true)
- {
- if (m_selDist.size() == 0) {
- DBG_FAILIF(!m_selDist.func().isValid(), ValueError,
- "Please specify either a distribution with parameter or a function.");
- } else if (static_cast<int>(m_selDist[0]) == CONSTANT) {
- DBG_FAILIF(m_selDist.size() < 2, ValueError, "At least one parameter is needed for constant selection coefficient.");
- } else if (static_cast<int>(m_selDist[0]) == GAMMA_DISTRIBUTION) {
- DBG_FAILIF(m_selDist.size() < 3, ValueError, "At least two parameters are needed for gamma distribution.");
- }
- }
-
-
- virtual ~MutSpaceSelector()
- {
- }
-
-
- /// HIDDEN Deep copy of a map selector
- virtual BaseOperator * clone() const
- {
- return new MutSpaceSelector(*this);
- }
-
-
- /** CPPONLY
- * calculate/return the fitness value, currently assuming diploid
- */
- virtual double indFitness(Population & pop, Individual * ind) const;
-
- /// HIDDEN
- string describe(bool format = true) const
- {
- (void)format; // avoid warning about unused parameter
- return "<simuPOP.MutSpaceSelector>" ;
- }
-
-
- /// CPPONLY
- bool apply(Population & pop) const;
-
- typedef std::pair<double, double> SelCoef;
-
-private:
- SelCoef getFitnessValue(size_t mutant) const;
-
-
- double randomSelAddFitness(GenoIterator it, GenoIterator it_end) const;
-
- double randomSelExpFitness(GenoIterator it, GenoIterator it_end) const;
-
- // extended models does not assume additivity (h != 0.5)
- double randomSelMulFitnessExt(GenoIterator it, GenoIterator it_end) const;
-
- double randomSelAddFitnessExt(GenoIterator it, GenoIterator it_end) const;
-
- double randomSelExpFitnessExt(GenoIterator it, GenoIterator it_end) const;
-
-private:
- ///
- floatListFunc m_selDist;
-
- int m_mode;
- ///
-#if TR1_SUPPORT == 0
- typedef std::map<unsigned int, SelCoef> SelMap;
- typedef std::map<unsigned int, int> MutCounter;
-#else
- // this is faster than std::map
- typedef std::tr1::unordered_map<size_t, SelCoef> SelMap;
- typedef std::tr1::unordered_map<size_t, size_t> MutCounter;
-#endif
- mutable SelMap m_selFactory;
- mutable vectoru m_newMutants;
- // whether or not all markers are additive.
- mutable bool m_additive;
-};
-
-
-/** This is an infite site mutation model in mutational space. The alleles
- * in the population is assumed to be locations of mutants. A mutation
- * rate is given that mutate alleles in 'regions'. If number of mutants
- * for an individual exceed the number of loci, 10 loci will be added
- * to everyone in the population.
- */
-class MutSpaceMutator : public BaseOperator
-{
-public:
- /** This operator accepts a list of ranges which is the 'real range' of
- * each chromosome. Mutation happens with muation rate \e rate and mutants
- * will be recorded to the population (instead of alleles). By default,
- * this mutator assumes an finite-allele model where all mutations are
- * allowed and if a mutant (allele 1) is mutated, it will be mutated to
- * allele 0 (back mutation). Alternatively (\e model = 2), an
- * infinite-sites mutation model can be used where mutations can happen
- * only at a new locus. Mutations happen at a locus with existing mutants
- * will be moved to a random locus without existing mutant. A warning
- * message will be printed if there is no vacant locus available. If a
- * valid \e output is given, mutants will be outputted in the format of
- * "gen mutant ind type" where type is 0 for forward (0->1), 1 for
- * backward (1->0), 2 for relocated mutations, and 3 for ignored mutation
- * because no vacent locus is available. The second mode has the
- * advantage that all mutants in the simulated population can be traced
- * to a single mutation event. If the regions are reasonably wide and
- * mutation rates are low, these two mutation models should yield
- * similar results.
- */
- MutSpaceMutator(double rate,
- // FIXME: we should not have ranges, because ranges are just
- // chromosomes, so, ranges=[1, 63000], [5000, 50000] shouldbe
- // loci=[63000, 450000], lociPos=range(63000) + range(5000, 500000)
- //
- // The problem is that you need to optimize the storage of loci positions
- // by saving ranges, not every locations in genoStru.h/cpp.
- //
- const intMatrix & ranges,
- int model = 1,
- 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()) :
- BaseOperator(output, begin, end, step, at, reps, subPops, infoFields),
- m_rate(rate), m_ranges(ranges), m_model(model)
- {
- const matrixi & rngs = m_ranges.elems();
-
- for (size_t i = 0; i < rngs.size(); ++i) {
- DBG_FAILIF(rngs[i].size() != 2, ValueError, "Ranges should have two elements");
- for (size_t j = i + 1; j < rngs.size(); ++j) {
- DBG_FAILIF(rngs[j].size() != 2, ValueError, "Ranges should have two elements");
- if (i == j)
- continue;
- if (rngs[i][0] >= rngs[j][0] && rngs[i][0] <= rngs[j][1])
- throw ValueError("Overlapping ranges are currently not supported because of potential conflict of mutant locations on different chromosomes.");
- if (rngs[i][1] >= rngs[j][0] && rngs[i][1] <= rngs[j][1])
- throw ValueError("Overlapping ranges are currently not supported because of potential conflict of mutant locations on different chromosomes.");
- }
- }
-#ifdef BINARYALLELE
- DBG_FAILIF(true, ValueError, "This operator does not work in binary allele type.");
-#endif
- }
-
-
- /// destructor.
- ~MutSpaceMutator()
- {
- }
-
-
- virtual bool apply(Population & pop) const;
-
- /// HIDDEN Deep copy of a \c MutSpaceMutator
- virtual BaseOperator * clone() const
- {
- return new MutSpaceMutator(*this);
- }
-
-
- /// HIDDEN
- string describe(bool format = true) const
- {
- (void)format; // avoid warning about unused parameter
- return "<simuPOP.MutSpaceMutator>";
- }
-
-
-private:
- size_t locateVacantLocus(Population & pop, size_t beg, size_t end, std::set<size_t> & mutants) const;
-
-private:
- const double m_rate;
-
- const intMatrix m_ranges;
-
- const int m_model;
-};
-
-
-/** This during mating operator recombine chromosomes, which records mutant
- * locations, using a fixed recombination rate (per base pair).
- */
- // FIXME:
- // This should jsut be the usual recombinator that is implemented for
- // mutational model. so this should be removed and be replaced with regular
- // recombinator.
- //
-class MutSpaceRecombinator : public GenoTransmitter
-{
-public:
- /** Create a Recombinator (a mendelian genotype transmitter with
- * recombination and gene conversion) that passes genotypes from parents
- * (or a parent in case of self-fertilization) to offspring. A
- * recombination \e rate in the unit of base pair is needed.
- */
- MutSpaceRecombinator(double rate, const intMatrix & ranges,
- 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())
- : GenoTransmitter(output, begin, end, step, at, reps, subPops, infoFields),
- m_rate(rate), m_ranges(ranges)
- {
- DBG_FAILIF(rate > 0.5 || rate < 0, ValueError, "Recombination rate should be between 0 and 0.5");
-#ifdef BINARYALLELE
- DBG_FAILIF(true, ValueError, "This operator does not work in binary allele type.");
-#endif
- }
-
-
- /// HIDDEN Deep copy of a Recombinator
- virtual BaseOperator * clone() const
- {
- return new MutSpaceRecombinator(*this);
- }
-
-
- virtual ~MutSpaceRecombinator()
- {
- }
-
-
- /// HIDDEN
- string describe(bool format = true) const
- {
- (void)format; // avoid warning about unused parameter
- return "<simuPOP.MutSpaceRecombinator>";
- }
-
-
- /** CPPONLY
- * Apply the Recombinator during mating
- */
- virtual bool applyDuringMating(Population & pop, Population & offPop,
- RawIndIterator offspring,
- Individual * dad, Individual * mom) const;
-
-private:
-#if TR1_SUPPORT == 0
- typedef std::map<unsigned int, int> MutCounter;
-#else
- // this is faster than std::map
- typedef std::tr1::unordered_map<unsigned int, int> MutCounter;
-#endif
- // use when m_rate = 0.5
- void transmitGenotype0(Population & pop, Population & offPop, const Individual & parent,
- size_t offIndex, int ploidy) const;
-
- // use when m_rate < 1e-4
- void transmitGenotype1(Population & pop, Population & offPop, const Individual & parent,
- size_t offIndex, int ploidy) const;
-
-private:
- /// recombination rate
- const double m_rate;
- const intMatrix m_ranges;
-};
-
-
/// ######################################
#else
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2012-03-02 15:53:00
|
Revision: 4410
http://simupop.svn.sourceforge.net/simupop/?rev=4410&view=rev
Author: simupop
Date: 2012-03-02 15:52:50 +0000 (Fri, 02 Mar 2012)
Log Message:
-----------
Clean up sandbox so that it only has the old functions for long allele modules
Modified Paths:
--------------
trunk/src/mutator.cpp
trunk/src/sandbox.cpp
trunk/src/sandbox.h
Modified: trunk/src/mutator.cpp
===================================================================
--- trunk/src/mutator.cpp 2012-03-02 15:48:18 UTC (rev 4409)
+++ trunk/src/mutator.cpp 2012-03-02 15:52:50 UTC (rev 4410)
@@ -499,4 +499,162 @@
}
+size_t MutSpaceMutator::locateVacantLocus(Population & /* pop */, size_t beg, size_t end, std::set<size_t> & mutants) const
+{
+ // FIXME: IGNORE this for now
+ // this get a random locations
+ size_t loc = getRNG().randInt(static_cast<ULONG>(end - beg)) + beg;
+
+ // see if it exists in existing mutants.
+ std::set<size_t>::iterator it= std::find(mutants.begin(), mutants.end(), loc);
+ if (it == mutants.end())
+ return loc;
+ //
+ // FIXME:
+ //
+ // assume mutants have 100 105 106 107
+ // and loc = 106,
+ // goes back and 104, and return
+ //
+ // look forward and backward
+ size_t loc1 = loc + 1;
+ std::set<size_t>::iterator it1(it);
+ ++it1;
+ for (; it1 != mutants.end() && loc1 != end; ++it1, ++loc1) {
+ if (*it1 != loc1)
+ return loc1;
+ }
+ size_t loc2 = loc - 1;
+ std::set<size_t>::reverse_iterator it2(it);
+ --it2;
+ for (; it2 != mutants.rend() && loc2 != beg; --it2, --loc2) {
+ if (*it2 != loc2)
+ return loc2;
+ }
+ // still cannot find
+ return 0;
}
+
+
+bool MutSpaceMutator::apply(Population & pop) const
+{
+ // FIXME:
+ //
+ const matrixi & ranges = m_ranges.elems();
+ vectoru width(ranges.size());
+
+ width[0] = ranges[0][1] - ranges[0][0];
+ for (size_t i = 1; i < width.size(); ++i)
+ width[i] = ranges[i][1] - ranges[i][0] + width[i - 1];
+ //
+ // FIXME: width[i] is the accumulative number of loci on each chromosome
+ //
+ // for three ranges(chromsomes) with loci: 500, 1000, 2000
+ // width[i] = 0, 500, 1500, 35000
+ //
+ // ploidy with = totalNumLoci()
+ //size_t ploidyWidth = width.back();
+ //
+ //size_t indWidth = pop.ploidy() * ploidyWidth;
+
+ size_t indWidth = pop.genoSize();
+ size_t ploidyWidth = pop.totNumLoci();
+
+ ostream * out = NULL;
+ if (!noOutput())
+ out = &getOstream(pop.dict());
+
+ // build a set of existing mutants
+ std::set<size_t> mutants;
+ //bool saturated = false;
+
+ subPopList subPops = applicableSubPops(pop);
+ subPopList::const_iterator sp = subPops.begin();
+ subPopList::const_iterator spEnd = subPops.end();
+ for (; sp != spEnd; ++sp) {
+ DBG_FAILIF(sp->isVirtual(), ValueError, "This operator does not support virtual subpopulation.");
+ for (size_t indIndex = 0; indIndex < pop.subPopSize(sp->subPop()); ++indIndex) {
+ size_t loc = 0;
+ while (true) {
+ // using a geometric distribution to determine mutants
+ loc += getRNG().randGeometric(m_rate);
+ if (loc > indWidth)
+ break;
+ Individual & ind = pop.individual(indIndex);
+ size_t p = (loc - 1) / ploidyWidth;
+ // chromosome and position on chromosome?
+ size_t mutLoc = (loc - 1) - p * ploidyWidth;
+ size_t ch = 0;
+ for (size_t reg = 0; reg < width.size(); ++reg) {
+ if (mutLoc < width[reg]) {
+ ch = reg;
+ break;
+ }
+ }
+ mutLoc += ranges[ch][0];
+ if (ch > 0)
+ mutLoc -= width[ch - 1];
+
+ // FIXME: as the first step, ignore model 2, so you do not have to figure out
+ // the function locateVacantLocus...
+ /*
+ // mutLoci is the location of the mutant, for individual ..., chromsome ...
+ if (m_model == 2) {
+ // under an infinite-site model
+ if (saturated) {
+ if (out)
+ (*out) << pop.gen() << '\t' << mutLoc << '\t' << indIndex
+ << "\t3\n";
+ continue;
+ }
+ bool ok = false;
+ //
+ // if the first time
+ if (mutants.empty()) {
+ // first try our luck...
+ // FIXME:
+ // Here we are checing all genotypes (mutant), if mutLoc exists.
+ // for the new module, yo uneed to check pop.mutBegin()??? (iterate through all
+ // loci with non-zero allele)
+ ok = find(pop.genoBegin(false), pop.genoEnd(false), ToAllele(mutLoc)) == pop.genoEnd(false);
+ if (!ok) {
+ std::set<size_t> existing(pop.genoBegin(false), pop.genoEnd(false));
+ mutants.swap(existing);
+ mutants.erase(0);
+ saturated = mutants.size() == ploidyWidth;
+ if (saturated)
+ cerr << "Failed to introduce new mutants at generation " << pop.gen() << " because all loci have existing mutants." << endl;
+ }
+ }
+ if (!ok && mutants.find(mutLoc) != mutants.end()) {
+
+ size_t newLoc = locateVacantLocus(pop, ranges[ch][0], ranges[ch][1], mutants);
+ // nothing is found
+ if (out)
+ (*out) << pop.gen() << '\t' << mutLoc << '\t' << indIndex
+ << (newLoc == 0 ? "\t3\n" : "\t2\n");
+ if (newLoc != 0)
+ mutLoc = newLoc;
+ else {
+ cerr << "Failed to introduce a new mutant at generation " << pop.gen() << " because all loci have existing mutants." << endl;
+ // ignore this mutation, and subsequent mutations...
+ saturated = true;
+ continue;
+ }
+ // if there is no existing mutant, new mutant is allowed
+ }
+ mutants.insert(mutLoc);
+ }
+ */
+ ind.setAllele(1, mutLoc, int(p), int(ch));
+ } // while
+ } // each individual
+ } // each subpopulation
+ if (out)
+ closeOstream();
+ return true;
+}
+
+
+
+}
Modified: trunk/src/sandbox.cpp
===================================================================
--- trunk/src/sandbox.cpp 2012-03-02 15:48:18 UTC (rev 4409)
+++ trunk/src/sandbox.cpp 2012-03-02 15:52:50 UTC (rev 4410)
@@ -29,9 +29,7 @@
namespace simuPOP {
namespace sandbox {
-#ifdef MUTANTALLELE
-// ########################### THE OLD VERSION FOR LONG ALLELE TYPE ####
-#else
+#ifdef LONGALLELE
bool RevertFixedSites::apply(Population & pop) const
{
@@ -351,9 +349,6 @@
bool MutSpaceMutator::apply(Population & pop) const
{
- #ifdef BINARYALLELE
- (void)pop; // avoid warning about unused parameter
- #else
const matrixi & ranges = m_ranges.elems();
vectoru width(ranges.size());
@@ -486,7 +481,6 @@
} // each subpopulation
if (out)
closeOstream();
- #endif
return true;
}
@@ -494,13 +488,6 @@
void MutSpaceRecombinator::transmitGenotype0(Population & pop, Population & offPop, const Individual & parent,
size_t offIndex, int ploidy) const
{
- #ifdef BINARYALLELE
- (void)offPop; // avoid warning about unused parameter
- (void)parent; // avoid warning about unused parameter
- (void)offIndex; // avoid warning about unused parameter
- (void)ploidy; // avoid warning about unused parameter
- (void)pop; // avoid warning about unused parameter
- #else
size_t nCh = parent.numChrom();
// count duplicates...
@@ -579,20 +566,12 @@
// fill the rest with 0.
std::fill(it, it_end, 0);
}
- #endif
}
void MutSpaceRecombinator::transmitGenotype1(Population & pop, Population & offPop, const Individual & parent,
size_t offIndex, int ploidy) const
{
- #ifdef BINARYALLELE
- (void)offPop; // avoid warning about unused parameter
- (void)parent; // avoid warning about unused parameter
- (void)offIndex; // avoid warning about unused parameter
- (void)ploidy; // avoid warning about unused parameter
- (void)pop; // avoid warning about unused parameter
- #else
const matrixi & ranges = m_ranges.elems();
for (size_t ch = 0; ch < parent.numChrom(); ++ch) {
@@ -694,7 +673,6 @@
// fill the rest with 0.
std::fill(it, it_end, 0);
}
- #endif
}
Modified: trunk/src/sandbox.h
===================================================================
--- trunk/src/sandbox.h 2012-03-02 15:48:18 UTC (rev 4409)
+++ trunk/src/sandbox.h 2012-03-02 15:52:50 UTC (rev 4410)
@@ -45,11 +45,8 @@
namespace simuPOP {
namespace sandbox {
-#ifdef MUTANTALLELE
+#ifdef LONGALLELE
-/// ######################################
-#else
-
/** This operator looks into a population in mutational space and revert a mutant
* to wildtype allele if it is fixed in the population. If a valid output is
* specifieid, fixed alleles will be outputed with a leading generation number.
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2012-03-02 16:02:07
|
Revision: 4411
http://simupop.svn.sourceforge.net/simupop/?rev=4411&view=rev
Author: simupop
Date: 2012-03-02 16:01:56 +0000 (Fri, 02 Mar 2012)
Log Message:
-----------
Rename MutSpaceMutator to FiniteSitesMutator
Modified Paths:
--------------
trunk/src/__init__.py
trunk/src/mutator.cpp
trunk/src/mutator.h
trunk/src/selector.cpp
trunk/src/selector.h
Modified: trunk/src/__init__.py
===================================================================
--- trunk/src/__init__.py 2012-03-02 15:52:50 UTC (rev 4410)
+++ trunk/src/__init__.py 2012-03-02 16:01:56 UTC (rev 4411)
@@ -175,11 +175,13 @@
'ContextMutator',
'KAlleleMutator',
'RevertFixedSites',
+ 'FiniteSitesMutator',
#
'MapSelector',
'MaSelector',
'MlSelector',
'PySelector',
+ 'RandomFitnessSelector',
#
'MaPenetrance',
'MapPenetrance',
Modified: trunk/src/mutator.cpp
===================================================================
--- trunk/src/mutator.cpp 2012-03-02 15:52:50 UTC (rev 4410)
+++ trunk/src/mutator.cpp 2012-03-02 16:01:56 UTC (rev 4411)
@@ -499,7 +499,7 @@
}
-size_t MutSpaceMutator::locateVacantLocus(Population & /* pop */, size_t beg, size_t end, std::set<size_t> & mutants) const
+size_t FiniteSitesMutator::locateVacantLocus(Population & /* pop */, size_t beg, size_t end, std::set<size_t> & mutants) const
{
// FIXME: IGNORE this for now
// this get a random locations
@@ -536,7 +536,7 @@
}
-bool MutSpaceMutator::apply(Population & pop) const
+bool FiniteSitesMutator::apply(Population & pop) const
{
// FIXME:
//
Modified: trunk/src/mutator.h
===================================================================
--- trunk/src/mutator.h 2012-03-02 15:52:50 UTC (rev 4410)
+++ trunk/src/mutator.h 2012-03-02 16:01:56 UTC (rev 4411)
@@ -706,7 +706,7 @@
* for an individual exceed the number of loci, 10 loci will be added
* to everyone in the population.
*/
-class MutSpaceMutator : public BaseOperator
+class FiniteSitesMutator : public BaseOperator
{
public:
/** This operator accepts a list of ranges which is the 'real range' of
@@ -728,7 +728,7 @@
* mutation rates are low, these two mutation models should yield
* similar results.
*/
- MutSpaceMutator(double rate,
+ FiniteSitesMutator(double rate,
// FIXME: we should not have ranges, because ranges are just
// chromosomes, so, ranges=[1, 63000], [5000, 50000] shouldbe
// loci=[63000, 450000], lociPos=range(63000) + range(5000, 500000)
@@ -766,17 +766,17 @@
/// destructor.
- ~MutSpaceMutator()
+ ~FiniteSitesMutator()
{
}
virtual bool apply(Population & pop) const;
- /// HIDDEN Deep copy of a \c MutSpaceMutator
+ /// HIDDEN Deep copy of a \c FiniteSitesMutator
virtual BaseOperator * clone() const
{
- return new MutSpaceMutator(*this);
+ return new FiniteSitesMutator(*this);
}
@@ -784,7 +784,7 @@
string describe(bool format = true) const
{
(void)format; // avoid warning about unused parameter
- return "<simuPOP.MutSpaceMutator>";
+ return "<simuPOP.FiniteSitesMutator>";
}
Modified: trunk/src/selector.cpp
===================================================================
--- trunk/src/selector.cpp 2012-03-02 15:52:50 UTC (rev 4410)
+++ trunk/src/selector.cpp 2012-03-02 16:01:56 UTC (rev 4411)
@@ -268,7 +268,6 @@
}
-#ifdef MUTANTALLELE
double RandomFitnessSelector::indFitness(Population & /* pop */, Individual * ind) const
{
if (m_mode == MULTIPLICATIVE) {
Modified: trunk/src/selector.h
===================================================================
--- trunk/src/selector.h 2012-03-02 15:52:50 UTC (rev 4410)
+++ trunk/src/selector.h 2012-03-02 16:01:56 UTC (rev 4411)
@@ -36,6 +36,16 @@
#include <numeric>
using std::min;
+#if TR1_SUPPORT == 0
+# include <map>
+#elif TR1_SUPPORT == 1
+# include <unordered_map>
+#else
+# include <tr1/unordered_map>
+#endif
+
+#include <set>
+
namespace simuPOP {
/** This class is the base class to all selectors, namely operators that
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: Bo P. <be...@gm...> - 2012-03-02 16:08:29
|
Hi, Panitee, Will these patches 1. Sandbox only works for the long allele module, with functions and operators untouched for the existing srv.py to work. 2. RevertFixedSites are now available for all modules. You should use this operator from the core for your new script. 3. MutSpaceMutator is renamed to FiniteSitesMutator in the core. It works only for mutant module right now. You should use this one for your new script. 4. MutSpaceSelector is renamed to RandomFitnessSelector in the core. It works only for mutant module right now. You should use this one for your new script. 5. MutSpaceRecombinator should not be used because the existing Recombinator should work. Please, 1. Update your srv.py to not using anything from the sandbox. 2. Update FiniteSitesMutator and RandomFitnessSelector to make them work for all modules (binary, short, long). Thanks, Bo On Fri, Mar 2, 2012 at 10:01 AM, <si...@us...> wrote: > Revision: 4411 > http://simupop.svn.sourceforge.net/simupop/?rev=4411&view=rev > Author: simupop > Date: 2012-03-02 16:01:56 +0000 (Fri, 02 Mar 2012) > Log Message: > ----------- > Rename MutSpaceMutator to FiniteSitesMutator > > Modified Paths: > -------------- > trunk/src/__init__.py > trunk/src/mutator.cpp > trunk/src/mutator.h > trunk/src/selector.cpp > trunk/src/selector.h > > Modified: trunk/src/__init__.py > =================================================================== > --- trunk/src/__init__.py 2012-03-02 15:52:50 UTC (rev 4410) > +++ trunk/src/__init__.py 2012-03-02 16:01:56 UTC (rev 4411) > @@ -175,11 +175,13 @@ > 'ContextMutator', > 'KAlleleMutator', > 'RevertFixedSites', > + 'FiniteSitesMutator', > # > 'MapSelector', > 'MaSelector', > 'MlSelector', > 'PySelector', > + 'RandomFitnessSelector', > # > 'MaPenetrance', > 'MapPenetrance', > > Modified: trunk/src/mutator.cpp > =================================================================== > --- trunk/src/mutator.cpp 2012-03-02 15:52:50 UTC (rev 4410) > +++ trunk/src/mutator.cpp 2012-03-02 16:01:56 UTC (rev 4411) > @@ -499,7 +499,7 @@ > } > > > -size_t MutSpaceMutator::locateVacantLocus(Population & /* pop */, size_t beg, size_t end, std::set<size_t> & mutants) const > +size_t FiniteSitesMutator::locateVacantLocus(Population & /* pop */, size_t beg, size_t end, std::set<size_t> & mutants) const > { > // FIXME: IGNORE this for now > // this get a random locations > @@ -536,7 +536,7 @@ > } > > > -bool MutSpaceMutator::apply(Population & pop) const > +bool FiniteSitesMutator::apply(Population & pop) const > { > // FIXME: > // > > Modified: trunk/src/mutator.h > =================================================================== > --- trunk/src/mutator.h 2012-03-02 15:52:50 UTC (rev 4410) > +++ trunk/src/mutator.h 2012-03-02 16:01:56 UTC (rev 4411) > @@ -706,7 +706,7 @@ > * for an individual exceed the number of loci, 10 loci will be added > * to everyone in the population. > */ > -class MutSpaceMutator : public BaseOperator > +class FiniteSitesMutator : public BaseOperator > { > public: > /** This operator accepts a list of ranges which is the 'real range' of > @@ -728,7 +728,7 @@ > * mutation rates are low, these two mutation models should yield > * similar results. > */ > - MutSpaceMutator(double rate, > + FiniteSitesMutator(double rate, > // FIXME: we should not have ranges, because ranges are just > // chromosomes, so, ranges=[1, 63000], [5000, 50000] shouldbe > // loci=[63000, 450000], lociPos=range(63000) + range(5000, 500000) > @@ -766,17 +766,17 @@ > > > /// destructor. > - ~MutSpaceMutator() > + ~FiniteSitesMutator() > { > } > > > virtual bool apply(Population & pop) const; > > - /// HIDDEN Deep copy of a \c MutSpaceMutator > + /// HIDDEN Deep copy of a \c FiniteSitesMutator > virtual BaseOperator * clone() const > { > - return new MutSpaceMutator(*this); > + return new FiniteSitesMutator(*this); > } > > > @@ -784,7 +784,7 @@ > string describe(bool format = true) const > { > (void)format; // avoid warning about unused parameter > - return "<simuPOP.MutSpaceMutator>"; > + return "<simuPOP.FiniteSitesMutator>"; > } > > > > Modified: trunk/src/selector.cpp > =================================================================== > --- trunk/src/selector.cpp 2012-03-02 15:52:50 UTC (rev 4410) > +++ trunk/src/selector.cpp 2012-03-02 16:01:56 UTC (rev 4411) > @@ -268,7 +268,6 @@ > } > > > -#ifdef MUTANTALLELE > double RandomFitnessSelector::indFitness(Population & /* pop */, Individual * ind) const > { > if (m_mode == MULTIPLICATIVE) { > > Modified: trunk/src/selector.h > =================================================================== > --- trunk/src/selector.h 2012-03-02 15:52:50 UTC (rev 4410) > +++ trunk/src/selector.h 2012-03-02 16:01:56 UTC (rev 4411) > @@ -36,6 +36,16 @@ > #include <numeric> > using std::min; > > +#if TR1_SUPPORT == 0 > +# include <map> > +#elif TR1_SUPPORT == 1 > +# include <unordered_map> > +#else > +# include <tr1/unordered_map> > +#endif > + > +#include <set> > + > namespace simuPOP { > > /** This class is the base class to all selectors, namely operators that > > This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. > > > ------------------------------------------------------------------------------ > Virtualization & Cloud Management Using Capacity Planning > Cloud computing makes use of virtualization - but cloud computing > also focuses on allowing computing to be delivered as a service. > http://www.accelacomm.com/jaw/sfnl/114/51521223/ > _______________________________________________ > simuPOP-devel mailing list > sim...@li... > https://lists.sourceforge.net/lists/listinfo/simupop-devel |
|
From: <ti...@us...> - 2012-03-02 21:22:23
|
Revision: 4418
http://simupop.svn.sourceforge.net/simupop/?rev=4418&view=rev
Author: tiwcpe8
Date: 2012-03-02 21:22:16 +0000 (Fri, 02 Mar 2012)
Log Message:
-----------
Add FixedSites in NumofSegSites class
Modified Paths:
--------------
trunk/src/stator.cpp
trunk/src/stator.h
Modified: trunk/src/stator.cpp
===================================================================
--- trunk/src/stator.cpp 2012-03-02 19:54:08 UTC (rev 4417)
+++ trunk/src/stator.cpp 2012-03-02 21:22:16 UTC (rev 4418)
@@ -670,9 +670,10 @@
: m_loci(loci), m_subPops(subPops), m_vars(), m_suffix(suffix)
{
const char * allowedVars[] = {
- numOfSegSites_String, numOfSegSites_sp_String, ""
+ numOfSegSites_String, numOfSegSites_sp_String,
+ numOfFixedSites_String, numOfFixedSites_sp_String,""
};
- const char * defaultVars[] = { numOfSegSites_String, "" };
+ const char * defaultVars[] = { numOfSegSites_String, numOfFixedSites_String,"" };
m_vars.obtainFrom(vars, allowedVars, defaultVars);
}
@@ -697,12 +698,14 @@
const vectoru & loci = m_loci.elems(&pop);
DBG_DO(DBG_STATOR, cerr << "Count number of segregating sites for " << loci.size() << " loci " << endl);
+ std::set<size_t> allFixedSites;
std::set<size_t> allSegSites;
// for each subpopulation.
subPopList subPops = m_subPops.expandFrom(pop);
subPopList::const_iterator sp = subPops.begin();
subPopList::const_iterator spEnd = subPops.end();
for (; sp != spEnd; ++sp) {
+ std::set<size_t> fixedSites;
std::set<size_t> segSites;
pop.activateVirtualSubPop(*sp);
@@ -710,6 +713,8 @@
#ifdef MUTANTALLELE
IndIterator ind = pop.indIterator(sp->subPop());
for (; ind.valid(); ++ind) {
+ std::set<size_t> indFixedSitesPrev;
+ std::set<size_t> indFixedSitesCurr;
GenoIterator it = ind->genoBegin();
GenoIterator it_end = ind->genoEnd();
compressed_vector<Allele>::index_array_type::iterator index_it = it.getIndexIterator();
@@ -717,42 +722,52 @@
compressed_vector<Allele>::value_array_type::iterator value_it = it.getValueIterator();
size_t indIndex = it.getIndex();
for (; index_it != index_it_end; ++index_it, ++value_it) {
- if (m_loci.allAvail()) {
- if (*value_it != 0) {
- segSites.insert(*index_it - indIndex);
+ for (size_t idx = 0; idx < loci.size(); ++idx) {
+ if (*index_it == indIndex + loci[idx] && *value_it != 0) {
+ segSites.insert(loci[idx]);
+ indFixedSitesCurr.insert(loci[idx]);
+ break;
}
- } else {
- for (size_t idx = 0; idx < loci.size(); ++idx) {
- if (*index_it == indIndex + loci[idx] && *value_it != 0) {
- segSites.insert(loci[idx]);
- break;
- }
- }
}
}
+ if (ind == pop.indIterator(sp->subPop()))
+ indFixedSitesPrev = indFixedSitesCurr;
+ else {
+ set_intersection(indFixedSitesCurr.begin(), indFixedSitesCurr.end(),
+ indFixedSitesPrev.begin(), indFixedSitesPrev.end(),
+ std::inserter(fixedSites, fixedSites.begin()));
+ indFixedSitesPrev = fixedSites;
+ }
}
#else
for (ssize_t idx = 0; idx < static_cast<ssize_t>(loci.size()); ++idx) {
size_t loc = loci[idx];
+ size_t siteCount = 0;
IndAlleleIterator a = pop.alleleIterator(loc, sp->subPop());
for (; a.valid(); ++a)
- if (*a != 0) {
- segSites.insert((ULONG)loc);
- break;
- }
+ if (*a != 0) siteCount++;
+ if (siteCount == pop.subPopSize(sp->subPop()) * pop.ploidy())
+ fixedSites.insert((ULONG)loc);
+ else if (siteCount > 0)
+ segSites.insert((ULONG)loc);
}
#endif
pop.deactivateVirtualSubPop(sp->subPop());
+ if (m_vars.contains(numOfFixedSites_sp_String))
+ pop.getVars().setVar(subPopVar_String(*sp, numOfFixedSites_String) + m_suffix, fixedSites.size());
if (m_vars.contains(numOfSegSites_sp_String))
pop.getVars().setVar(subPopVar_String(*sp, numOfSegSites_String) + m_suffix, segSites.size());
+ allFixedSites.insert(fixedSites.begin(), fixedSites.end());
allSegSites.insert(segSites.begin(), segSites.end());
}
// output whole population
+ if (m_vars.contains(numOfFixedSites_String))
+ pop.getVars().setVar(numOfFixedSites_String + m_suffix, allFixedSites.size());
if (m_vars.contains(numOfSegSites_String))
pop.getVars().setVar(numOfSegSites_String + m_suffix, allSegSites.size());
return true;
Modified: trunk/src/stator.h
===================================================================
--- trunk/src/stator.h 2012-03-02 19:54:08 UTC (rev 4417)
+++ trunk/src/stator.h 2012-03-02 21:22:16 UTC (rev 4418)
@@ -405,7 +405,9 @@
{
private:
#define numOfSegSites_String "numOfSegSites"
+#define numOfFixedSites_String "numOfFixedSites"
#define numOfSegSites_sp_String "numOfSegSites_sp"
+#define numOfFixedSites_sp_String "numOfFixedSites_sp"
public:
statNumOfSegSites(const lociList & loci, const subPopList & subPops,
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2012-03-05 01:41:11
|
Revision: 4422
http://simupop.svn.sourceforge.net/simupop/?rev=4422&view=rev
Author: simupop
Date: 2012-03-05 01:41:04 +0000 (Mon, 05 Mar 2012)
Log Message:
-----------
Lineage Patch 4: add m_lineagePtr to the Individual class and implement individul level lineage manipulation and access functions
Modified Paths:
--------------
trunk/src/individual.cpp
trunk/src/individual.h
trunk/src/simuPOP_cfg.h
Modified: trunk/src/individual.cpp
===================================================================
--- trunk/src/individual.cpp 2012-03-04 06:02:31 UTC (rev 4421)
+++ trunk/src/individual.cpp 2012-03-05 01:41:04 UTC (rev 4422)
@@ -40,6 +40,7 @@
m_flags = rhs.m_flags;
setGenoPtr(rhs.genoPtr());
setInfoPtr(rhs.infoPtr());
+ LINEAGE_EXPR(setLineagePtr(rhs.lineagePtr()));
// also copy genoStru pointer...
this->setGenoStruIdx(rhs.genoStruIdx());
return *this;
@@ -55,6 +56,7 @@
copy(rhs.genoBegin(), rhs.genoEnd(), genoBegin());
#endif
copy(rhs.infoBegin(), rhs.infoEnd(), infoBegin());
+ LINEAGE_EXPR(copy(rhs.lineageBegin(), rhs.lineageEnd(), lineageBegin()));
// also copy genoStru pointer...
this->setGenoStruIdx(rhs.genoStruIdx());
return *this;
@@ -81,6 +83,12 @@
if (*(m_genoPtr + i) != *(rhs.m_genoPtr + i))
return false;
+#ifdef LINEAGE
+ for (size_t i = 0, iEnd = genoSize(); i < iEnd; ++i)
+ if (*(m_lineagePtr + i) != *(rhs.m_lineagePtr + i))
+ return false;
+#endif
+
for (size_t i = 0, iEnd = infoSize(); i < iEnd; ++i)
if (*(m_infoPtr + i) != *(rhs.m_infoPtr + i)) {
DBG_DO(DBG_POPULATION, cerr << "Information field " << infoField(i) << " differ" << endl);
@@ -133,22 +141,22 @@
}
-UINT Individual::allele(size_t idx, ssize_t p, ssize_t chrom) const
+ULONG Individual::allele(size_t idx, ssize_t p, ssize_t chrom) const
{
DBG_FAILIF(p < 0 && chrom >= 0, ValueError,
"A valid ploidy index has to be specified if chrom is non-positive");
if (p < 0) {
CHECKRANGEGENOSIZE(idx);
- return static_cast<UINT>(*(m_genoPtr + idx));
+ return static_cast<ULONG>(*(m_genoPtr + idx));
} else if (chrom < 0) {
CHECKRANGEABSLOCUS(idx);
CHECKRANGEPLOIDY(static_cast<size_t>(p));
- return static_cast<UINT>(*(m_genoPtr + idx + p * totNumLoci()));
+ return static_cast<ULONG>(*(m_genoPtr + idx + p * totNumLoci()));
} else {
CHECKRANGELOCUS(chrom, idx);
CHECKRANGEPLOIDY(static_cast<size_t>(p));
CHECKRANGECHROM(static_cast<size_t>(chrom));
- return static_cast<UINT>(*(m_genoPtr + idx + p * totNumLoci() + chromBegin(chrom)));
+ return static_cast<ULONG>(*(m_genoPtr + idx + p * totNumLoci() + chromBegin(chrom)));
}
}
@@ -316,6 +324,72 @@
#endif
+#ifdef LINEAGE
+PyObject * Individual::lineage(const uintList & ply, const uintList & ch)
+{
+ DBG_WARNIF(true, "The returned object of function Individual.lineage() is a special "
+ "carray_lineage object that reflects the underlying genotype lineage of an individual. "
+ "It will become invalid once the population changes. Please use "
+ "list(ind.lineage()) if you would like to keep a copy of lineages");
+
+ size_t beginP = 0;
+ size_t endP = 0;
+ size_t beginCh = 0;
+ size_t endCh = 0;
+
+ if (ply.allAvail())
+ endP = ploidy();
+ else {
+ const vectoru & ploidys = ply.elems();
+ if (ploidys.empty()) {
+ Py_INCREF(Py_None);
+ return Py_None;
+ }
+ beginP = ploidys[0];
+ endP = ploidys[0];
+ CHECKRANGEPLOIDY(static_cast<size_t>(beginP));
+ for (size_t i = 1; i < ploidys.size(); ++i) {
+ CHECKRANGEPLOIDY(static_cast<size_t>(ploidys[i]));
+ if (beginP > ploidys[i])
+ beginP = ploidys[i];
+ if (endP < ploidys[i])
+ endP = ploidys[i];
+ }
+ ++endP;
+ }
+ if (ch.allAvail())
+ endCh = numChrom();
+ else {
+ const vectoru & chroms = ch.elems();
+ if (chroms.empty()) {
+ Py_INCREF(Py_None);
+ return Py_None;
+ }
+ beginCh = chroms[0];
+ endCh = chroms[0];
+ CHECKRANGECHROM(static_cast<size_t>(beginCh));
+ for (size_t i = 1; i < chroms.size(); ++i) {
+ CHECKRANGECHROM(static_cast<size_t>(chroms[i]));
+ if (beginCh > chroms[i])
+ beginCh = chroms[i];
+ if (endCh < chroms[i])
+ endCh = chroms[i];
+ }
+ ++endCh;
+ }
+
+ if (endP > beginP + 1) {
+ // has to be all chromosomes
+ DBG_FAILIF(beginCh != 0 || endCh != numChrom(), ValueError,
+ "If multiple ploidy are chosen, all chromosomes has to be chosen.");
+ return Lineage_Vec_As_NumArray(m_lineagePtr + beginP * totNumLoci(),
+ m_lineagePtr + endP * totNumLoci());
+ } else
+ return Lineage_Vec_As_NumArray(m_lineagePtr + beginP * totNumLoci() + chromBegin(beginCh),
+ m_lineagePtr + beginP * totNumLoci() + chromEnd(endCh - 1));
+}
+#endif
+
// Fix me: This one has to optimize
PyObject * Individual::genoAtLoci(const lociList & lociList)
{
@@ -324,7 +398,7 @@
if (isHaplodiploid() && sex() == MALE)
ply = 1;
- vector<UINT> alleles;
+ vector<ULONG> alleles;
if (lociList.allAvail()) {
@@ -407,6 +481,26 @@
}
}
+#ifdef LINEAGE
+void Individual::setLineage(long lineage, size_t idx, int p, int chrom)
+{
+ DBG_FAILIF(p < 0 && chrom >= 0, ValueError,
+ "A valid ploidy index has to be specified if chrom is non-positive");
+ if (p < 0) {
+ CHECKRANGEGENOSIZE(idx);
+ *(m_lineagePtr + idx) = lineage;
+ } else if (chrom < 0) {
+ CHECKRANGEABSLOCUS(idx);
+ CHECKRANGEPLOIDY(static_cast<size_t>(p));
+ *(m_lineagePtr + idx + p * totNumLoci()) = lineage;
+ } else {
+ CHECKRANGELOCUS(static_cast<size_t>(chrom), idx);
+ CHECKRANGEPLOIDY(static_cast<size_t>(p));
+ CHECKRANGECHROM(static_cast<size_t>(chrom));
+ *(m_lineagePtr + idx + p * totNumLoci() + chromBegin(chrom)) = lineage;
+ }
+}
+#endif
void Individual::setGenotype(const uintList & genoList, const uintList & ply, const uintList & ch)
{
@@ -472,13 +566,22 @@
if (swapContent) {
Allele tmp;
+ LINEAGE_EXPR(long tmpLineage);
for (size_t i = 0, iEnd = genoSize(); i < iEnd; i++) {
tmp = m_genoPtr[i];
m_genoPtr[i] = ind.m_genoPtr[i];
ind.m_genoPtr[i] = tmp;
}
+#ifdef LINEAGE
+ for (size_t i = 0, iEnd = genoSize(); i < iEnd; i++) {
+ tmpLineage = m_lineagePtr[i];
+ m_lineagePtr[i] = ind.m_lineagePtr[i];
+ ind.m_lineagePtr[i] = tmpLineage;
+ }
+#endif
} else {
std::swap(m_genoPtr, ind.m_genoPtr);
+ LINEAGE_EXPR(std::swap(m_lineagePtr, ind.m_lineagePtr));
}
}
Modified: trunk/src/individual.h
===================================================================
--- trunk/src/individual.h 2012-03-04 06:02:31 UTC (rev 4421)
+++ trunk/src/individual.h 2012-03-05 01:41:04 UTC (rev 4422)
@@ -109,6 +109,9 @@
Individual(const Individual & ind) :
GenoStruTrait(ind), m_flags(ind.m_flags),
m_genoPtr(ind.m_genoPtr),
+#ifdef LINEAGE
+ m_lineagePtr(ind.m_lineagePtr),
+#endif
m_infoPtr(ind.m_infoPtr)
{
}
@@ -126,6 +129,13 @@
m_genoPtr = pos;
}
+#ifdef LINEAGE
+ /// CPPONLY set lineage pointer
+ void setLineagePtr(LineageIterator pos)
+ {
+ m_lineagePtr = pos;
+ }
+#endif
/// CPPONLY set pointer to individual info
void setInfoPtr(InfoIterator pos)
@@ -150,6 +160,13 @@
return m_genoPtr;
}
+#ifdef LINEAGE
+ /// CPPONLY pointer to lineage
+ LineageIterator lineagePtr() const
+ {
+ return m_lineagePtr;
+ }
+#endif
/// CPPONLY
InfoIterator infoPtr() const
@@ -170,7 +187,7 @@
*
* <group>1-allele</group>
*/
- UINT allele(size_t idx, ssize_t ploidy = -1, ssize_t chrom = -1) const;
+ ULONG allele(size_t idx, ssize_t ploidy = -1, ssize_t chrom = -1) const;
/** return the name of <tt>allele(idx, ploidy, chrom)</tt>. If idx is
@@ -188,6 +205,18 @@
*/
void setAllele(Allele allele, size_t idx, int ploidy = -1, int chrom = -1);
+#ifdef LINEAGE
+ /** set lineage \e lineage to an ID, using its absolute index \e idx.
+ * If a ploidy \e ploidy and/or a chromosome indexes are given, \e idx is
+ * relative to the beginning of specified homologous copy of chromosomes
+ * (if \e chrom=-1) or the beginning of the specified homologous copy of
+ * specified chromosome (if \e chrom >= 0). <bf>This function is only
+ * available in modules with lineage allele type.</bf>
+ * <group>1-allele</group>
+ */
+ void setLineage(long lineage, size_t idx, int ploidy = -1, int chrom = -1);
+#endif
+
/** return an editable array (a \c carray object) that represents all
* alleles of an individual. If \e ploidy or \e chroms is given, only
* alleles on the specified chromosomes and homologous copy of chromosomes
@@ -201,6 +230,19 @@
mutantList mutants(const uintList & ploidy = uintList(), const uintList & chroms = uintList());
#endif
+
+#ifdef LINEAGE
+ /** return an editable array (a \c carray_lineage object) that represents the
+ * lineages of all alleles of an individual. If \e ploidy or \e chroms is
+ * given, only lineages on the specified chromosomes and homologous copy of
+ * chromosomes will be returned. If multiple chromosomes are specified,
+ * there should not be gaps between chromosomes. <bf>This function is only
+ * available in modules with lineage allele type.</bf>
+ * <group>2-genotype</group>
+ */
+ PyObject * lineage(const uintList & ploidy = uintList(), const uintList & chroms = uintList());
+#endif
+
/** CPPONLY
* Return a Python object with alleles at specified loci. This function
* is usually used to collect alleles to send to a user-provided function.
@@ -366,7 +408,20 @@
return m_genoPtr + genoSize();
}
+#ifdef LINEAGE
+ /// CPPONLY start of lineage
+ LineageIterator lineageBegin() const
+ {
+ return m_lineagePtr;
+ }
+ /// CPPONLY end of lineage
+ LineageIterator lineageEnd() const
+ {
+ return m_lineagePtr + genoSize();
+ }
+#endif
+
/// CPPONLY start of allele of the pth set of chromosome
GenoIterator genoBegin(size_t p) const
{
@@ -383,7 +438,23 @@
return m_genoPtr + (p + 1) * totNumLoci();
}
+#ifdef LINEAGE
+ /// CPPONLY start of lineage of the pth set of chromosome
+ LineageIterator lineageBegin(size_t p) const
+ {
+ CHECKRANGEPLOIDY(p);
+ return m_lineagePtr + p * totNumLoci();
+ }
+
+ /// CPPONLY end of lineage of the pth set of chromosome
+ LineageIterator lineageEnd(size_t p) const
+ {
+ CHECKRANGEPLOIDY(p);
+ return m_lineagePtr + (p + 1) * totNumLoci();
+ }
+#endif
+
/// CPPONLY start of allele of the pth set of chromosome, chrom ch
GenoIterator genoBegin(size_t p, size_t chrom) const
{
@@ -403,7 +474,27 @@
}
+#ifdef LINEAGE
+ /// CPPONLY start of lineage of the pth set of chromosome, chrom ch
+ LineageIterator lineageBegin(size_t p, size_t chrom) const
+ {
+ CHECKRANGEPLOIDY(p);
+ CHECKRANGECHROM(chrom);
+ return m_lineagePtr + p * totNumLoci() + chromBegin(chrom);
+ }
+
+
+ /// CPPONLY end of lineage of the pth set of chromosome
+ LineageIterator lineageEnd(size_t p, size_t chrom) const
+ {
+ CHECKRANGEPLOIDY(p);
+ CHECKRANGECHROM(chrom);
+ return m_lineagePtr + p * totNumLoci() + chromEnd(chrom);
+
+ }
+#endif
+
/// CPPONLY start of info
InfoIterator infoBegin() const
{
@@ -508,6 +599,11 @@
/// pointer to genotype.
GenoIterator m_genoPtr;
+#ifdef LINEAGE
+ /// pointer to lineage.
+ LineageIterator m_lineagePtr;
+#endif
+
/// pointer to info
InfoIterator m_infoPtr;
};
@@ -1157,7 +1253,275 @@
typedef CombinedAlleleIterator<RawIndIterator> IndAlleleIterator;
typedef CombinedAlleleIterator<ConstRawIndIterator> ConstIndAlleleIterator;
+#ifdef LINEAGE
+/* This is the lineageIterator version of the CombinedAlleleIterator. It cannot
+ * be derived from the previous template because of the use of lineage specific
+ * functions such as LineageBegin().
+ */
+template <typename T>
+class CombinedLineageIterator
+{
+public:
+ typedef std::forward_iterator_tag iterator_category;
+ typedef long value_type;
+ typedef long int difference_type;
+ typedef long & reference;
+ //typedef GenoIterator pointer;
+ CombinedLineageIterator()
+ {
+ }
+
+ CombinedLineageIterator(size_t shift, IndividualIterator<T> it,
+ LineageIterator ptr, LineageIterator ptrEnd, size_t size)
+ : m_useGappedIterator(true), m_valid(true), m_shift(shift),
+ m_ptr(ptr), m_ptrBegin(ptr), m_ptrEnd(ptrEnd), m_size(size), m_it(it),
+ // ignored
+ m_index(0), m_ploidy(0), m_chromType(0),
+ m_haplodiploid(false), m_p(0)
+ {
+ m_valid = m_ptr != m_ptrEnd;
+ m_ploidy = it->ploidy();
+ }
+
+
+ CombinedLineageIterator(size_t idx, IndividualIterator<T> it)
+ : m_useGappedIterator(false), m_valid(true), m_shift(),
+ m_ptr(), m_ptrBegin(), m_ptrEnd(), m_size(0), // belong to a previous one
+ m_it(it), m_index(idx), m_ploidy(0), m_chromType(0),
+ m_haplodiploid(false), m_p(0)
+ {
+ if (!it.valid()) {
+ m_valid = false;
+ return;
+ }
+
+ m_size = it->totNumLoci();
+ m_ploidy = it->ploidy();
+ m_haplodiploid = it->isHaplodiploid();
+ m_chromType = it->chromType(it->chromLocusPair(idx).first);
+ // we do not know anything about customized chromosome
+ // so we just assume it is autosome.
+ if (m_chromType == CUSTOMIZED)
+ m_chromType = AUTOSOME;
+ //
+ if (m_chromType == CHROMOSOME_Y) {
+ if (m_it->sex() == FEMALE) {
+ while (m_it.valid())
+ if ((++m_it)->sex() == MALE)
+ break;
+ m_valid = m_it.valid();
+ }
+ m_p = 1;
+ }
+ }
+
+
+ bool valid()
+ {
+ return m_valid;
+ }
+
+ IndividualIterator<T> individual()
+ {
+ if (m_useGappedIterator) {
+ int offset = (m_ptr - m_ptrBegin) / (m_size * m_ploidy);
+ return(m_it + offset);
+ } else {
+ return(m_it);
+ }
+ }
+
+ // this is the most important part!
+ long & operator *() const
+ {
+ if (m_useGappedIterator)
+ return *(m_ptr + m_shift);
+ else {
+ DBG_ASSERT(m_it.valid(), SystemError, "Cannot refer to an invalid individual iterator");
+ return *(m_it->lineageBegin() + m_index + m_p * m_size);
+ }
+ }
+
+
+ LineageIterator ptr()
+ {
+ if (m_useGappedIterator)
+ return m_ptr + m_shift;
+ else
+ return m_it->lineageBegin() + m_index + m_p * m_size;
+ }
+
+
+ void advance(IndividualIterator<T> & it, size_t & p, bool & valid)
+ {
+ DBG_ASSERT(valid, RuntimeError, "Can not advance invalid allele iterator");
+ if (m_chromType == AUTOSOME) {
+ ++p;
+ if (p == m_ploidy) {
+ p = 0;
+ ++it;
+ valid = it.valid();
+ }
+ } else if (m_chromType == CHROMOSOME_X) {
+ if (it->sex() == FEMALE) {
+ // X0 -> X1
+ if (p == 0)
+ ++p;
+ // X1 -> X0 of next ind (male or female)
+ else {
+ p = 0;
+ ++it;
+ valid = it.valid();
+ }
+ } else {
+ // male, no X1.
+ DBG_ASSERT(p == 0, SystemError,
+ "Male Individual only has the first homologous copy of chromosome X");
+ // next Individual, ploidy 0, sex does not matter.
+ ++it;
+ valid = it.valid();
+ }
+ } else if (m_chromType == CHROMOSOME_Y) {
+ DBG_ASSERT(it->sex() == MALE, SystemError,
+ "There is no chromosome Y for FEMALE Individuals");
+ while ((++it).valid())
+ if (it->sex() == MALE)
+ break;
+ p = 1;
+ valid = it.valid();
+ } else if (m_chromType == MITOCHONDRIAL) {
+ // only the first homologous copy is valid
+ DBG_ASSERT(p == 0, SystemError, "Only the first homologous copy of mitochondrial DNA can be iterated.");
+ ++ it;
+ valid = it.valid();
+ }
+ }
+
+
+ // return, then advance.
+ CombinedLineageIterator operator++(int)
+ {
+ // save current state
+ CombinedLineageIterator tmp(*this);
+
+ if (!valid())
+ return tmp;
+
+ if (m_useGappedIterator) {
+ m_ptr += m_size;
+ m_valid = m_ptr != m_ptrEnd;
+ } else {
+ DBG_ASSERT(m_it.valid(), SystemError, "Cannot refer to an invalid individual iterator");
+ advance(m_it, m_p, m_valid);
+ }
+ return tmp;
+ }
+
+
+ CombinedLineageIterator operator++()
+ {
+ if (!valid())
+ return *this;
+ if (m_useGappedIterator) {
+ m_ptr += m_size;
+ m_valid = m_ptr != m_ptrEnd;
+ } else {
+ DBG_ASSERT(m_it.valid(), SystemError, "Cannot refer to an invalid individual iterator");
+ advance(m_it, m_p, m_valid);
+ }
+ return *this;
+ }
+
+
+ CombinedLineageIterator & operator+=(difference_type diff)
+ {
+ if (!valid())
+ return *this;
+ if (m_useGappedIterator) {
+ m_ptr += diff * m_size;
+ if (m_ptr > m_ptrEnd)
+ m_ptr = m_ptrEnd;
+ m_valid = m_ptr != m_ptrEnd;
+ } else {
+ DBG_ASSERT(m_it.valid(), SystemError, "Cannot refer to an invalid individual iterator");
+ for (int i = 0; i < diff && m_valid; ++i)
+ advance(m_it, m_p, m_valid);
+ }
+ return *this;
+ }
+
+
+ CombinedLineageIterator operator+(difference_type diff)
+ {
+ CombinedLineageIterator tmp(*this);
+
+ if (!valid())
+ return tmp;
+
+ if (m_useGappedIterator) {
+ tmp.m_ptr += diff * m_size;
+ if (tmp.m_ptr > tmp.m_ptrEnd)
+ tmp.m_ptr = tmp.m_ptrEnd;
+ tmp.m_valid = tmp.m_ptr != tmp.m_ptrEnd;
+ } else {
+ DBG_ASSERT(m_it.valid(), SystemError, "Cannot refer to an invalid individual iterator");
+ for (int i = 0; i < diff && tmp.m_valid; ++i)
+ advance(tmp.m_it, tmp.m_p, tmp.m_valid);
+ }
+ return tmp;
+ }
+
+
+ bool operator!=(const CombinedLineageIterator & rhs)
+ {
+ if (m_useGappedIterator)
+ return m_ptr != rhs.m_ptr || m_shift != rhs.m_shift;
+ else {
+ //DBG_FAILIF(m_it.valid() && rhs.m_it.valid() &&
+ // (m_ploidy != rhs.m_ploidy || m_size != rhs.m_size
+ // || m_chromType != rhs.m_chromType), ValueError,
+ // "Iterator comparison fails");
+ return !(m_index == rhs.m_index && m_it == rhs.m_it &&
+ (m_p == rhs.m_p || !m_it.valid() || !rhs.m_it.valid()));
+ }
+ }
+
+
+private:
+ ///
+ bool m_useGappedIterator;
+ ///
+ bool m_valid;
+ //
+ size_t m_shift;
+ //
+ LineageIterator m_ptr;
+ //
+ LineageIterator m_ptrBegin;
+ //
+ LineageIterator m_ptrEnd;
+ // genosize
+ size_t m_size;
+ //
+ // The second iteration method
+ // Individual iterator
+ IndividualIterator<T> m_it;
+ // index of the locus
+ size_t m_index;
+ // overall ploidy
+ size_t m_ploidy;
+ // chromosome type
+ size_t m_chromType;
+ //
+ bool m_haplodiploid;
+ // current ploidy, used in individualiterator
+ size_t m_p;
+};
+
+typedef CombinedLineageIterator<RawIndIterator> IndLineageIterator;
+#endif
+
}
#ifndef SWIG
Modified: trunk/src/simuPOP_cfg.h
===================================================================
--- trunk/src/simuPOP_cfg.h 2012-03-04 06:02:31 UTC (rev 4421)
+++ trunk/src/simuPOP_cfg.h 2012-03-05 01:41:04 UTC (rev 4422)
@@ -122,6 +122,12 @@
#define toID(val) (static_cast<size_t>((val) + 0.5))
+#ifdef LINEAGE
+#define LINEAGE_EXPR(expr) expr
+#else
+#define LINEAGE_EXPR(expr)
+#endif
+
/// needed by the following typedefs
#include <vector>
using std::vector;
@@ -342,6 +348,8 @@
// if this is changed Info_Var_As_Numarray in utility.cpp also needs to be changed.
typedef std::vector<double>::iterator InfoIterator;
typedef std::vector<double>::const_iterator ConstInfoIterator;
+typedef std::vector<long>::iterator LineageIterator;
+typedef std::vector<long>::const_iterator ConstLineageIterator;
extern const size_t InvalidValue;
// FIXME: I need a type that is 32 or 64 bit long depending on platform
@@ -557,7 +565,9 @@
#define CHECKRANGEABSLOCUS(locus) DBG_FAILIF(locus >= totNumLoci(), IndexError, "absolute locus index (" + toStr(locus) + ") out of range of 0 ~ " + toStr(totNumLoci() - 1))
#define CHECKRANGEGENOSIZE(p) DBG_FAILIF(p >= genoSize(), IndexError, "locus index (" + toStr(p) + ") out of range of 0 ~ " + toStr(genoSize() - 1))
#define CHECKRANGESUBPOPMEMBER(ind, sp) DBG_FAILIF(subPopSize(sp) > 0 && ind >= subPopSize(sp), IndexError, "individual index (" + toStr(ind) + ") out of range 0 ~" + toStr(subPopSize(sp) - 1) + " in subpopulation " + toStr(sp))
-#define CHECKRANGEIND(ind) DBG_FAILIF(ind >= popSize(), IndexError, "individual index (" + toStr(ind) + ") out of range of 0 ~ " + toStr(popSize() - 1))
-#define CHECKRANGEINFO(ind) DBG_FAILIF(ind >= infoSize(), IndexError, "info index (" + toStr(ind) + ") out of range of 0 ~ " + toStr(infoSize() - 1))
+#define CHECKRANGEIND(ind) DBG_FAILIF(ind >= popSize(), IndexError, "individual index (" + toStr(ind) + ") " \
+ + (popSize() > 0 ? (" out of range of 0 ~ " + toStr(popSize() - 1)) : "invoked on a population without any individual."))
+#define CHECKRANGEINFO(ind) DBG_FAILIF(ind >= infoSize(), IndexError, "info index (" + toStr(ind) + ") " \
+ + (infoSize() > 0 ? (" out of range of 0 ~ " + toStr(infoSize() - 1)) : "invoked on a population without any information field."))
}
#endif
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2012-03-05 03:37:53
|
Revision: 4423
http://simupop.svn.sourceforge.net/simupop/?rev=4423&view=rev
Author: simupop
Date: 2012-03-05 03:37:47 +0000 (Mon, 05 Mar 2012)
Log Message:
-----------
Expose functions Individual.lineage and Individual.setLineage in other modules (although doing nothing there) so that scripts written in the lineage module could be executed directly in other modules
Modified Paths:
--------------
trunk/src/individual.cpp
trunk/src/individual.h
Modified: trunk/src/individual.cpp
===================================================================
--- trunk/src/individual.cpp 2012-03-05 01:41:04 UTC (rev 4422)
+++ trunk/src/individual.cpp 2012-03-05 03:37:47 UTC (rev 4423)
@@ -324,9 +324,9 @@
#endif
-#ifdef LINEAGE
PyObject * Individual::lineage(const uintList & ply, const uintList & ch)
{
+#ifdef LINEAGE
DBG_WARNIF(true, "The returned object of function Individual.lineage() is a special "
"carray_lineage object that reflects the underlying genotype lineage of an individual. "
"It will become invalid once the population changes. Please use "
@@ -387,8 +387,13 @@
} else
return Lineage_Vec_As_NumArray(m_lineagePtr + beginP * totNumLoci() + chromBegin(beginCh),
m_lineagePtr + beginP * totNumLoci() + chromEnd(endCh - 1));
+#else
+ (void) ply;
+ (void) ch;
+ Py_INCREF(Py_None);
+ return Py_None;
+#endif
}
-#endif
// Fix me: This one has to optimize
PyObject * Individual::genoAtLoci(const lociList & lociList)
@@ -481,9 +486,9 @@
}
}
-#ifdef LINEAGE
void Individual::setLineage(long lineage, size_t idx, int p, int chrom)
{
+#ifdef LINEAGE
DBG_FAILIF(p < 0 && chrom >= 0, ValueError,
"A valid ploidy index has to be specified if chrom is non-positive");
if (p < 0) {
@@ -499,8 +504,13 @@
CHECKRANGECHROM(static_cast<size_t>(chrom));
*(m_lineagePtr + idx + p * totNumLoci() + chromBegin(chrom)) = lineage;
}
+#else
+ (void) lineage;
+ (void) idx;
+ (void) p;
+ (void) chrom;
+#endif
}
-#endif
void Individual::setGenotype(const uintList & genoList, const uintList & ply, const uintList & ch)
{
Modified: trunk/src/individual.h
===================================================================
--- trunk/src/individual.h 2012-03-05 01:41:04 UTC (rev 4422)
+++ trunk/src/individual.h 2012-03-05 03:37:47 UTC (rev 4423)
@@ -205,17 +205,15 @@
*/
void setAllele(Allele allele, size_t idx, int ploidy = -1, int chrom = -1);
-#ifdef LINEAGE
/** set lineage \e lineage to an ID, using its absolute index \e idx.
* If a ploidy \e ploidy and/or a chromosome indexes are given, \e idx is
* relative to the beginning of specified homologous copy of chromosomes
* (if \e chrom=-1) or the beginning of the specified homologous copy of
- * specified chromosome (if \e chrom >= 0). <bf>This function is only
- * available in modules with lineage allele type.</bf>
+ * specified chromosome (if \e chrom >= 0). <bf>This function does nothing
+ * for modules without lineage information.</bf>
* <group>1-allele</group>
*/
void setLineage(long lineage, size_t idx, int ploidy = -1, int chrom = -1);
-#endif
/** return an editable array (a \c carray object) that represents all
* alleles of an individual. If \e ploidy or \e chroms is given, only
@@ -231,17 +229,15 @@
#endif
-#ifdef LINEAGE
/** return an editable array (a \c carray_lineage object) that represents the
* lineages of all alleles of an individual. If \e ploidy or \e chroms is
* given, only lineages on the specified chromosomes and homologous copy of
* chromosomes will be returned. If multiple chromosomes are specified,
- * there should not be gaps between chromosomes. <bf>This function is only
- * available in modules with lineage allele type.</bf>
+ * there should not be gaps between chromosomes. <bf>A \c None object will be
+ * returned for modules without lineage information.</bf>
* <group>2-genotype</group>
*/
PyObject * lineage(const uintList & ploidy = uintList(), const uintList & chroms = uintList());
-#endif
/** CPPONLY
* Return a Python object with alleles at specified loci. This function
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2012-03-05 03:52:07
|
Revision: 4424
http://simupop.svn.sourceforge.net/simupop/?rev=4424&view=rev
Author: simupop
Date: 2012-03-05 03:52:00 +0000 (Mon, 05 Mar 2012)
Log Message:
-----------
Add parameter lineageField to operator InitGenotype to optionally record lineage from the specified field
Modified Paths:
--------------
trunk/src/initializer.cpp
trunk/src/initializer.h
Modified: trunk/src/initializer.cpp
===================================================================
--- trunk/src/initializer.cpp 2012-03-05 03:37:47 UTC (rev 4423)
+++ trunk/src/initializer.cpp 2012-03-05 03:52:00 UTC (rev 4424)
@@ -163,10 +163,12 @@
const uintList & ploidy,
int begin, int end, int step, const intList & at,
const intList & reps, const subPopList & subPops,
- const stringList & infoFields)
+ const stringList & infoFields,
+ const string & lineageField)
: BaseOperator("", begin, end, step, at, reps, subPops, infoFields),
m_freq(freq), m_genotype(genotype.elems()), m_prop(prop),
- m_haplotypes(haplotypes.elems()), m_loci(loci), m_ploidy(ploidy)
+ m_haplotypes(haplotypes.elems()), m_loci(loci),
+ m_ploidy(ploidy), m_lineageField(lineageField)
{
for (size_t i = 0; i < m_freq.size(); ++i)
PARAM_FAILIF(fcmp_lt(m_freq[i], 0) || fcmp_gt(m_freq[i], 1), ValueError,
@@ -306,6 +308,25 @@
}
}
}
+#ifdef LINEAGE
+#pragma omp parallel if(numThreads() > 1)
+ {
+#ifdef _OPENMP
+ IndIterator it = pop.indIterator(sp->subPop(), omp_get_thread_num());
+#else
+ IndIterator it = pop.indIterator(sp->subPop());
+#endif
+ bool assignLineage = ! m_lineageField.empty();
+ int idIdx = assignLineage ? pop.infoIdx(m_lineageField) : 0;
+
+ for (; it.valid(); ++it) {
+ long lineage = assignLineage ? toID(it->info(idIdx)) : 0;
+ for (vectoru::iterator p = ploidy.begin(); p != ploidy.end(); ++p)
+ for (vectoru::const_iterator loc = loci.begin(); loc != loci.end(); ++loc)
+ it->setLineage(lineage, *loc, static_cast<int>(*p));
+ }
+ }
+#endif // LINEAGE
pop.deactivateVirtualSubPop(sp->subPop());
}
return true;
Modified: trunk/src/initializer.h
===================================================================
--- trunk/src/initializer.h 2012-03-05 03:37:47 UTC (rev 4423)
+++ trunk/src/initializer.h 2012-03-05 03:52:00 UTC (rev 4424)
@@ -184,7 +184,9 @@
* If the length of a haplotype is not enough to fill all loci, the
* haplotype will be reused. If a list (or a single) haplotypes are
* specified without \e freq or \e prop, they are used with equal
- * probability.
+ * probability. If \e lineageFrom is set for modules with lineage
+ * information, allele lineage (source) will be set to IDs in the
+ * specified information field (usually \e ind_id).
*
* In the last case, if a sequence of genotype is specified, it will be
* uesd repeatedly to initialize all alleles sequentially. This works
@@ -198,7 +200,8 @@
const uintList & ploidy = uintList(),
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 = vectorstr(),
+ const string & lineageField = "");
~InitGenotype()
@@ -231,6 +234,8 @@
//
const uintList m_ploidy;
+
+ const string m_lineageField;
};
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2012-03-05 20:40:24
|
Revision: 4426
http://simupop.svn.sourceforge.net/simupop/?rev=4426&view=rev
Author: simupop
Date: 2012-03-05 20:40:15 +0000 (Mon, 05 Mar 2012)
Log Message:
-----------
Rename Individual.setLineage to Individual.setAlleleLineage and add two functions Individual.alleleLineage and Individual.setLineage
Modified Paths:
--------------
trunk/src/individual.cpp
trunk/src/individual.h
trunk/src/initializer.cpp
Modified: trunk/src/individual.cpp
===================================================================
--- trunk/src/individual.cpp 2012-03-05 17:50:02 UTC (rev 4425)
+++ trunk/src/individual.cpp 2012-03-05 20:40:15 UTC (rev 4426)
@@ -161,6 +161,32 @@
}
+long Individual::alleleLineage(size_t idx, ssize_t p, ssize_t chrom) const
+{
+#ifdef LINEAGE
+ DBG_FAILIF(p < 0 && chrom >= 0, ValueError,
+ "A valid ploidy index has to be specified if chrom is non-positive");
+ if (p < 0) {
+ CHECKRANGEGENOSIZE(idx);
+ return *(m_lineagePtr + idx);
+ } else if (chrom < 0) {
+ CHECKRANGEABSLOCUS(idx);
+ CHECKRANGEPLOIDY(static_cast<size_t>(p));
+ return *(m_lineagePtr + idx + p * totNumLoci());
+ } else {
+ CHECKRANGELOCUS(chrom, idx);
+ CHECKRANGEPLOIDY(static_cast<size_t>(p));
+ CHECKRANGECHROM(static_cast<size_t>(chrom));
+ return *(m_lineagePtr + idx + p * totNumLoci() + chromBegin(chrom));
+ }
+#else
+ (void) idx;
+ (void) p;
+ (void) chrom;
+ return 0;
+#endif
+}
+
string Individual::alleleChar(size_t idx, ssize_t p, ssize_t chrom) const
{
DBG_FAILIF(p < 0 && chrom >= 0, ValueError,
@@ -486,7 +512,7 @@
}
}
-void Individual::setLineage(long lineage, size_t idx, int p, int chrom)
+void Individual::setAlleleLineage(long lineage, size_t idx, int p, int chrom)
{
#ifdef LINEAGE
DBG_FAILIF(p < 0 && chrom >= 0, ValueError,
@@ -567,6 +593,55 @@
}
+void Individual::setLineage(const uintList & lineageList, const uintList & ply, const uintList & ch)
+{
+#ifdef LINEAGE
+ const vectoru & lineage = lineageList.elems();
+
+ size_t sz = lineage.size();
+ size_t idx = 0;
+
+ vectoru ploidys = ply.elems();
+
+ if (ply.allAvail()) {
+ for (size_t i = 0; i < ploidy(); ++i)
+ ploidys.push_back(i);
+ } else {
+#ifndef OPTIMIZED
+ for (size_t i = 0; i < ploidys.size(); ++i) {
+ CHECKRANGEPLOIDY(static_cast<size_t>(ploidys[i]));
+ }
+#endif
+ }
+ vectoru chroms = ch.elems();
+ if (ch.allAvail()) {
+ for (size_t i = 0; i < numChrom(); ++i)
+ chroms.push_back(i);
+ } else {
+#ifndef OPTIMIZED
+ for (size_t i = 0; i < chroms.size(); ++i) {
+ CHECKRANGECHROM(static_cast<size_t>(chroms[i]));
+ }
+#endif
+ }
+ for (size_t i = 0; i < ploidys.size(); ++i) {
+ size_t p = ploidys[i];
+ for (size_t j = 0; j < chroms.size(); ++j) {
+ size_t chrom = chroms[j];
+ LineageIterator ptr = m_lineagePtr + p * totNumLoci() + chromBegin(chrom);
+
+ for (size_t i = 0; i < numLoci(chrom); i++, ++idx)
+ *(ptr + i) = static_cast<long>(lineage[idx % sz]);
+ }
+ }
+#else
+ (void) lineageList;
+ (void) ply;
+ (void) ch;
+#endif
+}
+
+
void Individual::swap(Individual & ind, bool swapContent)
{
if (genoStruIdx() != ind.genoStruIdx())
Modified: trunk/src/individual.h
===================================================================
--- trunk/src/individual.h 2012-03-05 17:50:02 UTC (rev 4425)
+++ trunk/src/individual.h 2012-03-05 20:40:15 UTC (rev 4426)
@@ -189,7 +189,6 @@
*/
ULONG allele(size_t idx, ssize_t ploidy = -1, ssize_t chrom = -1) const;
-
/** return the name of <tt>allele(idx, ploidy, chrom)</tt>. If idx is
* invalid (e.g. second homologus copy of chromosome Y), '_' is returned.
* <group>1-allele</group>
@@ -205,7 +204,19 @@
*/
void setAllele(Allele allele, size_t idx, int ploidy = -1, int chrom = -1);
- /** set lineage \e lineage to an ID, using its absolute index \e idx.
+ /** return the lineage of the allele at a locus, using its absolute index
+ * \e idx. If a ploidy \e ploidy and/or a chromosome indexes is given,
+ * \e idx is relative to the beginning of specified homologous copy of
+ * chromosomes (if \e chrom=-1) or the beginning of the specified
+ * homologous copy of specified chromosome (if \e chrom >= 0). <bf>
+ * This function returns 0 for modules without lineage information.</bf>
+ *
+ * <group>1-allele</group>
+ */
+ long alleleLineage(size_t idx, ssize_t ploidy = -1, ssize_t chrom = -1) const;
+
+
+ /** set lineage \e lineage to an allele, using its absolute index \e idx.
* If a ploidy \e ploidy and/or a chromosome indexes are given, \e idx is
* relative to the beginning of specified homologous copy of chromosomes
* (if \e chrom=-1) or the beginning of the specified homologous copy of
@@ -213,7 +224,7 @@
* for modules without lineage information.</bf>
* <group>1-allele</group>
*/
- void setLineage(long lineage, size_t idx, int ploidy = -1, int chrom = -1);
+ void setAlleleLineage(long lineage, size_t idx, int ploidy = -1, int chrom = -1);
/** return an editable array (a \c carray object) that represents all
* alleles of an individual. If \e ploidy or \e chroms is given, only
@@ -254,6 +265,17 @@
*/
void setGenotype(const uintList & geno, const uintList & ploidy = uintList(), const uintList & chroms = uintList());
+ /** Fill the lineage of an individual using a list of IDs \e lineage.
+ * If parameters \e ploidy and/or \e chroms are specified, lineages will
+ * be copied to only all or specified chromosomes on selected homologous
+ * copies of chromosomes. \c lineage will be reused if its length is less than
+ * number of allelic lineage to be filled. <bf>This function does nothing
+ * for modules without lineage information.</bf>
+ * <group>2-genotype</group>
+ */
+ void setLineage(const uintList & lineage, const uintList & ploidy = uintList(), const uintList & chroms = uintList());
+
+
/** return the sex of an individual, \c 1 for male and \c 2 for female.
* <group>3-sex</group>
*/
Modified: trunk/src/initializer.cpp
===================================================================
--- trunk/src/initializer.cpp 2012-03-05 17:50:02 UTC (rev 4425)
+++ trunk/src/initializer.cpp 2012-03-05 20:40:15 UTC (rev 4426)
@@ -323,7 +323,7 @@
long lineage = assignLineage ? toID(it->info(idIdx)) : 0;
for (vectoru::iterator p = ploidy.begin(); p != ploidy.end(); ++p)
for (vectoru::const_iterator loc = loci.begin(); loc != loci.end(); ++loc)
- it->setLineage(lineage, *loc, static_cast<int>(*p));
+ it->setAlleleLineage(lineage, *loc, static_cast<int>(*p));
}
}
#endif // LINEAGE
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2012-03-06 04:56:49
|
Revision: 4428
http://simupop.svn.sourceforge.net/simupop/?rev=4428&view=rev
Author: simupop
Date: 2012-03-06 04:56:42 +0000 (Tue, 06 Mar 2012)
Log Message:
-----------
Lineage patch 4: Add lineage support for mutation opeartors. Revise lineage interface for operator InitGenotype to make it more consistent with the mutation operators
Modified Paths:
--------------
trunk/src/initializer.cpp
trunk/src/initializer.h
trunk/src/mutator.cpp
trunk/src/mutator.h
Modified: trunk/src/initializer.cpp
===================================================================
--- trunk/src/initializer.cpp 2012-03-05 21:09:14 UTC (rev 4427)
+++ trunk/src/initializer.cpp 2012-03-06 04:56:42 UTC (rev 4428)
@@ -316,15 +316,16 @@
#else
IndIterator it = pop.indIterator(sp->subPop());
#endif
- bool assignLineage = ! m_lineageField.empty();
- int idIdx = assignLineage ? pop.infoIdx(m_lineageField) : 0;
+ if (pop.hasInfoField(m_lineageField)) {
+ int idIdx = pop.infoIdx(m_lineageField);
- for (; it.valid(); ++it) {
- long lineage = assignLineage ? toID(it->info(idIdx)) : 0;
- for (vectoru::iterator p = ploidy.begin(); p != ploidy.end(); ++p)
- for (vectoru::const_iterator loc = loci.begin(); loc != loci.end(); ++loc)
- it->setAlleleLineage(lineage, *loc, static_cast<int>(*p));
- }
+ for (; it.valid(); ++it) {
+ long lineage = toID(it->info(idIdx));
+ for (vectoru::iterator p = ploidy.begin(); p != ploidy.end(); ++p)
+ for (vectoru::const_iterator loc = loci.begin(); loc != loci.end(); ++loc)
+ it->setAlleleLineage(lineage, *loc, static_cast<int>(*p));
+ }
+ }
}
#endif // LINEAGE
pop.deactivateVirtualSubPop(sp->subPop());
Modified: trunk/src/initializer.h
===================================================================
--- trunk/src/initializer.h 2012-03-05 21:09:14 UTC (rev 4427)
+++ trunk/src/initializer.h 2012-03-06 04:56:42 UTC (rev 4428)
@@ -184,9 +184,9 @@
* If the length of a haplotype is not enough to fill all loci, the
* haplotype will be reused. If a list (or a single) haplotypes are
* specified without \e freq or \e prop, they are used with equal
- * probability. If \e lineageFrom is set for modules with lineage
- * information, allele lineage (source) will be set to IDs in the
- * specified information field (usually \e ind_id).
+ * probability. For modules with lineage allele type, if individual IDs
+ * are assigned to a field \e lineageField (default to \c ind_id), its
+ * value will be saved as the lineage of modified alleles.
*
* In the last case, if a sequence of genotype is specified, it will be
* uesd repeatedly to initialize all alleles sequentially. This works
@@ -201,7 +201,7 @@
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 string & lineageField = "");
+ const string & lineageField = "ind_id");
~InitGenotype()
Modified: trunk/src/mutator.cpp
===================================================================
--- trunk/src/mutator.cpp 2012-03-05 21:09:14 UTC (rev 4427)
+++ trunk/src/mutator.cpp 2012-03-06 04:56:42 UTC (rev 4428)
@@ -93,6 +93,11 @@
{
DBG_DO(DBG_MUTATOR, cerr << "Mutate replicate " << pop.rep() << endl);
+#ifdef LINEAGE
+ bool assignLineage = pop.hasInfoField(m_lineageField);
+ int lineageIdx = assignLineage ? pop.infoIdx(m_lineageField) : 0;
+#endif
+
// mapIn and mapOut
bool mapIn = !m_mapIn.empty() || m_mapIn.func().isValid();
vectoru const & mapInList = m_mapIn.elems();
@@ -142,12 +147,22 @@
size_t pos = bt.trialFirstSucc(i);
size_t lastPos = 0;
IndAlleleIterator ptr = pop.alleleIterator(locus, sp);
+ LINEAGE_EXPR(IndLineageIterator lineagePtr = pop.lineageIterator(locus, sp));
if (pos != Bernullitrials::npos) {
do {
+#ifdef LINEAGE
+ long lineage = 0;
+ if (assignLineage) {
+ lineagePtr += static_cast<IndLineageIterator::difference_type>(pos - lastPos);
+ lineage = toID(lineagePtr.individual()->info(lineageIdx));
+ }
+#endif
ptr += static_cast<IndAlleleIterator::difference_type>(pos - lastPos);
lastPos = pos;
if (!ptr.valid())
continue;
+ Allele alleleValue = *ptr;
+ (void) alleleValue; // suppress a warning for unused variable
DBG_DO(DBG_MUTATOR, cerr << "Allele " << int(*ptr) << " at locus " << locus);
if (mapIn) {
if (numMapInAllele > 0) {
@@ -171,7 +186,15 @@
static_cast<int>(*ptr)));
}
}
- DBG_DO(DBG_MUTATOR, cerr << " is mutated to " << int(*ptr) << endl);
+ DBG_DO(DBG_MUTATOR, cerr << " is mutated from ");
+ DBG_DO(DBG_MUTATOR, cerr << int(alleleValue) << " to " << int(*ptr) << endl);
+#ifdef LINEAGE
+ if (assignLineage && alleleValue != *ptr) {
+ DBG_DO(DBG_MUTATOR, cerr << "Lineage updated from " << *lineagePtr);
+ DBG_DO(DBG_MUTATOR, cerr << " to " << lineage << endl);
+ *lineagePtr = lineage;
+ }
+#endif
} while ( (pos = bt.trialNextSucc(i, pos)) != Bernullitrials::npos);
} // succ.any
}
@@ -188,9 +211,10 @@
const stringFunc & output,
int begin, int end, int step, const intList & at,
const intList & reps, const subPopList & subPops,
- const stringList & infoFields)
+ const stringList & infoFields,
+ const string & idField)
: BaseMutator(vectorf(1, 0), loci, mapIn, mapOut, 0, output, begin, end, step,
- at, reps, subPops, infoFields)
+ at, reps, subPops, infoFields, idField)
{
matrixf rateMatrix = rate.elems();
// step 0, determine mu
@@ -270,8 +294,11 @@
double incProb, UINT maxAllele, const floatListFunc & mutStep,
const uintListFunc & mapIn, const uintListFunc & mapOut, const stringFunc & output,
int begin, int end, int step, const intList & at,
- const intList & reps, const subPopList & subPops, const stringList & infoFields)
- : BaseMutator(rates, loci, mapIn, mapOut, 0, output, begin, end, step, at, reps, subPops, infoFields),
+ const intList & reps, const subPopList & subPops,
+ const stringList & infoFields,
+ const string & idField)
+ : BaseMutator(rates, loci, mapIn, mapOut, 0, output, begin, end,
+ step, at, reps, subPops, infoFields, idField),
m_incProb(incProb), m_maxAllele(maxAllele), m_mutStep(mutStep)
{
#ifdef BINARYALLELE
@@ -421,6 +448,11 @@
subPopList::const_iterator sp = subPops.begin();
subPopList::const_iterator spEnd = subPops.end();
+#ifdef LINEAGE
+ bool assignLineage = pop.hasInfoField(m_lineageField);
+ int lineageIdx = assignLineage ? pop.infoIdx(m_lineageField) : 0;
+#endif
+
for (; sp != spEnd; ++sp) {
if (sp->isVirtual())
pop.activateVirtualSubPop(*sp);
@@ -431,9 +463,14 @@
IndIterator ind = pop.indIterator(sp->subPop()) + static_cast<IndIterator::difference_type>(*it);
if (!ind.valid())
continue;
+ LINEAGE_EXPR(long lineage = assignLineage ? toID(ind->info(lineageIdx)) : 0);
for (size_t p = 0; p < m_ploidy.size(); ++p) {
if (m_loci.allAvail()) {
for (size_t i = 0; i < pop.totNumLoci(); ++i) {
+#ifdef LINEAGE
+ if (assignLineage && m_allele != ind->allele(i, static_cast<int>(m_ploidy[p])))
+ ind->setAlleleLineage(lineage, i, static_cast<int>(m_ploidy[p]));
+#endif
ind->setAllele(m_allele, i, static_cast<int>(m_ploidy[p]));
DBG_DO(DBG_MUTATOR,
cerr << "Mutate locus " << i << " at ploidy " << m_ploidy[p]
@@ -443,6 +480,10 @@
} else {
const vectoru & loci = m_loci.elems(&pop);
for (size_t i = 0; i < loci.size(); ++i) {
+#ifdef LINEAGE
+ if (assignLineage && m_allele != ind->allele(loci[i], static_cast<int>(m_ploidy[p])))
+ ind->setAlleleLineage(lineage, loci[i], static_cast<int>(m_ploidy[p]));
+#endif
ind->setAllele(m_allele, loci[i], static_cast<int>(m_ploidy[p]));
DBG_DO(DBG_MUTATOR,
cerr << "Mutate locus " << loci[i] << " at ploidy " << m_ploidy[p]
Modified: trunk/src/mutator.h
===================================================================
--- trunk/src/mutator.h 2012-03-05 21:09:14 UTC (rev 4427)
+++ trunk/src/mutator.h 2012-03-06 04:56:42 UTC (rev 4428)
@@ -86,6 +86,13 @@
* between a large number of alleles and advanced models such as random
* emission of alleles.
*
+ * If a valid information field is specified for parameter \e lineageField
+ * (default to \c ind_id) for modules with lineage allele type, the lineage
+ * of the mutated alleles will be the ID (stored in field \e linkageField)
+ * of individuals that harbor the mutated alleles. The lineage information
+ * will be transmitted along with the alleles so this feature allows you to
+ * track the source of mutants during evolution.
+ *
* Some mutation models are context dependent. Namely, how an allele
* mutates will depend on its adjecent alleles. Whereas most simuPOP
* mutators are context independent, some of them accept a parameter
@@ -101,10 +108,11 @@
int context = 0, 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 = vectorstr(),
+ const string & lineageField = "ind_id")
: BaseOperator(output, begin, end, step, at, reps, subPops, infoFields),
m_rates(rates.elems()), m_loci(loci), m_mapIn(mapIn), m_mapOut(mapOut),
- m_context(context * 2)
+ m_lineageField(lineageField), m_context(context * 2)
{
// NOTE: empty rates is allowed because a mutator might be
// used in a mixed mutator.
@@ -188,6 +196,8 @@
const uintListFunc m_mapOut;
+ const string m_lineageField;
+
// Be careful about this variable, which is not constant.
mutable vectoru m_context;
};
@@ -220,7 +230,7 @@
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 = vectorstr(), const string & lineageField = "ind_id");
/// destructor.
~MatrixMutator()
@@ -273,9 +283,10 @@
const uintListFunc & mapIn = uintListFunc(), const uintListFunc & mapOut = uintListFunc(),
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 intList & reps = intList(), const subPopList & subPops = subPopList(),
+ const stringList & infoFields = vectorstr(), const string & lineageField = "ind_id")
: BaseMutator(rates, loci, mapIn, mapOut, 0, output, begin, end, step, at,
- reps, subPops, infoFields), m_k(k)
+ reps, subPops, infoFields, lineageField), m_k(k)
{
#ifndef BINARYALLELE
if (m_k > 1 && static_cast<ULONG>(m_k - 1) > ModuleMaxAllele)
@@ -354,7 +365,7 @@
const uintListFunc & mapIn = uintListFunc(), const uintListFunc & mapOut = uintListFunc(), 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 = vectorstr(), const string & lineageField = "ind_id");
~StepwiseMutator()
{
@@ -415,8 +426,10 @@
PyObject * func = NULL, int context = 0, const uintListFunc & mapIn = uintListFunc(),
const uintListFunc & mapOut = uintListFunc(), 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())
- : BaseMutator(rates, loci, mapIn, mapOut, context, output, begin, end, step, at, reps, subPops, infoFields),
+ const intList & reps = intList(), const subPopList & subPops = subPopList(),
+ const stringList & infoFields = vectorstr(), const string & lineageField = "ind_id")
+ : BaseMutator(rates, loci, mapIn, mapOut, context, output, begin, end,
+ step, at, reps, subPops, infoFields, lineageField),
m_func(func)
{
DBG_ASSERT(m_func.isValid(), ValueError,
@@ -469,8 +482,10 @@
const uintListFunc & mapIn = uintListFunc(), const uintListFunc & mapOut = uintListFunc(),
int context = 0, 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())
- : BaseMutator(rates, loci, mapIn, mapOut, context, output, begin, end, step, at, reps, subPops, infoFields),
+ const intList & reps = intList(), const subPopList & subPops = subPopList(),
+ const stringList & infoFields = vectorstr(), const string & lineageField = "ind_id")
+ : BaseMutator(rates, loci, mapIn, mapOut, context, output, begin, end,
+ step, at, reps, subPops, infoFields, lineageField),
m_mutators(mutators), m_sampler()
{
DBG_FAILIF(m_mutators.size() != prob.size(), ValueError,
@@ -539,8 +554,10 @@
const uintListFunc & mapIn = uintListFunc(), const uintListFunc & mapOut = uintListFunc(),
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())
- : BaseMutator(rates, loci, mapIn, mapOut, 0, output, begin, end, step, at, reps, subPops, infoFields),
+ const intList & reps = intList(), const subPopList & subPops = subPopList(),
+ const stringList & infoFields = vectorstr(), const string & lineageField = "ind_id")
+ : BaseMutator(rates, loci, mapIn, mapOut, 0, output, begin, end,
+ step, at, reps, subPops, infoFields, lineageField),
m_mutators(mutators), m_contexts(contexts.elems())
{
if (m_contexts.size() != 0) {
@@ -605,9 +622,10 @@
const uintList & inds = vectoru(), const stringFunc & output = ">",
int begin = 0, int end = -1, int step = 1, const intList & at = vectori(),
const intList & reps = intList(), const subPopList & subPops = 0,
- const stringList & infoFields = vectorstr())
+ const stringList & infoFields = vectorstr(), const string & lineageField = "ind_id")
: BaseOperator(output, begin, end, step, at, reps, subPops, infoFields),
- m_loci(loci), m_allele(allele), m_ploidy(ploidy.elems()), m_inds(inds.elems())
+ m_lineageField(lineageField), m_loci(loci), m_allele(allele),
+ m_ploidy(ploidy.elems()), m_inds(inds.elems())
{
}
@@ -638,6 +656,8 @@
private:
/// applicable loci.
+ const string m_lineageField;
+
lociList m_loci;
Allele m_allele;
vectoru m_ploidy;
@@ -741,8 +761,8 @@
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()) :
- BaseOperator(output, begin, end, step, at, reps, subPops, infoFields),
+ const stringList & infoFields = vectorstr(), const string & lineageField = "ind_it") :
+ BaseOperator(output, begin, end, step, at, reps, subPops, infoFields, lineageField),
m_rate(rate), m_ranges(ranges), m_model(model)
{
const matrixi & rngs = m_ranges.elems();
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2012-03-07 03:27:22
|
Revision: 4435
http://simupop.svn.sourceforge.net/simupop/?rev=4435&view=rev
Author: simupop
Date: 2012-03-07 03:27:16 +0000 (Wed, 07 Mar 2012)
Log Message:
-----------
Lineage path 5: Add lineage support for genotype transmitters
Modified Paths:
--------------
trunk/src/transmitter.cpp
trunk/src/transmitter.h
Modified: trunk/src/transmitter.cpp
===================================================================
--- trunk/src/transmitter.cpp 2012-03-07 03:25:23 UTC (rev 4434)
+++ trunk/src/transmitter.cpp 2012-03-07 03:27:16 UTC (rev 4435)
@@ -43,10 +43,10 @@
void GenoTransmitter::initialize(const Individual & ind) const
{
- m_hasCustomizedChroms = !ind.customizedChroms().empty();
+ m_hasCustomizedChroms = !ind.customizedChroms().empty() || ind.mitochondrial() > 0;
m_lociToCopy.clear();
for (size_t ch = 0; ch < ind.numChrom(); ++ch)
- if (ind.chromType(ch) == CUSTOMIZED)
+ if (ind.chromType(ch) == CUSTOMIZED || ind.chromType(ch) == MITOCHONDRIAL)
m_lociToCopy.push_back(0);
else
m_lociToCopy.push_back(ind.numLoci(ch));
@@ -71,6 +71,10 @@
ind.genoBegin(ploidy) + m_chromIdx[chrom + 1], 0);
# endif
#endif
+#ifdef LINEAGE
+ fill(ind.lineageBegin(ploidy) + m_chromIdx[chrom],
+ ind.lineageBegin(ploidy) + m_chromIdx[chrom + 1], 0);
+#endif
}
@@ -96,6 +100,12 @@
offspring.genoBegin(ploidy) + m_chromIdx[chrom]);
# endif
#endif
+
+#ifdef LINEAGE
+ copy(parent.lineageBegin(parPloidy) + m_chromIdx[chrom],
+ parent.lineageBegin(parPloidy) + m_chromIdx[chrom + 1],
+ offspring.lineageBegin(ploidy) + m_chromIdx[chrom]);
+#endif
}
@@ -112,6 +122,12 @@
GenoIterator par = parent.genoBegin(parPloidy, ch);
GenoIterator off = offspring.genoBegin(ploidy, ch);
+#ifdef LINEAGE
+ LineageIterator parLineage = parent.lineageBegin(parPloidy, ch);
+ LineageIterator offLineage = offspring.lineageBegin(ploidy, ch);
+ LineageIterator lineage_end = parent.lineageEnd(parPloidy, ch);
+#endif
+
#ifdef BINARYALLELE
copyGenotype(par, off, m_lociToCopy[ch]);
#else
@@ -122,6 +138,7 @@
copy(par, par_end, off);
# endif
#endif
+ LINEAGE_EXPR(copy(parLineage, lineage_end, offLineage));
}
} else { // easy
#ifdef BINARYALLELE
@@ -134,6 +151,7 @@
copy(parent.genoBegin(parPloidy), parent.genoEnd(parPloidy), offspring.genoBegin(ploidy));
# endif
#endif
+ LINEAGE_EXPR(copy(parent.lineageBegin(parPloidy), parent.lineageEnd(parPloidy), offspring.lineageBegin(ploidy)));
}
}
@@ -166,6 +184,13 @@
size_t ch = chroms[i];
GenoIterator par = parent->genoBegin(p, ch);
GenoIterator off = offspring->genoBegin(p, ch);
+
+#ifdef LINEAGE
+ LineageIterator parLineage = parent->lineageBegin(p, ch);
+ LineageIterator offLineage = offspring->lineageBegin(p, ch);
+ LineageIterator lineage_end = parent->lineageEnd(p, ch);
+#endif
+
#ifdef BINARYALLELE
copyGenotype(par, off, offspring->numLoci(ch));
#else
@@ -176,6 +201,7 @@
copy(par, par_end, off);
# endif
#endif
+ LINEAGE_EXPR(copy(parLineage, lineage_end, offLineage));
}
}
} else if (m_hasCustomizedChroms) {
@@ -185,6 +211,13 @@
continue;
GenoIterator par = parent->genoBegin(p, ch);
GenoIterator off = offspring->genoBegin(p, ch);
+
+#ifdef LINEAGE
+ LineageIterator parLineage = parent->lineageBegin(p, ch);
+ LineageIterator offLineage = offspring->lineageBegin(p, ch);
+ LineageIterator lineage_end = parent->lineageEnd(p, ch);
+#endif
+
#ifdef BINARYALLELE
copyGenotype(par, off, m_lociToCopy[ch]);
#else
@@ -195,6 +228,7 @@
copy(par, par_end, off);
# endif
#endif
+ LINEAGE_EXPR(copy(parLineage, lineage_end, offLineage));
}
}
} else { // easy
@@ -208,6 +242,7 @@
copy(parent->genoBegin(), parent->genoEnd(), offspring->genoBegin());
# endif
#endif
+ LINEAGE_EXPR(copy(parent->lineageBegin(), parent->lineageEnd(), offspring->lineageBegin()));
}
// for clone transmitter, sex is also transmitted
offspring->setSex(parent->sex());
@@ -229,6 +264,7 @@
GenoTransmitter::initialize(ind);
m_chromX = ind.chromX();
m_chromY = ind.chromY();
+ m_mitochondrial = ind.mitochondrial();
m_numChrom = ind.numChrom();
}
@@ -246,11 +282,22 @@
if (m_chromX < 0 && m_chromY < 0 && !m_hasCustomizedChroms) {
// pointer to parental, and offspring chromosome copies
GenoIterator par[2];
- GenoIterator off;
+ GenoIterator off;
+#ifdef LINEAGE
+ LineageIterator parLineage[2];
+ LineageIterator offLineage;
+#endif
+
//
par[0] = parent.genoBegin(0);
par[1] = parent.genoBegin(1);
off = offspring.genoBegin(ploidy);
+#ifdef LINEAGE
+ parLineage[0] = parent.lineageBegin(0);
+ parLineage[1] = parent.lineageBegin(1);
+ offLineage = offspring.lineageBegin(ploidy);
+#endif
+
//
//
// 1. try to copy in blocks,
@@ -278,10 +325,15 @@
size_t length = parEnd - parBegin;
//
// the easiest case, try to get some speed up...
- if (length == 1)
+ if (length == 1) {
off[parBegin] = par[parPloidy][parBegin];
- else
+ LINEAGE_EXPR(offLineage[parBegin] = parLineage[parPloidy][parBegin]);
+ } else {
copyGenotype(par[parPloidy] + parBegin, off + parBegin, length);
+ LINEAGE_EXPR(copy(parLineage[parPloidy] + parBegin,
+ parLineage[parPloidy] + parBegin + length,
+ offLineage + parBegin));
+ }
//
if (ch != m_numChrom - 1)
parPloidy = nextParPloidy;
@@ -299,7 +351,8 @@
if ((ploidy == 0 && ch == m_chromY) || // maternal, Y chromosome
(ploidy == 1 &&
((ch == m_chromX && offspring.sex() == MALE) ||
- (ch == m_chromY && offspring.sex() == FEMALE)))) {
+ (ch == m_chromY && offspring.sex() == FEMALE) ||
+ (ch == m_mitochondrial)))) {
clearChromosome(offspring, ploidy, ch);
continue;
}
@@ -440,6 +493,12 @@
size_t src = getRNG().randInt(static_cast<ULONG>(m_mitoChroms.size()));
GenoIterator par = parent->genoBegin(0, m_mitoChroms[src]);
GenoIterator off = offspring->genoBegin(0, *it);
+#ifdef LINEAGE
+ LineageIterator parLineage = parent->lineageBegin(0, m_mitoChroms[src]);
+ LineageIterator offLineage = offspring->lineageBegin(0, *it);
+ LineageIterator lineage_end = parent->lineageEnd(0, m_mitoChroms[src]);
+#endif
+
#ifdef BINARYALLELE
copyGenotype(par, off, m_numLoci);
#else
@@ -450,6 +509,7 @@
copy(par, par_end, off);
# endif
#endif
+ LINEAGE_EXPR(copy(parLineage, lineage_end, offLineage));
for (size_t p = 1; p < pldy; ++p)
clearChromosome(*offspring, 1, static_cast<int>(*it));
}
@@ -465,7 +525,7 @@
:
GenoTransmitter(output, begin, end, step, at, reps, subPops, infoFields),
m_intensity(intensity), m_rates(rates.elems()), m_loci(loci),
- m_recBeforeLoci(0), m_convMode(convMode.elems()), m_chromX(-1), m_chromY(-1),
+ m_recBeforeLoci(0), m_convMode(convMode.elems()), m_chromX(-1), m_chromY(-1), m_mitochondrial(-1),
m_customizedBegin(-1), m_customizedEnd(-1), m_algorithm(0), m_debugOutput(NULL),
#ifdef _OPENMP
m_bt(numThreads(), getRNG())
@@ -539,9 +599,16 @@
m_chromX = ind.chromX();
m_chromY = ind.chromY();
+ m_mitochondrial = ind.mitochondrial();
if (!ind.customizedChroms().empty()) {
m_customizedBegin = static_cast<int>(ind.chromBegin(ind.customizedChroms()[0]));
m_customizedEnd = static_cast<int>(ind.chromEnd(ind.customizedChroms().back()));
+ if (m_mitochondrial > 0)
+ // mitochondrial must be the last chromosome before customized
+ m_customizedBegin -= 1;
+ } else if (m_mitochondrial > 0) {
+ m_customizedBegin = m_mitochondrial;
+ m_customizedEnd = m_mitochondrial;
}
// prepare m_bt
vectorf vecP;
@@ -571,9 +638,9 @@
if (chBegin == chEnd)
continue;
- if (ind.chromType(ch) == CUSTOMIZED) {
+ if (ind.chromType(ch) == CUSTOMIZED || ind.chromType(ch) == MITOCHONDRIAL) {
// recombine before customized chromosome.
- if (ind.numChrom() != ch + 1 && ind.chromType(ch + 1) != CUSTOMIZED) {
+ if (ind.numChrom() != ch + 1 && (ind.chromType(ch + 1) != CUSTOMIZED || ind.chromType(ch + 1) != MITOCHONDRIAL)) {
m_recBeforeLoci.push_back(chEnd);
vecP.push_back(0.5);
}
@@ -677,10 +744,18 @@
// use which copy of chromosome
GenoIterator cp[2], off;
+#ifdef LINEAGE
+ LineageIterator lineagep[2], lineageOff;
+#endif
cp[0] = parent.genoBegin(0);
cp[1] = parent.genoBegin(1);
off = offspring.genoBegin(ploidy);
+#ifdef LINEAGE
+ lineagep[0] = parent.lineageBegin(0);
+ lineagep[1] = parent.lineageBegin(1);
+ lineageOff = offspring.lineageBegin(ploidy);
+#endif
// handling of sex chromosomes, by specifying chromsome
// ranges with specified ploidy.
@@ -737,9 +812,11 @@
for (size_t gt = 0, bl = 0; gt < gtEnd; ++gt, --convCount) {
// do not copy genotype in the ignored region.
if ((ignoreBegin < 0 || gt < static_cast<size_t>(ignoreBegin) || gt >= static_cast<size_t>(ignoreEnd)) &&
- (m_customizedBegin < 0 || gt < static_cast<size_t>(m_customizedBegin) || gt >= static_cast<size_t>(m_customizedEnd)))
+ (m_customizedBegin < 0 || gt < static_cast<size_t>(m_customizedBegin) || gt >= static_cast<size_t>(m_customizedEnd))) {
// copy
off[gt] = cp[curCp][gt];
+ LINEAGE_EXPR(lineageOff[gt] = lineagep[curCp][gt]);
+ }
// look ahead
if (convCount == 0) { // conversion ...
if (forceFirstBegin > 0 && gt + 1 >= static_cast<size_t>(forceFirstBegin)
@@ -807,8 +884,10 @@
simuPOP::copy(cp[curCp] + gt, cp[curCp] + gt + m_recBeforeLoci[pos], off + gt);
gt = m_recBeforeLoci[pos];
# else
- for (; gt < m_recBeforeLoci[pos]; ++gt)
+ for (; gt < m_recBeforeLoci[pos]; ++gt) {
off[gt] = cp[curCp][gt];
+ LINEAGE_EXPR(lineageOff[gt] = lineagep[curCp][gt]);
+ }
# endif
curCp = (curCp + 1) % 2;
if (m_debugOutput)
@@ -831,8 +910,10 @@
simuPOP::copy(cp[curCp] + gt, cp[curCp] + gt + convEnd, off + gt);
gt = convEnd;
# else
- for (; gt < convEnd; ++gt)
+ for (; gt < convEnd; ++gt) {
off[gt] = cp[curCp][gt];
+ LINEAGE_EXPR(lineageOff[gt] = lineagep[curCp][gt]);
+ }
# endif
curCp = (curCp + 1) % 2;
if (m_debugOutput)
@@ -846,8 +927,10 @@
simuPOP::copy(cp[curCp] + gt, cp[curCp] + gt + gtEnd, off + gt);
gt = gtEnd;
# else
- for (; gt < gtEnd; ++gt)
+ for (; gt < gtEnd; ++gt) {
off[gt] = cp[curCp][gt];
+ LINEAGE_EXPR(lineageOff[gt] = lineagep[curCp][gt]);
+ }
# endif
curCp = (curCp + 1) % 2;
if (m_debugOutput)
@@ -872,8 +955,10 @@
simuPOP::copy(cp[curCp] + gt, cp[curCp] + gt + convEnd, off + gt);
gt = convEnd;
# else
- for (; gt < convEnd; ++gt)
+ for (; gt < convEnd; ++gt) {
off[gt] = cp[curCp][gt];
+ LINEAGE_EXPR(lineageOff[gt] = lineagep[curCp][gt]);
+ }
# endif
curCp = (curCp + 1) % 2;
if (m_debugOutput)
@@ -884,8 +969,10 @@
simuPOP::copy(cp[curCp] + gt, cp[curCp] + gt + gtEnd, off + gt);
gt = gtEnd;
# else
- for (; gt < gtEnd; ++gt)
+ for (; gt < gtEnd; ++gt) {
off[gt] = cp[curCp][gt];
+ LINEAGE_EXPR(lineageOff[gt] = lineagep[curCp][gt]);
+ }
# endif
#else
size_t gt = 0, gtEnd = 0;
@@ -897,6 +984,9 @@
// first piece
gtEnd = m_recBeforeLoci[pos];
copyGenotype(cp[curCp] + gt, off + gt, m_recBeforeLoci[pos] - gt);
+ // for binary module, this is not needed, but it is kept here anyway, in case some
+ // people would like to implement lineage feature for binary modules
+ LINEAGE_EXPR(copy(lineagep[curCp] + gt, lineagep[curCp] + m_recBeforeLoci[pos], lineageOff + gt));
gt = gtEnd;
curCp = (curCp + 1) % 2;
if (m_debugOutput)
@@ -913,6 +1003,8 @@
convEnd = gt + convCount;
if (convEnd < gtEnd) {
copyGenotype(cp[curCp] + gt, off + gt, convCount);
+ // not used for binary module
+ LINEAGE_EXPR(copy(lineagep[curCp] + gt, lineagep[curCp] + gt + convCount, lineageOff + gt));
gt = convEnd;
curCp = (curCp + 1) % 2;
if (m_debugOutput)
@@ -923,6 +1015,8 @@
}
// copy from the end of conversion to the next recombination point
copyGenotype(cp[curCp] + gt, off + gt, m_recBeforeLoci[pos] - gt);
+ // not used for binary module
+ LINEAGE_EXPR(copy(lineagep[curCp] + gt, lineagep[curCp] + m_recBeforeLoci[pos], lineageOff + gt));
gt = gtEnd;
curCp = (curCp + 1) % 2;
if (m_debugOutput)
@@ -943,6 +1037,8 @@
convEnd = gt + convCount;
if (convEnd < gtEnd) {
copyGenotype(cp[curCp] + gt, off + gt, convCount);
+ // not used for binary module
+ LINEAGE_EXPR(copy(lineagep[curCp] + gt, lineagep[curCp] + gt + convCount, lineageOff + gt));
gt = convEnd;
curCp = (curCp + 1) % 2;
if (m_debugOutput)
@@ -950,6 +1046,8 @@
}
}
copyGenotype(cp[curCp] + gt, off + gt, gtEnd - gt);
+ // not used for binary module
+ LINEAGE_EXPR(copy(lineagep[curCp] + gt, lineagep[curCp] + gtEnd, lineageOff + gt));
#endif
}
if (m_debugOutput)
Modified: trunk/src/transmitter.h
===================================================================
--- trunk/src/transmitter.h 2012-03-07 03:25:23 UTC (rev 4434)
+++ trunk/src/transmitter.h 2012-03-07 03:27:16 UTC (rev 4435)
@@ -65,14 +65,17 @@
/** Clear (set alleles to zero) chromosome \e chrom on the \e ploidy-th
* homologous set of chromosomes of individual \e ind. It is equivalent to
- * <tt>ind.setGenotype([0], ploidy, chrom)</tt>.
+ * <tt>ind.setGenotype([0], ploidy, chrom)</tt>, except that it also clears
+ * allele lineage if it is executed in a module with lineage allele type.
*/
void clearChromosome(const Individual & ind, int ploidy, size_t chrom) const;
/** Transmit chromosome \e chrom on the \e parPloidy set of homologous
* chromosomes from \e parent to the \e ploidy set of homologous
* chromosomes of \e offspring. It is equivalent to
- * <tt>offspring.setGenotype(parent.genotype(parPloidy, chrom), polidy, chrom)</tt>.
+ * <tt>offspring.setGenotype(parent.genotype(parPloidy, chrom), polidy, chrom)</tt>,
+ * except that it also copies allelic lineage when it is executed in a
+ * module with lineage allele type.
*/
void copyChromosome(const Individual & parent, int parPloidy,
Individual & offspring, int ploidy, size_t chrom) const;
@@ -80,7 +83,9 @@
/** Transmit the \e parPloidy set of homologous chromosomes from \e parent
* to the \e ploidy set of homologous chromosomes of \e offspring.
* Customized chromosomes are not copied. It is equivalent to
- * <tt>offspring.setGenotype(parent.genotype(parPloidy), ploidy)</tt>.
+ * <tt>offspring.setGenotype(parent.genotype(parPloidy), ploidy)</tt>,
+ * except that it also copies allelic lineage when it is executed in a
+ * module with lineage allele type.
*/
void copyChromosomes(const Individual & parent, int parPloidy,
Individual & offspring, int ploidy) const;
@@ -149,6 +154,8 @@
* is ignored. This operator by default copies genotypes on all
* autosome and sex chromosomes (excluding customized chromosomes), unless
* a parameter \e chroms is used to specify which chromosomes to copy.
+ * This operator also copies allelic lineage when it is executed in a
+ * module with lineage allele type.
*/
CloneGenoTransmitter(const stringFunc & output = "", const uintList & chroms = uintList(),
int begin = 0, int end = -1, int step = 1, const intList & at = vectori(),
@@ -203,13 +210,15 @@
/** Create a Mendelian genotype transmitter (a during-mating operator) that
* transmits genotypes from parents to offspring following Mendel's laws.
* Autosomes and sex chromosomes are handled but customized chromosomes
- * are ignored. Parameters \e subPops and \e infoFields are ignored.
+ * are ignored. Parameters \e subPops and \e infoFields are ignored. This
+ * operator also copies allelic lineage when it is executed in a module
+ * with lineage allele type.
*/
MendelianGenoTransmitter(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()) :
GenoTransmitter(output, begin, end, step, at, reps, subPops, infoFields),
- m_chromX(-1), m_chromY(-1), m_numChrom(0)
+ m_chromX(-1), m_chromY(-1), m_mitochondrial(-1), m_numChrom(0)
{
}
@@ -264,6 +273,8 @@
mutable int m_chromY;
+ mutable int m_mitochondrial;
+
mutable size_t m_numChrom;
};
@@ -279,7 +290,8 @@
/** Create a self-fertilization genotype transmitter that transmits
* genotypes of a parent to an offspring through self-fertilization.
* Cutsomized chromosomes are not handled. Parameters \e subPops and
- * \e infoFields are ignored.
+ * \e infoFields are ignored. This operator also copies allelic lineage
+ * when it is executed in a module with lineage allele type.
*/
SelfingGenoTransmitter(const stringFunc & output = "", int begin = 0, int end = -1, int step = 1, const intList & at = vectori(),
const intList & reps = intList(), const subPopList & subPops = subPopList(),
@@ -327,7 +339,8 @@
/** Create a haplodiploid genotype transmitter (during-mating operator)
* that transmit parental genotypes from parents to offspring in a
* haplodiploid population. Parameters \e subPops and \e infoFields
- * are ignored.
+ * are ignored. This operator also copies allelic lineage when it is
+ * executed in a module with lineage allele type.
*/
HaplodiploidGenoTransmitter(const stringFunc & output = "", int begin = 0, int end = -1, int step = 1, const intList & at = vectori(),
const intList & reps = intList(), const subPopList & subPops = subPopList(),
@@ -387,6 +400,8 @@
* mitochondrial chromosomes. These chromosomes should have the same
* length and the same number of loci. This operator transmits these
* chromosomes randomly from the female parent to offspring of both sexes.
+ * It also copies allelic lineage when it is executed in a module with
+ * lineage allele type.
*/
MitochondrialGenoTransmitter(const stringFunc & output = "",
const uintList & chroms = uintList(),
@@ -546,6 +561,9 @@
* have two lines of output. Note that individual IDs need to be set
* (using a \c IdTagger operator) before this Recombinator is applied.
*
+ * In addition to genotypes, this operator also copies alleleic lineage if
+ * it is executed in a module with lineage allele type.
+ *
* \note conversion tract length is usually short, and is estimated to be
* between 337 and 456 bp, with overall range between maybe 50 - 2500
* bp. This is usually not enough to convert, for example, two adjacent
@@ -631,6 +649,7 @@
// locataion of special chromosomes
mutable int m_chromX;
mutable int m_chromY;
+ mutable int m_mitochondrial;
mutable int m_customizedBegin;
mutable int m_customizedEnd;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2012-03-07 05:10:52
|
Revision: 4436
http://simupop.svn.sourceforge.net/simupop/?rev=4436&view=rev
Author: simupop
Date: 2012-03-07 05:10:45 +0000 (Wed, 07 Mar 2012)
Log Message:
-----------
Lineage patch 6: update the population class with all the lineage operations. This patch needs more testing and documentation.
Modified Paths:
--------------
trunk/src/population.cpp
trunk/src/population.h
Modified: trunk/src/population.cpp
===================================================================
--- trunk/src/population.cpp 2012-03-07 03:27:16 UTC (rev 4435)
+++ trunk/src/population.cpp 2012-03-07 05:10:45 UTC (rev 4436)
@@ -56,6 +56,9 @@
m_subPopIndex(size.elems().size() + 1),
m_vspSplitter(NULL),
m_genotype(0),
+#ifdef LINEAGE
+ m_lineage(0),
+#endif
m_info(0),
m_inds(0),
m_ancestralGens(ancGen),
@@ -113,6 +116,9 @@
m_subPopIndex(rhs.m_subPopIndex),
m_vspSplitter(NULL),
m_genotype(0),
+#ifdef LINEAGE
+ m_lineage(0),
+#endif
m_info(0),
m_inds(0),
m_ancestralGens(rhs.m_ancestralGens),
@@ -128,6 +134,7 @@
try {
m_inds.resize(rhs.m_popSize);
m_genotype.resize(m_popSize * genoSize());
+ LINEAGE_EXPR(m_lineage.resize(m_popSize * genoSize()));
// have 0 length for mpi/non-head node
m_info.resize(rhs.m_popSize * infoSize());
} catch (...) {
@@ -147,7 +154,13 @@
size_t infoStep = infoSize();
size_t step = genoSize();
GenoIterator ptr = m_genotype.begin();
+#ifdef LINEAGE
+ LineageIterator lineagePtr = m_lineage.begin();
+ for (size_t i = 0; i < m_popSize; ++i, ptr += step, lineagePtr += step, infoPtr += infoStep) {
+ m_inds[i].setLineagePtr(lineagePtr);
+#else
for (size_t i = 0; i < m_popSize; ++i, ptr += step, infoPtr += infoStep) {
+#endif
m_inds[i].setGenoPtr(ptr);
m_inds[i].setInfoPtr(infoPtr);
m_inds[i].copyFrom(rhs.m_inds[i]);
@@ -166,6 +179,11 @@
const vector<Individual> & rinds = rp.m_inds;
GenoIterator lg = lp.m_genotype.begin();
ConstGenoIterator rg = rp.m_genotype.begin();
+#ifdef LINEAGE
+ LineageIterator llin = lp.m_lineage.begin();
+ ConstLineageIterator rlin = rp.m_lineage.begin();
+#endif
+
InfoIterator li = lp.m_info.begin();
ConstInfoIterator ri = rp.m_info.begin();
@@ -174,6 +192,7 @@
for (size_t i = 0; i < ps; ++i) {
linds[i].setGenoPtr(rinds[i].genoPtr() - rg + lg);
linds[i].setInfoPtr(rinds[i].infoPtr() - ri + li);
+ LINEAGE_EXPR(linds[i].setLineagePtr(rinds[i].lineagePtr() - rlin + llin));
}
}
} catch (...) {
@@ -200,6 +219,7 @@
pop.m_subPopSize.swap(m_subPopSize);
pop.m_subPopNames.swap(m_subPopNames);
pop.m_genotype.swap(m_genotype);
+ LINEAGE_EXPR(pop.m_lineage.swap(m_lineage));
pop.m_info.swap(m_info);
pop.m_inds.swap(m_inds);
std::swap(pop.m_indOrdered, m_indOrdered);
@@ -533,7 +553,7 @@
size_t ct = chromType(chromLocusPair(locus).first);
// this is a complex case
if (hasActivatedVirtualSubPop() || !indOrdered()
- || (ct != AUTOSOME && ct != CUSTOMIZED) || isHaplodiploid())
+ || (ct != AUTOSOME && ct != CUSTOMIZED && ct != MITOCHONDRIAL) || isHaplodiploid())
// this is a complex case
return IndAlleleIterator(locus, indIterator(subPop));
else
@@ -544,7 +564,65 @@
totNumLoci());
}
+#ifdef LINEAGE
+/// CPPONLY allele begin, for given subPop
+IndLineageIterator Population::lineageIterator(size_t locus, size_t subPop)
+{
+ CHECKRANGEABSLOCUS(locus);
+ CHECKRANGESUBPOP(subPop);
+
+ size_t ct = chromType(chromLocusPair(locus).first);
+ // this is a complex case
+ if (hasActivatedVirtualSubPop() || !indOrdered()
+ || (ct != AUTOSOME && ct != CUSTOMIZED && ct != MITOCHONDRIAL) || isHaplodiploid())
+ // this is a complex case
+ return IndLineageIterator(locus, indIterator(subPop));
+ else
+ // this is a complex case
+ return IndLineageIterator(locus,
+ indIterator(subPop),
+ m_lineage.begin() + m_subPopIndex[subPop] * genoSize(),
+ m_lineage.begin() + m_subPopIndex[subPop + 1] * genoSize(),
+ totNumLoci());
+}
+#endif // LINEAGE
+
+PyObject * Population::lineage(vspID subPopID)
+{
+#ifdef LINEAGE
+ DBG_WARNIF(true, "The returned object of function Population.lineage() is a special "
+ "carray_lineage object that reflects the underlying lineage of a "
+ "population's alleles. It will become invalid once the population changes. "
+ "Please use list(pop.lineage()) if you would like to keep a copy of lineages");
+
+ vspID vsp = subPopID.resolve(*this);
+
+ DBG_FAILIF(vsp.isVirtual(), ValueError,
+ "Function genotype currently does not support virtual subpopulation");
+ DBG_FAILIF(hasActivatedVirtualSubPop(), ValueError,
+ "This operation is not allowed when there is an activated virtual subpopulation");
+
+ syncIndPointers();
+ if (!vsp.valid()) {
+ // directly expose values. Do not copy data over.
+ return Lineage_Vec_As_NumArray(m_lineage.begin(), m_lineage.end());
+ } else {
+ size_t subPop = vsp.subPop();
+ CHECKRANGESUBPOP(subPop);
+ // directly expose values. Do not copy data over.
+ return Lineage_Vec_As_NumArray(lineageBegin(subPop, true), lineageEnd(subPop, true));
+ }
+ Py_INCREF(Py_None);
+ return Py_None;
+#else
+ (void) subPopID;
+ Py_INCREF(Py_None);
+ return Py_None;
+#endif
+}
+
+
PyObject * Population::genotype(vspID subPopID)
{
@@ -612,6 +690,48 @@
}
+void Population::setLineage(const uintList & lineageList, vspID subPopID)
+{
+#ifdef LINEAGE
+ const vectoru & lineage = lineageList.elems();
+
+ vspID subPop = subPopID.resolve(*this);
+
+ syncIndPointers();
+ if (!subPop.valid()) {
+ LineageIterator ptr = m_lineage.begin();
+ size_t sz = lineage.size();
+ for (size_t i = 0; i < popSize() * genoSize(); ++i)
+ *(ptr++) = static_cast<long>(lineage[i % sz]);
+ return;
+ }
+
+ DBG_FAILIF(hasActivatedVirtualSubPop(), ValueError,
+ "This operation is not allowed when there is an activated virtual subpopulation");
+
+ size_t sp = subPop.subPop();
+ CHECKRANGESUBPOP(sp);
+
+ size_t sz = lineage.size();
+ if (!subPop.isVirtual()) {
+ LineageIterator ptr = lineageBegin(sp, true);
+ for (size_t i = 0; i < subPopSize(sp) * genoSize(); ++i)
+ *(ptr++) = static_cast<long>(lineage[i % sz]);
+ } else {
+ activateVirtualSubPop(subPop);
+ IndIterator it = indIterator(sp);
+ size_t i = 0;
+ for (; it.valid(); ++it)
+ for (LineageIterator git = it->lineageBegin(); git != it->lineageEnd(); ++git, ++i)
+ *git = static_cast<long>(lineage[i % sz]);
+ deactivateVirtualSubPop(subPop.subPop());
+ }
+#else
+ (void) lineageList;
+ (void) subPopID;
+#endif
+}
+
void Population::validate(const string & msg) const
{
#ifdef OPTIMIZED
@@ -658,6 +778,7 @@
if (step != 0 && m_popSize > MaxIndexSize / step)
throw RuntimeError("Population size times number of loci exceed maximum index size.");
m_genotype.resize(m_popSize * step);
+ LINEAGE_EXPR(m_lineage.resize(m_popSize * step));
m_info.resize(m_popSize * is);
m_inds.resize(m_popSize);
} catch (...) {
@@ -678,7 +799,13 @@
// reset individual pointers
InfoIterator infoPtr = m_info.begin();
GenoIterator ptr = m_genotype.begin();
+#ifdef LINEAGE
+ LineageIterator lineagePtr = m_lineage.begin();
+ for (size_t i = 0; i < m_popSize; ++i, ptr += step, infoPtr += is, lineagePtr += step) {
+ m_inds[i].setLineagePtr(lineagePtr);
+#else
for (size_t i = 0; i < m_popSize; ++i, ptr += step, infoPtr += is) {
+#endif
m_inds[i].setGenoPtr(ptr);
m_inds[i].setInfoPtr(infoPtr);
m_inds[i].setGenoStruIdx(genoStruIdx());
@@ -713,14 +840,22 @@
for (int depth = ancestralGens(); depth >= 0; --depth) {
useAncestralGen(depth);
- if (oldSize != newSize)
+ if (oldSize != newSize) {
m_genotype.resize(newSize * popSize());
+ LINEAGE_EXPR(m_lineage.resize(newSize * popSize()));
+ }
if (oldInfoSize != newInfoSize)
m_info.resize(newInfoSize * popSize());
// reset structure
InfoIterator infoPtr = m_info.begin();
GenoIterator ptr = m_genotype.begin();
+#ifdef LINEAGE
+ LineageIterator lineagePtr = m_lineage.begin();
+ for (size_t i = 0; i < m_popSize; ++i, ptr += newSize, infoPtr += newInfoSize, lineagePtr += newSize) {
+ m_inds[i].setLineagePtr(lineagePtr);
+#else
for (size_t i = 0; i < m_popSize; ++i, ptr += newSize, infoPtr += newInfoSize) {
+#endif
m_inds[i].setGenoStruIdx(stru);
m_inds[i].setGenoPtr(ptr);
m_inds[i].setInfoPtr(infoPtr);
@@ -841,6 +976,7 @@
#else
vectora newGenotype(genoSize() * newPopSize);
#endif
+ LINEAGE_EXPR(vectori newLineage(genoSize() * newPopSize));
vectorf newInfo(newPopSize * infoSize());
vector<Individual> newInds(newPopSize);
@@ -849,17 +985,23 @@
size_t step = genoSize();
size_t infoStep = infoSize();
GenoIterator ptr = newGenotype.begin();
+#ifdef LINEAGE
+ LineageIterator lineagePtr = newLineage.begin();
+ for (size_t i = 0; i < newPopSize; ++i, ptr += step, ++it, lineagePtr += step, infoPtr += infoStep) {
+ newInds[i].setLineagePtr(lineagePtr);
+#else
for (size_t i = 0; i < newPopSize; ++i, ptr += step, ++it, infoPtr += infoStep) {
+#endif
newInds[i].setGenoStruIdx(genoStruIdx());
newInds[i].setGenoPtr(ptr);
newInds[i].setInfoPtr(infoPtr);
newInds[i].copyFrom(*it); // copy everything, with info value
}
-
// now, switch!
m_genotype.swap(newGenotype);
m_info.swap(newInfo);
m_inds.swap(newInds);
+ LINEAGE_EXPR(m_lineage.swap(newLineage));
m_popSize = newPopSize;
setIndOrdered(true);
#ifdef MUTANTALLELE
@@ -958,6 +1100,10 @@
RawIndIterator newInd = m_inds.begin();
GenoIterator oldPtr = m_genotype.begin();
GenoIterator newPtr = m_genotype.begin();
+#ifdef LINEAGE
+ LineageIterator oldLineagePtr = m_lineage.begin();
+ LineageIterator newLineagePtr = m_lineage.begin();
+#endif
InfoIterator oldInfoPtr = m_info.begin();
InfoIterator newInfoPtr = m_info.begin();
@@ -976,6 +1122,7 @@
oldInd += spSize;
oldPtr += step * spSize;
oldInfoPtr += infoStep * spSize;
+ LINEAGE_EXPR(oldLineagePtr += step * spSize);
} else if (subPops.overlap(sp)) {
// partial removal
//
@@ -1000,14 +1147,17 @@
copy(oldPtr, oldPtr + step, newPtr);
#endif
copy(oldInfoPtr, oldInfoPtr + infoStep, newInfoPtr);
+ LINEAGE_EXPR(copy(oldLineagePtr, oldLineagePtr + step, newLineagePtr));
}
++newInd;
newPtr += step;
newInfoPtr += infoStep;
+ LINEAGE_EXPR(newLineagePtr += step);
}
++oldInd;
oldPtr += step;
oldInfoPtr += infoStep;
+ LINEAGE_EXPR(oldLineagePtr += step);
}
//
new_size.push_back(newSize);
@@ -1027,6 +1177,7 @@
#endif
copy(oldInd, oldInd + spSize, newInd);
copy(oldInfoPtr, oldInfoPtr + infoStep * spSize, newInfoPtr);
+ LINEAGE_EXPR(copy(oldLineagePtr, oldLineagePtr + step * spSize, newLineagePtr));
}
newInd += spSize;
newPtr += step * spSize;
@@ -1034,18 +1185,27 @@
oldInd += spSize;
oldPtr += step * spSize;
oldInfoPtr += infoStep * spSize;
+ LINEAGE_EXPR(newLineagePtr += step * spSize);
+ LINEAGE_EXPR(oldLineagePtr += step * spSize);
}
}
//
m_inds.erase(newInd, m_inds.end());
m_genotype.erase(newPtr, m_genotype.end());
m_info.erase(newInfoPtr, m_info.end());
+ LINEAGE_EXPR(m_lineage.erase(newLineagePtr, m_lineage.end()));
m_popSize = std::accumulate(new_size.begin(), new_size.end(), size_t(0));
setSubPopStru(new_size, new_spNames);
//
InfoIterator infoPtr = m_info.begin();
GenoIterator ptr = m_genotype.begin();
+#ifdef LINEAGE
+ LineageIterator lineagePtr = m_lineage.begin();
+ for (size_t i = 0; i < m_popSize; ++i, ptr += step, infoPtr += infoStep, lineagePtr += step) {
+ m_inds[i].setLineagePtr(lineagePtr);
+#else
for (size_t i = 0; i < m_popSize; ++i, ptr += step, infoPtr += infoStep) {
+#endif
m_inds[i].setGenoPtr(ptr);
m_inds[i].setInfoPtr(infoPtr);
}
@@ -1065,6 +1225,10 @@
InfoIterator newInfoPtr = m_info.begin();
GenoIterator oldPtr = m_genotype.begin();
GenoIterator newPtr = m_genotype.begin();
+#ifdef LINEAGE
+ LineageIterator oldLineagePtr = m_lineage.begin();
+ LineageIterator newLineagePtr = m_lineage.begin();
+#endif
//
for (size_t sp = 0; sp < numSubPop(); ++sp) {
size_t newSize = 0;
@@ -1081,14 +1245,17 @@
copy(oldPtr, oldPtr + step, newPtr);
#endif
copy(oldInfoPtr, oldInfoPtr + infoStep, newInfoPtr);
+ LINEAGE_EXPR(copy(oldLineagePtr, oldLineagePtr + step, newLineagePtr));
}
++newInd;
newPtr += step;
newInfoPtr += infoStep;
+ LINEAGE_EXPR(newLineagePtr += step);
}
++oldInd;
oldPtr += step;
oldInfoPtr += infoStep;
+ LINEAGE_EXPR(oldLineagePtr += step);
}
new_size[sp] = newSize;
}
@@ -1096,6 +1263,7 @@
m_inds.erase(newInd, m_inds.end());
m_genotype.erase(newPtr, m_genotype.end());
m_info.erase(newInfoPtr, m_info.end());
+ LINEAGE_EXPR(m_lineage.erase(newLineagePtr, m_lineage.end()));
m_popSize = std::accumulate(new_size.begin(), new_size.end(), size_t(0));
setSubPopStru(new_size, m_subPopNames);
//
@@ -1105,6 +1273,12 @@
m_inds[i].setGenoPtr(ptr);
m_inds[i].setInfoPtr(infoPtr);
}
+#ifdef LINEAGE
+ LineageIterator lineagePtr = m_lineage.begin();
+ for (size_t i = 0; i < m_popSize; ++i, lineagePtr += step) {
+ m_inds[i].setLineagePtr(lineagePtr);
+ }
+#endif
}
@@ -1268,6 +1442,10 @@
vectora new_genotype;
new_genotype.reserve(step * popSize());
#endif
+#ifdef LINEAGE
+ vectori new_lineage;
+ new_lineage.reserve(step * popSize());
+#endif
new_inds.reserve(popSize());
new_info.reserve(infoStep * popSize());
@@ -1278,6 +1456,7 @@
// do not remove.
new_inds.insert(new_inds.end(), rawIndBegin(src), rawIndEnd(src));
new_genotype.insert(new_genotype.end(), genoBegin(src, true), genoEnd(src, true));
+ LINEAGE_EXPR(new_lineage.insert(new_lineage.end(), lineageBegin(src, true), lineageEnd(src, true)));
if (infoStep > 0)
new_info.insert(new_info.end(),
@@ -1292,6 +1471,7 @@
m_inds.swap(new_inds);
m_genotype.swap(new_genotype);
m_info.swap(new_info);
+ LINEAGE_EXPR(m_lineage.swap(new_lineage));
setSubPopStru(new_size, new_names);
//
InfoIterator infoPtr = m_info.begin();
@@ -1300,6 +1480,12 @@
m_inds[i].setGenoPtr(ptr);
m_inds[i].setInfoPtr(infoPtr);
}
+#ifdef LINEAGE
+ LineageIterator lineagePtr = m_lineage.begin();
+ for (size_t i = 0; i < m_popSize; ++i, lineagePtr += step) {
+ m_inds[i].setLineagePtr(lineagePtr);
+ }
+#endif
return sps[0];
}
@@ -1328,6 +1514,10 @@
#endif
// append pop2 chromosomes to the first one
GenoIterator ptr = newGenotype.begin();
+#ifdef LINEAGE
+ vectori newLineage(genoSize() * m_popSize);
+ LineageIterator lineagePtr = newLineage.begin();
+#endif
size_t pEnd = ploidy();
for (size_t i = 0; i < m_popSize; ++i) {
@@ -1336,15 +1526,25 @@
GenoIterator ptr1 = m_inds[i].genoPtr();
GenoIterator ptr2 = pop.m_inds[i].genoPtr();
m_inds[i].setGenoPtr(ptr);
+#ifdef LINEAGE
+ LineageIterator linPtr1 = m_inds[i].lineagePtr();
+ LineageIterator linPtr2 = pop.m_inds[i].lineagePtr();
+ m_inds[i].setLineagePtr(lineagePtr);
+#endif
for (size_t p = 0; p < pEnd; ++p) {
- for (size_t j = 0; j < numLoci1; ++j)
+ for (size_t j = 0; j < numLoci1; ++j) {
*(ptr++) = *(ptr1++);
- for (size_t j = 0; j < numLoci2; ++j)
- *(ptr++) = *(ptr2++);
+ LINEAGE_EXPR(*(lineagePtr++) = *(linPtr1++));
+ }
+ for (size_t j = 0; j < numLoci2; ++j) {
+ *(ptr++) = *(ptr2++);
+ LINEAGE_EXPR(*(lineagePtr++) = *(linPtr2++));
+ }
}
}
m_genotype.swap(newGenotype);
+ LINEAGE_EXPR(m_lineage.swap(newLineage));
#ifdef MUTANTALLELE
// compressed_vectora must be setGenoPtr after swap
ptr = m_genotype.begin();
@@ -1379,6 +1579,7 @@
m_inds.insert(m_inds.end(), pop.m_inds.begin(), pop.m_inds.end());
m_genotype.insert(m_genotype.end(), pop.m_genotype.begin(), pop.m_genotype.end());
m_info.insert(m_info.end(), pop.m_info.begin(), pop.m_info.end());
+ LINEAGE_EXPR(m_lineage.insert(m_lineage.end(), pop.m_lineage.begin(), pop.m_lineage.end()));
// iterators ready
InfoIterator infoPtr = m_info.begin();
size_t step = genoSize();
@@ -1390,6 +1591,12 @@
m_inds[i].setGenoPtr(ptr);
m_inds[i].setInfoPtr(infoPtr);
}
+#ifdef LINEAGE
+ LineageIterator lineagePtr = m_lineage.begin();
+ for (size_t i = 0; i < m_popSize; ++i, lineagePtr += step) {
+ m_inds[i].setLineagePtr(lineagePtr);
+ }
+#endif
// rebuild index
m_subPopIndex.resize(numSubPop() + 1);
size_t j = 1;
@@ -1439,6 +1646,10 @@
#endif
// merge chromosome by chromosome
GenoIterator ptr = newGenotype.begin();
+#ifdef LINEAGE
+ vectori newLineage(genoSize() * m_popSize);
+ LineageIterator lineagePtr = newLineage.begin();
+#endif
size_t pEnd = ploidy();
size_t newSize = totNumLoci();
for (size_t i = 0; i < m_popSize; ++i) {
@@ -1448,16 +1659,27 @@
GenoIterator ptr2 = pop.m_inds[i].genoPtr();
// new genotype
m_inds[i].setGenoPtr(ptr);
+#ifdef LINEAGE
+ m_inds[i].setLineagePtr(lineagePtr);
+ LineageIterator lineagePtr1 = m_inds[i].lineagePtr();
+ LineageIterator lineagePtr2 = pop.m_inds[i].lineagePtr();
+#endif
// copy each allele
for (size_t p = 0; p < pEnd; ++p) {
- for (size_t i = 0; i < size1; ++i)
+ for (size_t i = 0; i < size1; ++i) {
ptr[indexes1[i]] = *(ptr1++);
- for (size_t i = 0; i < size2; ++i)
+ LINEAGE_EXPR(lineagePtr[indexes1[i]] = *(lineagePtr1++));
+ }
+ for (size_t i = 0; i < size2; ++i) {
ptr[indexes2[i]] = *(ptr2++);
+ LINEAGE_EXPR(lineagePtr[indexes2[i]] = *(lineagePtr2++));
+ }
ptr += newSize;
+ LINEAGE_EXPR(lineagePtr += newSize);
}
}
m_genotype.swap(newGenotype);
+ LINEAGE_EXPR(m_lineage.swap(newLineage));
#ifdef MUTANTALLELE
// compressed_vectora must be setGenoPtr after swap
ptr = m_genotype.begin();
@@ -1499,6 +1721,10 @@
// copy data over
GenoIterator newPtr = newGenotype.begin();
+#ifdef LINEAGE
+ vectori newLineage(newPopGenoSize, 0);
+ LineageIterator newLineagePtr = newLineage.begin();
+#endif
size_t pEnd = ploidy();
size_t gap = totNumLoci() - oldNumLoci;
@@ -1509,15 +1735,21 @@
GenoIterator oldPtr = m_inds[i].genoPtr();
// new genotype
m_inds[i].setGenoPtr(newPtr);
+ LINEAGE_EXPR(LineageIterator oldLineagePtr = m_inds[i].lineagePtr());
+ LINEAGE_EXPR(m_inds[i].setLineagePtr(newLineagePtr));
// copy each chromosome
for (size_t p = 0; p < pEnd; ++p) {
- for (size_t i = 0; i < oldNumLoci; ++i)
+ for (size_t i = 0; i < oldNumLoci; ++i) {
*(newPtr++) = *(oldPtr++);
+ LINEAGE_EXPR(*(newLineagePtr++) = *(oldLineagePtr++));
+ }
newPtr += gap;
+ LINEAGE_EXPR(newLineagePtr += gap);
}
}
m_genotype.swap(newGenotype);
+ LINEAGE_EXPR(m_lineage.swap(newLineage));
#ifdef MUTANTALLELE
// compressed_vectora must be setGenoPtr after swap
GenoIterator ptr = m_genotype.begin();
@@ -1571,6 +1803,10 @@
#endif
// copy data over
GenoIterator newPtr = newGenotype.begin();
+#ifdef LINEAGE
+ vectori newLineage(newPopGenoSize, 0);
+ LineageIterator newLineagePtr = newLineage.begin();
+#endif
size_t pEnd = ploidy();
for (size_t i = 0; i < m_popSize; ++i) {
// set new geno structure
@@ -1578,15 +1814,21 @@
GenoIterator oldPtr = m_inds[i].genoPtr();
// new genotype
m_inds[i].setGenoPtr(newPtr);
+ LINEAGE_EXPR(LineageIterator oldLineagePtr = m_inds[i].lineagePtr());
+ LINEAGE_EXPR(m_inds[i].setLineagePtr(newLineagePtr));
// copy each chromosome
for (size_t p = 0; p < pEnd; ++p) {
vectoru::iterator loc = loci.begin();
- for (; loc != loci.end(); ++loc)
+ for (; loc != loci.end(); ++loc) {
newPtr[*loc] = *(oldPtr++);
+ LINEAGE_EXPR(newLineagePtr[*loc] = *(oldLineagePtr++));
+ }
newPtr += totNumLoci();
+ LINEAGE_EXPR(newLineagePtr += totNumLoci());
}
}
m_genotype.swap(newGenotype);
+ LINEAGE_EXPR(m_lineage.swap(newLineage));
#ifdef MUTANTALLELE
// compressed_vectora must be setGenoPtr after swap
GenoIterator ptr = m_genotype.begin();
@@ -1631,6 +1873,13 @@
// set pointers
newInds[i].setInfoPtr(infoPtr);
}
+#ifdef LINEAGE
+ vectori newLineage(genoSize() * newPopSize);
+ LineageIterator lineagePtr = newLineage.begin();
+ for (size_t i = 0; i < newPopSize; ++i, lineagePtr += step) {
+ newInds[i].setLineagePtr(lineagePtr);
+ }
+#endif
// copy stuff over
size_t startSP = 0;
for (size_t sp = 0; sp < numSubPop(); ++sp) {
@@ -1648,6 +1897,7 @@
m_genotype.swap(newGenotype);
m_info.swap(newInfo);
m_inds.swap(newInds);
+ LINEAGE_EXPR(m_lineage.swap(newLineage));
m_popSize = newPopSize;
setIndOrdered(true);
m_subPopSize = newSubPopSizes;
@@ -1694,6 +1944,7 @@
#else
vectora new_genotype;
#endif
+ LINEAGE_EXPR(vectori new_lineage);
vectorf new_info;
if (rearrange) {
@@ -1710,10 +1961,12 @@
new_inds.resize(sz);
new_genotype.resize(sz * step);
new_info.resize(sz * infoStep);
+ LINEAGE_EXPR(new_lineage.resize(sz * step));
//
RawIndIterator newInd = new_inds.begin();
GenoIterator newPtr = new_genotype.begin();
InfoIterator newInfoPtr = new_info.begin();
+ LINEAGE_EXPR(LineageIterator newLineagePtr = new_lineage.begin());
it = subPops.begin();
for (; it != itEnd; ++it) {
@@ -1727,9 +1980,11 @@
copy(oldInd->genoBegin(), oldInd->genoEnd(), newPtr);
#endif
copy(oldInd->infoBegin(), oldInd->infoEnd(), newInfoPtr);
+ LINEAGE_EXPR(copy(oldInd->lineageBegin(), oldInd->lineageEnd(), newLineagePtr));
++newInd;
newPtr += step;
newInfoPtr += infoStep;
+ LINEAGE_EXPR(newLineagePtr += step);
}
deactivateVirtualSubPop(it->subPop());
}
@@ -1743,15 +1998,18 @@
new_inds.resize(sz);
new_genotype.resize(sz * step);
new_info.resize(sz * infoStep);
+ LINEAGE_EXPR(new_lineage.resize(sz * step));
//
RawIndIterator newInd = new_inds.begin();
GenoIterator newPtr = new_genotype.begin();
InfoIterator newInfoPtr = new_info.begin();
+ LINEAGE_EXPR(LineageIterator newLineagePtr = new_lineage.begin());
//
ConstRawIndIterator oldInd = m_inds.begin();
ConstGenoIterator oldPtr = m_genotype.begin();
ConstInfoIterator oldInfoPtr = m_info.begin();
+ LINEAGE_EXPR(ConstLineageIterator oldLineagePtr = m_lineage.begin());
for (size_t sp = 0; sp < numSubPop(); ++sp) {
size_t spSize = subPopSize(sp);
if (subPops.contains(sp)) {
@@ -1767,13 +2025,16 @@
copy(oldPtr, oldPtr + step * spSize, newPtr);
#endif
copy(oldInfoPtr, oldInfoPtr + infoStep * spSize, newInfoPtr);
+ LINEAGE_EXPR(copy(oldLineagePtr, oldLineagePtr + step * spSize, newLineagePtr));
oldInd += spSize;
oldPtr += step * spSize;
oldInfoPtr += infoStep * spSize;
+ LINEAGE_EXPR(oldLineagePtr += step * spSize);
newInd += spSize;
newPtr += step * spSize;
newInfoPtr += infoStep * spSize;
+ LINEAGE_EXPR(newLineagePtr += step * spSize);
} else if (subPops.overlap(sp)) {
// partial copy
//
@@ -1797,13 +2058,16 @@
copy(oldPtr, oldPtr + step, newPtr);
#endif
copy(oldInfoPtr, oldInfoPtr + infoStep, newInfoPtr);
+ LINEAGE_EXPR(copy(oldLineagePtr, oldLineagePtr + step, newLineagePtr));
++newInd;
newPtr += step;
newInfoPtr += infoStep;
+ LINEAGE_EXPR(newLineagePtr += step);
}
++oldInd;
oldPtr += step;
oldInfoPtr += infoStep;
+ LINEAGE_EXPR(oldLineagePtr += step);
}
//
new_size.push_back(newSize);
@@ -1814,6 +2078,7 @@
oldInd += spSize;
oldPtr += step * spSize;
oldInfoPtr += infoStep * spSize;
+ LINEAGE_EXPR(oldLineagePtr += step * spSize);
}
}
}
@@ -1821,10 +2086,12 @@
new_inds.resize(sz);
new_genotype.resize(sz * step);
new_info.resize(sz * infoStep);
+ LINEAGE_EXPR(new_lineage.resize(sz * step));
//
pop.m_inds.swap(new_inds);
pop.m_genotype.swap(new_genotype);
pop.m_info.swap(new_info);
+ LINEAGE_EXPR(pop.m_lineage.swap(new_lineage));
pop.m_popSize = sz;
pop.setSubPopStru(new_size, new_spNames);
//
@@ -1834,7 +2101,12 @@
pop.m_inds[i].setGenoPtr(ptr);
pop.m_inds[i].setInfoPtr(infoPtr);
}
-
+#ifdef LINEAGE
+ LineageIterator lineagePtr = pop.m_lineage.begin();
+ for (size_t i = 0; i < pop.m_popSize; ++i, lineagePtr += step) {
+ pop.m_inds[i].setLineagePtr(lineagePtr);
+ }
+#endif
return pop;
}
@@ -1854,6 +2126,7 @@
size_t infoStep = infoSize();
ConstRawIndIterator oldInd = m_inds.begin();
ConstGenoIterator oldPtr = m_genotype.begin();
+ LINEAGE_EXPR(ConstLineageIterator oldLineagePtr = m_lineage.begin());
ConstInfoIterator oldInfoPtr = m_info.begin();
size_t sz = 0;
@@ -1869,10 +2142,12 @@
#else
vectora new_genotype(sz * step);
#endif
+ LINEAGE_EXPR(vectori new_lineage(sz * step));
vectorf new_info(sz * infoStep);
RawIndIterator newInd = new_inds.begin();
GenoIterator newPtr = new_genotype.begin();
+ LINEAGE_EXPR(LineageIterator newLineagePtr = new_lineage.begin());
InfoIterator newInfoPtr = new_info.begin();
//
for (size_t sp = 0; sp < numSubPop(); ++sp) {
@@ -1889,13 +2164,16 @@
copy(oldPtr, oldPtr + step, newPtr);
#endif
copy(oldInfoPtr, oldInfoPtr + infoStep, newInfoPtr);
+ LINEAGE_EXPR(copy(oldLineagePtr, oldLineagePtr + step, newLineagePtr));
++newInd;
newPtr += step;
newInfoPtr += infoStep;
+ LINEAGE_EXPR(newLineagePtr += step);
}
++oldInd;
oldPtr += step;
oldInfoPtr += infoStep;
+ LINEAGE_EXPR(oldLineagePtr += step);
}
new_size.push_back(newSize);
}
@@ -1903,6 +2181,7 @@
pop.m_inds.swap(new_inds);
pop.m_genotype.swap(new_genotype);
pop.m_info.swap(new_info);
+ LINEAGE_EXPR(pop.m_lineage.swap(new_lineage));
pop.m_popSize = std::accumulate(new_size.begin(), new_size.end(), size_t(0));
pop.setSubPopStru(new_size, m_subPopNames);
//
@@ -1912,6 +2191,12 @@
pop.m_inds[i].setGenoPtr(ptr);
pop.m_inds[i].setInfoPtr(infoPtr);
}
+#ifdef LINEAGE
+ LineageIterator lineagePtr = pop.m_lineage.begin();
+ for (size_t i = 0; i < pop.m_popSize; ++i, lineagePtr += step) {
+ pop.m_inds[i].setLineagePtr(lineagePtr);
+ }
+#endif
return pop;
}
@@ -2160,6 +2445,10 @@
vectora new_genotype;
new_genotype.reserve(size * step);
#endif
+#ifdef LINEAGE
+ vectori new_lineage;
+ new_lineage.reserve(size * step);
+#endif
vectorf new_info;
new_inds.reserve(size);
@@ -2168,21 +2457,24 @@
if (!removeInd) {
new_inds.insert(new_inds.end(), m_inds.begin(), m_inds.end());
// handle genotype
- if (!removeLoci)
+ if (!removeLoci) {
new_genotype.insert(new_genotype.end(), m_genotype.begin(), m_genotype.end());
- else {
+ LINEAGE_EXPR(new_lineage.insert(new_lineage.end(), m_lineage.begin(), m_lineage.end()));
+ } else {
ConstRawIndIterator it = rawIndBegin();
ConstRawIndIterator it_end = rawIndEnd();
for (; it != it_end; ++it) {
GenoIterator ptr = it->genoBegin();
+ LINEAGE_EXPR(LineageIterator lineagePtr = it->lineageBegin());
for (size_t p = 0; p < pEnd; p += pStep) {
- for (lociPtr = new_loci.begin(); lociPtr != lociEnd; ++lociPtr)
+ for (lociPtr = new_loci.begin(); lociPtr != lociEnd; ++lociPtr) {
#ifdef MUTANTALLELE
new_genotype.push_back(newIdx++, *(ptr + *lociPtr + p));
#else
new_genotype.push_back(*(ptr + *lociPtr + p));
#endif
-
+ LINEAGE_EXPR(new_lineage.push_back(*(lineagePtr + *lociPtr + p)));
+ }
}
}
}
@@ -2208,11 +2500,14 @@
vectoru::iterator it_end = idx.end();
for (; it != it_end; ++it) {
new_inds.push_back(m_inds[*it]);
- if (!removeLoci)
+ if (!removeLoci) {
new_genotype.insert(new_genotype.end(),
indGenoBegin(*it), indGenoEnd(*it));
- else {
+ LINEAGE_EXPR(new_lineage.insert(new_lineage.end(),
+ indLineageBegin(*it), indLineageEnd(*it)));
+ } else {
GenoIterator ptr = indGenoBegin(*it);
+ LINEAGE_EXPR(LineageIterator lineagePtr = indLineageBegin(*it));
for (size_t p = 0; p < pEnd; p += pStep) {
for (lociPtr = new_loci.begin(); lociPtr != lociEnd; ++lociPtr)
#ifdef MUTANTALLELE
@@ -2220,6 +2515,7 @@
#else
new_genotype.push_back(*(ptr + *lociPtr + p));
#endif
+ LINEAGE_EXPR(new_lineage.push_back(*(lineagePtr + *lociPtr + p)));
}
}
if (!removeInfo)
@@ -2250,6 +2546,12 @@
new_inds[i].setGenoPtr(ptr);
new_inds[i].setInfoPtr(infoPtr);
}
+#ifdef LINEAGE
+ vectori::iterator lineagePtr = new_lineage.begin();
+ for (size_t i = 0; i < size; ++i, lineagePtr += step) {
+ new_inds[i].setLineagePtr(lineagePtr);
+ }
+#endif
// the arrays are ready, are they?
DBG_ASSERT(new_inds.size() == size && (new_genotype.size() == size * step)
&& (new_info.size() == size * infoStep), SystemError,
@@ -2262,6 +2564,7 @@
pop.m_inds.swap(new_inds);
pop.m_genotype.swap(new_genotype);
pop.m_info.swap(new_info);
+ LINEAGE_EXPR(pop.m_lineage.swap(new_lineage));
#ifdef MUTANTALLELE
// compressed_vectora must be setGenoPtr after swap
GenoIterator ptr = pop.m_genotype.begin();
@@ -2275,6 +2578,7 @@
pd.m_subPopSize.swap(spSizes);
pd.m_genotype.swap(new_genotype);
pd.m_info.swap(new_info);
+ LINEAGE_EXPR(pd.m_lineage.swap(new_lineage));
pd.m_inds.swap(new_inds);
#ifdef MUTANTALLELE
// compressed_vectora must be setGenoPtr after swap
@@ -2338,22 +2642,32 @@
#endif
// copy data over
GenoIterator newPtr = newGenotype.begin();
+#ifdef LINEAGE
+ vectori newLineage(genoSize() * m_popSize);
+ LineageIterator newLineagePtr = newLineage.begin();
+#endif
size_t pEnd = ploidy();
for (size_t i = 0; i < m_popSize; ++i) {
// set new geno structure
m_inds[i].setGenoStruIdx(genoStruIdx());
GenoIterator oldPtr = m_inds[i].genoPtr();
+ LINEAGE_EXPR(LineageIterator oldLineagePtr = m_inds[i].lineagePtr());
// new genotype
m_inds[i].setGenoPtr(newPtr);
+ LINEAGE_EXPR(m_inds[i].setLineagePtr(newLineagePtr));
for (size_t p = 0; p < pEnd; ++p) {
vectoru::iterator loc = kept.begin();
- for (; loc != kept.end(); ++loc)
+ for (; loc != kept.end(); ++loc) {
// this line needs ordered kept array
- *(newPtr++) = oldPtr[*loc]; //assginGenotype
+ *(newPtr++) = oldPtr[*loc]; //assignGenotype
+ LINEAGE_EXPR(*(newLineagePtr++) = oldLineagePtr[*loc]); //assignLineage
+ }
oldPtr += oldTotNumLoci;
+ LINEAGE_EXPR(oldLineagePtr += oldTotNumLoci);
}
}
m_genotype.swap(newGenotype);
+ LINEAGE_EXPR(m_lineage.swap(newLineage));
#ifdef MUTANTALLELE
// compressed_vectora must be setGenoPtr after swap
GenoIterator ptr = m_genotype.begin();
@@ -2557,6 +2871,7 @@
m_subPopIndex.swap(rhs.m_subPopIndex);
std::swap(m_vspSplitter, rhs.m_vspSplitter);
m_genotype.swap(rhs.m_genotype);
+ LINEAGE_EXPR(m_lineage.swap(rhs.m_lineage));
m_info.swap(rhs.m_info);
m_inds.swap(rhs.m_inds);
std::swap(m_indOrdered, rhs.m_indOrdered);
@@ -2827,6 +3142,7 @@
pd1.m_subPopSize.swap(pd.m_subPopSize);
pd1.m_subPopNames.swap(pd.m_subPopNames);
pd1.m_genotype.swap(pd.m_genotype);
+ LINEAGE_EXPR(pd1.m_lineage.swap(pd.m_lineage));
pd1.m_info.swap(pd.m_info);
pd1.m_inds.swap(pd.m_inds);
std::swap(pd1.m_indOrdered, pd.m_indOrdered);
@@ -2884,7 +3200,7 @@
//template<class Archive>
-void Population::save(boost::archive::text_oarchive & ar, const unsigned int /* version */) const
+void Population::save(boost::archive::text_oarchive & ar, const unsigned int version) const
{
// deep adjustment: everyone in order
const_cast<Population *>(this)->syncIndPointers();
@@ -2918,6 +3234,17 @@
#else
ar & m_genotype;
#endif
+ if (version > 0) {
+#ifdef LINEAGE
+ bool has_lineage = true;
+ DBG_DO(DBG_POPULATION, cerr << "Handling lineage" << endl);
+ ar & has_lineage;
+ ar & m_lineage;
+#else
+ bool has_lineage = false;
+ ar & has_lineage;
+#endif
+ }
DBG_DO(DBG_POPULATION, cerr << "Handling information" << endl);
ar & m_info;
DBG_DO(DBG_POPULATION, cerr << "Handling Individuals" << endl);
@@ -2949,6 +3276,17 @@
#else
ar & m_genotype;
#endif
+ if (version > 0) {
+#ifdef LINEAGE
+ bool has_lineage = true;
+ ar & has_lineage;
+ ar & m_lineage;
+#else
+ bool has_lineage = false;
+ ar & has_lineage;
+#endif
+ }
+
ar & m_info;
ar & m_inds;
}
@@ -2963,7 +3301,7 @@
//template<class Archive>
-void Population::load(boost::archive::text_iarchive & ar, const unsigned int /* version */)
+void Population::load(boost::archive::text_iarchive & ar, const unsigned int version)
{
size_t ma;
ar & ma;
@@ -3060,7 +3398,23 @@
# endif
}
#endif
-
+ if (version > 0) {
+ bool has_lineage;
+ ar & has_lineage;
+#ifdef LINEAGE
+ if (has_lineage) {
+ DBG_DO(DBG_POPULATION, cerr << "Handling lineage" << endl);
+ ar & m_lineage;
+ } else {
+ m_lineage.resize(m_popSize * genoSize());
+ }
+#else
+ if (has_lineage) {
+ vectori lineage;
+ ar & lineage;
+ }
+#endif
+ }
DBG_DO(DBG_POPULATION, cerr << "Handling info" << endl);
ar & m_info;
@@ -3097,6 +3451,14 @@
m_inds[i].setGenoPtr(ptr);
m_inds[i].setInfoPtr(infoPtr);
}
+
+#ifdef LINEAGE
+ LineageIterator lineagePtr = m_lineage.begin();
+ for(size_t i = 0; i < m_popSize; ++i, lineagePtr += step) {
+ m_inds[i].setLineagePtr(lineagePtr);
+ }
+#endif
+
m_ancestralGens = 0;
m_ancestralPops.clear();
@@ -3189,6 +3551,23 @@
# endif
}
#endif
+ if (version > 0) {
+ bool has_lineage;
+ ar & has_lineage;
+#ifdef LINEAGE
+ if (has_lineage) {
+ DBG_DO(DBG_POPULATION, cerr << "Handling lineage" << endl);
+ ar & pd.m_lineage;
+ } else {
+ pd.m_lineage.resize(pd.m_genotype.size());
+ }
+#else
+ if (has_lineage) {
+ vectori lineage;
+ ar & lineage;
+ }
+#endif
+ }
ar & pd.m_info;
ar & pd.m_inds;
// set pointer after copy this thing again (push_back)
@@ -3206,6 +3585,12 @@
// set new genoStructure
inds[i].setGenoStruIdx(genoStruIdx());
}
+#ifdef LINEAGE
+ lineagePtr = p.m_lineage.begin();
+ for (size_t i = 0; i < ps; ++i, lineagePtr += step) {
+ inds[i].setLineagePtr(lineagePtr);
+ }
+#endif
}
// load vars from string
@@ -3351,6 +3736,11 @@
vectora tmpGenotype(m_popSize * genoSize());
vectora::iterator it = tmpGenotype.begin();
#endif
+#ifdef LINEAGE
+ vectori tmpLineage(m_popSize * genoSize());
+ vectori::iterator lineagePtr = tmpLineage.begin();
+#endif
+
vectorf tmpInfo(m_popSize * infoSize());
vectorf::iterator infoPtr = tmpInfo.begin();
@@ -3365,9 +3755,11 @@
copy(ind->genoBegin(), ind->genoEnd(), it);
# endif
#endif
-
+ LINEAGE_EXPR(copy(ind->lineageBegin(), ind->lineageEnd(), lineagePtr));
ind->setGenoPtr(it);
+ LINEAGE_EXPR(ind->setLineagePtr(lineagePtr));
it += sz;
+ LINEAGE_EXPR(lineagePtr += sz);
copy(ind->infoBegin(), ind->infoEnd(), infoPtr);
ind->setInfoPtr(infoPtr);
infoPtr += is;
@@ -3375,6 +3767,7 @@
// discard original genotype
const_cast<Population *>(this)->m_genotype.swap(tmpGenotype);
const_cast<Population *>(this)->m_info.swap(tmpInfo);
+ LINEAGE_EXPR(const_cast<Population *>(this)->m_lineage.swap(tmpLineage));
}
setIndOrdered(true);
}
Modified: trunk/src/population.h
===================================================================
--- trunk/src/population.h 2012-03-07 03:27:16 UTC (rev 4435)
+++ trunk/src/population.h 2012-03-07 05:10:45 UTC (rev 4436)
@@ -285,6 +285,9 @@
m_subPopNames.swap(rhs.m_subPopNames);
m_subPopIndex.swap(rhs.m_subPopIndex);
m_genotype.swap(rhs.m_genotype);
+#ifdef LINEAGE
+ m_lineage.swap(rhs.m_lineage);
+#endif
m_info.swap(rhs.m_info);
m_inds.swap(rhs.m_inds);
std::swap(m_ancestralGens, rhs.m_ancestralGens);
@@ -832,6 +835,10 @@
/// CPPONLY allele begin, for given subPop
IndAlleleIterator alleleIterator(size_t locus, size_t subPop);
+#ifdef LINEAGE
+ /// CPPONLY lineage begin, for given subPop
+ IndLineageIterator lineageIterator(size_t locus, size_t subPop);
+#endif
/// CPPONLY allele iterator, go through all allels one by one, without subPop info
/**
@@ -861,7 +868,38 @@
return m_genotype.end();
}
+#ifdef LINEAGE
+ /// CPPONLY allele iterator, go through all lineages one by one, without subPop info
+ /**
+ if order, in order
+ otherwise, do not even respect subpopulation structure
+ */
+ LineageIterator lineageBegin(bool order)
+ {
+ DBG_FAILIF(hasActivatedVirtualSubPop(), ValueError,
+ "This function is not valid with an activated virtual subpopulation");
+
+ if (order && !indOrdered())
+ syncIndPointers();
+
+ return m_lineage.begin();
+ }
+
+
+ /// CPPONLY lineage iterator
+ LineageIterator lineageEnd(bool order)
+ {
+ DBG_FAILIF(hasActivatedVirtualSubPop(), ValueError,
+ "This function is not valid with an activated virtual subpopulation");
+ if (order && !indOrdered())
+ syncIndPointers();
+
+ return m_lineage.end();
+ }
+
+#endif // LINEAGE
+
/// CPPONLY allele iterator, go through all allels one by one in a subpopulation
/**
if order, keep order
@@ -890,7 +928,38 @@
return m_genotype.begin() + m_subPopIndex[subPop + 1] * genoSize();
}
+#ifdef LINEAGE
+ /// CPPONLY lineage iterator, go through all lineages one by one in a subpopulation
+ /**
+ if order, keep order
+ if not order, respect subpopulation structure
+ */
+ LineageIterator lineageBegin(size_t subPop, bool order)
+ {
+ DBG_FAILIF(hasActivatedVirtualSubPop(), ValueError,
+ "This function is not valid with an activated virtual subpopulation");
+ CHECKRANGESUBPOP(subPop);
+
+ syncIndPointers(order);
+
+ return m_lineage.begin() + m_subPopIndex[subPop] * genoSize();
+ }
+
+
+ /// CPPONLY lineage iterator in a subpopulation.
+ LineageIterator lineageEnd(size_t subPop, bool order)
+ {
+ DBG_FAILIF(hasActivatedVirtualSubPop(), ValueError,
+ "This function is not valid with an activated virtual subpopulation");
+ CHECKRANGESUBPOP(subPop);
+ syncIndPointers(order);
+
+ return m_lineage.begin() + m_subPopIndex[subPop + 1] * genoSize();
+ }
+
+#endif // LINEAGE
+
/// CPPONLY genoIterator --- beginning of individual ind.
GenoIterator indGenoBegin(size_t ind) const
{
@@ -906,7 +975,25 @@
return m_inds[ind].genoEnd();
}
+#ifdef LINEAGE
+ /// CPPONLY lineageIterator --- beginning of individual ind.
+ LineageIterator indLineageBegin(size_t ind) const
+ {
+ CHECKRANGEIND(ind);
+ return m_inds[ind].lineageBegin();
+ }
+
+
+ /// CPPONLY lineageIterator -- end of individual ind.
+ LineageIterator indLineageEnd(size_t ind) const
+ {
+ CHECKRANGEIND(ind);
+ return m_inds[ind].lineageEnd();
+ }
+
+#endif
+
/// CPPONLY genoIterator --- beginning of individual ind.
GenoIterator indGenoBegin(size_t ind, size_t subPop) const
{
@@ -926,7 +1013,29 @@
return m_inds[ subPopBegin(subPop) + ind].genoEnd();
}
+#ifdef LINEAGE
+ /// CPPONLY lineageIterator --- beginning of individual ind.
+ LineageIterator indLineageBegin(size_t ind, size_t subPop) const
+ {
+ CHECKRANGESUBPOP(subPop);
+ CHECKRANGESUBPOPMEMBER(ind, subPop);
+
+ return m_inds[ subPopBegin(subPop) + ind].lineageBegin();
+ }
+
+
+ /// CPPONLY genoIterator -- end of individual ind.
+ LineageIterator indLineageEnd(size_t ind, size_t subPop) const
+ {
+ CHECKRANGESUBPOP(subPop);
+ CHECKRANGESUBPOPMEMBER(ind, subPop);
+
+ return m_inds[ subPopBegin(subPop) + ind].lineageEnd();
+ }
+
+#endif
+
/** Return an editable array of the genotype of all individuals in
* a population (if <tt>subPop=[]</tt>, default), or individuals in a
* subpopulation \e subPop. Virtual subpopulation is unsupported.
@@ -934,6 +1043,14 @@
*/
PyObject * genotype(vspID subPop = vspID());
+ /** Return an editable array of the lineage of alleles for all individuals in
+ * a population (if <tt>subPop=[]</tt>, default), or individuals in a
+ * subpopulation \e subPop. Virtual subpopulation is unsupported. <bf>
+ * This function returns \c None for modules without lineage information.</bf>
+ * <group>5-genotype</group>
+ */
+ PyObject * lineage(vspID subPop = vspID());
+
/** Fill the genotype of all individuals in a population (if
* <tt>subPop=[]</tt>) or in a (virtual) subpopulation \e subPop (if
* <tt>subPop=sp</tt> or <tt>(sp, vsp)</tt>) using a list of alleles
@@ -943,6 +1060,16 @@
*/
void setGenotype(const uintList & geno, vspID subPop = vspID());
+ /** Fill the lineage of all individuals in a population (if
+ * <tt>subPop=[]</tt>) or in a (virtual) subpopulation \e subPop (if
+ * <tt>subPop=sp</tt> or <tt>(sp, vsp)</tt>) using a list of IDs
+ * \e lineage. \e lineage will be reused if its length is less than
+ * <tt>subPopSize(subPop)*totNumLoci()*ploidy()</tt>. This function
+ * returns directly for modules without lineage information.
+ * <group>5-genotype</group>
+ */
+ void setLineage(const uintList & geno, vspID subPop = vspID());
+
//@}
/** @name utility functions, set subpopulation, save and load etc.
@@ -1505,6 +1632,10 @@
vectora m_genotype;
#endif
+#ifdef LINEAGE
+ vectori m_lineage;
+#endif
+
/// information
/// only in head node
vectorf m_info;
@@ -1529,6 +1660,11 @@
#else
vectora m_genotype;
#endif
+
+#ifdef LINEAGE
+ vectori m_lineage;
+#endif
+
vectorf m_info;
vector<Individual> m_inds;
bool m_indOrdered;
@@ -1604,7 +1740,8 @@
#ifndef SWIG
# ifndef _NO_SERIALIZATION_
// version 0: base (reset for version 1.0)
-BOOST_CLASS_VERSION(simuPOP::Population, 0)
+// version 1: with lineage information for lineage-aware modules
+BOOST_CLASS_VERSION(simuPOP::Population, 1)
# endif
#endif
#endif
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <si...@us...> - 2012-03-09 04:32:16
|
Revision: 4449
http://simupop.svn.sourceforge.net/simupop/?rev=4449&view=rev
Author: simupop
Date: 2012-03-09 04:32:09 +0000 (Fri, 09 Mar 2012)
Log Message:
-----------
Windows VC++ compatibility patch for the lineage modules.
Modified Paths:
--------------
trunk/src/individual.h
trunk/src/initializer.cpp
trunk/src/mutator.cpp
Modified: trunk/src/individual.h
===================================================================
--- trunk/src/individual.h 2012-03-09 04:02:22 UTC (rev 4448)
+++ trunk/src/individual.h 2012-03-09 04:32:09 UTC (rev 4449)
@@ -1346,7 +1346,10 @@
// NOTE: this iterator is used only when indOrdered() is set to true for
// the whole population (see Population.lineageIterator()). It is therefore
// possible to get the index of individual from index of m_ptr.
- int offset = (m_ptr - m_ptrBegin) / (m_size * m_ploidy);
+ //
+ // There is a conversion from size_t to long (difference_type), a possible
+ // loss of data
+ difference_type offset = static_cast<difference_type>((m_ptr - m_ptrBegin) / (m_size * m_ploidy));
return(m_it + offset);
} else {
return(m_it);
Modified: trunk/src/initializer.cpp
===================================================================
--- trunk/src/initializer.cpp 2012-03-09 04:02:22 UTC (rev 4448)
+++ trunk/src/initializer.cpp 2012-03-09 04:32:09 UTC (rev 4449)
@@ -317,10 +317,10 @@
IndIterator it = pop.indIterator(sp->subPop());
#endif
if (pop.hasInfoField(m_lineageField)) {
- int idIdx = pop.infoIdx(m_lineageField);
+ size_t idIdx = pop.infoIdx(m_lineageField);
for (; it.valid(); ++it) {
- long lineage = toID(it->info(idIdx));
+ long lineage = static_cast<long>(toID(it->info(idIdx)));
for (vectoru::iterator p = ploidy.begin(); p != ploidy.end(); ++p)
for (vectoru::const_iterator loc = loci.begin(); loc != loci.end(); ++loc)
it->setAlleleLineage(lineage, *loc, static_cast<int>(*p));
Modified: trunk/src/mutator.cpp
===================================================================
--- trunk/src/mutator.cpp 2012-03-09 04:02:22 UTC (rev 4448)
+++ trunk/src/mutator.cpp 2012-03-09 04:32:09 UTC (rev 4449)
@@ -95,7 +95,7 @@
#ifdef LINEAGE
bool assignLineage = pop.hasInfoField(m_lineageField);
- int lineageIdx = assignLineage ? pop.infoIdx(m_lineageField) : 0;
+ size_t lineageIdx = assignLineage ? pop.infoIdx(m_lineageField) : 0;
#endif
// mapIn and mapOut
@@ -154,7 +154,7 @@
long lineage = 0;
if (assignLineage) {
lineagePtr += static_cast<IndLineageIterator::difference_type>(pos - lastPos);
- lineage = toID(lineagePtr.individual()->info(lineageIdx));
+ lineage = static_cast<long>(toID(lineagePtr.individual()->info(lineageIdx)));
}
#endif
ptr += static_cast<IndAlleleIterator::difference_type>(pos - lastPos);
@@ -450,7 +450,7 @@
#ifdef LINEAGE
bool assignLineage = pop.hasInfoField(m_lineageField);
- int lineageIdx = assignLineage ? pop.infoIdx(m_lineageField) : 0;
+ size_t lineageIdx = assignLineage ? pop.infoIdx(m_lineageField) : 0;
#endif
for (; sp != spEnd; ++sp) {
@@ -463,7 +463,7 @@
IndIterator ind = pop.indIterator(sp->subPop()) + static_cast<IndIterator::difference_type>(*it);
if (!ind.valid())
continue;
- LINEAGE_EXPR(long lineage = assignLineage ? toID(ind->info(lineageIdx)) : 0);
+ LINEAGE_EXPR(long lineage = assignLineage ? static_cast<long>(toID(ind->info(lineageIdx))) : 0);
for (size_t p = 0; p < m_ploidy.size(); ++p) {
if (m_loci.allAvail()) {
for (size_t i = 0; i < pop.totNumLoci(); ++i) {
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
Thanks for helping keep SourceForge clean.
X