text:
# Test of PhyloNetworklm on trees

###############################################################################
## Caudata dataset - shared paths matrix
###############################################################################
@testset "phyloNetworklm: Caudata Dataset" begin
## Export "caudata" dataset (from geiger)
phy = readTopology(joinpath(@__DIR__, "..", "examples", "caudata_tree.txt"));

V = sharedPathMatrix(phy);
VR = DataFrame!(CSV.File(joinpath(@__DIR__, "..", "examples", "caudata_shared_paths.txt")))
VR = convert(Matrix, VR);

# Tips
@test V[:Tips] ≈ VR[1:197, 1:197]

# Internal nodes
@test V[:InternalNodes] ≈ VR[-V.internalNodeNumbers .+ 196, -V.internalNodeNumbers .+ 196]

# Tips Nodes
@test V[:TipsNodes] ≈ VR[1:197, -V.internalNodeNumbers .+ 196]

### R Code to get those results
# library(geiger)
# ## Load data caudata (salamanders)
# data("caudata")
# ## Save tree
# write.tree(caudata$phy, file = "caudata_tree.txt", append = FALSE,
#            digits = 10, tree.names = FALSE)
#
# ## Times shared
# V <- node.depth.edgelength(caudata$phy)
# prac <- mrca(caudata$phy, full = TRUE)
# V <- matrix(V[prac], dim(prac))
# write.table(V,
#             file = "caudata_shared_paths.txt",
#             sep = ",", row.names = FALSE,
#             col.names = TRUE)

###############################################################################
## Caudata dataset - BM
###############################################################################

## Export "caudata" dataset (from geiger)
phy = readTopology(joinpath(@__DIR__, "..", "examples", "caudata_tree.txt"));
dat = DataFrame!(CSV.File(joinpath(@__DIR__, "..", "examples", "caudata_trait.txt")));

## Fit a BM
fitBM = phyloNetworklm(@formula(trait ~ 1), dat, phy)

# Tests against results obtained with geiger::fitContinuous or phylolm::phylolm
@test loglikelihood(fitBM) ≈ -78.9611507833 atol=1e-10
@test dof(fitBM) ≈ 2.0 atol=1e-10
@test aic(fitBM) ≈ 161.9223015666 atol=1e-10
@test aicc(fitBM) ≈ 161.9841572367 atol=1e-10
@test coef(fitBM) ≈ [4.6789989001] atol=1e-8
@test vcov(fitBM) ≈ [0.1093144100] atol=1e-10
@test nobs(fitBM) ≈ 197.0 atol=1e-10
@test sum(residuals(fitBM)) ≈ -115.5767321312 atol=1e-8
@test dof_residual(fitBM) ≈ 196.0 atol=1e-10
@test sigma2_estim(fitBM) ≈ 0.0029452097 atol=1e-10
@test stderror(fitBM) ≈ [0.3306272978] atol=1e-10
@test confint(fitBM)[1] ≈ 4.0269551772 atol=1e-10
@test confint(fitBM)[2] ≈ 5.3310426231 atol=1e-10
tmp = predict(fitBM);
# below: print for when the test was broken
#@show length(tmp) # fixit: tests fail here
#@show tmp[1:6], tmp[190:end] # 4.679 except for last 5: 9.31707, 9.33832, 9.33052, 9.34356, 9.31707
#println("are they all 4.6789989001?")
# next: looks random. sometimes passes, most times fails
@test predict(fitBM) ≈ [4.6789989001 for i in 1:197] atol=1e-8
# @test_approx_eq_eps deviance(fitBM)  deviance(fitbis)
# @test_approx_eq_eps nulldeviance(fitBM)  nulldeviance(fitbis)
# @test_approx_eq_eps nullloglikelihood(fitBM)  nullloglikelihood(fitbis)
# @test_approx_eq_eps r2(fitBM)  r2(fitbis)
# @test_approx_eq_eps adjr2(fitBM)  adjr2(fitbis)
# @test_approx_eq_eps bic(fitBM)  bic(fitbis)
# @test_approx_eq_eps mu_estim(fitBM)  mu_estim(fitbis)

### Ancestral state reconstruction (with Rphylopars)
anc = (@test_logs (:warn, r"^These prediction intervals show uncertainty in ancestral values") ancestralStateReconstruction(fitBM));
ancR = DataFrame!(CSV.File(joinpath(@__DIR__, "..", "examples", "caudata_Rphylopars.txt")))

## Expectations
expe = expectations(anc)
expeR = ancR[!,:trait]
# Matching tips ?
tipsR = expeR[expe[197:393, :nodeNumber]]
tipsJulia = expe[197:393, :condExpectation]
@test tipsR ≈ tipsJulia
# Matching nodes ?
nodesR = expeR[-expe[1:196, :nodeNumber] .+ 196]
nodesJulia = expe[1:196, :condExpectation]
# below: print for when the test was broken
#@show nodesR[1:6],    nodesR[190:end]
#@show nodesJulia[1:6],nodesJulia[190:end]
@test isapprox(nodesR, nodesJulia)

## Variances
vars = diag(anc.variances_nodes)
# Rphylopars
varsR = ancR[!,:var]
# Matching nodes ?
nodesR = varsR[-expe[1:196, :nodeNumber] .+ 196]
@test nodesR ≈ vars atol=1e-3 ## RK: Small tol !!

### Ancestral state reconstruction (with Phytools)
ancRt = DataFrame!(CSV.File(joinpath(@__DIR__, "..", "examples", "caudata_Phytools.txt")));

## Expectations
expe = expectations(anc)
expeRt = ancRt[!,:trait]
# Matching nodes ?
nodesRt = expeRt[-expe[1:196, :nodeNumber] .+ (196 - 197)]
nodesJulia = expe[1:196, :condExpectation]
# below: print for when the test was broken
#@show nodesRt[1:6],   nodesRt[190:end]
#@show nodesJulia[1:6],nodesJulia[190:end]
@test isapprox(nodesRt, nodesJulia)

## Variances
vars = diag(anc.variances_nodes)
# Rphylopars
varsRt = ancRt[!,:var]
# Matching nodes ?
nodesRt = varsRt[-expe[1:196, :nodeNumber] .+ (196 - 197)]
@test nodesRt ≈ vars atol=2e-3 ## RK: Small tol !!

### Comparison between Rphylopars and Phytools:
@test nodesRt ≈ nodesR atol=0.003 ## RK: Small tol !!

