poonlab/kaphi

Kernel-embedded ABC-SMC for phylodynamic inference

kaphi's Issues

Failing to set C igraph attributes

When calling EAN(graph, attribute, edge id), R crashes. This seems to be related to the general lack of support for attributes in C igraph:

Error at cattributes.c:2385 :Unknown attribute, Invalid value
Abort trap: 6

Pass an igraph object from R to C

Writing a simple bit of C code test.c and an R script test.R that will create an igraph object from a phylogeny and then pass it to C. Working off the rinterface.c code from rigraph as a template. There seem to be some versioning issues with Rosemary's fork of igraph and function calls in rinterface.c that I've temporarily resolved by stripping out the offending functions in the latter.

In coalescent example, particles are not being resampled

When I ran the script in examples/example-coalescent.R, the sampler ran for over 50 iterations but did not resample particles when the effective sample size had gone below the set tolerance (50). Need to dig into why.

Problems running Kaphi tests

After recursively cloning the repo and installing the required R packages, I ran runTestSuite.sh and got the following:

fabusard@Nagisa:~/git/Kaphi$ sh runTestSuite.sh

R version 3.2.3 (2015-12-10) -- "Wooden Christmas-Tree"
Copyright (C) 2015 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> require(RUnit, quietly=TRUE)
> 
> testsuite.Kaphi <- defineTestSuite('Kaphi',
+     dirs=file.path(getwd(), 'tests'),
+     testFileRegexp='^test_.+\\.R',
+     testFuncRegexp='^test.+'
+ )
> testResult <- runTestSuite(testsuite.Kaphi)


Executing test function test.issue.22  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "parse.input.tree"
In addition: Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE,  :
  there is no package called ‘Kaphi’
 done successfully.



Executing test function test.kernel.label  ... essentially unlabelled
Timing stopped at: 0 0 0 
Error in func() : could not find function "tree.kernel"
 done successfully.



Executing test function test.kernel.normalized  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "tree.kernel"
 done successfully.



Executing test function test.kernel.trivial  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "tree.kernel"
 done successfully.



Executing test function test.kernel.unnormalized  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "tree.kernel"
 done successfully.



Executing test function test.rescale.tree  ... Timing stopped at: 0 0 0.001 
Error in func() : could not find function "rescale.tree"
 done successfully.

Error in eval(expr, envir, enclos) : could not find function "parse.ode"
In addition: Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE,  :
  there is no package called ‘Kaphi’


Executing test function test.distance  ... Error in eval(expr, envir = parent.frame()) : 
  could not find function "distance"
In addition: Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE,  :
  there is no package called ‘Kaphi’
Timing stopped at: 0 0 0 
Error in func() : could not find function "tree.kernel"
 done successfully.



Executing test function test.epsilon.obj.func  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "epsilon.obj.func"
 done successfully.



Executing test function test.ess  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "ess"
 done successfully.



Executing test function test.initialize.smc  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "load.config"
 done successfully.



Executing test function test.next.epsilon  ... Timing stopped at: 0 0 0.001 
Error in func() : could not find function "ess"
 done successfully.



Executing test function test.perturb.particles  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "load.config"
 done successfully.



Executing test function test.resample.particles  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "resample.particles"
 done successfully.



Executing test function test.run.smc  ... Timing stopped at: 0 0 0.001 
Error in DEACTIVATED("Used for debugging; produces a lot of output!") : 
  Used for debugging; produces a lot of output!
 done successfully.



Executing test function test.simulate.trees  ... Timing stopped at: 0 0 0 
Error in func() : object 'const.coalescent' not found
 done successfully.



Executing test function test.load.config  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "load.config"
In addition: Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE,  :
  there is no package called ‘Kaphi’
 done successfully.



Executing test function test.prior.density  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "load.config"
 done successfully.



Executing test function test.proposal.density  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "load.config"
 done successfully.



Executing test function test.propose  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "load.config"
 done successfully.



Executing test function test.sample.priors  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "load.config"
 done successfully.



Executing test function test.set.model  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "load.config"
 done successfully.



Executing test function test.get.tip.heights  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "get.tip.heights"
In addition: Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE,  :
  there is no package called ‘Kaphi’
 done successfully.



Executing test function test.init.workspace  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "init.workspace"
 done successfully.



