Skip to content

Commit

Permalink
Streamline timeseries printing.
Browse files Browse the repository at this point in the history
  • Loading branch information
@trvrb
trvrb committed Feb 22, 2014
1 parent a6026f0 commit fd3e07d
Show file tree
Hide file tree
Showing 4 changed files with 90 additions and 87 deletions.
100 changes: 57 additions & 43 deletions HostPopulation.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,19 +8,23 @@ public class HostPopulation {

// fields
private int deme;
private String name;
private int cases;
private List<Host> susceptibles = new ArrayList<Host>();
private List<Host> infecteds = new ArrayList<Host>();
private List<Host> recovereds = new ArrayList<Host>(); // this is the transcendental class, immune to all forms of virus
private double diversity;
private double tmrca;
private double netau;
private double serialInterval;
private double antigenicDiversity;

// construct population, using Virus v as initial infection
public HostPopulation(int d) {

// basic parameters
deme = d;
name = Parameters.demeNames[deme];
int initialR = 0;
if (Parameters.transcendental) {
initialR = (int) ((double) Parameters.initialNs[deme] * Parameters.initialPrT);
Expand Down Expand Up @@ -61,6 +65,7 @@ public HostPopulation(int d, boolean checkpoint) {
if (checkpoint == true) {

deme = d;
name = Parameters.demeNames[deme];

try {
BufferedReader in = new BufferedReader(new FileReader("out.hosts"));
Expand Down Expand Up @@ -182,14 +187,22 @@ public int getCases() {

public double getDiversity() {
return diversity;
}

public double getNetau() {
return netau;
}

public double getTmrca() {
return tmrca;
}

public double getSerialInterval() {
return serialInterval;
}

public double getNetau() {
return netau;
public double getAntigenicDiversity() {
return antigenicDiversity;
}

public void removeSusceptible(int i) {
Expand Down Expand Up @@ -453,54 +466,55 @@ public void makeTrunk() {
}

public void updateDiversity() {

diversity = 0.0;
tmrca = 0.0;
int sampleCount = 0;
for (int i = 0; i < Parameters.diversitySamplingCount; i++) {
Virus vA = getRandomInfection();
Virus vB = getRandomInfection();
if (vA != null && vB != null) {
double dist = vA.distance(vB);
diversity += dist;
if (dist > tmrca) {
tmrca = dist;
}
sampleCount += 1;
}
}
if (sampleCount > 0) {
tmrca = 0.0;
antigenicDiversity = 0.0;
netau = 0.0;
serialInterval = 0.0;

if (getI()>1) {

double coalCount = 0.0;
double coalOpp = 0.0;
double coalWindow = Parameters.netauWindow / 365.0;
int sampleCount = Parameters.diversitySamplingCount;

for (int i = 0; i < sampleCount; i++) {
Virus vA = getRandomInfection();
Virus vB = getRandomInfection();
if (vA != null && vB != null) {
double dist = vA.distance(vB);
diversity += dist;
if (dist > tmrca) {
tmrca = dist;
}
antigenicDiversity += vA.antigenicDistance(vB);
coalOpp += coalWindow;
coalCount += vA.coalescence(vB, coalWindow);
serialInterval += vA.serialInterval();
}
}

diversity /= (double) sampleCount;
}
tmrca /= 2.0;
}

public void updateNetau() {
double coalCount = 0.0;
double coalOpp = 0.0;
double coalWindow = Parameters.netauWindow / 365.0;
int sampleCount = 0;
for (int i = 0; i < Parameters.netauSamplingCount; i++) {
Virus vA = getRandomInfection();
Virus vB = getRandomInfection();
if (vA != null && vB != null) {
coalOpp += coalWindow;
coalCount += vA.coalescence(vB, coalWindow);
sampleCount += 1;
}
}
if (sampleCount > 0) {
tmrca /= 2.0;
antigenicDiversity /= (double) sampleCount;
netau = coalOpp / coalCount;
serialInterval /= (double) sampleCount;

}
}


}

public void printState(PrintStream stream) {
if (Parameters.day > Parameters.burnin) {
updateDiversity();
updateNetau();
stream.printf("\t%.4f\t%.4f\t%.4f\t%d\t%d\t%d\t%d\t%d", getDiversity(), getTmrca(), getNetau(), getN(), getS(), getI(), getR(), getCases());
}
updateDiversity();
stream.printf("\t%.4f\t%.4f\t%.4f\t%.5f\t%.4f\t%d\t%d\t%d\t%d\t%d", getDiversity(), getTmrca(), getNetau(), getSerialInterval(), getAntigenicDiversity(), getN(), getS(), getI(), getR(), getCases());
}

public void printHeader(PrintStream stream) {
stream.printf("\t%sDiv\t%sTmrca\t%sNetau\t%sSerialInterval\t%sAntigenicDiv\t%sN\t%sS\t%sI\t%sR\t%sCases", name, name, name, name, name, name, name, name, name, name);
}

// reset population to factory condition
public void reset() {

Expand Down
10 changes: 3 additions & 7 deletions Parameters.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,10 +21,8 @@ public class Parameters {
public static int tipSamplesPerDeme = 1000;
public static boolean tipSamplingProportional = true; // whether to sample proportional to prevalance
public static double treeProportion = 0.1; // proportion of tips to use in tree reconstruction
public static int diversitySamplingCount = 1000; // how many samples to draw to calculate diversity
public static int netauSamplingCount = 1000; // how many samples to draw to calculate Ne*tau
public static int netauWindow = 100; // window in days to calculate Ne*tau
public static int serialIntervalSamplingCount = 1000; // how many samples to draw to calculate serial interval
public static int diversitySamplingCount = 1000; // how many samples to draw to calculate diversity, Ne*tau, serial interval
public static int netauWindow = 100; // window in days to calculate Ne*tau
public static boolean repeatSim = true; // repeat simulation until endDay is reached?
public static boolean immunityReconstruction = false; // whether to print immunity reconstruction to out.immunity
public static boolean memoryProfiling = false; // requires -javaagent:classmexer.jar to run
Expand Down Expand Up @@ -111,10 +109,8 @@ public static void load() {
tipSamplesPerDeme = (int) map.get("tipSamplesPerDeme");
tipSamplingProportional = (boolean) map.get("tipSamplingProportional");
treeProportion = (double) map.get("treeProportion");
diversitySamplingCount = (int) map.get("diversitySamplingCount");
netauSamplingCount = (int) map.get("netauSamplingCount");
diversitySamplingCount = (int) map.get("diversitySamplingCount");
netauWindow = (int) map.get("netauWindow");
serialIntervalSamplingCount = (int) map.get("serialIntervalSamplingCount");
repeatSim = (boolean) map.get("repeatSim");
immunityReconstruction = (boolean) map.get("immunityReconstruction");
memoryProfiling = (boolean) map.get("memoryProfiling");
Expand Down
59 changes: 25 additions & 34 deletions Simulation.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ public class Simulation {
private double diversity;
private double tmrca;
private double netau;
private double serialInterval;
private double antigenicDiversity;

// constructor
Expand Down Expand Up @@ -87,10 +88,14 @@ public double getTmrca() {
return tmrca;
}

public double getSerialInterval() {
return serialInterval;
}

public double getAntigenicDiversity() {
return antigenicDiversity;
}

// proportional to infecteds in each deme
public int getRandomDeme() {
int n = Random.nextInt(0,getN()-1);
Expand Down Expand Up @@ -155,19 +160,7 @@ public double getAverageRisk(Phenotype p) {
return averageRisk;

}

public double getSerialInterval() {
double interval = 0.0;
if (getI()>0) {
for (int i = 0; i < Parameters.serialIntervalSamplingCount; i++) {
Virus v = getRandomInfection();
interval += v.serialInterval();
}
interval /= (double) Parameters.serialIntervalSamplingCount;
}
return interval;
}


public void printImmunity() {

try {
Expand Down Expand Up @@ -249,8 +242,8 @@ public void printState() {
public void printHeader(PrintStream stream) {
stream.print("date\tdiv\ttmrca\tnetau\tserialInterval\tantigenicDiv\ttotalN\ttotalS\ttotalI\ttotalR\ttotalCases");
for (int i = 0; i < Parameters.demeCount; i++) {
String name = Parameters.demeNames[i];
stream.printf("\t%sDiv\t%sTmrca\t%sNetau\t%sN\t%sS\t%sI\t%sR\t%sCases", name, name, name, name, name, name, name, name);
HostPopulation hp = demes.get(i);
hp.printHeader(stream);
}
stream.println();
}
Expand All @@ -267,10 +260,18 @@ public void printState(PrintStream stream) {
}

public void updateDiversity() {

diversity = 0.0;
tmrca = 0.0;
antigenicDiversity = 0.0;
netau = 0.0;
serialInterval = 0.0;

double coalCount = 0.0;
double coalOpp = 0.0;
double coalWindow = Parameters.netauWindow / 365.0;
int sampleCount = Parameters.diversitySamplingCount;

for (int i = 0; i < sampleCount; i++) {
Virus vA = getRandomInfection();
Virus vB = getRandomInfection();
Expand All @@ -281,29 +282,20 @@ public void updateDiversity() {
tmrca = dist;
}
antigenicDiversity += vA.antigenicDistance(vB);
coalOpp += coalWindow;
coalCount += vA.coalescence(vB, coalWindow);
serialInterval += vA.serialInterval();
}
}

diversity /= (double) sampleCount;
tmrca /= 2.0;
antigenicDiversity /= (double) sampleCount;
}

public void updateNetau() {
double coalCount = 0.0;
double coalOpp = 0.0;
double coalWindow = Parameters.netauWindow / 365.0;
int sampleCount = Parameters.netauSamplingCount;
for (int i = 0; i < sampleCount; i++) {
Virus vA = getRandomInfection();
Virus vB = getRandomInfection();
if (vA != null && vB != null) {
coalOpp += coalWindow;
coalCount += vA.coalescence(vB, coalWindow);
}
}
netau = coalOpp / coalCount;
}

serialInterval /= (double) sampleCount;

}

public void resetCases() {
for (int i = 0; i < Parameters.demeCount; i++) {
HostPopulation hp = demes.get(i);
Expand Down Expand Up @@ -345,7 +337,6 @@ public void run() {

if (Parameters.day % Parameters.printStep == 0) {
updateDiversity();
updateNetau();
printState();
printState(seriesStream);
resetCases();
Expand Down
8 changes: 5 additions & 3 deletions parameters.yml
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,7 @@
# This contains all the parameters used in the model.
# These values correspond to defaults.
# Leaving an entry out is fine, the value will remain at its default.

# simulation parameters
burnin: 0 # days to wait before logging output
endDay: 5000 # number of days to simulate
Expand All @@ -6,10 +10,8 @@ tipSamplingRate: 0.0002 # store X samples per deme per day
tipSamplesPerDeme: 2000 # cap number of samples per deme
tipSamplingProportional: true # whether to sample proportional to prevalence
treeProportion: 0.1 # proportion of tips to use in tree reconstruction
diversitySamplingCount: 1000 # how many samples to draw to calculate diversity
netauSamplingCount: 1000 # how many samples to draw to calculate Ne*ta
diversitySamplingCount: 1000 # how many samples to draw to calculate diversity, netau and serial interval
netauWindow: 100 # window in days to calculate Ne*tau
serialIntervalSamplingCount: 1000 # how many samples to draw to calculate serial interval
repeatSim: true # repeat simulation until endDay is reached?
immunityReconstruction: false # whether to print immunity reconstruction to out.immunity
memoryProfiling: false # requires -javaagent:classmexer.jar to run
Expand Down

0 comments on commit fd3e07d

Please sign in to comment.