### R script to get the above values:
# library(geiger)
#
# ## Load data caudata (salamanders)
# data("caudata")
#
# ## Save tree and data
# write.tree(caudata$phy, file = "caudata_tree.txt", append = FALSE,
#            digits = 10, tree.names = FALSE)
#
# write.table(data.frame(tipsNames = names(caudata$dat),
#                        trait = unname(caudata$dat)),
#             file = "caudata_trait.txt",
#             sep = ",", row.names = FALSE)
#
# ## Fit using Geiger
# fitgeiger <- fitContinuous(caudata$phy, caudata$dat, model = "BM")
#
# ## Fit using phylolm
# library(phylolm)
# fitphylolm <- phylolm(trait ~ 1,
#                       data.frame(trait = caudata$dat),
#                       caudata$phy, model = "BM")
#
# ## Fit using Rphylopars
# library(Rphylopars)
# fitphylopars <- phylopars(data.frame(species = names(caudata$dat),
#                           trait = unname(caudata$dat)),
#                           caudata$phy,
#                           pheno_error = FALSE,
#                           pheno_correlated = FALSE,
#                           REML = FALSE)
#
# # Save results of Rphylopars for ancestral trait reconstruction
# write.table(data.frame(trait = unname(fitphylopars$anc_recon),
#                        var = unname(fitphylopars$anc_var)),
#             file = "caudata_Rphylopars.txt",
#             sep = ",", row.names = FALSE)
#
# ## Ancestral State reconstruction using phytools
# library(phytools)
# fitphytools <- fastAnc(caudata$phy, caudata$dat, vars = TRUE)
#
# # Save results of Rphylopars for ancestral trait reconstruction
# write.table(data.frame(trait = unname(fitphytools$ace),
#                        var = unname(fitphytools$var)),
#             file = "caudata_Phytools.txt",
#             sep = ",", row.names = FALSE)
#
# ## Quantities to compare
# sprintf("%.10f", fitgeiger$opt$ln) # log likelihood
# sprintf("%.10f", fitphylolm$logLik)
# sprintf("%.10f", fitphylopars$logLik)
# sprintf("%.10f", fitgeiger$opt$aic) # aic
# sprintf("%.10f", fitphylolm$aic)
# sprintf("%.10f", fitgeiger$opt$aicc) # aicc
# sprintf("%.10f", fitphylolm$coefficients) # coef
# sprintf("%.10f", fitgeiger$opt$z0)
# sprintf("%.10f", fitphylopars$mu)
# sprintf("%.10f", fitphylolm$vcov) # vcov
# sprintf("%.10f", fitphylolm$n) # nobs
# sprintf("%.10f", sum(fitphylolm$residuals)) # residuals (sum)
# sprintf("%.10f", fitphylolm$n - fitphylolm$d) # df residuals
# sprintf("%.10f", fitphylolm$sigma2) # sigma 2
# sprintf("%.10f", fitgeiger$opt$sigsq)
# sprintf("%.10f", fitphylopars$pars$phylocov)
# sprintf("%.10f", summary(fitphylolm)$coefficients[2]) # std error
# sprintf("%.10f", summary(fitphylolm)$df) # df
# sprintf("%.10f", fitgeiger$opt$k)
# sprintf("%.10f", coef(fitphylolm) + summary(fitphylolm)$coefficients[2] * qt(0.025, 196))
# sprintf("%.10f", coef(fitphylolm) + summary(fitphylolm)$coefficients[2] * qt(0.975, 196))
# sprintf("%.10f", predict(fitphylolm)) # df

###############################################################################
## Caudata dataset - Pagel's lambda
###############################################################################

## Fit Pagel's lambda
fitLambda = (@test_logs (:info, r"^Maximum lambda value") match_mode=:any phyloNetworklm(@formula(trait ~ 1), dat, phy, model = "lambda"));

@test lambda_estim(fitLambda) ≈ 0.9193 atol=1e-4 # Due to convergence issues, tolerance is lower.
@test loglikelihood(fitLambda) ≈ -51.684379 atol=1e-6
@test dof(fitLambda) ≈ 3.0 atol=1e-10
@test aic(fitLambda) ≈ 109.368759 atol=1e-6
@test aicc(fitLambda) ≈ 109.493111 atol=1e-6
@test coef(fitLambda) ≈ [4.66893] atol=1e-5
@test vcov(fitLambda) ≈ [0.05111] atol=1e-5
@test nobs(fitLambda) ≈ 197.0 atol=1e-10
# below: print for when the test was broken
#@show sum(residuals(fitLambda)) # locally: -115.91591040367894
#println("is this -113.594?")
@test isapprox(sum(residuals(fitLambda)), -113.594, atol=1e-2) ## Low Tolerance !!
@test dof_residual(fitLambda) ≈ 196.0 atol=1e-10 ## Correct Definition ?
@test sigma2_estim(fitLambda) ≈ 0.0014756 atol=1e-7
@test stderror(fitLambda) ≈ [0.22608] atol=1e-5
@test confint(fitLambda)[1] ≈ 4.2230 atol=1e-4
@test confint(fitLambda)[2] ≈ 5.114 atol=1e-3
tmp = predict(fitLambda);
# below: print for when the test was broken
#@show length(tmp)
#@show tmp[1:6], tmp[190:end] # all 4.66893 except for last 5: 8.42676, 8.44585, etc.
#println("are they all 4.66893?")
# next: looks random. sometimes passes, most times fails
@test predict(fitLambda) ≈ [4.66893 for i in 1:197] atol=1e-5 norm=x->norm(x,Inf)

### R script to get the above values:
# library(geiger)
#
# ## Load data caudata (salamanders)
# data("caudata")
#
# ## Fit using Geiger
# fitgeiger <- fitContinuous(caudata$phy, caudata$dat, model = "lambda")
#
# ## Fit using phylolm
# library(phylolm)
# fitphylolm <- phylolm(trait ~ 1,
#                       data.frame(trait = caudata$dat),
#                       caudata$phy, model = "lambda",
#                       starting.value = 0.9)
#
# ## Fit using Rphylopars
# library(Rphylopars)
# fitphylopars <- phylopars(data.frame(species = names(caudata$dat),
#                           trait = unname(caudata$dat)),
#                           caudata$phy,
#                           model = "lambda",
#                           model_par_start = 0.9,
#                           pheno_error = FALSE,
#                           pheno_correlated = FALSE,
#                           REML = FALSE)
#
# ## Quantities to compare
# sprintf("%.10f", fitgeiger$opt$ln) # log likelihood
# sprintf("%.10f", fitphylolm$logLik)
# sprintf("%.10f", fitphylopars$logLik)
# sprintf("%.10f", fitgeiger$opt$aic) # aic
# sprintf("%.10f", fitphylolm$aic)
# sprintf("%.10f", fitgeiger$opt$aicc) # aicc
# sprintf("%.10f", fitphylolm$coefficients) # coef
# sprintf("%.10f", fitgeiger$opt$z0)
# sprintf("%.10f", fitphylopars$mu)
# sprintf("%.10f", fitphylolm$vcov) # vcov
# sprintf("%.10f", fitphylolm$n) # nobs
# sprintf("%.10f", sum(fitphylolm$residuals)) # residuals (sum)
# sprintf("%.10f", fitphylolm$n - fitphylolm$d) # df residuals
# sprintf("%.10f", fitphylolm$sigma2) # sigma 2
# sprintf("%.10f", fitgeiger$opt$sigsq)
# sprintf("%.10f", fitphylopars$pars$phylocov)
# sprintf("%.10f", summary(fitphylolm)$coefficients[2]) # std error
# sprintf("%.10f", summary(fitphylolm)$df) # df
# sprintf("%.10f", fitgeiger$opt$k)
# sprintf("%.10f", coef(fitphylolm) + summary(fitphylolm)$coefficients[2] * qt(0.025, 196))
# sprintf("%.10f", coef(fitphylolm) + summary(fitphylolm)$coefficients[2] * qt(0.975, 196))
# sprintf("%.10f", predict(fitphylolm)) # df