Executing test function test.parse.input.tree  ... Timing stopped at: 0 0 0.001 
Error in func() : could not find function "parse.input.tree"
 done successfully.



Executing test function test.parse.labels  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "parse.labels"
 done successfully.



Executing test function test.rescale.tree  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "parse.input.tree"
 done successfully.

Error in eval(expr, envir, enclos) : could not find function "parse.ode"
In addition: Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE,  :
  there is no package called ‘Kaphi’


Executing test function test.colless  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "colless"
In addition: Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE,  :
  there is no package called ‘Kaphi’
 done successfully.



Executing test function test.nLTT  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "nLTT"
 done successfully.



Executing test function test.sackin  ... Timing stopped at: 0 0 0 
Error in func() : could not find function "sackin"
 done successfully.

> printTextProtocol(testResult)
RUNIT TEST PROTOCOL -- Tue May  2 10:09:00 2017 
*********************************************** 
Number of test functions: 30 
Number of deactivated test functions: 1 
Number of errors: 30 
Number of failures: 0 

 
1 Test Suite : 
Kaphi - 30 test functions, 30 errors, 0 failures
ERROR in test.issue.22: Error in func() : could not find function "parse.input.tree"
ERROR in test.kernel.label: Error in func() : could not find function "tree.kernel"
ERROR in test.kernel.normalized: Error in func() : could not find function "tree.kernel"
ERROR in test.kernel.trivial: Error in func() : could not find function "tree.kernel"
ERROR in test.kernel.unnormalized: Error in func() : could not find function "tree.kernel"
ERROR in test.rescale.tree: Error in func() : could not find function "rescale.tree"
ERROR in /home/fabusard/git/Kaphi/tests/test_simcoal.R: Error while sourcing  /home/fabusard/git/Kaphi/tests/test_simcoal.R : Error in eval(expr, envir, enclos) : could not find function "parse.ode"
ERROR in test.distance: Error in func() : could not find function "tree.kernel"
ERROR in test.epsilon.obj.func: Error in func() : could not find function "epsilon.obj.func"
ERROR in test.ess: Error in func() : could not find function "ess"
ERROR in test.initialize.smc: Error in func() : could not find function "load.config"
ERROR in test.next.epsilon: Error in func() : could not find function "ess"
ERROR in test.perturb.particles: Error in func() : could not find function "load.config"
ERROR in test.resample.particles: Error in func() : could not find function "resample.particles"
DEACTIVATED test.run.smc: Used for debugging; produces a lot of output!
ERROR in test.simulate.trees: Error in func() : object 'const.coalescent' not found
ERROR in test.load.config: Error in func() : could not find function "load.config"
ERROR in test.prior.density: Error in func() : could not find function "load.config"
ERROR in test.proposal.density: Error in func() : could not find function "load.config"
ERROR in test.propose: Error in func() : could not find function "load.config"
ERROR in test.sample.priors: Error in func() : could not find function "load.config"
ERROR in test.set.model: Error in func() : could not find function "load.config"
ERROR in test.get.tip.heights: Error in func() : could not find function "get.tip.heights"
ERROR in test.init.workspace: Error in func() : could not find function "init.workspace"
ERROR in test.parse.input.tree: Error in func() : could not find function "parse.input.tree"
ERROR in test.parse.labels: Error in func() : could not find function "parse.labels"
ERROR in test.rescale.tree: Error in func() : could not find function "parse.input.tree"
ERROR in /home/fabusard/git/Kaphi/tests/test_solve_ODE.R: Error while sourcing  /home/fabusard/git/Kaphi/tests/test_solve_ODE.R : Error in eval(expr, envir, enclos) : could not find function "parse.ode"
ERROR in test.colless: Error in func() : could not find function "colless"
ERROR in test.nLTT: Error in func() : could not find function "nLTT"
ERROR in test.sackin: Error in func() : could not find function "sackin"



