Restarting simulations in SLiM engine (community.allScriptBlocks.active property)

[jmurga]

Hi all, first, I want to congratulate you on such a fantastic project you're working on.

I've been using Stdpopsim to automatize sweep simulations given the time (t) at which beneficial mutation arises, the final frequency (f), and the selection coefficient (s) range that can produce such a sweep. In my approach, I solved the Otto and Whitelock 2013 Eq. 5 to approximate the time of fixation or the desired final frequency to obtain proper t and s combinations producing the expected sweep. Because SLiM + Stdpopsim restarts the simulation as much as needed, such an approximation works well (I've tested constant population size and OutOfAfrica_2T12 demography models).

I was using SLiM 4.0.1, and everything worked fine. Today, I upgraded SLiM to version 4.2.1, and the simulations don't work. I dumped the SLiM recipe and tested it in different environments in addition to trying to debug it in SLiMgui. Nonetheless, it does not work on SLiM 4.2.1 but still works on 4.0.1.

The error I obtained is the following one:

ERROR (EidosPropertySignature::CheckAssignedValue): value cannot be type logical for read-write property active.  
Error on script line 25, character 33 (inside runtime script block):  
    sb[sb.type == "late"].active = F;  
                                 ^

I checked the SLiM manual, and now the community.allScriptBlocks.active property is defined as an integer. The problem should be related to assigning a boolean value instead of an integer.

I changed the line sb[sb.type == "late"].active = F; to sb[sb.type == "late"].active = 0;` and tested in SLiM and SLiMgui, and it works now as expected. There is no error, and I tracked the mutation frequency over time using the SLiMgui, and the beneficial mutation reached the proper frequency.

How to reproduce the problem:
I am attaching a piece of code to simulate a sweep with constant population size and uniform recombination rate in a 1.2MB region. Timing, selection coefficient, and initial and final frequency are defined to produce an ongoing sweep that started ~1300 generations ago and reaches a final frequency in the range [0.8,0.9].

import stdpopsim
import msprime
import numpy as np
from contextlib import redirect_stdout

slim_script = "/home/jmurgamoreno/test.slim"
# selection coefficient, time to beneficial, end time to partial sweep if needed, burning, iteration number
(_s, _t, _t_end, _f_i, _f_t, _burn_in) = (
    0.0119,
    1384,
    0,
    0.0001,
    0.8,
    0.15000000000000002,
)


# Retrieve species information
species = stdpopsim.get_species("HomSap")
model = stdpopsim.PiecewiseConstantSize(species.population_size)
model.generation_time = species.generation_time
# Force to sample to pop0 is PiecewiseConstantSize selected
sample = {"pop_0": 50}

contig = species.get_contig("chr2", mutation_rate=model.mutation_rate, right=1.2e6)

# Change manually to test og Flex-sweep
contig.mutation_rate = np.random.uniform(2 * 1e-9, 5.2 * 1e-8)
contig.recombination_map = msprime.RateMap(
    position=contig.recombination_map.position,
    rate=np.array([np.random.exponential(1e-8)]),
)

# adaptive_position = np.random.randint(1.2e6*0.25,1.2e6*0.75)
adaptive_position = sum(contig.original_coordinates[1:]) / 2

contig.add_single_site(id="hard", coordinate=adaptive_position)


extended_events = stdpopsim.ext.selective_sweep(
    single_site_id="hard",
    population=list(sample.keys())[0],
    selection_coeff=_s,
    min_freq_at_end=_f_t,
    mutation_generation_ago=_t,
    end_generation_ago=_t_end,
)

if _f_t != 1:
    end_condition = stdpopsim.ext.ConditionOnAlleleFrequency(
        start_time=_t_end,
        end_time=_t_end,
        single_site_id="hard",
        population=list(sample.keys())[0],
        op="<",
        allele_frequency=_f_t + 0.1,
    )

    extended_events.append(end_condition)

engine = stdpopsim.get_engine("slim")


with open(slim_script, "w") as f:
    with redirect_stdout(f):
        ts_sweep = engine.simulate(
            model,
            contig,
            sample,
            extended_events=extended_events,
            slim_scaling_factor=1,
            slim_burn_in=_burn_in,
            slim_script=True,
            verbosity=0,
        )

ts_sweep = engine.simulate(
    model,
    contig,
    sample,
    extended_events=extended_events,
    slim_scaling_factor=1,
    slim_burn_in=_burn_in,
)

Suggested fixes:
The community.allScriptBlocks.active property seems unable to deal with boolean values now. Using 0 values as False solves the problem. So, changing line 219 on slim_engine.py is a potential solution in case I am not doing anything wrong.

Best,

Jesús Murga-Moreno

[nspope]

Thanks for the very thorough bug report @jmurga!

@petrelharp do you see any potential issues with the solution above? e.g. do we need to handle back-compatibility with earlier slim versions?

[petrelharp]

Looking at the help (in SLiMgui), the documentation says:

If this evaluates to logical F (i.e., is equal to 0), the script block is inactive and will not be called.  The value of active for all registered script blocks is reset to -1 at the beginning of each tick, prior to script events being called, thus activating all blocks (except callbacks associated with a species that is not active in that tick, which are deactivated as part of the deactivation of the species).  Any integer value other than -1 may be used instead of -1 to represent that a block is active; for example, active may be used as a counter to make a block execute a fixed number of times in each tick.  This value is not cached by SLiM; if it is changed, the new value takes effect immediately.  For example, a callback might be activated and inactivated repeatedly during a single tick.

... so, it looks like the type got changed to allow more behavior (in particular, using -1 for something). This might actually be a bit of unintended backwards incompatibility? FYI, @bhaller.

But - I don't see any problem with this solution, no. Thanks very much for the report!

[bhaller]

Looking at the help (in SLiMgui), the documentation says:

If this evaluates to logical F (i.e., is equal to 0), the script block is inactive and will not be called.  The value of active for all registered script blocks is reset to -1 at the beginning of each tick, prior to script events being called, thus activating all blocks (except callbacks associated with a species that is not active in that tick, which are deactivated as part of the deactivation of the species).  Any integer value other than -1 may be used instead of -1 to represent that a block is active; for example, active may be used as a counter to make a block execute a fixed number of times in each tick.  This value is not cached by SLiM; if it is changed, the new value takes effect immediately.  For example, a callback might be activated and inactivated repeatedly during a single tick.

... so, it looks like the type got changed to allow more behavior (in particular, using -1 for something). This might actually be a bit of unintended backwards incompatibility? FYI, @bhaller.

But - I don't see any problem with this solution, no. Thanks very much for the report!

From the 4.2 release notes:

Assignment into a property of an object must now match the type of the property; type promotion will no longer occur. Existing code that now errors should use asInteger(), asFloat(), etc. as needed.

It broke a few patterns like this; but active has actually always been type integer, not type logical. Tightening up implicit type coercion in Eidos a bit. Sorry!