###############################################################################
## Caudata dataset - BM with shifts
###############################################################################

## Export "caudata" dataset (from geiger)
phy = readTopology(joinpath(@__DIR__, "..", "examples", "caudata_tree.txt"));
dat = DataFrame!(CSV.File(joinpath(@__DIR__, "..", "examples", "caudata_trait.txt")));

## Add some shifts in the model
df_shift = regressorShift(phy.edge[[98, 326, 287]], phy)
dat = innerjoin(dat, df_shift, on=:tipNames)
## Fit a BM
fitBM = phyloNetworklm(@formula(trait ~ shift_98 + shift_326 + shift_287), dat, phy)

# Tests against results obtained with geiger::fitContinuous or phylolm::phylolm
@test loglikelihood(fitBM) ≈ -76.1541605207 atol=1e-10
@test dof(fitBM) ≈ 5.0 atol=1e-10
@test aic(fitBM) ≈ 162.3083210414 atol=1e-10
@test coef(fitBM) ≈ [4.8773804547 -0.6664678924 -1.0692950215 -0.0159349201]' atol=1e-10
vcovR = [ 0.1165148753  -0.0431446679 -0.0305707092 0.0000000000
         -0.0431446679 0.3212796435  0.0340202694  -0.0000000000
         -0.0305707092 0.0340202694 0.2472613252  -0.0236868278
          0.0000000000  -0.0000000000 -0.0236868278 0.1263794810]
@test vcov(fitBM) ≈ vcovR atol=1e-8
@test nobs(fitBM) ≈ 197.0 atol=1e-10
@test sum(residuals(fitBM)) ≈ -10.1752334428 atol=1e-10
@test dof_residual(fitBM) ≈ 193.0 atol=1e-10
@test sigma2_estim(fitBM) ≈ 0.0028624636 atol=1e-10
@test stderror(fitBM) ≈ [0.3413427535 0.5668153522 0.4972537835 0.3554989184]' atol=1e-10
@test confint(fitBM)[:,1] ≈ [4.2041393297 -1.7844157659 -2.0500444097 -0.7170966977]' atol=1e-10
@test confint(fitBM)[:,2] ≈ [5.5506215797 0.4514799811 -0.0885456333 0.6852268574]' atol=1e-10
predictR = [4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 4.2109125624, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.7921505131, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 3.8080854332, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547, 4.8773804547]
tmp  = predict(fitBM)
tmp2 = predictR[fitBM.model.ind]
# below: print for when the test was broken
#@show tmp[1:6],  tmp[190:end]
#@show tmp2[1:6],tmp2[190:end] # the last 5 values are different
@test isapprox(predict(fitBM), predictR[fitBM.model.ind], atol=1e-8)

# ## R code to get those results
# library(geiger)
# ## Load data caudata (salamanders)
# data("caudata")
# ## Re-order with the tree
# dat <- caudata$dat[match(caudata$phy$tip.label, names(caudata$dat))]
# ## Incidence matrix
# Tm = PhylogeneticEM::incidence.matrix(caudata$phy) + 0
# ## Fit using phylolm
# library(phylolm)
# fitphylolm <- phylolm(trait ~ shift_37 + shift_105 + shift_222,
#                       data.frame(trait = dat,
#                                  shift_37 = Tm[, 37],
#                                  shift_105 = Tm[, 105],
#                                  shift_222 = Tm[, 222]),
#                       caudata$phy, model = "BM")
# ## Quantities to compare
# sprintf("%.10f", fitphylolm$logLik)
# sprintf("%.10f", fitphylolm$aic)
# sprintf("%.10f", fitphylolm$coefficients) # coef
# sprintf("%.10f", fitphylolm$vcov) # vcov
# sprintf("%.10f", fitphylolm$n) # nobs
# sprintf("%.10f", sum(fitphylolm$residuals)) # residuals (sum)
# sprintf("%.10f", fitphylolm$n - fitphylolm$d) # df residuals
# sprintf("%.10f", fitphylolm$sigma2) # sigma 2
# sprintf("%.10f", summary(fitphylolm)$coefficients[,2]) # std error
# sprintf("%.10f", summary(fitphylolm)$df) # df
# sprintf("%.10f", coef(fitphylolm) + summary(fitphylolm)$coefficients[,2] * qt(0.025, 193))
# sprintf("%.10f", coef(fitphylolm) + summary(fitphylolm)$coefficients[,2] * qt(0.975, 193))
# sprintf("%.10f", predict(fitphylolm))

end

###############################################################################
## Lizard dataset - BM
###############################################################################

@testset "phyloNetworklm: Lizard Dataset" begin

## Export "lizard" dataset (Mahler et al 2013)
phy = readTopology(joinpath(@__DIR__, "..", "examples", "lizard_tree.txt"));
dat = DataFrame!(CSV.File(joinpath(@__DIR__, "..", "examples", "lizard_trait.txt")));
categorical!(dat, :region)

## Fit a BM
fitBM = phyloNetworklm(@formula(AVG_SVL ~ AVG_ltoe_IV + AVG_lfing_IV * region), dat, phy)