Details 
*************************** 
Test Suite: Kaphi 
Test function regexp: ^test.+ 
Test file regexp: ^test_.+\.R 
Involved directory: 
/home/fabusard/git/Kaphi/tests 
--------------------------- 
Test file: /home/fabusard/git/Kaphi/tests/test_kernel.R 
test.issue.22: ERROR !! 
Error in func() : could not find function "parse.input.tree"
test.kernel.label: ERROR !! 
Error in func() : could not find function "tree.kernel"
test.kernel.normalized: ERROR !! 
Error in func() : could not find function "tree.kernel"
test.kernel.trivial: ERROR !! 
Error in func() : could not find function "tree.kernel"
test.kernel.unnormalized: ERROR !! 
Error in func() : could not find function "tree.kernel"
test.rescale.tree: ERROR !! 
Error in func() : could not find function "rescale.tree"
--------------------------- 
Test file: /home/fabusard/git/Kaphi/tests/test_simcoal.R 
/home/fabusard/git/Kaphi/tests/test_simcoal.R: ERROR !! 
Error while sourcing  /home/fabusard/git/Kaphi/tests/test_simcoal.R : Error in eval(expr, envir, enclos) : could not find function "parse.ode"
--------------------------- 
Test file: /home/fabusard/git/Kaphi/tests/test_smc_abc.R 
test.distance: ERROR !! 
Error in func() : could not find function "tree.kernel"
test.epsilon.obj.func: ERROR !! 
Error in func() : could not find function "epsilon.obj.func"
test.ess: ERROR !! 
Error in func() : could not find function "ess"
test.initialize.smc: ERROR !! 
Error in func() : could not find function "load.config"
test.next.epsilon: ERROR !! 
Error in func() : could not find function "ess"
test.perturb.particles: ERROR !! 
Error in func() : could not find function "load.config"
test.resample.particles: ERROR !! 
Error in func() : could not find function "resample.particles"
test.run.smc : DEACTIVATED, Used for debugging; produces a lot of output!
test.simulate.trees: ERROR !! 
Error in func() : object 'const.coalescent' not found
--------------------------- 
Test file: /home/fabusard/git/Kaphi/tests/test_smc_config.R 
test.load.config: ERROR !! 
Error in func() : could not find function "load.config"
test.prior.density: ERROR !! 
Error in func() : could not find function "load.config"
test.proposal.density: ERROR !! 
Error in func() : could not find function "load.config"
test.propose: ERROR !! 
Error in func() : could not find function "load.config"
test.sample.priors: ERROR !! 
Error in func() : could not find function "load.config"
test.set.model: ERROR !! 
Error in func() : could not find function "load.config"
--------------------------- 
Test file: /home/fabusard/git/Kaphi/tests/test_smc_workspace.R 
test.get.tip.heights: ERROR !! 
Error in func() : could not find function "get.tip.heights"
test.init.workspace: ERROR !! 
Error in func() : could not find function "init.workspace"
test.parse.input.tree: ERROR !! 
Error in func() : could not find function "parse.input.tree"
test.parse.labels: ERROR !! 
Error in func() : could not find function "parse.labels"
test.rescale.tree: ERROR !! 
Error in func() : could not find function "parse.input.tree"
--------------------------- 
Test file: /home/fabusard/git/Kaphi/tests/test_solve_ODE.R 
/home/fabusard/git/Kaphi/tests/test_solve_ODE.R: ERROR !! 
Error while sourcing  /home/fabusard/git/Kaphi/tests/test_solve_ODE.R : Error in eval(expr, envir, enclos) : could not find function "parse.ode"
--------------------------- 
Test file: /home/fabusard/git/Kaphi/tests/test_tree_stats.R 
test.colless: ERROR !! 
Error in func() : could not find function "colless"
test.nLTT: ERROR !! 
Error in func() : could not find function "nLTT"
test.sackin: ERROR !! 
Error in func() : could not find function "sackin"
> 

Refactor deSolve:ode(method="rk4")

The corresponding function is call_rk4. When we do this, call the function in both R and from C and check that the results match.

Merge labeled kernel

I accidentally pushed a merge of the labeled kernel.

You should delete my last push and merge my labeled kernel yourself to have a cleaner repo.

Write unit tests for main SMC loop

Finished porting code from netabc to Kaphi.

compiling igraph from source

No instructions provided on developer website as far as I can tell. We need to install the following dependencies before compiling from source:

• bison
• automake -- to modify Makefile.am (has to be version 1.14?!?)
• autoconf

.simulate.tree loop terminates early

Some example output:

190 221 3267.577 3294.234 31 coalesce
191 221 3294.234 3304.434 30 coalesce
192 221 3304.434 3335.753 29 coalesce
Error in sample.path(a.state, b.state, 0, z$edge.length[u], qm) : 
  The rate matrix is not valid. The rates must sum to 0 zero within each row.-5.91058801335229000.00183720687095107-0.0027569008437571405.908750806481330.002756900843757140

