Life-long sexual monogamy

20 views
Skip to first unread message

Kamolphat Atsawawaranunt

unread,
May 8, 2025, 12:28:22 AMMay 8
to slim-discuss
Hi Ben,

I have been trying to model life-long sexual monogamous species. Individuals find new mates if the partner dies. I think I have successfully made it run with the help of the individual.tag and the sortBy function. In the manual, it says that "SLiM does not - for reasons of speed - guarantee a random order to the individuals in the subpopulation". When the sortBy function is used, is the order of the individuals "random" when the "tag" is of the same value? Here I applied sortBy to the array of individuals after I have used the sample function.

My SLiM code is run in SLiM 4.3 and is as follows:

initialize() {
initializeSLiMModelType("nonWF");
initializeSex("A");
defineConstant("K", 500); // carrying capacity
initializeMutationType("m1", 0.5, "f", 0.0);
m1.convertToSubstitution = T;
initializeGenomicElementType("g1", m1, 1.0);
initializeGenomicElement(g1, 0, 99999);
initializeMutationRate(1e-7);
initializeRecombinationRate(1e-8);
}

reproduction() {
// Find females and males
female = which(p1.individuals.sex == "F");
male = which(p1.individuals.sex == "M");
npairs = min(c(length(female), length(male)));
fparents = p1.sampleIndividuals(length(female), replace = F, sex = "F");
mparents = p1.sampleIndividuals(length(male), replace = F, sex = "M");
fparents = sortBy(fparents, "tag", ascending = F); // Sort by tags
mparents = sortBy(mparents, "tag", ascending = F); // Sort by tags
// draw monogamous pairs and generate litters
for (i in seqLen(npairs))
{
parent1 = fparents[i];
parent2 = mparents[i];
parent1.tag = i + 1;
parent2.tag = i + 1;
// To alter the number of offsprings by the cached fitness
offspring = p1.addCrossed(parent1, parent2, count = rpois(1, 2.7));
// New offsprings have tag of zero
offspring.tag = 0;
}
// disable this callback for this cycle
self.active = 0;
}

1 early() {
sim.addSubpop("p1", 10);
p1.individuals.tag = 0; // Add empty tags
}

early() {
p1.fitnessScaling = K / p1.individualCount;
}

late() {
// Get parent tags
parentTag = p1.individuals.tag[which(p1.individuals.tag != 0)];
// For each parent tag
for (t in unique(parentTag)){
// if only one parent is left, reset the tag to zero
if (length(which(p1.individuals.tag == t)) < 2){
ptag = p1.individuals[which(p1.individuals.tag == t)];
ptag.tag = 0;
}
}
}

3000 late() {
sim.outputFixedMutations(); }


Yours sincerely,
A (Kamolphat Atsawawaranunt)

Thank you so much and I will try to move to SLiM 5.0 soon! :)

Ben Haller

unread,
May 8, 2025, 6:56:20 AMMay 8
to Kamolphat Atsawawaranunt, slim-discuss
Hi A!

Interesting algorithm, I like the idea.  It seems like you need to treat the tag==0 individuals separately in your reproduction callback, though, right?  There might not be the same number of tag==0 males and tag==0 females, especially after you sample (but even before you sample, actually), so that means your sorted vectors won't actually correspond one-to-one with each other correctly, doesn't it?  Perhaps I'm missing something, but I don't see how your algorithm handles that problem.

Once you've separated out the tag==0 individuals for separate treatment, the remainder – the mated pairs – should just have one male and one female for any given tag value, so in that context your question doesn't matter any more, right?  It would only matter for the tag==0 individuals (of which there will often be more than one male and more than one female); but since (it seems to me) you need to treat them separately anyway, in a different way, I think your question doesn't matter there either, perhaps?

But to answer your question: no, Eidos makes no guarantee about the randomness or otherwise of a sorted vector, in the case where the sort key is identical between more than one element.  There is no guarantee as to whether the sort is "stable" or not, in computer science terms; but even if it isn't stable, the order of elements with identical sort keys will not be random in any reliable way.  (They certainly aren't deliberately randomized with draws from a random number generator; they might be scrambled in some manner by the underlying algorithm, but that is not reliably random and would probably exhibit sort-order artifacts.)

A few notes on your code, while I'm here:

- You do "which(p1.individuals.sex == "F")" to select the females, but seem to use that only to count the number of females (and likewise for the males); and then when you want to draw females you call p1.sampleIndividuals() to do it.  It would probably be much more efficient to get the vector of females initially with "females = p1.subsetIndividuals(sex='F')", and then scramble their order with "fparents = sample(females, length(females), replace='F')".