# Tests against results obtained with geiger::fitContinuous or phylolm::phylolm
@test loglikelihood(fitBM) ≈ 105.17337853473711 atol=1e-10
@test dof(fitBM) ≈ 10.0 atol=1e-10
@test aic(fitBM) ≈ -190.3467570695 atol=1e-10
# @test aicc(fitBM)
@test coef(fitBM) ≈ [2.7925712673 -0.2010704391 0.9832555589 -0.1021226296 -0.3703658712 0.1557471731 0.0374549036 0.1805667675 -0.0495767233]' atol=1e-10
vcovR =  [0.0200086273  -0.0136717540 0.0084815090  -0.0093192029 -0.0114417825 -0.0113346813 0.0041102304  0.0053787287  0.0050521693
          -0.0136717540 0.0185396965  -0.0169114682 0.0020645005  0.0036352899  0.0026227856 -0.0012281620 -0.0018838231 -0.0014800242
          0.0084815090  -0.0169114682 0.0174647413  0.0016403871  0.0005624488  0.0015003301 -0.0005714235 -0.0003201257 -0.0005423354
          -0.0093192029 0.0020645005  0.0016403871  0.0167953394  0.0078012534  0.0086329399 -0.0080782771 -0.0037333495 -0.0039327836
          -0.0114417825 0.0036352899  0.0005624488  0.0078012534  0.0490482083  0.0092203882 -0.0033670465 -0.0191567265 -0.0040068947
          -0.0113346813 0.0026227856  0.0015003301  0.0086329399  0.0092203882  0.0331395502 -0.0037513830 -0.0041592743 -0.0146108207
          0.0041102304  -0.0012281620 -0.0005714235 -0.0080782771 -0.0033670465 -0.0037513830 0.0045172675  0.0018165174  0.0020857846
          0.0053787287  -0.0018838231 -0.0003201257 -0.0037333495 -0.0191567265 -0.0041592743 0.0018165174  0.0093292284  0.0020427637
          0.0050521693  -0.0014800242 -0.0005423354 -0.0039327836 -0.0040068947 -0.0146108207 0.0020857846  0.0020427637  0.0074817942]
@test vcov(fitBM) ≈ vcovR atol=1e-8
@test nobs(fitBM) ≈ 100.0 atol=1e-10
# below: print for when the test was broken (with @test_skip)
#@show sum(residuals(fitBM)) # looks random, e.g. 1.6091413520064477, or 0.8338189090359597
#println("is this equal to 0.6352899255?") # sometimes NO, yet the test passes below!!
@test sum(residuals(fitBM)) ≈ 0.6352899255 atol=1e-10
@test dof_residual(fitBM) ≈ 91.0 atol=1e-10
@test sigma2_estim(fitBM) ≈ 0.0003025014 atol=1e-10
@test stderror(fitBM) ≈ [0.1414518551,0.1361605540,0.1321542330,0.1295968341,0.2214683008,0.1820427154,0.0672106202,0.0965879311,0.0864973651] atol=1e-10
@test confint(fitBM)[:,1] ≈ [2.5115945339,-0.4715366529,0.7207474097,-0.3595508202,-0.8102854443,-0.2058583178,-0.0960507369,-0.0112932922,-0.2213931131] atol=1e-10 norm=x->norm(x,Inf)
@test confint(fitBM)[:,2] ≈ [3.0735480006,0.0693957746,1.2457637082,0.1553055609,0.0695537019,0.5173526640,0.1709605441,0.3724268272,0.1222396666] atol=1e-10
# @test_approx_eq_eps predict(fitBM)[fitbis.model.ind] predict(fitbis)
# @test_approx_eq_eps deviance(fitBM)  deviance(fitbis)
# @test_approx_eq_eps nulldeviance(fitBM)  nulldeviance(fitbis)
# @test_approx_eq_eps nullloglikelihood(fitBM)  nullloglikelihood(fitbis)
# @test_approx_eq_eps r2(fitBM)  r2(fitbis)
# @test_approx_eq_eps adjr2(fitBM)  adjr2(fitbis)
# @test_approx_eq_eps bic(fitBM)  bic(fitbis)
# @test_approx_eq_eps mu_estim(fitBM)  mu_estim(fitbis)

### R script to get the above values
# ## Data
# dat <- read.csv(file = "GA_Anolis_traits.csv")
# geo <- read.csv(file = "GA_Anolis_biogeography.csv")
#
# dat <- merge(dat, geo, by = "species")
# # keep only toe and hand length
# dat <- dat[, c("species", "AVG.SVL", "AVG.ltoe.IV", "AVG.lfing.IV", "region")]
# colnames(dat)[1] <- "tipsNames"
#
# write.table(dat,
#             file = "lizard_trait.txt",
#             sep = ",", row.names = FALSE)
#
# ## Tree
# phy <- read.tree(file = "GA_Anolis_MCC.tre")
#
# write.tree(phy, file = "lizards_tree.txt", append = FALSE,
#            digits = 10, tree.names = FALSE)
#
# rownames(dat) <- dat$tipsNames
# dat <- dat[, -1]
# dat$region <- as.factor(dat$region)
#
# ## Fit
# fitphylolm <- phylolm(AVG.SVL ~ 1 + AVG.ltoe.IV + AVG.lfing.IV * region, dat, phy, model = "BM")
# ## Quantities to compare
# sprintf("%.10f", fitphylolm$logLik)
# sprintf("%.10f", summary(fitphylolm)$df) # df
# sprintf("%.10f", fitphylolm$aic)
# sprintf("%.10f", fitphylolm$coefficients) # coef
# matrix(sprintf("%.10f", fitphylolm$vcov), 9, 9) # vcov
# sprintf("%.10f", fitphylolm$n) # nobs
# sprintf("%.10f", sum(fitphylolm$residuals)) # residuals (sum)
# sprintf("%.10f", fitphylolm$n - fitphylolm$d) # df residuals
# sprintf("%.10f", fitphylolm$sigma2) # sigma 2
# sprintf("%.10f", summary(fitphylolm)$coefficients[,2]) # std error
# sprintf("%.10f", coef(fitphylolm) + summary(fitphylolm)$coefficients[, 2] * qt(0.025, 91))
# sprintf("%.10f", coef(fitphylolm) + summary(fitphylolm)$coefficients[, 2] * qt(0.975, 91))

###############################################################################
## Lizard dataset - lambda
###############################################################################

## Fit lambda
fitLambda = (@test_logs (:info, r"^Maximum lambda value") match_mode=:any phyloNetworklm(@formula(AVG_SVL ~ AVG_ltoe_IV + AVG_lfing_IV * region), dat, phy, model = "lambda"))