The first three lines are produced by some debug statements that I inserted into the routine that loops over events (sampling tips or coalescent events). There should be 220 events, but for some reason the program is exiting after the 192nd event. As a result, the edge matrix of the new tree is not completely populated and is most likely the reason for downstream bugs, e.g., "Error in sample.path()"

tree simulation code is locking up

Need to do some profiling of R code to determine what the issue is.
Encountered while running batch processing loop in rong-perelson.R script.

Port bisection solver for updating epsilon from C to R

Make ECctmc package optional

Currently, installation of Kaphi is dependent on the R package ECctmc. This package should be made optional as it is only necessary for certain projects.

Implement standard speciation models

This should include the following models:

• Yule model
• birth-death models
• binary state speciation extinction, BiSSE
• MuSSE
• QuaSSE
• GeoSSE
• BiSS-ness
• ClaSSE

Map edges to ape:phylo edges

We are mapping state transitions to branches of the simulated tree, and returning these transitions as a list of matrices keyed by node numbers. However, when we convert the tree into an ape:phylo object then the node numbering is scrambled.

The phylo tip labels retain the original node numbering for tips. For example, using the attached R object:

> as.integer(phylo$tip.label)
  [1]  79  92  55  13  46  23   8  93  38  78   4  85  31  44  56  88   1  45  64  81  10  94  70
 [24]  83  29  77  35  76  57  28  50  18  59  20  36  86  91  17  26 100  47  71  22  58   2  87
 [47]  43  98  80   5  68  75  27  32  84   7  25  48  30  51  53  54  99  63  16  67  82  96  41
 [70]  15  42  66  69  49  89  12  19  21   9  39  52  97  14  65  95  37  61  24  40  60  74   6
 [93]  72  33  73  34  90  11   3  62

phylo node 1 corresponds to z node 79. This means that z$edge.length[79] should map to phylo$edge.length[x] as follows:

> phylo$edge.length[which(phylo$edge[,2]==1)]
[1] 190.3292
> z$edge.length[79]
[1] 190.3292

So this should be straight-forward for terminal nodes, but we need to map the internal nodes as well.

zphylo.RData.zip

Validate ABC-SMC estimation for compartmental (SIR) models

Just like kernel-ABC validation

Problem simulating tree from RP model

Error in sample.path(a.state, b.state, 0, z$edge.length[u], qm) : 
  The process cannot start in an absorbing state if the endpoints are different.

The corresponding line in sample.path() is:

 if (all(Q[a, ] == 0) & a != b) {
        stop("The process cannot start in an absorbing state if the endpoints are different.", Q)
    }

Refactor R code so that non-user functions have dot prefix

For example, make.fgy() should be called .make.fgy().
Also remove these files from the package metadata so they don't get exposed by default.

Kernel normalization problem

This code produces a ws$dists matrix with negative values:

require(Kaphi)
require(yaml)
    config <- load.config('tests/fixtures/coalescent.yaml')
    config <- set.model(config, const.coalescent)
    nparticle <- config$nparticle
    theta <- c(Ne.tau=100)
    set.seed(100)
    obs.tree <- const.coalescent(theta, nsim=1, n.tips=20)[[1]]
    ws <- init.workspace(obs.tree, config)
ws
config <- ws$config
ws <- initialize.smc(ws)
ws$epsilon
ws$accept
ws$alive
next.epsilon(ws)
epsilon.obj.func(ws, 1)
epsilon.obj.func(ws, 0.1)
ws$dists

Output:

            [,1]         [,2]        [,3]        [,4]        [,5]         [,6]
[1,]  0.09855444 -0.003547483  0.01952083  0.26282476 -0.08976043 -0.358024602
[2,]  0.18498320 -0.019951910 -0.05274277 -0.20163216 -0.03755821 -0.313424967
[3,] -0.08907720 -0.003324091  0.10716642  0.03362616 -0.07638166  0.019548475
[4,] -0.01809160 -0.111153807 -0.12428295 -0.03709038  0.36552203 -0.009851596
[5,]  0.09247920 -0.074805804 -0.04914617 -0.39641339 -0.03768756 -0.066284821
[6,]  0.02441562  0.020498423  0.29638671 -0.04224123 -0.33942571  0.041807836
            [,7]        [,8]        [,9]       [,10]