- But since the next thing you do is sort by tag, the scrambling is unnecessary anyway, if you exclude the tag==0 individuals from the process first.  (Which I think you need to do anyway, as I discussed above, because there won't be an equal number of male and female tag==0 individuals, as far as I can see.)

- The way that you call "p1.individuals.tag" over and over is quite inefficient; every time you do that, Eidos has to fetch p1, then fetch every individual from it (p1.individuals), then fetch every tag value from those individuals (p1.individuals.tag).  If you can fetch that once and cache it in a local variable, it's much more efficient.  Of course this only matters if you need to increase the performance of your code.

- You lean on "which()" more than I'd recommend for efficiency.  For example, "length(which(p1.individuals.tag == t))" to get the number of individuals with tag t is better done as "sum(p1.individuals.tag == t)"; and "p1.individuals[which(p1.individuals.tag == t)]" is better done as "p1.subsetIndividuals(tag=t)".  (In fact, the count can be obtained with "size(p1.subsetIndividuals(tag=t))" too, which might be even better/faster for your particular algorithm since then you've already got the subset of the individuals that you need in hand, after getting the size.)

- Beyond the comments above, the unmatched parents algorithm can perhaps be optimized further, if it's a bottleneck (profile in SLiMgui to find out).  One simple thing is that you should pass "preserveOrder=F" to the call to unique(); that makes it much faster, and you don't need the order preserved for your algorithm.  That converts it from an O(N^2) algorithm to an O(N log N) algorithm.  But even faster, for large N, would probably be to fetch p1.individuals.tag once and then use tabulate() to determine which tag values have a count of 1; that is, in effect, a radix sort and so would be O(N), I think.  For large N if there are many unmatched parents then the lookup of the particular parents that are unmatched might still make things slow; I think you could speed that part of the algorithm up using match() to give you the indices of the individuals that have the unmatched tag values, making that part of the algorithm, also, O(N).  So if you plan to scale this model up to large N, an approach like that would perhaps be important.  You could also use match() instead of sortBy() to find the matching male mate for each tag != 0 female, again converting from an O(N log N) algorithm to an O(N) algorithm.  But if performance is not a problem for you so far, or if profiling doesn't show these algorithms to be bottlenecks, then don't worry about it.  :->

It'd be nice to add a recipe to the manual demonstrating how to model life-long monogamy like this; this comes up as a question from time to time.  I'm feeling inspired to do that by the approach you've taken here, which has given me some ideas as you can see :->.  Would it be OK with you if I adapted your code into a recipe and added it to the manual, with a credit to you for the idea?  :->

Good luck, and happy modeling!

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University



Kamolphat Atsawawaranunt wrote on 5/8/25 12:28 AM:
--
SLiM forward genetic simulation: http://messerlab.org/slim/
---
You received this message because you are subscribed to the Google Groups "slim-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to slim-discuss...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/slim-discuss/a6a0c864-c833-4024-bcd7-96fc4ecc333dn%40googlegroups.com.

Ben Haller

unread,
May 8, 2025, 7:00:35 PMMay 8
to slim-discuss
A quick followup on this thread.  I've written a new recipe for the manual, demonstrating life-long monogamous mating.  It's been pushed to GitHub just now, and will appear in the next version of SLiM, whenever that happens.  It ended up being fairly different from A's model in the original post, and uses match() and tabulate() for efficiency as I suggested previously; but some ideas did carry over in spirit, at least.  So I wrote "Thanks to Kamolphat Atsawawaranunt for the inspiration behind this recipe."  :->  Anyone who wants a copy of the recipe and the write-up can email me off-list (do NOT reply to the whole list, please!).  Thanks A!

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

Message has been deleted

Kamolphat Atsawawaranunt

unread,
May 9, 2025, 4:11:20 PMMay 9
to slim-discuss
Hi Ben,

Thank you so much for your reply and all the suggestions and the tips for optimisation. I am new to SLiM (I attended the SLiM workshop in Auckland a month or so ago) and your suggestions for optimisation has already improved the speed greatly. And yes, you are very welcome to adapt my code into a recipe. I just noticed that you have already added a new recipe to the github page as I was writing this reply so I thought I would write it here anyways in case there are some ideas. The recipe looks great with the minimum and maximum reproduction age and also mean fecundity. I will definitely adopt some of the principles you used in the recipe to my own codes. I tried running the recipe (on SLiM 4.3, adding "A" to initalizeSex()) and played around with a couple of parameters and ran into some problems when I alter FECUN to values >= 1.4 and R_AGE_M and R_AGE_F to zero. The error says:

Simulation Runtime Error
ERROR (SubsetEidosValue): out-of-range index -1 used with the '[]' operator.
This error has invalidated the simulation; it cannot be run further. Once the script is fixed, you can recycle the simulation and try again.

Regarding how I treated the tag==0 individuals separately in my original codes, I treated by mating the paired individuals from previous generations first as organised by the sortBy() function and only looped until "npairs" individuals. That way, the loop only reaches the minimum number of male or female individuals, with a few remainders which will not have a partner and not mated. The reason I did it this way was to avoid having to treat the tag==0 individuals separately, hoping to avoid having two for-loops. Following your suggestion, I created a list of tag>0 individuals and sorted them, and created a list of tag==0 individuals, which I then concatenate to one long list to avoid two for-loops (or if-statements). Currently, I am having some trouble selecting for the tag>0 individuals efficiently. I currently use three lines of codes per sex (first six lines of the reproduction callback) and I felt that it can be more efficient. I could not use the subsetIndividuals() function to select individuals with tag>0. Is there a way to make this more efficiently?