# Tests against results obtained with geiger::fitContinuous or phylolm::phylolm
@test lambda_estim(fitLambda) ≈ 0.9982715594 atol=1e-5
@test loglikelihood(fitLambda) ≈ 105.1769275564 atol=1e-8
@test dof(fitLambda) ≈ 11.0 atol=1e-10
@test aic(fitLambda) ≈ -188.3538551128 atol=1e-8
# @test aicc(fitBM)
@test coef(fitLambda) ≈ [2.7940573420 -0.2066584606 0.9897083949 -0.1004840950 -0.3677991157 0.1576743022 0.0367633665 0.1792502383 -0.0505291142]' atol=1e-5
vcovR =  [0.0200251600  -0.0137474015 0.0085637021  -0.0092973836 -0.0114259722 -0.0113056243 0.0041037877  0.0053740100  0.0050429112
          -0.0137474015 0.0186885224  -0.0170645512 0.0020509207  0.0036334103 0.0026066694  -0.0012237488 -0.0018826836 -0.0014743137
          0.0085637021  -0.0170645512 0.0176200143  0.0016494733  0.0005604169 0.0015116125  -0.0005735573 -0.0003192320 -0.0005457420
          -0.0092973836 0.0020509207  0.0016494733  0.0167461876  0.0077885115 0.0086173037  -0.0080563819 -0.0037287856 -0.0039275469
          -0.0114259722 0.0036334103  0.0005604169  0.0077885115  0.0490092393 0.0092036032  -0.0033631662 -0.0191657329 -0.0040017905
          -0.0113056243 0.0026066694  0.0015116125  0.0086173037  0.0092036032 0.0330248707  -0.0037465110 -0.0041543671 -0.0145663751
          0.0041037877  -0.0012237488 -0.0005735573 -0.0080563819 -0.0033631662 -0.0037465110 0.0045042057  0.0018142470  0.0020823721
          0.0053740100  -0.0018826836 -0.0003192320 -0.0037287856 -0.0191657329 -0.0041543671 0.0018142470  0.0093334212  0.0020404652
          0.0050429112  -0.0014743137 -0.0005457420 -0.0039275469 -0.0040017905 -0.0145663751 0.0020823721  0.0020404652  0.0074600880]
@test vcov(fitLambda) ≈ vcovR atol=3e-7
@test nobs(fitLambda) ≈ 100.0 atol=1e-10
# below: print for when the test was broken (with @test_skip)
#@show sum(residuals(fitLambda)) # looks random, eg 0.033126277561337916 or 0.644941961666333
#println("is this equal to 0.6369008979?")
@test sum(residuals(fitLambda)) ≈ 0.6369008979 atol=1e-5
@test dof_residual(fitLambda) ≈ 91.0 atol=1e-10
@test sigma2_estim(fitLambda) ≈ 0.0003009914 atol=1e-9
@test stderror(fitLambda) ≈ [0.1415102824,0.1367059706,0.1327404019,0.1294070617,0.2213803048,0.1817274626,0.0671133793,0.0966096332,0.0863718011] atol=1e-6
@test confint(fitLambda)[:,1] ≈ [2.5129645499,-0.4782080775,0.7260358930,-0.3575353260,-0.8075438955,-0.2033049779,-0.0965491169,-0.0126529301,-0.2220960868] atol=1e-5
@test confint(fitLambda)[:,2] ≈ [3.0751501341,0.0648911562,1.2533808968,0.1565671360,0.0719456641,0.5186535822,0.1700758500,0.3711534067,0.1210378584] atol=1e-5
# @test_approx_eq_eps predict(fitLambda)[fitbis.model.ind] predict(fitbis)
# @test_approx_eq_eps deviance(fitLambda)  deviance(fitbis)
# @test_approx_eq_eps nulldeviance(fitLambda)  nulldeviance(fitbis)
# @test_approx_eq_eps nullloglikelihood(fitLambda)  nullloglikelihood(fitbis)
# @test_approx_eq_eps r2(fitLambda)  r2(fitbis)
# @test_approx_eq_eps adjr2(fitLambda)  adjr2(fitbis)
# @test_approx_eq_eps bic(fitLambda)  bic(fitbis)
# @test_approx_eq_eps mu_estim(fitLambda)  mu_estim(fitbis)

### R script to get the above values
# ## Data
# dat <- read.csv(file = "GA_Anolis_traits.csv")
# geo <- read.csv(file = "GA_Anolis_biogeography.csv")
#
# dat <- merge(dat, geo, by = "species")
# # keep only toe and hand length
# dat <- dat[, c("species", "AVG.SVL", "AVG.ltoe.IV", "AVG.lfing.IV", "region")]
# colnames(dat)[1] <- "tipsNames"
#
# write.table(dat,
#             file = "lizard_trait.txt",
#             sep = ",", row.names = FALSE)
#
# ## Tree
# phy <- read.tree(file = "GA_Anolis_MCC.tre")
#
# write.tree(phy, file = "lizards_tree.txt", append = FALSE,
#            digits = 10, tree.names = FALSE)
#
# rownames(dat) <- dat$tipsNames
# dat <- dat[, -1]
# dat$region <- as.factor(dat$region)
#
# ## Fit
# fitphylolm <- phylolm(AVG.SVL ~ 1 + AVG.ltoe.IV + AVG.lfing.IV * region, dat, phy, model = "lambda")
# ## Quantities to compare
# sprintf("%.10f", fitphylolm$logLik)
# sprintf("%.10f", summary(fitphylolm)$df) # df
# sprintf("%.10f", fitphylolm$aic)
# sprintf("%.10f", fitphylolm$coefficients) # coef
# matrix(sprintf("%.10f", fitphylolm$vcov), 9, 9) # vcov
# sprintf("%.10f", fitphylolm$n) # nobs
# sprintf("%.10f", sum(fitphylolm$residuals)) # residuals (sum)
# sprintf("%.10f", fitphylolm$n - fitphylolm$d) # df residuals
# sprintf("%.10f", fitphylolm$sigma2) # sigma 2
# sprintf("%.10f", summary(fitphylolm)$coefficients[,2]) # std error
# sprintf("%.10f", coef(fitphylolm) + summary(fitphylolm)$coefficients[, 2] * qt(0.025, 91))
# sprintf("%.10f", coef(fitphylolm) + summary(fitphylolm)$coefficients[, 2] * qt(0.975, 91))

end


text:
using PyCall

@pyimport __builtin__
@pyimport Bio.SeqIO as SeqIO
@pyimport Bio.Seq as Seq
@pyimport Bio.Alphabet as IUPAC

function translate_sixframes(filenames, format="fastq")
	for filename in filenames
		file = open(filename)
		output_file = open(filename*".faa", "w")

		for record in SeqIO.parse(file, format)
			nt_sequence = record[:seq]
			sequence_id = record[:id]

			write(output_file, ">"*sequence_id*".0\n"*__builtin__.str(pycall(nt_sequence[:translate],PyAny, table=11))*"\n")

			frame_1_str = __builtin__.str(nt_sequence)[2:end]
			frame_1_seq = Seq.Seq(frame_1_str, IUPAC.generic_dna)
			write(output_file, ">"*sequence_id*".1\n"*__builtin__.str(pycall(frame_1_seq[:translate],PyAny, table=11))*"\n")

			frame_2_str = __builtin__.str(nt_sequence)[3:end]
			frame_2_seq = Seq.Seq(frame_2_str, IUPAC.generic_dna)
			write(output_file, ">"*sequence_id*".2\n"*__builtin__.str(pycall(frame_2_seq[:translate],PyAny, table=11))*"\n")

			rc_nt_sequence = pycall(nt_sequence[:reverse_complement], PyAny)
			write(output_file, ">"*sequence_id*".0rc\n"*__builtin__.str(pycall(rc_nt_sequence[:translate],PyAny, table=11))*"\n")

			frame_1rc_str = __builtin__.str(rc_nt_sequence)[2:end]
			frame_1rc_seq = Seq.Seq(frame_1rc_str, IUPAC.generic_dna)
			write(output_file, ">"*sequence_id*".1rc\n"*__builtin__.str(pycall(frame_1rc_seq[:translate],PyAny, table=11))*"\n")

			frame_2rc_str = __builtin__.str(rc_nt_sequence)[3:end]
			frame_2rc_seq = Seq.Seq(frame_2rc_str, IUPAC.generic_dna)
			write(output_file, ">"*sequence_id*".2rc\n"*__builtin__.str(pycall(frame_2rc_seq[:translate],PyAny, table=11))*"\n")

		end

	end
