Merge branch 'master' of https://github.com/BEAST2-Dev/bdsky

examples/HCV_oup_bdsky.R

@@ -0,0 +1,41 @@
+# This script plots OU-BDSKY plot
+
+#  origin_post is a posterior vector of origins
+#  r0 is a data table of posterior samples of r0 vectors, one row per sample
+#  time_grid is a vector of times to evaluate the skyline at
+bdsky_post <- function(origin_post, r0, time_grid) {
+
+	r0_time_gridded <- list()
+
+	n <- ncol(r0)
+
+	for (s in 1:length(origin_post)) {
+
+		origin <- origin_post[s]
+		r0_vec <- r0[s,]
+
+		ind <- pmax(1,n - floor(time_grid / origin * n))
+		r0_time_gridded[[s]] <- r0_vec[ind]			
+	}
+
+	return (r0_time_gridded)
+}
+
+setwd("~/Git/bdsky/examples/")
+
+lf <- read.table("HCV_oup_40_1447852063188.log", sep="\t", header=T)
+
+origin_post <- lf$orig_root
+
+r0_subset <- lf[grepl("R0", names(lf))]
+
+time_grid <- 1:400
+
+bdskypost <- bdsky_post(origin_post, r0_subset, time_grid)
+
+pdf("hcv_oubdsky.pdf")
+plot(time_grid,bdskypost[[950]],type='S', xlab="Time (years before present)", ylab="R0",col=rgb(0,0,1,0.1))
+for (s in 20:200*50) {
+	lines(time_grid, bdskypost[[s]], type='S',col=rgb(0,0,1,0.1))
+}
+dev.off()

examples/testOUPrior.xml

@@ -0,0 +1,31 @@
+<beast version='2.0'
+       namespace='beast.core:beast.evolution.speciation:beast.core.util:beast.evolution.nuc:beast.evolution.operators:beast.evolution.sitemodel:beast.math.distributions:beast.evolution.branchratemodel:beast.evolution.likelihood:beast.core.parameter'>
+
+    <!-- Test XML for the Ornstein-Uhlenbeck prior -->
+
+
+    <run spec="MCMC" id="mcmc" chainLength="100000"> <!--autoOptimize="true"-->
+        <state>
+            <input name='stateNode' idref='x'/>
+        </state>
+
+        <distribution spec="CompoundDistribution" id="posterior">
+            <distribution spec="OUPrior">
+                <parameter spec='RealParameter' name='x' id='x' dimension="10" value="1 1 1 1 1 1 1 1 1 1"/>
+                <parameter spec='RealParameter' name='times' id='t' dimension="10" value="0 1 2 3 4 5 6 7 8 9"/>
+                <parameter spec='RealParameter' name='mean' id='mu' dimension="1" value="3"/>
+                <parameter spec='RealParameter' name='sigma' id='sigma' dimension="1" value="1"/>
+                <parameter spec='RealParameter' name='nu' id='nu' dimension="1" value="0.1"/>
+                <distribution spec='Exponential' name="x0Prior" mean="0.1"/>
+            </distribution>
+
+        </distribution>
+
+        <operator id='scaler' spec='ScaleOperator' scaleFactor=".75" weight="1" parameter="@x"/>
+
+        <logger logEvery="10" fileName="ouprior.log">
+            <log idref="x"/>
+        </logger>
+    </run>
+
+</beast>

src/bdsky/BDSSkylineSegment.java

@@ -0,0 +1,35 @@
+package bdsky;
+
+/**
+ * A piecewise constant segment of a skyline.
+ */
+public class BDSSkylineSegment extends SkylineSegment {
+
+    public BDSSkylineSegment(double lambda, double mu, double psi, double r, double t1, double t2) {
+
+        super(t1, t2, new double[]{lambda, mu, psi, r});
+    }
+
+    /**
+     * @return the birth rate per unit time.
+     */
+    public double lambda() { return value[0]; };
+
+    /**
+     * @return the death rate per unit time.
+     */
+    public double mu() { return value[1]; };
+
+    /**
+     * @return the sampling rate per unit time.
+     */
+    public double psi() { return value[2]; };
+
+    /**
+     * @return the removal probability, i.e. the probability that sampling causes recovery/removal/death.
+     */
+    public double r() { return value[3]; };
+
+    // TODO fold rho sampling events into the Skyline
+    // public boolean hasRho();
+}