[1,] -0.11590842 -0.08407187 -0.05446778  0.03628705
[2,] -0.01963804 -0.03468434 -1.32059338  0.11950563
[3,] -0.03013048 -0.26036020 -0.02060306 -0.06903098
[4,] -0.02251815  0.01507980 -0.25298418 -0.35906854
[5,] -0.46878691 -0.00202518 -0.07799568 -0.02806073
[6,] -0.72675567 -0.06024982 -0.05166268 -1.10345042

Something is very wrong here!

Judy problem

Error in dyn.load(file, DLLpath = DLLpath, ...) : 
  unable to load shared object '/Library/Frameworks/R.framework/Versions/3.2/Resources/library/Kaphi/libs/Kaphi.so':
  dlopen(/Library/Frameworks/R.framework/Versions/3.2/Resources/library/Kaphi/libs/Kaphi.so, 6): Symbol not found: _JudyLFreeArray
  Referenced from: /Library/Frameworks/R.framework/Versions/3.2/Resources/library/Kaphi/libs/Kaphi.so
  Expected in: flat namespace
 in /Library/Frameworks/R.framework/Versions/3.2/Resources/library/Kaphi/libs/Kaphi.so
Error: loading failed

I just fixed the same type of problem where R couldn't find the symbol for _gsl_rng_default, which was resolved by adding -L/usr/local/lib -lgsl to Makevars. I probably just need to add the same thing for the Judy library.

Possible for cherries to have label mismatch

Currently when computing tree kernel at node pair n1 and n2, if left and right children of both nodes are tips then these subtrees rooted at n1 and n2 are a match. However, this situation may arise:

 A  B    B  A
 \  /    \  /
  n1      n2

This should still be counted as a match. Ladderizing the tree does not rotate cherries with respect to tip labels.

Can't compile igraph

Current objective is to implement and call the kernel function from R:

• Pass tree strings as SEXP objects from R
• parse tree strings using bison and flex functionality (see issue #9)
• calculate kernel on resulting igraph objects

Currently unable to compile code due to this error:

is2882:Kaphi art$ R CMD INSTALL pkg
* installing to library ‘/Library/Frameworks/R.framework/Versions/3.3/Resources/library’
* installing *source* package ‘Kaphi’ ...
** libs
clang -dynamiclib -Wl,-headerpad_max_install_names -undefined dynamic_lookup -single_module -multiply_defined suppress -L/Library/Frameworks/R.framework/Resources/lib -L/usr/local/lib -o Kaphi.so kernel.o newick_lexer.o newick_parser.o newick_parser.tab.o rinterface.o rinterface_extra.o tree.o treekernel.o treestats.o util.o -Ligraph/src -ligraph -F/Library/Frameworks/R.framework/.. -framework R -Wl,-framework -Wl,CoreFoundation
ld: library not found for -ligraph
clang: error: linker command failed with exit code 1 (use -v to see invocation)
make: *** [Kaphi.so] Error 1
ERROR: compilation failed for package ‘Kaphi’
* removing ‘/Library/Frameworks/R.framework/Versions/3.3/Resources/library/Kaphi’
* restoring previous ‘/Library/Frameworks/R.framework/Versions/3.3/Resources/library/Kaphi’

There is no libigraph available in igraph subdirectory because I haven't been able to get it to compile.

allow Kaphi to fit model to more than one tree

For example, if the user has a sample of trees from a posterior distribution, rather than an ML point estimate

Linear interpolation for not-yet-sampled lineages

The functionsimcoal.R:solve.A.mx() attempts a numerical solution of the ODE defined for the time differential of the A matrix, which represents the number of lineages per deme (columns) over coalescent (reverse) time (rows).

There seems to be a problem when one or more tips share the same sample height. In the following example, there are 10 tips (sampled lineages), of which 5 were sampled at height 0 and 5 were sampled at height 5.

The numerical solution for A seems to make sense:

> odA[,2]
  [1] 10.0000000  9.9990524  9.9980434  9.9969767  9.9958335  9.9946077  9.9912967
  [8]  9.9857463  9.9798173  9.9734765  9.9666921  9.9594249  9.9515829  9.9431616
 [15]  9.9341489  9.9244512  9.9140338  9.9028100  9.8907590  9.8778023  9.8638435
 [22]  9.8488504  9.8327174  9.8152404  9.7964088  9.7760982  9.7541963  9.7306515
 [29]  9.7053315  9.6780195  9.6486798  9.6171166  9.5831978  9.5467410  9.5072452
 [36]  9.4647858  9.4191915  9.3703111  9.3178661  9.2616768  9.2015377  9.1372304
 [43]  9.0685395  8.9952256  8.9170925  8.8331706  8.7438679  8.6490141  8.5484277
 [50]  8.4419071  8.3293246  8.2105541  8.0855163  7.9541241  7.8164045  7.6723116
 [57]  7.5205391  7.3624673  7.1982683  7.0281699  6.8524534  6.6714669  6.4855980
 [64]  6.2953147  6.1011332  5.9035997  5.7033039  5.4990122  5.2932205  5.0865718
 [71]  4.8797999  4.6736230  4.4687199  4.2657660  4.0654267  3.8683308  3.6750242
 [78]  3.4860542  3.3002343  3.1196935  2.9448147  2.7759163  2.6132292  2.4569670
 [85]  2.3072419  2.1641459  2.0277004  1.8978861  1.7746300  1.6568001  1.5454055
 [92]  1.4402775  1.3412391  1.2480833  1.1606064  1.0785600  1.0017248  0.9298505
 [99]  0.8626970  0.8000171

However, the not.sampled.yet function returns the following values on the same timeline:

> t(sapply(haxis, function(h) c(h, not.sampled.yet(h))))
             [,1] [,2]
  [1,]  0.0000000    8
  [2,]  0.9090909    5
  [3,]  1.8181818    5
  [4,]  2.7272727    5
  [5,]  3.6363636    5
  [6,]  4.5454545    5
  [7,]  5.4545455    0
  [8,]  6.3636364    0
  [9,]  7.2727273    0
 [10,]  8.1818182    0
 [11,]  9.0909091    0
 [12,] 10.0000000    0

Something silly is happening at h=0. The not.sampled.yet function results from the following linear interpolation:

    nsy.index <- approxfun(
        x=sort(jitter(sorted.sample.heights,
            factor=max(1e-6, sorted.sample.heights[length(sorted.sample.heights)]/1e6))),
        y=1:n, method='constant', rule=2
    )
    not.sampled.yet <- function(h) {
        cumul.sns[nsy.index(h), ]
    }

Need R header files (R.h) to compile C extensions for R

Unfortunately MacPorts doesn't have R-devel. I'm going to have to compile R from source 😱

const.coalescent is producing trees with nonsense tip labels

Encountered this while debugging issue #22
For example:

> write.tree(sim.tree)
[1] "((((((1.11022302462516e-16:0.02524922802,1.11022302462516e-16:0.02524922802):0.05033006789,1.11022302462516e-16:0.0755792959):0.05444047653,(1.11022302462516e-16:0.02979646736,1.11022302462516e-16:0.02979646736):0.1002233051):0.08255408061,((1.11022302462516e-16:0.04722236556,1.11022302462516e-16:0.04722236556):0.01307631053,1.11022302462516e-16:0.06029867609):0.152275177):0.2675979671,((((((1.11022302462516e-16:0.02267684476,1.11022302462516e-16:0.02267684476):0.03520747157,1.11022302462516e-16:0.05788431632):0.01945058017,1.11022302462516e-16:0.07733489649):0.0300051726,1.11022302462516e-16:0.1073400691):0.07328217883,1.11022302462516e-16:0.1806222479):0.1200019126,1.11022302462516e-16:0.3006241606):0.1795476596):2.64717087,(((1.11022302462516e-16:0.057328283,1.11022302462516e-16:0.057328283):0.02427477499,0:0.08160305799):0.1849555688,(1.11022302462516e-16:0.06503798585,0:0.06503798585):0.201520641):2.860784063);"

Minimize rinterface code but keep acknowledgement and license

rinterface.c is a huge file. Cut out functions that we aren't going to use and check that the package still compiles.

Store deme migration events in coalescent tree as internal nodes

Currently, migrations between demes (state transitions) occur as events along branches of a coalescent tree, but there is no explicit representation of this in the output. We can modify the simulation code to return an igraph object instead of a phylo object so that the state transitions are represented by internal nodes of degree 2. It should then be straight-forward to convert the igraph object into a phylo.

Write unit tests for R_Kaphi_parse_newick

This should probably be run 🕶️😏 at C level 😎 .

confirm that branch lengths and topology in igraph object is correctly passed from R to C

implement this in a unit test suite

Parallelize SMC routines

Using function mclapply from R package parallel. Good place to start may be simulation of trees.

Tree kernel requires exact label matches

In src/treestats.c, node labels (strings) are converted into integers using a deterministic heap sort and search. Thereafter, labels are compared as integer values within the kernel code. However, we may want match labels as patterns. For example, if a tip has multiple states A and B, it should match another tip in state A just as it matches tips in state B.

problem cloning igraph with fresh clone of repo

When cloning the Kaphi repo to the lab iMac, we ran into this error message:

is2914:git art$ git clone --recursive https://github.com/PoonLab/Kaphi.git
Cloning into 'Kaphi'...
remote: Counting objects: 989, done.
remote: Total 989 (delta 0), reused 0 (delta 0), pack-reused 989
Receiving objects: 100% (989/989), 464.65 KiB | 0 bytes/s, done.
Resolving deltas: 100% (653/653), done.
Submodule 'igraph' (git://github.com/rmcclosk/igraph.git) registered for path 'pkg/src/igraph'
Cloning into '/Users/art/git/Kaphi/pkg/src/igraph'...
Fetched in submodule path 'pkg/src/igraph', but it did not contain 6055bff266ae84fcf4e48fb7a505fd97dec3911f. Direct fetching of that commit failed.

Unit testing coalescent simulation

Using undampened SI model to simulate trees as exponential coalescent. The resulting trees should have coalescent intervals that concord with model expectations.

• Use regression to assess model fit.
• Compare rate parameter estimate to tree simulation setting
• Look at tree imbalance statistics?

Mapping `phylo` to `z` is failing

Error in .simulate.ode.tree(sample.times, sample.states, fgy, solve.QAL,  : 
  (list) object cannot be coerced to type 'integer'

The code is failing at this step:

index <- as.integer(izx[names(ipx)])

where ipx is an inverted associative list with concatenated tip labels as keys and phylo node indices as values, and izx is also an inverted associative list with what should be the same keys and z node indices as values.

Incorporate tree shape statistics into SMC-ABC

We need to evaluate the kernel method against tree shape statistics and to provide the option to use any combination of such statistics as a similarity measure. It would be nice to use R's model specification format as in glm for this, i.e.: ~ colless + cherries + kernel

Problem with coalescent tree simulation for multiple demes

Error in t(Q) %*% mstates[is.extant, ] : non-conformable arguments

plot.smc.config should show multiple prior distributions

Prompt user to hit to view next plot to cycle through each prior distribution

Write Makefile to build newick_parser.h from .l and .y sources

The current file was built within the netabc repo that includes autogen.sh and Makefile.ac scripts. In order to make Kaphi portable across platforms, we need to be able to recompile these sources. In addition, manually compiling newick_parser.y generates newick_parser.tab.h and newick_parser.tab.c files (see #9) that conflict with newick_parser.h because of duplication of symbols when trying to build the Kaphi package. We need a Makefile that will build these sources and then clean up files that are no longer needed.

Pass tree as igraph object from R to tree.c:height()

Write manual .Rd pages

Here is a current list of R functions - those with a star next to them are meant to be internal functions and should probably be prefixed with a . (see #48)

• treekernel.R
  • rescale.tree*
  • preprocess.tree*
  • to.newick*
  • utk (stands for unlabeled tree kernel)
  • tree.kernel
  • nLTT (normalized lineages through time metric)
  • sackin (the following functions are tree shape statistics)
  • colless
  • cophenetic
  • ladder.length
  • IL.nodes
  • tree.width
  • max.delta.width
  • n.cherries
  • prop.unbalanced
  • avg.unbalance
  • pybus.gamma
  • internal.terminal.ratio
• solveODE.R
  • parse.ode
  • get.times*
  • solve.ode
• smcWorkspace.R
  • parse.input.tree
  • get.tip.heights*
  • get.node.heights*
  • parse.labels
  • init.workspace
  • print.smc.workspace
• smcConfig.R
  • load.config
  • set.model
  • sample.priors
  • prior.density
  • propose
  • proposal.density
  • print.smc.config
  • plot.smc.config
• smcABC.R (Mathias)
  • simulate.trees
  • distance
  • initialize.smc
  • ess*
  • epsilon.obj.func*
  • next.epsilon*
  • resample.particles*
  • perturb.particles*
  • run.smc
• simcoal.R (Tammy)
  • init.fgy*
  • init.QAL.solver*
  • solve.A.mx*
  • get.event.times*
  • update.mstates*
  • sample.path*
  • coalesce.lineages*
  • get.terminals*
  • invert.list*
  • simulate.ode.tree
• ape-wrappers.R (in progress)
  • [ ] const.coalescent
• rcolgem-wrappers.R (in progress)
  • [ ] rcolgem.SI

How are we going to implement prior distributions?

We need a nice way for the user to be able to specify prior distributions on the model parameters. I like the concept of reading in a YAML file that will contain all SMC control parameters, kernel tuning parameters, but do we want the priors to be specified in the same file?

Implement skyline coalescent simulator

This would make a nice basis for #42

Check if nLTT statistic actually works

I was in the middle of writing unit tests but I don't think it's complete yet.

R crashes when I run tests/test_kernel.R

I tried to build this package on my mac. I had to edit a few of the files since i installed igraph to /usr/local on my machine. I also had to remove the last parameter from igraph_neighborhood and igraph_neighborhood_size.

After I installed the Kaphi R pkg and ran the test.kernel.trivial() function in tests/test_kernel.R, R crashes with:
*** caught segfault ***
address 0x7fc7e5802c50, cause 'memory not mapped'

Traceback:
1: .Call("R_Kaphi_kernel", nwk1, nwk2, lambda, sigma, rho, normalize, PACKAGE = "Kaphi")
2: tree.kernel(t1, t1, lambda = 0.1, sigma = 2, rho = 1, normalize = 1, rescale.mode = "NONE")
3: test.kernel.trivial()

I wonder if it's an igraph version mismatch; I have version 0.7.1 of the igraph c library.

uniroot fails in run.smc

Reproducible with unit test test.run.smc.

f() values at end points not of opposite sign

This happens because after one iteration of the main loop in run.smc, ws$epsilon is set to 0.4227564and when we callnext.epsilon(ws)withlower=0andupper=ws$epsilon, we get the following values for f.lowerandf.upper`:

> epsilon.obj.func(ws, 0)
[1] -1
> epsilon.obj.func(ws, ws$epsilon)
[1] -10.1716

Port/refactor wrapper function for deSolve

deSolve is the R package for interfacing with the FORTRAN library ODEPACK.
rcolgem/phydynR uses this library for the numerical solution of ODE systems.

Generate newick_parser.h from .y and .l files

Okay, if we're going to go the route of passing Newick tree strings from R into C, then we're going to have to get bison and flex working. Specifically, we need to be able to compile newick_parser.h from newick_parser.y and newick_lexer.l. So far I haven't been able to get this to work on a Mac; see rmcclosk/thesis#37.

Unit test failure

ERROR in test.run.smc: Error in if (is.na(f.upper)) stop("f.upper = f(upper) is NA") : 
  argument is of length zero

Coalescent events are lower in tree than sample times

I'm currently dealing with a tree where sample times are:

> sort(table(sample.times))
sample.times
1075 2292 2626 1927  223  771 1409 1714 2079 2444 2991   10  466 
   8    8    9   10   11   20   20   20   20   20   20   21   22 

The problem is that the numerical solution of the A matrix returned by solve.A.mx has a very steep decline in the Ts deme (actively infected cells):

> A.mx$mat[1:20,]
             V         L           Ts
 [1,]  88.0000   9.00000 112.00000000
 [2,] 179.6151  16.86399   9.14353328
 [3,] 179.3197  24.32650   0.05054836
 [4,] 178.9634  31.43167   0.00000000
 [5,] 178.6122  38.19613   0.00000000
 [6,] 178.2655  44.63573   0.00000000

As expected, the first row corresponds to the observed sample states:

> apply(sample.states, 2, sum)
  V   L  Ts 
 88   9 112

R yaml.load_file fails to parse 1e-5 as float in scientific notation

Instead this gets parsed as a string and throws an exception down the line. We may be able to get around this by implementing a custom handler in smcConfig.R:load.config()