end

ext: jl

max_stars_repo_name: ianpgm/apfamet

max_issues_repo_name: ianpgm/apfamet

max_forks_repo_name: ianpgm/apfamet


text:
""" Tensor networks via the operad of undirected wiring diagrams.
"""
module TensorNetworks
export RelationDiagram, @tensor_network, parse_tensor_network,
  contract_tensor_network, @contract_tensors_with, gen_tensor_notation

using MLStyle: @match

using ...CategoricalAlgebra.CSets
using ...WiringDiagrams.UndirectedWiringDiagrams
import ...WiringDiagrams: oapply
using ...Programs.RelationalPrograms: RelationDiagram

# Evaluation
############

""" Contract tensors according to UWD.

This method simply wraps [`contract_tensor_network`](@ref), which also presents
an interface closer to the standard tensor algebra packages.
"""
oapply(d::UndirectedWiringDiagram, tensors::AbstractVector{<:AbstractArray}) =
  contract_tensor_network(d, tensors)

""" Contract a tensor network all at once.

Performs a generalized contraction of a list of tensors of arbitrary rank. In
the binary case, this function includes matrix multiplication, tensor product,
and trace as special cases. In general, it implements the "generalized matrix
multiplication" described in [Spivak et al's "Pixel
Arrays"](https://arxiv.org/abs/1609.00061).

Be warned that this method is inefficient for contracting networks with more
than two junctions. Specifically, if the network has ``n`` junctions with
dimensions ``k₁, …, kₙ``, then this function performs exactly ``k₁ ⋯ kₙ``
multiplications in the base ring (or rig). For example, a binary matrix multiply
of square matrices performs ``k³`` multiplications. A ternary matrix multiply
performs ``k⁴`` multiplications, whereas a sequence of two binary multiplies
performs ``2k³``. Thus, when there are more than two junctions, it is usually
better to schedule a computation using [`schedule`](@ref) or another package,
possibly via [`@contract_tensors_with`](@ref).

In addition to the method taking an undirected wiring diagram, this function has
methods similar to the `ncon` ("Network CONtractor") function in:

- the [NCON package](https://arxiv.org/abs/1402.0939) for MATLAB
- the Julia package [TensorOperations.jl](https://github.com/Jutho/TensorOperations.jl)

except that here the outer junctions are represented explictly by a third
argument rather than implicitly by using negative numbers in the second
argument.
"""
function contract_tensor_network(d::UndirectedWiringDiagram,
                                 tensors::AbstractVector{<:AbstractArray})
  @assert nboxes(d) == length(tensors)
  juncs = [junction(d, ports(d, b)) for b in boxes(d)]
  j_out = junction(d, ports(d, outer=true), outer=true)
  contract_tensor_network(tensors, juncs, j_out)
end

function contract_tensor_network(tensors::AbstractVector{<:AbstractArray{T}},
                                 juncs::AbstractVector, j_out) where T
  # Handle important binary case with specialized code.
  if length(tensors) == 2 && length(juncs) == 2
    return contract_tensor_network(Tuple(tensors), Tuple(juncs), j_out)
  end

  jsizes = Tuple(infer_junction_sizes(tensors, juncs, j_out))
  juncs, j_out = map(Tuple, juncs), Tuple(j_out)
  C = zeros(T, Tuple(jsizes[j] for j in j_out))
  for index in CartesianIndices(jsizes)
    x = one(T)
    for (A, junc) in zip(tensors, juncs)
      x *= A[(index[j] for j in junc)...]
    end
    C[(index[j] for j in j_out)...] += x
  end
  return C
end

function contract_tensor_network( # Binary case.
    (A, B)::Tuple{<:AbstractArray{T},<:AbstractArray{T}},
    (jA, jB), j_out) where T
  jsizes = Tuple(infer_junction_sizes((A, B), (jA, jB), j_out))
  jA, jB, j_out = Tuple(jA), Tuple(jB), Tuple(j_out)
  C = zeros(T, Tuple(jsizes[j] for j in j_out))
  for index in CartesianIndices(jsizes)
    C[(index[j] for j in j_out)...] +=
      A[(index[j] for j in jA)...] * B[(index[j] for j in jB)...]
  end
  return C
end

function infer_junction_sizes(tensors, juncs, j_out)
  @assert length(tensors) == length(juncs)
  njunc = maximum(Iterators.flatten((Iterators.flatten(juncs), j_out)))
  jsizes = fill(-1, njunc)
  for (A, junc) in zip(tensors, juncs)
    for (i, j) in enumerate(junc)
      if jsizes[j] == -1
        jsizes[j] = size(A, i)
      else
        @assert jsizes[j] == size(A, i)
      end
    end
  end
  @assert all(s >= 0 for s in jsizes)
  jsizes
end

# Parsing and code generation
#############################

""" Construct an undirected wiring diagram using tensor notation.

The syntax is compatible with the Einstein-style tensor notation used by tensor
algebra packages like
[TensorOperations.jl](https://github.com/Jutho/TensorOperations.jl) and
[TensorCast.jl](https://github.com/mcabbott/TensorCast.jl). For example, the
wiring diagram for the composite of two matrices (or two binary relations) is
constructed by

```julia
@tensor_network (i,j,k) C[i,k] := A[i,j] * B[j,k]
```

The leading context, or list of variables, may be omitted, in which case it is
inferred from the variables used in the tensor expression. So, in this example,
an equivalent macro call is

```julia
@tensor_network C[i,k] := A[i,j] * B[j,k]
```

See also: [`gen_tensor_notation`](@ref), the "inverse" to this macro.
"""
macro tensor_network(exprs...)
  :(parse_tensor_network($((QuoteNode(expr) for expr in exprs)...)))
end

