PoonLab / Kaphi

Kernel-embedded ABC-SMC for phylodynamic inference
GNU Affero General Public License v3.0
4 stars 2 forks source link

ws$sim.trees contains integers, not tree objects after run.smc is finished #76

Closed MathiasRenaud closed 7 years ago

MathiasRenaud commented 7 years ago

run.smc with 10 particles results in:

> ws$sim.trees
[[1]]
[[1]][[1]]
[1] 1

[[2]]
[[2]][[1]]
[1] 2

[[3]]
[[3]][[1]]
[1] 3

[[4]]
[[4]][[1]]
[1] 4

[[5]]
[[5]][[1]]
[1] 5

[[6]]
[[6]][[1]]
[1] 6

[[7]]
[[7]][[1]]
[1] 7

[[8]]
[[8]][[1]]
[1] 8

[[9]]
[[9]][[1]]
[1] 9

[[10]]
[[10]][[1]]
[1] 10

ws$sim.trees should contain the simulated trees instead of integers.

MathiasRenaud commented 7 years ago

Currently tracing run.smc to solve issue

MathiasRenaud commented 7 years ago

Running initialize.smc with the same config as above gives 50 trees in ws$sim.trees (5 samples * 10 particles).

ArtPoon commented 7 years ago

@MathiasRenaud noticed that in smcABC.R we are needlessly replicating samples:

        for (j in 1:config$nsample) {
            # retain sim.trees in case we revert to previous particle
            new.trees <- simulate.trees(ws, new.particle, model=model)
            new.dists[,i] <- sapply(new.trees, function(sim.tree) {
                distance(ws$obs.tree, sim.tree, config)
            })
        }

simulate.trees should already be simulating config$nsample trees per particle, so there should be no need to loop over nsample. Fix and see if this explains this issue.

MathiasRenaud commented 7 years ago

I have figured out that manually running each line of run.smc works and updates the workspace to store the simulated trees, particles, weights, new.weights, dists. However, calling the function directly gives the same result but does not update the items in workspace.

.epsilon.obj.func is called by .next.epsilon to solve for the new epsilon. The function takes the arguments ws & epsilon and returns the value of (.ess(ws$new.weights) - config$quality * .ess(ws$weights)) but not ws.

Currently working on a way to pass both the value of the new epsilon and the updated workspace back to .next.epsilon. The tricky part is that .epsilon.obj.func is called within uniroot, so the return value must numerical to be usable by uniroot.

MathiasRenaud commented 7 years ago

Here is the use in next.epsilon:

res <- uniroot(function(x) .epsilon.obj.func(ws, x), lower=0,
        upper=ws$epsilon, tol=config$step.tolerance, maxiter=1E6)
    root <- res$root
ArtPoon commented 7 years ago

run.smc is not returning workspace object ws, instead it returns list object results. Try returning both as a list containing the two objects and see if ws contains expected data, i.e., newick tree strings.

MathiasRenaud commented 7 years ago

run.smc now returns ret, where ret is a named list of workspace and result. Workspace now stores all expected data from run.smc.

    # pack ws and result into one list to be returned
    ret <- list(workspace=ws, result=result)

    return (ret)