split.slim

initialize() {
  defineConstant('tmu', 1e-8);
  defineConstant('nmu', 1e-8);
  defineConstant('rbp', 1e-8);
  defineConstant('N', 1000);
  defineConstant('alpha', 0.01);
  defineConstant('nrep', 1);
      
  defineConstant("seed", getSeed());
  defineConstant("data_dir", '../data/sims/');
  defineConstant("region_length", 50e6);

  // Stabilizing selection variance
  defineConstant("Vs", 1);

  // The SLiM uses a fitness offset, where fitness is offset + z_i where z_i =
  // Σ α_i g_i.  In contrast, most GSS models don't do this. A fitness offset
  // of 1 buffers the differences in breeding values, e.g. z_1 = 0.01 and z_2 =
  // 0.03 relative fitness are 1 and 1/3 --- fairly large fitness differences,
  // versus 1 and 0.98 if we have a fitness offset of one.
  defineConstant("fitness_offset", 0.0);

  initializeMutationRate(nmu + tmu);
  initializeMutationType("m1", 0.5, "s", "if (runif(1) < 0.5) -alpha; else alpha;");
  initializeMutationType("m2", 0.5, "f", 0.0);
  initializeGenomicElementType("g1", c(m1, m2), c(tmu, nmu));
  initializeGenomicElement(g1, 0, region_length-1);
  initializeRecombinationRate(rbp);
  m1.convertToSubstitution = T;
  m1.mutationStackPolicy = "f";
  m2.convertToSubstitution = T;
  m2.mutationStackPolicy = "f";

  // initialize output file headers
  bname = ("gss_burnin_" + N + "N_" + rbp + "rbp_" +
            alpha + "alpha_" + nmu + "nmu_" +
            tmu + "tmu_" + nrep);
  defineConstant("basename", bname);
}


fitness(m1) {
    // fitness is assigned based on phenotypic value through global
    // fitness callback; this just multiplies each fitness by 1.
    return 1.0;
}

fitness(NULL) {
  return 1.0;
}

s3 fitness(NULL) {
    return fitness_offset + dnorm(individual.tagF, mean=0.0, sd=sqrt(Vs));
}


// --- blocks
1 early() {
    burnin = 10*N;
    end = burnin + 100;
    sim.rescheduleScriptBlock(s1, start=burnin-1, end=end);
    sim.rescheduleScriptBlock(s2, start=end, end=end);
    sim.rescheduleScriptBlock(s3, start=burnin-1, end=end);

}

1 late() {
	sim.readFromPopulationFile("/Users/vinceb/tmp/gss_burnin_1000N_1e-08rbp_0.01alpha_1e-08nmu_1e-08tmu_12_fullsims.txt");
    sim.addSubpopSplit("p2", 100, p1);
    sim.addSubpopSplit("p3", 100, p1);
    p1.setSubpopulationSize(0);
     
}
    


s1 early() {
    inds = sim.subpopulations.individuals;
    phenotypes = inds.sumOfMutationsOfType(m1);
    fixed_trait = sum(sim.substitutions.selectionCoeff);
    inds.tagF = phenotypes;
    print("generation: " + sim.generation);
}

	
s2 late() {
  sim.simulationFinished();
}