""" Parse an undirected wiring diagram from a tensor expression.

For more information, see the corresponding macro [`@tensor_network`](@ref).
"""
function parse_tensor_network(context::Expr, expr::Expr)
  all_vars = @match context begin
    Expr(:tuple, args...) => collect(Symbol, args)
    Expr(:vect, args...) => collect(Symbol, args)
  end
  parse_tensor_network(expr, all_vars=all_vars)
end

function parse_tensor_network(expr::Expr; all_vars=nothing)
  # Parse tensor expression.
  (outer_name, outer_vars), body = @match expr begin
    Expr(:(=), outer, body) => (parse_tensor_term(outer), body)
    Expr(:(:=), outer, body) => (parse_tensor_term(outer), body)
    _ => error("Tensor expression $expr must be an assignment, either = or :=")
  end
  names_and_vars = map(parse_tensor_term, @match body begin
    Expr(:call, :(*), args...) => args
    1 => [] # No terms
    arg => [arg] # One term
  end)

  # Check compatibility of used variables with declared variables.
  used_vars = unique!(reduce(vcat, ([[outer_vars]; last.(names_and_vars)])))
  if isnothing(all_vars)
    all_vars = sort!(used_vars)
  else
    used_vars ⊆ all_vars || error("One of variables $used_vars is not declared")
  end

  # Construct the undirected wiring diagram.
  d = RelationDiagram{Symbol}(length(outer_vars))
  add_junctions!(d, length(all_vars), variable=all_vars)
  set_junction!(d, ports(d, outer=true),
                incident(d, outer_vars, :variable), outer=true)
  for (name, vars) in names_and_vars
    box = add_box!(d, length(vars), name=name)
    set_junction!(d, ports(d, box), incident(d, vars, :variable))
  end
  return d
end

function parse_tensor_term(expr)
  @match expr begin
    Expr(:ref, name::Symbol, args...) => (name, collect(Symbol, args))
    name::Symbol => (name, Symbol[]) # Scalar
    _ => error("Invalid syntax in term $expr in tensor expression")
  end
end

""" Contract a tensor network using a macro from another package.

This macro takes two arguments: an undirected wiring diagram and a macro call
supporting the tensor contraction notation that is de facto standard among Julia
tensor packages. For example, to evaluate a tensor network using the `@tullio`
macro, use:

```julia
A = @contract_tensors_with @tullio diagram
```

The following macros should work:

- `@tensor` from
  [TensorOperations.jl](https://github.com/Jutho/TensorOperations.jl).
- `@tullio` from [Tullio.jl](https://github.com/mcabbott/Tullio.jl)
- `@einsum` from [Einsum.jl](https://github.com/ahwillia/Einsum.jl)
- `@ein` from [OMEinsum.jl](https://github.com/under-Peter/OMEinsum.jl)

However, the macros `@cast` and `@reduce` from
[TensorCast.jl](https://github.com/mcabbott/TensorCast.jl) will *not* work
because they do not support implicit summation.
"""
macro contract_tensors_with(diagram, macro_expr)
  # XXX: We cannot use `GeneralizedGenerated.mk_function` here because packages
  # like Tullio generate code with type parameters, which GG does not allow.
  compile_expr = :(gen_tensor_notation($(esc(diagram)),
    assign_op=:(:=), assign_name=gensym("out")))
  Expr(:call, esc(:eval),
       :(_contract_tensors_with($compile_expr, $(QuoteNode(macro_expr)))))
end
function _contract_tensors_with(tensor_expr, macro_expr)
  @match macro_expr begin
    Expr(:macrocall, args...) => Expr(:macrocall, args..., tensor_expr)
    _ => error("Expression $macro_expr is not a macro call")
  end
end

""" Generate Julia expression in tensor notation from undirected wiring diagram.

This function is used to implement [`@contract_tensors_with`](@ref) but may be
useful in its own right.

See also [`@tensor_network`](@ref), the "inverse" to this function.
"""
function gen_tensor_notation(d::UndirectedWiringDiagram;
    assign_op::Symbol=:(=), assign_name::Symbol=:out)
  vars = j -> subpart(d, j, :variable)
  outer_vars = vars(junction(d, ports(d, outer=true), outer=true))
  terms = map(boxes(d)) do box
    ref_expr(subpart(d, box, :name), vars(junction(d, ports(d, box))))
  end
  lhs = ref_expr(assign_name, outer_vars)
  rhs = if isempty(terms); 1
    elseif length(terms) == 1; first(terms)
    else Expr(:call, :(*), terms...) end
  Expr(assign_op, lhs, rhs)
end

function ref_expr(name::Symbol, vars)
  isempty(vars) ? name : Expr(:ref, name, vars...)
end

end


text:
using InteractionNetworkModels, Plots, Distributions
using StatsPlots, Plots.Measures, StatsBase
using Distances, StructuredDistances

V = 50
Sᵐ = [rand(1:V,rand(2:5)) for i in 1:10]
d = FastEditDistance(FastLCS(100), 100)
γ = 6.0
K_inner, K_outer = (DimensionRange(1,10), DimensionRange(1,25))
# K_inner, K_outer = (Inf,Inf)
model = SIS(
    Sᵐ, γ,
    d,
    1:V,
    K_inner, K_outer)

d(Sᵐ,Sᵐ)

mcmc_sampler = SisMcmcInsertDelete(
    ν_ed=1, ν_td=1, β=0.7, len_dist=TrGeometric(0.9,1,K_inner.u),
    lag=30, burn_in=1000,
    K=200
)

@time test = mcmc_sampler(
    model, desired_samples=50000, burn_in=0, lag=1)

plot(test)
summaryplot(test)
test.sample

@time mcmc_out = mcmc_sampler(
    model,
    desired_samples=100,
    lag=500,
    burn_in=10000
    )

plot(mcmc_out)
summaryplot(mcmc_out)

data = mcmc_out.sample
S_prior = SIS(Sᵐ, 0.1, model.dist, model.V, model.K_inner, model.K_outer)
γ_prior = Uniform(0.5,10.0)

posterior = SisPosterior(data, S_prior, γ_prior)

posterior_sampler = SisIexInsertDelete(
    mcmc_sampler,
    ν_ed=1, ν_td=1, len_dist=TrGeometric(0.8, 1, K_inner.u),
    ε=0.2, β=0.7,
    K=200,
    α=0.0,
    desired_samples=50, burn_in=0, lag=1
)

# Mode Conditional
S_init, ind = sample_frechet_mean(posterior.data, posterior.dist)
d(S_init, Sᵐ)
initialiser = InitRandEdit(2)

@time posterior_out = posterior_sampler(
    posterior,
    S_init=S_init,
    γ_init=4.0,
    # aux_init_at_prev=true,
    desired_samples=100
);

posterior_out
posterior_out.S_sample
plot(posterior_out, Sᵐ)