src/bdsky/MultiSkyline.java

@@ -0,0 +1,135 @@
+package bdsky;
+
+import beast.core.CalculationNode;
+import beast.core.Input;
+
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.Collections;
+import java.util.List;
+
+/**
+ * A multiskyline made up of simple skylines.
+ * This skyline will have a number of segments equal to the union of the number of unique boundaries in the daughter skylines.
+ */
+public class MultiSkyline extends CalculationNode implements Skyline {
+
+    public Input<List<SimpleSkyline>> skylineInput = new Input<>("skyline", "the simple skylines making up this multiple skyline", new ArrayList<>());
+
+    public MultiSkyline(SimpleSkyline... skyline) {
+        try {
+            initByName("skyline", Arrays.asList(skyline));
+        } catch (Exception e) {
+            throw new RuntimeException(e);
+        }
+    }
+
+    @Override
+    public void initAndValidate() throws Exception {}
+
+    @Override
+    public List<SkylineSegment> getSegments() {
+
+        List<SimpleSkyline> skylines = skylineInput.get();
+
+
+        List<Boundary2> boundaries = new ArrayList<>();
+
+        for (int i = 0; i < skylines.size(); i++) {
+            Skyline skyline = skylines.get(i);
+            List<SkylineSegment> segments = skyline.getSegments();
+            for (int j = 0; j < segments.size(); j++) {
+                boundaries.add(new Boundary2(j, segments.get(j).t1, skyline, i));
+            }
+        }
+        Collections.sort(boundaries, (o1, o2) -> Double.compare(o1.time, o2.time));
+
+        System.out.println(boundaries);
+
+        int[] index = new int[skylines.size()];
+
+        List<SkylineSegment> segments = new ArrayList<>();
+
+
+        System.out.println("Boundaries.size = " + boundaries.size());
+
+        int i = 0;
+        double start = boundaries.get(0).time;
+        while (i < boundaries.size()) {
+
+            int j = i + 1;
+
+            double end = Double.POSITIVE_INFINITY;
+            if (j != boundaries.size()) {
+                Boundary2 boundary = boundaries.get(j);
+                end = boundary.time;
+
+                while (j < boundaries.size() && end == start) {
+                    j += 1;
+                    if (j == boundaries.size()) {
+                        end = Double.POSITIVE_INFINITY;
+                    } else {
+                        end = boundaries.get(j).time;
+                    }
+                }
+                System.out.println("next end = " + boundaries.get(j));
+            }
+
+            double[] value = new double[index.length];
+            for (int k = 0; k < index.length; k++) {
+
+                int ind = index[k];
+
+                value[k] = skylines.get(k).getValues().get(ind)[0];
+            }
+
+            SkylineSegment segment = new SkylineSegment(start, end, value);
+            segments.add(segment);
+            System.out.println("Added segment: " + segment);
+
+            if (j != boundaries.size()) {
+                index[boundaries.get(j).skylineIndex] += 1;
+                System.out.println("incremented index for skyline " + boundaries.get(j).skylineIndex);
+            }
+            i = j;
+            start = end;
+        }
+
+        return segments;
+    }
+
+    @Override
+    public int getDimension() {
+
+        int dim = 0;
+        for (SimpleSkyline skyline : skylineInput.get()) {
+
+            dim += skyline.getDimension();
+        }
+        return dim;
+    }
+
+    class Boundary2 {
+        // the index
+        int index;
+
+        // time of the boundary
+        double time;
+
+        // the skyline the boundary is in
+        Skyline skyline;
+
+        int skylineIndex;
+
+        Boundary2(int index, double time, Skyline skyline, int skylineIndex) {
+            this.index = index;
+            this.time = time;
+            this.skyline = skyline;
+            this.skylineIndex = skylineIndex;
+        }
+
+        public String toString() {
+            return "skyline[" + skylineIndex + "].time(" + index + ")=" + time;
+        }
+    }
+}