It does appear that the line with the unique() function is indeed a bottleneck, so I have followed your suggestion to use tabulate() and it sped up the codes significantly. However, I am still unsure if I have done this most efficiently (I am still using one which statement in that late() callback).

My current codes:

initialize() {
initializeSLiMModelType("nonWF");
initializeSex("A");
defineConstant("K", 500); // carrying capacity
initializeMutationType("m1", 0.5, "f", 0.0);
m1.convertToSubstitution = T;
initializeGenomicElementType("g1", m1, 1.0);
initializeGenomicElement(g1, 0, 99999);
initializeMutationRate(1e-7);
initializeRecombinationRate(1e-8);
}

reproduction() {
// Get all males and females
females = p1.subsetIndividuals(sex='F');
males = p1.subsetIndividuals(sex='M');
// tag!=0 individuals
femalesnotzero = females[which(females.tag > 0)];
malesnotzero = males[which(males.tag > 0)];
femalesnotzero = sortBy(femalesnotzero, "tag");
malesnotzero = sortBy(malesnotzero, "tag");
// tag==0 individuals
femaleszero = p1.subsetIndividuals(sex='F', tag = 0);
maleszero = p1.subsetIndividuals(sex='M', tag = 0);
// Shuffle the tag==0 individuals
femaleszero = sample(females, size(femaleszero), replace=F);
maleszero = sample(males, size(maleszero), replace=F);
// Concatenate the tag>0 individuals and tag==0 individuals
fparents = c(femalesnotzero, femaleszero);
mparents = c(malesnotzero, maleszero);
// Get minimum number of males or females for looping
npairs = min(c(size(females), size(males)));

// draw monogamous pairs and generate litters
for (i in seqLen(npairs))
{
parent1 = fparents[i];
parent2 = mparents[i];
parent1.tag = i + 1;
parent2.tag = i + 1;
offspring = p1.addCrossed(parent1, parent2, count = rpois(1, 2.0));

offspring.tag = 0;
}
// disable this callback for this cycle
self.active = 0;
}

1 first() {

sim.addSubpop("p1", 10);
p1.individuals.tag = 0; // Add empty tags
}

early() {
p1.fitnessScaling = K / p1.individualCount;
}

late() {
// Reset the tags
// Get tags
ptags = p1.individuals.tag;
// tabulate tags
ptags_tab = tabulate(ptags);
// find the ones that are unmatched, i.e., partner dies
parentTag = which(ptags_tab == 1);
// get index
unmatchedp = match(parentTag, ptags);
// subset widowed individuals
ptag1 = p1.individuals[unmatchedp];
// Reset tag to zero
ptag1.tag = 0;
}

3000 late() {
sim.outputFixedMutations(); }



Thank you so much for all your suggestions. They really improved the performance significantly. I will try to adopt the rpois() function outside of the loop in my next iterations of codes.

Yours sincerely,
A

Ben Haller

unread,
May 9, 2025, 4:57:07 PMMay 9
to Kamolphat Atsawawaranunt, slim-discuss
Hi A!

You found a bug in the recipe, thanks!  The line that reads:

    tag_counts = tabulate(tags, maxbin=size(mated_individuals)-1);

should instead simply read:

    tag_counts = tabulate(tags);

I was imposing an incorrect max bin on the tabulation because I thought about the problem a little bit wrong.  :-O  I've pushed the new recipe to GitHub.

Regarding your question: "Currently, I am having some trouble selecting for the tag>0 individuals efficiently."  Your code now reads:


    females = p1.subsetIndividuals(sex='F');
    femalesnotzero = females[which(females.tag > 0)];

You are correct that subsetIndividuals() does not currently support minTag= and maxTag= options; that would be a nice addition.  The best you can do for now is:


    females = p1.subsetIndividuals(sex='F');
    femalesnotzero = females[females.tag > 0];

As I wrote in my earlier comments, it seems like you have a tendency to overuse the which() function.  By using which(), you are converting from a logical vector (the result of females.tag > 0) to an integer vector of indices where the logical vector was true.  But there's no need to do that; you can simply subset females with the logical vector directly, which is faster.  You might review the "Eidos Overview" worksheet, which covers this material in the workshop.  :->  It's still a little slower than it would be if subsetIndividuals() supported min/max tag values; if you decide you'd like to see that feature, please open a new issue on GitHub for it.  :->

Regarding the rest of your questions, I'd suggest that you review the new recipe in the manual.  I will email you a PDF of the writeup for it, off-list, momentarily.


Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University


Kamolphat Atsawawaranunt wrote on 5/8/25 8:07 PM:
Reply all
Reply to author
Forward
0 new messages