# Gamma conditional

S_fix = Sᵐ
posterior_out = posterior_sampler(
    posterior,
    S_fix,
    desired_samples=500, lag=1, burn_in=0,
    γ_init=4.6
)

plot(posterior_out)

# Joint

using BenchmarkTools

@time posterior_out = posterior_sampler(
    posterior,
    desired_samples=500, lag=2, burn_in=1000,
    S_init = S_init, γ_init=6.0
)

posterior_out
posterior_out.S_sample
plot(posterior_out, Sᵐ)

# Predictives
# -----------

# Distances to mode
dist_pred = DistanceToModePredictive(posterior_out)
outer_dim_pred = OuterDimensionPredictive(posterior_out)

pred_sample = draw_sample(mcmc_sampler, dist_pred, n_samples=500, n_reps=length(data))
obs_dists = map(x->d(x,model.mode), data)
density(pred_sample, alpha=0.1, label=nothing, lw=2)
density!(obs_dists, c=:red, lw=2, label="Observed")

# Outer dimension (length)
pred_sample = draw_sample(mcmc_sampler, outer_dim_pred, n_samples=500, n_reps=length(data))
obs_outer_dim = length.(data)
density(pred_sample, alpha=0.1, label=nothing, lw=2)
density!(obs_outer_dim, c=:red, lw=2, label="Observed")

# Mean inner length
inner_pred = MeanInnerDimensionPredictive(posterior_out)
pred_sample = draw_sample(mcmc_sampler, inner_pred, n_samples=10, n_reps=50)

# Missing entry predictive
data_test = draw_sample(
    mcmc_sampler,
    model,
    desired_samples=1,
    lag=1,
    burn_in=10000
    )
S_test = data_test[1]
Sᵐ
posterior_pred = pred_missing(S_test, (1,1), posterior_out)
get_prediction(posterior_pred)

model_pred = pred_missing(S_test, (1,1), model)
get_prediction(model_pred)

# Checking flip proposal
x = [1,3,4,5,6,7]
P, vmap, vmap_inv = get_informed_proposal_matrix(posterior, 0.0)
flip_informed_excl!(x,[1],P)
x
out = Int[]
for i in 1:1000
    flip_informed_excl!(x,[1],P)
    push!(out, x[1])
end

counts(out)


text:
#function to be called when calculating wigner 3j-symbols
#17/12/2015

function tfunction3j(tt,aa,bb,cc,alpha,beta)
return xoft3j = factorial(tt)*factorial(cc-bb+tt+alpha)*factorial(cc-aa+tt-beta)*factorial(aa+bb-cc-tt)*factorial(aa-tt-alpha)*factorial(bb-tt+beta)
end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
o1 = UncertainValue(Normal, 0, 0.5)
o2 = UncertainValue(Normal, 2, 0.3)
o3 = UncertainValue(Uniform, 0, 4)
o4 = UncertainValue(Uniform, rand(100))

#####################
# UncertainDataset
####################

# Iteration
D = UncertainDataset([o1, o2, o3])
@test length(D) == 3
@test length([x for x in D]) == 3

# Indexing
@test D[1] isa AbstractUncertainValue
@test D[end] isa AbstractUncertainValue
@test D[1:end] isa AbstractVector{<:AbstractUncertainValue}
@test D[[1, 2]] isa AbstractVector{<:AbstractUncertainValue}
@test D[:] isa AbstractVector{<:AbstractUncertainValue}

########################
# UncertainIndexDataset
########################

# Construction
UV = UncertainIndexDataset(D.values)
@test UV isa UncertainIndexDataset

# Iteration
@test length(UV) == 3
@test length([x for x in UV]) == 3

# Indexing
@test UV[1] isa AbstractUncertainValue
@test UV[end] isa AbstractUncertainValue
@test UV[1:end] isa AbstractVector{<:AbstractUncertainValue}
@test UV[[1, 2]] isa AbstractVector{<:AbstractUncertainValue}
@test UV[:] isa AbstractVector{<:AbstractUncertainValue}

########################
# UncertainValueDataset
########################

# Construction
UV = UncertainValueDataset(D.values)
@test UV isa UncertainValueDataset

# Iteration
@test length(UV) == 3
@test length([x for x in UV]) == 3

# Indexing
@test UV[1] isa AbstractUncertainValue
@test UV[end] isa AbstractUncertainValue
@test UV[1:end] isa AbstractVector{<:AbstractUncertainValue}
@test UV[[1, 2]] isa AbstractVector{<:AbstractUncertainValue}
@test UV[:] isa AbstractVector{<:AbstractUncertainValue}


text:
AdjMat = Matrix{Bool}
NeighborList = Vector{Vector{Int}}

function n_nodes(graph::NeighborList)
    return length(graph)
end

function n_nodes(graph::AdjMat)
    return size(graph,1)
end

function get_neighbors(graph::NeighborList, idx::Int)::Set{Int}
    # assume graph is represented in list format
    neighb = Set(graph[idx])
    return push!(neighb,idx)
end

function get_neighbors(graph::AdjMat, idx::Int)::Set{Int}
    # assume graph is represented in Adjacent matrix format
    neighb = Set([i for (i, islink) in enumerate(graph[idx,:]) if islink])
    return push!(neighb,idx)
end

function get_connected(graph, idx)::Set{Int}
    connected = Set([idx])
    to_visit = Stack{Int}()
    push!(to_visit,idx)
    while ~ isempty(to_visit)
        i = pop!(to_visit)
        push!(connected, i)
        neighb = get_neighbors(graph, i)
        for k in neighb
            if ! (k in connected)
            push!(to_visit, k)
            end
        end
    end
    return connected
end

function connected_components(graph)
    components = []
    visited = Set{Int}([])
    for i in 1:n_nodes(graph)
        if (i in visited)
            continue
        end
        comp = get_connected(graph, i)
        push!(components, comp)
        union!(visited, comp)
    end
    return components
end


text:
convert_point(p::Vector, m::Matrix) = (([p, 1]*m)/m[3,3])[1:2]

M(s) = Diagonal([1, 1, 1/s]);

V(Θ) = [cos(Θ)  -sin(Θ) 0;
        sin(Θ)  cos(Θ)  0;
        0       0       1]

P(Δx, Δy) =[1   0   Δx;
            0   1   Δy;
            0   0   1]

size = 300;
models = Set{Matrix}();

function create_image(mvp::Matrix)
    mvp[3, 3] > size && return
    for u in [M(1/2)*V(π/6)*P(0.2, 0.2)]
        create_image(mvp*u)
    end
    push!(models, mvp)
end

function draw()
    for m in models

    end
end

max_stars_repo_name: sclereid/collection

max_issues_repo_name: sclereid/collection

max_forks_repo_name: sclereid/collection