text:
## traits.jl : phylolm for within-species variation, continuous traits

@testset "phylolm: within-species variation, star" begin

#= simulation of data used below (test more reproducible when using fixed data)
# simulate traits at the species-level, then
# repeat identical trait values across m individuals
function simTraits(net, m, paramsprocess)
    sim = simulate(net, paramsprocess)
    trait = sim[:Tips] # simple vector now
    return repeat(trait, inner=m)
end

n = 4; m = 3
starnet = readTopology(PhyloNetworks.startree_newick(n))
Random.seed!(6591)
trait1 = simTraits(starnet, m, ParamsBM(2, 0.5)) # var: 0.5
trait2 = simTraits(starnet, m, ParamsBM(-2, 1))

phylo_noise_var = 2; meas_noise_var = 1
phylo_noise = simTraits(starnet, m, ParamsBM(0, phylo_noise_var))
meas_noise = randn(n*m) * sqrt(meas_noise_var)
trait3 = 10 .+ 2*trait1 + phylo_noise + meas_noise
print(round.(trait1, digits=4)) # etc
labels = repeat(starnet.names, inner=m)
df = DataFrame(trait1=trait1, trait2=trait2, trait3=trait3,tipNames=labels)
=#

n = 4; m = 3
starnet = readTopology(PhyloNetworks.startree_newick(n))
netnames = ["t1","t2","t3","t4"] # = tipLabels(starnet) # was generated to have length n
df = DataFrame(
    species = repeat(netnames, inner=m),
    trait1 = [2.8564, missing, 2.8564, 2.8457, 2.8457, 2.8457, 0.4197, 0.4197, 0.4197, 2.2359, 2.2359, 2.2359],
    trait2 = [-2.0935, -2.0935, -2.0935, 0.4955, 0.4955, 0.4955, -2.1977, -2.1977, -2.1977, -2.618, -2.618, -2.618],
    trait3 = [missing, 15.7869, 14.0615, 14.3547, 13.3721, 15.7062, 9.2764, 8.8553, 8.7627, 15.0298, 15.8258, 15.3248]
    # missing replaced 14.0467
)
Y = df[!,:trait3] # nm vector
X = fill(1.0, (n,2)); X[:,2] = df[1:m:n*m,:trait1] # nx2 matrix
#= reduced data: one row per species, extra columns for SD and n of trait 3
gdf = groupby(df, :species)
df_r = combine(gdf, :trait1 => (x -> mean(skipmissing(x))) => :trait1,
                    :trait2 => (x -> mean(skipmissing(x))) => :trait2,
                    :trait3 => (x -> mean(skipmissing(x))) => :trait3,
                    :trait3 => (x ->  std(skipmissing(x))) => :trait3_sd,
                    :trait3 => (x -> sum(.!ismissing.(x))) => :trait3_n)
=#
df_r = DataFrame(
  species = ["t1","t2","t3","t4"], trait1 = [2.8564,2.8457,.4197,2.2359],
  trait2 = [-2.0935,0.4955,-2.1977,-2.618],
  trait3 = [14.0615, 14.477666666666666, 8.9648, 15.393466666666667],
  #        [14.9242,... if we include ind 2
  trait3_sd = [0,1.1718985891848046,.2737966581242371,.4024181076111425],
  # [1.2200420402592675,... if we include ind 2
  trait3_n = [1,3,3,3]) # [2,3,3,3]

#= star tree: R code to check
df = data.frame(species = rep(c("t1","t2","t3","t4"), each=3),
  trait1 = c(2.8564,NA,2.8564,2.8457,2.8457,2.8457,.4197,.4197,.4197,2.2359,2.2359,2.2359),
  trait2 = c(-2.0935,-2.0935,-2.0935,.4955,.4955,.4955,-2.1977,-2.1977,-2.1977, -2.618,-2.618,-2.618),
  trait3 = c(NA,15.7869,14.0615,14.3547,13.3721,15.7062,9.2764,8.8553,8.7627,15.0298,15.8258,15.3248))
mR = lmer(trait3 ~ trait1 + (1|species), df)
fixef(mR) # 8.457020 2.296978
print(mR, digits=7, ranef.comp="Var") # variances: 2.1216068,0.5332865
logLik(mR) # -13.37294
logLik(mR, REML=FALSE) # -14.57265
vcov(mR) # matrix(c(3.111307,-1.219935,-1.219935,0.5913808), nrow=2)
mR = lmer(trait3 ~ trait1 + (1|species), df, REML=FALSE)
fixef(mR) # 8.439909 2.318488
print(mR, digits=7, ranef.comp="Var") # variances: 0.9470427,0.5299540
logLik(mR) # -14.30235
vcov(mR) # matrix(c(1.5237486,-0.6002803,-0.6002803,0.2941625), nrow=2)
=#

#= Alternatively: Julia code to check
using MixedModels
mm1 = fit(MixedModel, @formula(trait3 ~ trait1 + (1|species)), df, REML=true) # compare with m1
mm3 = fit(MixedModel, @formula(trait3 ~ trait1 + (1|species)), df) # compare with m3
fixef(mm1) # fixed-effect param estimates: 8.45702,2.29698
VarCorr(mm1) # estimated variance-components: species-variance=2.121607, residual-variance=0.533287
# `objective(m)` returns -2 * (log-likelihood of model m)
objective(mm1)/(-2) # -13.3729434, `loglikelihood` not available for models fit by REML
loglikelihood(mm3) # -14.3023458
vcov(mm1) # var-cov matrix for fixef coeffs: [3.11131 -1.21993; -1.21993 0.591381]
=#
m1 = phylolm(@formula(trait3 ~ trait1), df, starnet; reml=true,
      tipnames=:species, withinspecies_var=true)
m2 = phylolm(@formula(trait3 ~ trait1), df_r, starnet; reml=true,
      tipnames=:species, withinspecies_var=true, y_mean_std=true)
@test m1.model.reml && m2.model.reml
@test coef(m1) ≈ [8.457020,2.296978] rtol=1e-5 # fixef(mR)
@test coef(m2) ≈ [8.457020,2.296978] rtol=1e-5
@test sigma2_phylo(m1) ≈ 2.1216068 rtol=1e-5 # print(mR, digits=7, ranef.comp="Var")
@test sigma2_phylo(m2) ≈ 2.1216068 rtol=1e-5
@test sigma2_within(m1) ≈ 0.5332865 rtol=1e-5
@test sigma2_within(m2) ≈ 0.5332865 rtol=1e-5
@test loglikelihood(m1) ≈ -13.37294 rtol=1e-5
@test loglikelihood(m2) ≈ -13.37294 rtol=1e-5
@test vcov(m1) ≈ [3.111307 -1.219935; -1.219935 0.5913808] rtol=1e-5
@test vcov(m2) ≈ [3.111307 -1.219935; -1.219935 0.5913808] rtol=1e-5
m3 = phylolm(@formula(trait3 ~ trait1), df_r, starnet; reml=false,
      tipnames=:species, withinspecies_var=true, y_mean_std=true)
@test !m3.model.reml
@test coef(m3) ≈ [8.439909,2.318488] rtol=1e-5
@test sigma2_phylo(m3) ≈ 0.9470427 rtol=1e-5
@test sigma2_within(m3) ≈ 0.5299540 rtol=1e-5
@test loglikelihood(m3) ≈ -14.30235 rtol=1e-5
@test vcov(m3) ≈ [1.5237486 -0.6002803; -0.6002803 0.2941625] rtol=1e-5

end

@testset "phylolm: binary tree, no withinspecies var" begin

#= Rcode to generate the newick string and dataset:
library(phytools)
set.seed(1)

tree <- pbtree(n=5,scale=1) # topology and branch lengths
tree.newick <- write.tree(tree) # newick format
# tree.newick <- "((t1:0.8121974445,(t2:0.4806586387,t3:0.4806586387):0.3315388057):0.1878025555,(t4:0.1206907041,t5:0.1206907041):0.8793092959);"

X <- fastBM(tree,sig2=2,nsim=2); colnames(X) <- c("x1","x2") # non-intercept predictors
# X <- matrix(c(-0.5503706,-0.7503593,-0.7599224,-0.4217977,-0.9839226,-1.3102488,0.4544195,0.2151272,1.0384023,1.3138847),nrow=5)
# colnames(X) <- c("x1","x2"); rownames(X) <- c("t1","t2","t3","t4","t5")

bsperr <- fastBM(tree) # phylogenetic variation, true BM variance-rate is 1
y <- cbind(rep(1,5),X)%*%(1:3)+bsperr # response, true beta vector is c(1,2,3)
# y <- matrix(c(-4.1836578,0.6258014,-0.3070109,2.6120026,2.2391044),nrow=5)
# rownames(y) <- c("t1","t2","t3","t4","t5")

df <- data.frame(y,X,species=tree$tip.label)
=#
net = readTopology("((t1:0.8121974445,(t2:0.4806586387,t3:0.4806586387):0.3315388057):0.1878025555,(t4:0.1206907041,t5:0.1206907041):0.8793092959);")
df = DataFrame(
      y = [-4.18366,0.625801,-0.307011,2.612,2.2391],
      x1 = [-0.550371,-0.750359,-0.759922,-0.421798,-0.983923],
      x2 = [-1.31025,0.45442,0.215127,1.0384,1.31388],
      species = net.names # ["t1","t2","t3","t4","t5"]
)

#= Rcode to check model fit:
library(nlme)
m1 <- gls(y~x1+x2,data=df,
        correlation=corBrownian(1,tree,form=~species),method="REML")
coef(m1) # coefficients estimates: c(0.6469652,2.0420889,2.8285257)
sigma(m1)^2 # reml BM variance-rate estimate: 0.0438973
logLik(m1) # restricted log-likelihood: -0.07529961
vcov(m1) # cov mat for coef estimates
# matrix(c(0.030948901,0.024782246,0.003569513,0.024782246,0.039001092,0.009137043,0.003569513,0.009137043,0.013366014),nrow=3)
summary(m1)$tTable[,"t-value"] # t-values for coef estimates: c(3.677548,10.340374,24.465786)
summary(m1)$tTable[,"p-value"] # p-values for coef estimates: c(0.066635016,0.009223303,0.001666460)
coef(m1) + qt(p=0.025,df=5-3)*sqrt(diag(vcov(m1))) # lower limit 95%-CI: c(-0.1099703,1.1923712,2.3310897)
coef(m1) + qt(p=0.975,df=5-3)*sqrt(diag(vcov(m1))) # upper limit 95%-CI: c(1.403901,2.891807,3.325962)

m2 <- gls(y~x1+x2,data=df,
        correlation=corBrownian(1,tree,form=~species),method="ML")
coef(m2) # coefficients estimates: c(0.6469652,2.0420889,2.8285257)
sigma(m2)^2 # ml BM variance-rate estimate: 0.01755892
logLik(m2) # log-likelihood: 3.933531
vcov(m2) # cov mat for coef estimates, this evaluates to the same value as vcov(m1)
# Note: vcov(gls(...,method="REML)) == vcov(gls(...,method="ML"))
# The same holds for t/p-values and CIs for the coef estimates
summary(m2)$tTable[,"t-value"] # t-values for coef estimates: c(3.677548,10.340374,24.465786)
summary(m2)$tTable[,"p-value"] # p-values for coef estimates: c(0.066635016,0.009223303,0.001666460)
coef(m2) + qt(p=0.025,df=5-3)*sqrt(diag(vcov(m2))) # lower limit 95%-CI: c(-0.1099703,1.1923712,2.3310897)
coef(m2) + qt(p=0.975,df=5-3)*sqrt(diag(vcov(m2))) # upper limit 95%-CI: c(1.403901,2.891807,3.325962)
=#
m1 = phylolm(@formula(y~x1+x2),df,net; tipnames=:species, reml=true)
m2 = phylolm(@formula(y~x1+x2),df,net; tipnames=:species, reml=false)
@test m1.model.reml
@test coef(m1) ≈ [0.6469652,2.0420889,2.8285257] rtol=1e-5
@test sigma2_phylo(m1) ≈ 0.0438973 rtol=1e-4
@test isnothing(sigma2_within(m1))
@test loglikelihood(m1) ≈ -0.07529961 rtol=1e-3
@test vcov(m1) ≈ [0.030948901 0.024782246 0.003569513;0.024782246 0.039001092 0.009137043;0.003569513 0.009137043 0.013366014] rtol=1e-4
@test coeftable(m1).cols[coeftable(m1).teststatcol] ≈ [3.677548,10.340374,24.465786] rtol=1e-4
@test coeftable(m1).cols[coeftable(m1).pvalcol] ≈ [0.066635016,0.009223303,0.001666460] rtol=1e-4
@test coeftable(m1).cols[findall(coeftable(m1).colnms .== "Lower 95%")[1]] ≈ [-0.1099703,1.1923712,2.3310897] rtol=1e-5
@test coeftable(m1).cols[findall(coeftable(m1).colnms .== "Upper 95%")[1]] ≈ [1.403901,2.891807,3.325962] rtol=1e-5
@test !m2.model.reml
@test coef(m2) ≈ [0.6469652,2.0420889,2.8285257] rtol=1e-5
@test sigma2_phylo(m2) ≈ 0.01755892 rtol=1e-4
@test isnothing(sigma2_within(m2))
@test loglikelihood(m2) ≈ 3.933531 rtol=1e-4
@test vcov(m2) ≈ [0.030948901 0.024782246 0.003569513;0.024782246 0.039001092 0.009137043;0.003569513 0.009137043 0.013366014] rtol=1e-4
@test coeftable(m2).cols[coeftable(m2).teststatcol] ≈ [3.677548,10.340374,24.465786] rtol=1e-4
@test coeftable(m2).cols[coeftable(m2).pvalcol] ≈ [0.066635016,0.009223303,0.001666460] rtol=1e-4
@test coeftable(m2).cols[findall(coeftable(m2).colnms .== "Lower 95%")[1]] ≈ [-0.1099703,1.1923712,2.3310897] rtol=1e-5
@test coeftable(m2).cols[findall(coeftable(m2).colnms .== "Upper 95%")[1]] ≈ [1.403901,2.891807,3.325962] rtol=1e-5

end

@testset "phylolm: binary tree, withinspecies var" begin

#= NOTE:
This testset DOES NOT check for similarity of the variance-components estimates
from 'pgls.SEy' (phytools) and 'phylolm' (PhyloNetworks). In fact,
'pgls.SEy' and 'phylolm' solve different problems!
ML/REML variance-component estimates are obtained by running 'phylolm'.
The code for 'pgls.SEy' is adapted to calculate coef estimates, cov matrix for
coef estimates, t/p-values for coef estimates, and full/restricted-loglikelihood,
given these variance-components estimates. These values are then compared against
those extracted from the output of 'phylolm'.
=#

#= Rcode to generate the newick string and dataset:
library(phytools)
set.seed(2)
m <- 5 # sample-size per taxa
n <- 5 # no. of taxa
p <- 3 # no. of predictors (including intercept)

tree <- pbtree(n=n,scale=1) # topology and branch lengths
tree.newick <- write.tree(tree) # newick format
# tree.newick <- "(((t4:0.2537636499,t5:0.2537636499):0.2103870459,t3:0.4641506959):0.5358493041,(t1:0.4807642475,t2:0.4807642475):0.5192357525);"

X <- fastBM(tree,sig2=2,nsim=2); colnames(X) <- c("x1","x2") # non-intercept predictors
# X <- matrix(c(0.7894387,0.3285286,2.1495131,-1.4935665,-2.3199935,2.8184268,0.4740687,2.5801004,1.9967379,-0.4121874),nrow=5)
# colnames(X) <- c("x1","x2"); rownames(X) <- c("t4","t5","t3","t1","t2")

bsperr <- fastBM(tree) # phylogenetic variation, true BM variance-rate is 1
wspvar <- 0.1 # within-species variance is 0.1
msrerr <- rnorm(n=n,sd=sqrt(wspvar/m)) # msr error in species-level mean responses
y_sd <- sqrt(wspvar*rchisq(n=n,df=m-1)/(m-1)) # sd in individual responses per species
# y_sd <- c(0.2463003,0.3236629,0.2458547,0.4866844,0.3434582);
y <- cbind(rep(1,5),X)%*%(1:3)+bsperr+msrerr # response, true beta vector is c(1,2,3)
# y <- matrix(c(11.399108,3.216645,13.648011,4.851454,-4.922803),nrow=5); rownames(y) <- c("t4","t5","t3","t1","t2")

df <- data.frame(y,y_sd,X,species=tree$tip.label)
=#
net = readTopology("(((t4:0.2537636499,t5:0.2537636499):0.2103870459,t3:0.4641506959):0.5358493041,(t1:0.4807642475,t2:0.4807642475):0.5192357525);")
df = DataFrame(
      y=[11.399108,3.216645,13.648011,4.851454,-4.922803],
      y_sd=[0.2463003,0.3236629,0.2458547,0.4866844,0.3434582],
      y_n=fill(5,5),
      x1=[0.7894387,0.3285286,2.1495131,-1.4935665,-2.3199935],
      x2=[2.8184268,0.4740687,2.5801004,1.9967379,-0.4121874],
      species=net.names
)

#= R/Julia code to check model fit:
(1) Julia code:
## To extract reml/ml msrerr variance estimates fitted by phylolm
m1 |> sigma2_within |> x -> round(x,sigdigits=6) # reml est msrerr var: 0.116721
m2 |> sigma2_within |> x -> round(x,sigdigits=6) # ml est msrerr var: 0.125628
m1 |> sigma2_phylo |> x -> round(x,sigdigits=6) # reml est BM var: 0.120746
m2 |> sigma2_phylo |> x -> round(x,sigdigits=6) # ml est BM var: 0.000399783

(2) R code:
## To check phylolm REML fit
library(nlme)
T1 <- tree
wspvar1 <- 0.116721 # reml est msrerr var of indiv-lvl rsps
bspvar1 <- 0.120746 # reml est BM var
se1 <- setNames(rep(sqrt(wspvar1/m),n),T1$tip.label) # msrerr sd of species-lvl mean rsp
T1$edge.length <- T1$edge.length*bspvar1 # scale all edges by est BM var
ii <- sapply(1:Ntip(T1),function(x,e) which(e==x),e=T1$edge[,2]) # indices of pendant edges
# extend pendant edges by msrerr sd
T1$edge.length[ii] <- T1$edge.length[ii]+se1[T1$tip.label]^2
covmat1 <- vcv(T1) # extract est vars of species-lvl mean rsp
m1 <- gls(y~x1+x2,data=cbind(df,vf=diag(covmat1)),
          correlation=corBrownian(1,T1,form=~species),
          method="REML",
          weights=varFixed(~vf))

coef(m1) # reml coefficient estimates: c(1.079839,1.976719,3.217391)

Xp <- model.matrix(m1,df) # predictor matrix
vcov1 <- solve(t(Xp)%*%solve(covmat1)%*%Xp) # cov mat for coef estimates
# matrix(c(0.09386431,0.02273458,-0.02602937,0.02273458,0.02172123,-0.01133032,-0.02602937,-0.01133032,0.01584198),nrow=3)

teststat1 <- coef(m1)/sqrt(diag(vcov1)) # test stat for coef estimates: c(3.524591,13.412285,25.562252)
pval1 <- sapply(X=pt(teststat1,n-p),FUN=function(p) 2*min(p,1-p)) # pval for coef est:
# c(0.071921276,0.005513044,0.001526885)
lowerci1 <- coef(m1) + qt(p=0.025,df=n-p)*sqrt(diag(vcov1)) # lower limit 95%-CI: c(-0.2383769,1.3425890,2.6758379)
upperci1 <- coef(m1) + qt(p=0.975,df=n-p)*sqrt(diag(vcov1)) # upper limit 95%-CI: c(2.398055,2.610850,3.758944)

RSS <- sum((m-1)*(df$y_sd^2)) # residual sum-of-squares wrt to the species means
logLik(m1) # species-lvl cond restricted-ll: -3.26788
sigm1 <- sigma(m1) # this is not bspvar1!, but rather the best "scaling" for covmat1
# indiv-lvl cond restricted ll
# "+ (n-p)*(2*log(sigm1)-sigm1^2+1)/2" un-scales the species-lvl rll returned by logLik
# "- n*(m-1)*(log(wspvar1)+log(2*pi))/2 - (n*log(m)+RSS/wspvar1)/2" corrects to indiv-lvl rll
rll.species <- logLik(m1)
rll.indiv <- (rll.species
              + (n-p)*(2*log(sigm1)-sigm1^2+1)/2
              - n*(m-1)*(log(wspvar1)+log(2*pi))/2 - (n*log(m)+RSS/wspvar1)/2
              )
round(rll.indiv,digits=6) # indiv-lvl rll of remles: -14.14184

## To check phylolm ML fit
T2 <- tree
wspvar2 <- 0.125628 # ml est msrerr var of indiv-lvl rsps
bspvar2 <- 0.000399783 # ml est BM var
se2 <- setNames(rep(sqrt(wspvar2/m),n),T2$tip.label)
T2$edge.length <- T2$edge.length*bspvar2
ii <- sapply(1:Ntip(T2),function(x,e) which(e==x),e=T2$edge[,2])
T2$edge.length[ii] <- T2$edge.length[ii]+se2[T2$tip.label]^2
covmat2 <- vcv(T2)
m2 <- gls(y~x1+x2,data=cbind(df,vf=diag(covmat2)),
          correlation=corBrownian(1,T2,form=~species),
          method="ML",
          weights=varFixed(~vf))

coef(m2) # ml coefficient estimates: c(0.9767352,1.9155142,3.2661862)

Xp <- model.matrix(m2,df) # predictor matrix
vcov2 <- solve(t(Xp)%*%solve(covmat2)%*%Xp) # cov mat for coef estimates
# matrix(c(0.019118171,0.004924625,-0.008977541,0.004924625,0.003577654,-0.003028413,-0.008977541,-0.003028413,0.005793475),nrow=3)

teststat2 <- coef(m2)/sqrt(diag(vcov2)) # test stat for coef est: c(7.064049,32.024784,42.911270)
pval2 <- sapply(X=pt(teststat2,n-p),FUN=function(p) 2*min(p,1-p)) # pval for coef est:
# c(0.0194568147,0.0009736278,0.0005426298)
lowerci2 <- coef(m2) + qt(p=0.025,df=n-p)*sqrt(diag(vcov2)) # lower limit 95%-CI: c(0.381814,1.658158,2.938690)
upperci2 <- coef(m2) + qt(p=0.975,df=n-p)*sqrt(diag(vcov2)) # upper limit 95%-CI: c(1.571656,2.172871,3.593682)

ll.species <- logLik(m2) # species-lvl cond ll: 1.415653
sigm2 <- sigma(m2) # this is not bspvar2!, but rather the best "scaling" for covmat2
RSS <- sum((m-1)*(df$y_sd^2)) # residual sum-of-squares wrt to the species means
# indiv-lvl cond ll
# "+ n*(2*log(sigm2)-sigm2^2+1)/2" un-scales the species-lvl ll returned by logLik
# "- n*(m-1)*(log(wspvar2)+log(2*pi))/2 - (n*log(m)+RSS/wspvar2)/2" corrects to indiv-lvl ll
ll.indiv <- (ll.species
             + n*(2*log(sigm2)-sigm2^2+1)/2
             - n*(m-1)*(log(wspvar2)+log(2*pi))/2 - (n*log(m)+RSS/wspvar2)/2
             )
round(ll.indiv,digits=6) # indiv-lvl ll of mles: -9.582357
=#
m1 = phylolm(@formula(y~x1+x2),df,net;
      tipnames=:species, reml=true, withinspecies_var=true, y_mean_std=true)
m2 = phylolm(@formula(y~x1+x2),df,net;
      tipnames=:species, reml=false, withinspecies_var=true, y_mean_std=true)

@test coef(m1) ≈ [1.079839,1.976719,3.217391] rtol=1e-4
@test coef(m2) ≈ [0.9767352,1.9155142,3.2661862] rtol=1e-4
@test loglikelihood(m1) ≈ -14.14184 rtol=1e-4
@test loglikelihood(m2) ≈ -9.582357 rtol=1e-4
@test vcov(m1) ≈ [0.09386431 0.02273458 -0.02602937;0.02273458 0.02172123 -0.01133032;-0.02602937 -0.01133032 0.01584198] rtol=1e-4
@test vcov(m2) ≈ [0.019118171 0.004924625 -0.008977541;0.004924625 0.003577654 -0.003028413;-0.008977541 -0.003028413 0.005793475] rtol=1e-4
@test coeftable(m1).cols[coeftable(m1).teststatcol] ≈ [3.524591,13.412285,25.562252] rtol=1e-4
@test coeftable(m2).cols[coeftable(m2).teststatcol] ≈ [7.064049,32.024784,42.911270] rtol=1e-4
@test coeftable(m1).cols[coeftable(m1).pvalcol] ≈ [0.071921276,0.005513044,0.001526885] rtol=1e-4
@test coeftable(m2).cols[coeftable(m2).pvalcol] ≈ [0.0194568147,0.0009736278,0.0005426298] rtol=1e-4
@test coeftable(m1).cols[findall(coeftable(m1).colnms .== "Lower 95%")[1]] ≈ [-0.2383769,1.3425890,2.6758379] rtol=1e-4
@test coeftable(m2).cols[findall(coeftable(m2).colnms .== "Lower 95%")[1]] ≈ [0.381814,1.658158,2.938690] rtol=1e-4
@test coeftable(m1).cols[findall(coeftable(m1).colnms .== "Upper 95%")[1]] ≈ [2.398055,2.610850,3.758944] rtol=1e-4
@test coeftable(m2).cols[findall(coeftable(m2).colnms .== "Upper 95%")[1]] ≈ [1.571656,2.172871,3.593682] rtol=1e-4

# missing y (but y_n/SD present); missing/infinite sample size; missing/infinite y_SD
df = df[[1,2,3,4],:] # delete t2
allowmissing!(df, [:y,:y_n,:y_sd])
# m0 = phylolm(@formula(y~1),df,net;  tipnames=:species, withinspecies_var=true, y_mean_std=true)
df1 = DataFrame(df) # makes a copy
push!(df1, (missing,0.,0,0.,0.,"t2"))
m1 = phylolm(@formula(y~1),df1,net; tipnames=:species, withinspecies_var=true, y_mean_std=true)
@test coef(m1) ≈ [7.8116494901242195] atol=1e-6 # same as m0
@test sigma2_within(m1) ≈ 0.11568962074173368 atol=1e-6
@test nobs(m1) == 20        # 4 species x 5 indiv/species
@test dof_residual(m1) == 3
# sample sizes y_n: 1, 0, infinite
@test dof(m1) == 3 # intercept, s2phylo, s2within
df1[5,:] .= (0.,0.,1,0.,0.,"t2") # then almost same within-sp variation
m1 = phylolm(@formula(y~1),df1,net; tipnames=:species, withinspecies_var=true, y_mean_std=true)
@test sigma2_within(m1) ≈ 0.11569 atol=1e-4
@test nobs(m1) == 21
@test dof_residual(m1) == 4
df1[5,:] .= (0.,0.,0,0.,0.,"t2") # some species have 0 or <0 num_individuals, column y_n
@test_throws ErrorException phylolm(@formula(y~1),df1,net; tipnames=:species, withinspecies_var=true, y_mean_std=true)
df1[!,:y_n] = Float64.(df1[!,:y_n])
df1[5,:] .= (0.,0.,Inf,0.,0.,"t2")#some species have infinite num_individuals, column y_n
@test_throws ErrorException phylolm(@formula(y~1),df1,net; tipnames=:species, withinspecies_var=true, y_mean_std=true)
df1[5,:] .= (0.,missing,1,0.,0.,"t2") # some SD values are missing, column y_sd
@test_throws ErrorException phylolm(@formula(y~1),df1,net; tipnames=:species, withinspecies_var=true, y_mean_std=true)
df1[5,:] .= (0.0,Inf,1,0.,0.,"t2")   # some SD values are infinite, column y_sd
@test_throws ErrorException phylolm(@formula(y~1),df1,net; tipnames=:species, withinspecies_var=true, y_mean_std=true)

end

@testset "phylolm: within-species var, network h=1" begin

#= Simulation of data used below:
(1) Individual-level data
function simTraits(net, m, paramsprocess)
      sim = simulate(net, paramsprocess)
      trait = sim[:Tips]
      return repeat(trait, inner=m)
end
n = 6; m = 3 # 6 species, 3 individuals per species
net = readTopology("((((D:0.4,C:0.4):4.8,((A:0.8,B:0.8):2.2)#H1:2.2::0.7):4.0,(#H1:0::0.3,E:3.0):6.2):2.0,O:11.2);");
Random.seed!(18480224);
trait1 = simTraits(net, m, ParamsBM(2, 0.5)) # simulate a BM with mean 2 and variance 0.5 on net
trait2 = simTraits(net, m, ParamsBM(-2, 1.0)) # simulate a BM with mean -2 and variance 1.0 on net
phylo_noise_var = 0.1 # BM variance-rate
meas_noise_var = 0.01 # individual-level msrerr variance
phylo_noise = simTraits(net, m, ParamsBM(0, phylo_noise_var))
meas_noise = randn(n*m)*sqrt(meas_noise_var)
trait3 = 10 .+ 2 * trait1 + phylo_noise + meas_noise
labels = repeat(names(vcv(net)), inner=m)
df = DataFrame(trait1=trait1, trait2=trait2, trait3=trait3, tipNames=labels)

(2) Species-level data (extra columns for SD and n of trait3)
gdf = groupby(df, :species)
df_r = combine(gdf, :trait1 => (x -> mean(x)) => :trait1,
                    :trait2 => (x -> mean(x)) => :trait2,
                    :trait3 => (x -> mean(x)) => :trait3,
                    :trait3 => (x -> std(x)) => :trait3_sd,
                    :trait3 => (x -> length(x)) => :trait3_n)
=#
n = 6; m = 3
net = readTopology("((((D:0.4,C:0.4):4.8,((A:0.8,B:0.8):2.2)#H1:2.2::0.7):4.0,(#H1:0::0.3,E:3.0):6.2):2.0,O:11.2);");
df = DataFrame(
      species = repeat(["D","C","A","B","E","O"],inner=m),
      trait1 = [4.08298,4.08298,4.08298,3.10782,3.10782,3.10782,2.17078,2.17078,2.17078,1.87333,1.87333,1.87333,2.8445,
                2.8445,2.8445,5.88204,5.88204,5.88204],
      trait2 = [-7.34186,-7.34186,-7.34186,-7.45085,-7.45085,-7.45085,-3.32538,-3.32538,-3.32538,-4.26472,-4.26472,
                -4.26472,-5.96857,-5.96857,-5.96857,-1.99388,-1.99388,-1.99388],
      trait3 = [18.8101,18.934,18.9438,17.0687,17.0639,17.0732,14.4818,14.1112,14.2817,13.0842,12.9562,12.9019,15.4373,
                15.4075,15.4317,24.2249,24.1449,24.1302]
)
df_r = DataFrame(
      species = ["D","C","A","B","E","O"],
      trait1 = [4.08298,3.10782,2.17078,1.87333,2.8445,5.88204],
      trait2 = [-7.34186,-7.45085,-3.32538,-4.26472,-5.96857,-1.99388],
      trait3 = [18.895966666666666,17.0686,14.291566666666668,12.980766666666666,15.4255,24.166666666666668],
      trait3_sd = [0.07452397824414288,0.004650806381693211,0.18549690922851855,0.09360001780626576,0.01583792915756317,0.05096433393397276],
      trait3_n = [3, 3, 3, 3, 3, 3]
)

m1 = phylolm(@formula(trait3 ~ trait1), df, net; reml=true,
      tipnames=:species, withinspecies_var=true)
m2 = phylolm(@formula(trait3 ~ trait1), df_r, net; reml=true,
      tipnames=:species, withinspecies_var=true, y_mean_std=true)
m3 = phylolm(@formula(trait3 ~ trait1), df, net; reml=false,
      tipnames=:species, withinspecies_var=true)
m4 = phylolm(@formula(trait3 ~ trait1), df_r, net; reml=false,
      tipnames=:species, withinspecies_var=true, y_mean_std=true)

@testset "agreement across data input" begin # nested test set
#= - individual-level data vs when the species-level data: see m1, m2, m3, m4
- permuted data rows: see m1permute, m2permute
- w/o withinspecies var: when a species has missing entries vs when the species
  is completely absent from the data (see m3missingentry, m4missingrow). =#
m1permute = phylolm(@formula(trait3 ~ trait1), df[[1,6,11,17,16,18,8,5,9,3,12,7,13,10,2,14,4,15],:], # permute predictor rows
                  net; reml=true,tipnames=:species, withinspecies_var=true)
m2permute = phylolm(@formula(trait3 ~ trait1), df_r[[5,2,3,1,6,4],:], net; reml=true, # permute predictor rows
                  tipnames=:species, withinspecies_var=true, y_mean_std=true)
@test coef(m1) ≈ [9.65347,2.30357] rtol=1e-4
@test coef(m1permute) ≈ [9.65347,2.30357] rtol=1e-4
@test coef(m2) ≈ [9.65347,2.30357] rtol=1e-4
@test coef(m2permute) ≈ [9.65347,2.30357] rtol=1e-4
@test sigma2_phylo(m1) ≈ 0.156188 rtol=1e-4
@test sigma2_phylo(m1permute) ≈ 0.156188 rtol=1e-4
@test sigma2_phylo(m2) ≈ 0.156188 rtol=1e-4
@test sigma2_phylo(m2permute) ≈ 0.156188 rtol=1e-4
@test sigma2_within(m1) ≈ 0.008634 rtol=1e-4
@test sigma2_within(m1permute) ≈ 0.008634 rtol=1e-4
@test sigma2_within(m2) ≈ 0.008634 rtol=1e-4
@test sigma2_within(m2permute) ≈ 0.008634 rtol=1e-4
@test loglikelihood(m1) ≈ 1.944626 rtol=1e-4
@test loglikelihood(m1permute) ≈ 1.944626 rtol=1e-4
@test loglikelihood(m2) ≈ 1.944626 rtol=1e-4
@test loglikelihood(m2permute) ≈ 1.944626 rtol=1e-4
@test stderror(m1) ≈ [1.3065987324433421,0.27616258597477233] atol=1e-6
@test stderror(m1permute) ≈ [1.3065987324433421,0.27616258597477233] atol=1e-6
@test stderror(m2) ≈ [1.3065987324433421,0.27616258597477233] atol=1e-6
@test stderror(m2permute) ≈ [1.3065987324433421,0.27616258597477233] atol=1e-6
@test nobs(m1) == 18 # 6 species, 18 ind
@test nobs(m1permute) == 18
@test residuals(m1) ≈ [-0.16295769587309603,0.2560302141911026,-0.36246083625543507,-0.9880623366773218,-0.78049231427113,0.9634719640663754] atol=1e-6
# indexin(df[[1,6,11,17,16,18,8,5,9,3,12,7,13,10,2,14,4,15],:].species |> unique, df.species |> unique) == [1,2,4,6,3,5]
@test residuals(m1permute) ≈ [-0.16295769587309603,0.2560302141911026,-0.36246083625543507,
                              -0.9880623366773218,-0.78049231427113,0.9634719640663754][[1,2,4,6,3,5]] atol=1e-6
@test residuals(m2) ≈ [-0.16295769587309603,0.2560302141911026,-0.36246083625543507,-0.9880623366773218,-0.78049231427113,0.9634719640663754] atol=1e-6
@test residuals(m2permute) ≈ [-0.16295769587309603,0.2560302141911026,-0.36246083625543507,
                              -0.9880623366773218,-0.78049231427113,0.9634719640663754][[5,2,3,1,6,4]] atol=1e-6
@test dof_residual(m1) == 4
@test dof_residual(m1permute) == 4
@test dof(m1) == 4

dfmissingentry = allowmissing(df); dfmissingentry.trait1[1:3] .= missing # remove entire trait1 col for species D
m3missingentry = phylolm(@formula(trait3 ~ trait1), dfmissingentry, net; reml=false, # not using data from species D
                        tipnames=:species, withinspecies_var=true)
df_rmissingrow = copy(df_r); df_rmissingrow = df_rmissingrow[2:6,:] # remove entire row for species D
m4missingrow = phylolm(@formula(trait3 ~ trait1), df_rmissingrow, net; reml=false, # not using data from species D
                        tipnames=:species, withinspecies_var=true, y_mean_std=true)
@test coef(m3) ≈ [9.63523,2.30832] rtol=1e-4
@test coef(m4) ≈ [9.63523,2.30832] rtol=1e-4
@test coef(m3missingentry) ≈ coef(m4missingrow) rtol=1e-4
@test sigma2_phylo(m3) ≈ 0.102255 rtol=1e-4
@test sigma2_phylo(m4) ≈ 0.102255 rtol=1e-4
@test sigma2_phylo(m3missingentry) ≈ sigma2_phylo(m4missingrow) rtol=1e-4
@test sigma2_within(m3) ≈ 0.008677 rtol=1e-4
@test sigma2_within(m4) ≈ 0.008677 rtol=1e-4
@test sigma2_within(m3missingentry) ≈ sigma2_within(m4missingrow) rtol=1e-4
@test loglikelihood(m3) ≈ 1.876606 rtol=1e-4
@test loglikelihood(m4) ≈ 1.876606 rtol=1e-4
@test loglikelihood(m3missingentry) ≈ loglikelihood(m4missingrow) rtol=1e-4
@test stderror(m3) ≈ [1.0619360781577734, 0.22496955609230126] atol=1e-6
@test stderror(m4) ≈ [1.0619360781577734, 0.22496955609230126] atol=1e-6
@test stderror(m3missingentry) ≈ stderror(m4missingrow) atol=1e-6
@test residuals(m3missingentry) ≈ residuals(m4missingrow) atol=1e-6
end # agreement test subset

@testset "model comparison & likelihood ratio test" begin
m3null = phylolm(@formula(trait3 ~ 1), df, net; reml=false, tipnames=:species, withinspecies_var=true)
m3full = phylolm(@formula(trait3 ~ trait1 + trait2), df, net; reml=false, tipnames=:species, withinspecies_var=true)
@test StatsModels.isnested(m3null, m3full)
tab = (@test_logs lrtest(m3null, m3, m3full))
@test all(isapprox.(tab.deviance, (13.720216379785523,-3.75321128763398,-4.243488375445074), atol=1e-6))
@test tab.dof == (3, 4, 5)
@test all(isapprox.(tab.pval[2:end], (2.9135155795020905e-5,0.4838037279625203), atol=1e-6))
m1null = phylolm(@formula(trait3 ~ 1), df, net; tipnames=:species, withinspecies_var=true) # REML
@test !(@test_logs (:error, r"fitted with ML") StatsModels.isnested(m1null, m1)) # REML, different FEs
@test_logs (:error, r"same criterion") (@test_throws ArgumentError lrtest(m3null, m1)) # ML and REML
m2w = phylolm(@formula(trait3 ~ trait1), df_r, net; tipnames=:species)
m5w = phylolm(@formula(trait3 ~ trait1), df[[1,2,4,7,13,18],:], net; tipnames=:species, withinspecies_var=true)
@test !(@test_logs (:error, r"same number of obs") StatsModels.isnested(m2, m2w))
@test !(@test_logs (:error, r"same response") StatsModels.isnested(m5w,m2w))
m6w = (@test_logs (:info,r"^Maximum lambda") phylolm(@formula(trait3 ~ trait1), df_r, net; tipnames=:species, model="lambda"))
tab = lrtest(m2w, m6w) # both REML, but same predictors
@test tab.pval[2] ≈ 0.03928341265297505 atol=1e-6
@test !PhyloNetworks.isnested(PhyloNetworks.PagelLambda(0.1),PhyloNetworks.BM())
@test !PhyloNetworks.isnested(PhyloNetworks.PagelLambda(),PhyloNetworks.ScalingHybrid())
@test !PhyloNetworks.isnested(PhyloNetworks.ScalingHybrid(2.0),PhyloNetworks.PagelLambda())
@test_throws ArgumentError ftest(m3null, m3, m3full) # not the same Y after transformation
end # lrt test subset

@testset "equivalence with Pagel's lambda on expanded net" begin
# find a given named tip or node
function getNode(name::String, net::HybridNetwork)
    i = findfirst(n -> n.name == name, net.node)
    isnothing(i) && error("Node $(name) was not found in the network.")
    return net.node[i]
end
# new column with tip ids, initialized to species names
dfbig = deepcopy(df)
insertcols!(dfbig, 1, :speciesIds => dfbig[!,:species])
# add tips with zero length branches on the network
netbig = deepcopy(net)
for i in 1:nrow(dfbig)
    sp_name = dfbig[i, :species]; tip_name = sp_name * string(i)
    PhyloNetworks.addleaf!(netbig, getNode(sp_name, netbig), tip_name, 0.0)
    dfbig[i, :speciesIds] = tip_name
end
# (((((D1:0.0,D2:0.0,D3:0.0)D:0.4,(C4:0.0,C5:0.0,C6:0.0)C:0.4):4.8,(((A7:0.0,A8:0.0,A9:0.0)A:0.8,(B10:0.0,B11:0.0,B12:0.0)B:0.8):2.2)#H1:2.2::0.7):4.0,(#H1:0.0::0.3,(E13:0.0,E14:0.0,E15:0.0)E:3.0):6.2):2.0,(O16:0.0,O17:0.0,O18:0.0)O:11.2);
# fit lambda on bigger network
m1big = (@test_logs (:info, r"^M") phylolm(@formula(trait3 ~ trait1), dfbig, netbig, model="lambda"; tipnames=:speciesIds))
# equivalent phylo and within variances
net_height = maximum(PhyloNetworks.getHeights(net));
sig_phy_2 = sigma2_phylo(m1big) * lambda_estim(m1big);
sig_err_2 = sigma2_phylo(m1big) * (1 - lambda_estim(m1big)) * net_height;
@test sigma2_phylo(m1) ≈ sig_phy_2 atol=1e-7
@test sigma2_within(m1) ≈ sig_err_2 atol=1e-9
# equivalent lambda value; equivalent coefficients and other
w_var_lambda = sigma2_phylo(m1) / (sigma2_within(m1) / net_height + sigma2_phylo(m1));
@test w_var_lambda ≈ lambda_estim(m1big) atol=1e-8
@test coef(m1) ≈ coef(m1big) atol=1e-7
@test mu_phylo(m1) ≈ mu_phylo(m1big) atol=1e-7
@test vcov(m1) ≈ vcov(m1big) atol=1e-6
@test stderror(m1) ≈ stderror(m1big) atol=1e-6
@test nobs(m1) ≈ nobs(m1big) atol=1e-10
@test dof(m1) ≈ dof(m1big) atol=1e-10
@test deviance(m1) ≈ deviance(m1big) atol=1e-10
@test loglikelihood(m1) ≈ loglikelihood(m1big) atol=1e-10
# what differs: residuals, response, predict (not same length)
# also differs: dof_residuals (18-2=16 vs 6-2=4), and so predint
end

allowmissing!(df,  [:trait3]); df[4:6,:trait3] .= missing # species C missing
allowmissing!(df_r,[:trait3]); df_r[2,:trait3] = missing  # to check imputation

@testset "ancestral state prediction, intercept only" begin
m1 = phylolm(@formula(trait3 ~ 1), df_r, net; tipnames=:species, withinspecies_var=true, y_mean_std=true)
ar1 = (@test_logs (:warn, r"^T") ancestralStateReconstruction(m1))
# ar.NodeNumbers[8] == 2 (looking at node #2), m1.model.V.tipNames[indexin([2],m1.model.V.tipNumbers)[1]] == "C" (looking at tip "C")
@test ar1.traits_nodes[8] ≈ 18.74416393519304 rtol=1e-5 # masked sampled C_bar was 17.0686
@test predint(ar1)[8,:] ≈ [15.24005506417728,22.2482728062088] rtol=1e-5
# on dataframe with model passed as keyword args. must be individual data.
ar2 = (@test_logs (:warn, r"^T") ancestralStateReconstruction(df[!,[:species,:trait3]], net; tipnames=:species, withinspecies_var=true))
@test ar2.traits_nodes ≈ ar1.traits_nodes rtol=1e-5
@test predint(ar2) ≈ predint(ar1) rtol=1e-5
# When withinspecies_var=true, predicted values at the tips are part of
# "traits_nodes", not "traits_tips", and differ from the observed sample means.
@test length(ar1.traits_tips) == 0
@test length(ar2.traits_tips) == 0
@test length(ar1.traits_nodes) == 13
@test length(ar2.traits_nodes) == 13
@test ar1.traits_nodes[9] ≈ 18.895327175656757 # for tip D: observed 18.896
@test predint(ar1)[9,:] ≈ [18.73255168713768,19.058102664175834] # narrow
end

@testset "ancestral state prediction, more than intercept" begin
m3 = phylolm(@formula(trait3 ~ trait1 + trait2), df[[1,6,11,17,16,18,8,5,9,3,12,7,13,10,2,14,4,15],:], net; tipnames=:species, withinspecies_var=true)
X_n = [m3.model.X;m3.model.X[1:3,:]] # 8x3 Array
ar3 = (@test_logs (:warn, r"^T") ancestralStateReconstruction(m3, X_n))
m4 = phylolm(@formula(trait3 ~ trait1 + trait2), df_r[[1,4,6,3,2,5],:], net; tipnames=:species, withinspecies_var=true, y_mean_std=true)
ar4 = (@test_logs (:warn, r"^T") ancestralStateReconstruction(m4, X_n))
@test ar3.NodeNumbers == ar4.NodeNumbers
@test ar3.traits_nodes ≈ ar4.traits_nodes rtol=1e-5
@test ar3.variances_nodes ≈ ar4.variances_nodes rtol=1e-5
@test size(expectations(ar4)) == (13, 2)
@test expectationsPlot(ar4)[8,2] == "24.84*"
@test expectationsPlot(ar3)[8,2] == "24.84*"
@test predint(ar3)[13,:] ≈ [15.173280800793783,15.677786825808212] rtol=1e-5 # narrow at tip with data
@test predint(ar3)[7,:] ≈ [10.668220499837211,15.322037065478693] rtol=1e-5  # wide at a root
p3 = predintPlot(ar3)
p4 = predintPlot(ar4)
@test p3[!,:nodeNumber] == p4[!,:nodeNumber]
@test p3[13,2] == "[15.17, 15.68]"
@test p4[13,2] == "[15.17, 15.68]"
@test p3[7,2] == "[10.67, 15.32]"
@test p4[7,2] == "[10.67, 15.32]"
end # test subset

end # test set: withinspecies_var on network h=1


text:
"""
    @sync ex

Run expression `ex` and synchronize the GPU afterwards. This is a CPU-friendly
synchronization, i.e. it performs a blocking synchronization without increasing CPU load. As
such, this operation is preferred over implicit synchronization (e.g. when performing a
memory copy) for high-performance applications.

It is also useful for timing code that executes asynchronously.
"""
macro sync(ex)
    quote
        local e = CuEvent(EVENT_BLOCKING_SYNC | EVENT_DISABLE_TIMING)
        local ret = $(esc(ex))
        record(e)
        synchronize(e)
        ret
    end
end

ext: jl


text:
module SqState
    using QuantumStateBase
    using QuantumStateDistributions
    using LinearAlgebra

    using Fetch, DataDeps, JLD2

    using CUDA
    using Flux
    using Flux.Data: DataLoader
    using Zygote, ChainRulesCore

    using NeuralOperators

    function __init__()
        register(DataDep(
            "SqState",
            """Data and models for SqState.""",
            "https://drive.google.com/file/d/1UzaPBpTuhxvmyUWnoOupEr3cRZUMc-0-/view?usp=sharing",
            "4f426e53a018ae0692fb656b5fa63712157f0889bb3b2cf7de2a169764f9a00e",
            fetch_method=gdownload,
            post_fetch_method=unpack
        ))
    end

    include("utils.jl")
    include("gen_data.jl")
    include("preprocess.jl")
    include("model.jl")
end


text:
grad_abs(x::T) where {T<:Real} = sign(x)
grad_abs(z::T) where {T<:Complex} = exp(-im * angle(z))

# NAME, forward, backward (for real)
__NON_HOLOMOPHICS__ = [
    (:RealOp, x->:(real.($x))),
    (:ImagOp, x->:(imag.($x))),
    (:Conj, x->:(conj.($x))),
    (:Abs, x->:(abs.($x))),
    (:Abs2, x->:(abs.($x))),
    (:Arg, x->:(angle.($x))),
    (:Sign, x->:(sign.($x)))
]

for (OP, FUNC) in __NON_HOLOMOPHICS__

    @eval begin
        mutable struct $OP{T, O} <: AbstractBlock
            input::O
            $OP{T, O}() where {T, O} = new{T, O}()
        end

        $OP(::Type{T}, ::Type{O}) where {T, O} = $OP{T, O}()
        $OP(::Type{T};nbatch=1) where T = $OP(T, nbatch > 1 ? Matrix{T} : Vector{T})
        $OP(;nbatch=1) = $OP(Float64, nbatch=nbatch)

        function forward(op::$OP{T, O}, input::O) where {T, O}
            op.input = input
            return $(FUNC(:(input)))
        end
    end

end

function backward(op::RealOp{T, O}, grad::O) where {T <: Real, O}
    return grad
end

function backward(op::RealOp{T, O}, grad::O) where {T <: Complex, O}
    return real(grad)
end

function backward(op::ImagOp{T, O}, grad::O) where {T <: Real, O}
    return zeros(grad)
end

function backward(op::ImagOp{T, O}, grad::O) where {T <: Complex, O}
    return -im .* real(grad)
end

function backward(op::Conj{T, O}, grad::O) where {T <: Real, O}
    return grad
end

function backward(op::Conj{T, O}, grad::O) where {T <: Complex, O}
    return conj(grad)
end

function backward(op::Abs{T, O}, grad::O) where {T <: Real, O}
    return grad .* sign(op.input)
end

function backward(op::Abs{T, O}, grad::O) where {T <: Complex, O}
    return real(grad) .* exp.(-im .* op.input)
end

function backward(op::Abs2{T, O}, grad::O) where {T <: Real, O}
    return 2 .* grad .* op.input
end

function backward(op::Abs2{T, O}, grad::O) where {T <: Complex, O}
    return 2 .* real(grad) .* conj(op.input)
end

function backward(op::Arg{T, O}, grad::O) where {T <: Real, O}
    error("Do not support real number")
end

function backward(op::Arg{T, O}, grad::O) where {T <: Complex, O}
    return -im ./ op.input .* real(grad)
end

function backward(op::Sign{T, O}, grad::O) where {T <: Real, O}
    throw(ErrorException("you cannot use Sign for real gradient"))
end

function backward(op::Sign{T, O}, grad::O) where {T <: Complex, O}
    s = exp.(im .* angle(op.input))
    conj(s) ./ abs.(op.input) .* im .* imag.(s .* grad)
end


text:
using Test, CopEnt
using LinearAlgebra: det

@testset "Entropy of univraite normal" begin
    for σ in [1, 2, 3]
        x = randn(1_000) * σ
        ent = 1 / 2 * log(2 * π * ℯ * σ^2)
        @test entropy_knn(x) ≈ ent atol=0.1
    end
end

@testset "Entropy of bivariate normal" begin
    for σ1 in [1, 2, 3], σ2 in [1, 2, 3]
        x = randn(2, 1_000) .* [σ1, σ2]
        ent = 1 / 2 * log(det(2 * π * ℯ * [σ1^2 0; 0 σ2^2]))
        @test entropy_knn(x) ≈ ent atol=0.1
    end
end

@testset "Checking against pycopent" begin
    x = hcat(
        [-0.89395472, -1.48611224],
        [ 0.82243227, -1.5848612 ],
        [-1.26209653, -0.24610746],
        [-0.20907981,  1.51772724],
        [-0.59596898, -0.24640368],
        [-0.28688517,  0.6399274 ],
        [-1.32899876, -1.13641078],
        [-0.36206995, -0.42210696],
        [ 0.21879637,  1.21598259],
        [-0.48603486, -1.84245597],
    )

    cpy = hcat(
        [0.3, 0.3],
        [1. , 0.2],
        [0.2, 0.7],
        [0.8, 1. ],
        [0.4, 0.6],
        [0.7, 0.8],
        [0.1, 0.4],
        [0.6, 0.5],
        [0.9, 0.9],
        [0.5, 0.1],
    )
    entpy_euclidean = 2.884534529955401
    entpy_chebyshev = 2.9414966444360235
    copentpy_euclidean = 0.7167856845750248
    copentpy_chebyshev = 0.6760921005096859

    @test empirical_copula(x) == cpy
    @test entropy_knn(x; dist=Euclidean()) ≈ entpy_euclidean
    @test entropy_knn(x; dist=Chebyshev()) ≈ entpy_chebyshev
    @test copula_entropy(x; dist=Euclidean()) ≈ copentpy_euclidean
    @test copula_entropy(x; dist=Chebyshev()) ≈ copentpy_chebyshev
end


text:
import brainflow

brainflow.enable_dev_board_logger()

params = brainflow.BrainFlowInputParams()
board_shim = brainflow.BoardShim(Integer(brainflow.SYNTHETIC_BOARD), params)
sampling_rate = brainflow.get_sampling_rate(Integer(brainflow.SYNTHETIC_BOARD))
nfft = brainflow.get_nearest_power_of_two(sampling_rate)

brainflow.prepare_session(board_shim)
brainflow.start_stream(board_shim)
sleep(5)
brainflow.stop_stream(board_shim)
data = brainflow.get_board_data(board_shim)
brainflow.release_session(board_shim)

eeg_channels = brainflow.get_eeg_channels(Integer(brainflow.SYNTHETIC_BOARD))

bands = brainflow.get_avg_band_powers(data, eeg_channels, sampling_rate, true)
feature_vector = vcat(bands[1], bands[2])

# calc concentration
concentration = brainflow.MLModel(Integer(brainflow.CONCENTRATION), Integer(brainflow.REGRESSION))
brainflow.prepare(concentration)
print(brainflow.predict(feature_vector, concentration))
brainflow.release(concentration)

# calc relaxation
model_params = brainflow.BrainFlowModelParams(Integer(brainflow.RELAXATION), Integer(brainflow.REGRESSION))
relaxation = brainflow.MLModel(model_params)
brainflow.prepare(relaxation)
print(brainflow.predict(feature_vector, relaxation))
brainflow.release(relaxation)


text:
export cofactor_det

"""
`cofactor(A::Matrix{T})` computes the determinant of `A` using
cofactor expansion (which can be slow).
The return type of this method is a number of type `T`.
The entires in `A` can be polynomials and that won't work with
Julia's `det`.
"""
function cofactor_det(A::Matrix{T}) where T
    r,c = size(A)
    @assert r==c "Matrix must be square"
    @assert r>0 "Matrix cannot be 0-by-0"
    if r==1
        return A[1,1]
    end

    total = 0 * A[1,1]
    sign = -1
    for i=1:r
        a = A[i,1]
        sign = -sign
        if a==0
            continue
        end
        ii = vcat(1:i-1, i+1:r)
        jj = 2:r
        B = A[ii,jj]
        total += sign * a * cofactor_det(B)
    end

    return total
end


text:
# This file calculates scattering from the Helmholtz equation via the adaptive spectral method.
# Δu + k^2u = 0,
# u|Γ = 0,
# u^i = e^{im k x⋅d},
# u = u^i + u^s.
# The normal derivative ∂u/∂n of the entire wave is calculated on the sound-soft boundaries.
# The scattered wave is calculated by convolving [∂u/∂n] with the fundamental solution.
# Then, the total wave is obtained by summing the incident and the scattered waves.

# For periodic domains, the convolution property can be leveraged to make the bandwidth more precise.

using ApproxFun, SingularIntegralEquations
include("Scatteraux.jl")

K = logspace(1,3,25)
ts = Float64[]

for k in K

    ω = 2π
    d = (1,-1)
    d = d[1]/hypot(d[1],d[2]),d[2]/hypot(d[1],d[2])
    ui = (x,y) ->  exp(im*k*(d⋅(x,y)))

    # The Helmholtz Green's function, split into singular and nonsingular pieces.
    g1 = z -> -besselj0(k*abs(2sin(log(z)/2im)))/2
    g2 = z -> abs(z-1) < 5eps() ? -(log(k/2)+γ)/2/π + im/4 : im/4*hankelh1(0,k*abs(2sin(log(z)/2im))) - g1(z).*logabs(abs(2sin(log(z)/2im)))/π

    #g1 = z -> -besselj0(k*sqrt(2*(1-(z+inv(z))/2)))/2
    #g2 = z -> abs(z-1) < 5eps() ? -(log(k/2)+γ)/2/π + im/4 : im/4*hankelh1(0,k*sqrt(2*(1-(z+inv(z))/2))) - g1(z).*logabs(sqrt(2*(1-(z+inv(z))/2)))/π

    # A single circular domain.

    dof = round(Int,2.15(k+30))

    dom = Circle()
    sp = Space(dom)
    uiΓ = Fun(t->ui(real(t),imag(t)),sp,dof)
    ⨍ = DefiniteLineIntegral(sp)
    SI = SingularIntegral(sp,0)

    G1 = Fun(g1,sp,dof)
    G2 = Fun(g2,sp,dof)
    SingleLayer = Convolution(SI,G1) + Convolution(⨍,G2)

    tic()
    ∂u∂n = SingleLayer\uiΓ
    t = toc()
    push!(ts,t)
    println("The length of ∂u∂n is: ",ncoefficients(∂u∂n))
    #us = (x,y) ->  -logkernel(g1,∂u∂n,complex(x,y))-linesum(g2,∂u∂n,complex(x,y))

end

#=
using Plots
pyplot()
plot(K,ts;xscale=:log10,yscale=:log10,label="Time")
plot!(K,0.000001K.^2;xscale=:log10,yscale=:log10,label="\$ \\mathcal{O}(k^2) \$")
xlabel!("\$ k \$")
ylabel!("Execution Time (s)")
savefig("time.pdf")

using Plots
pyplot()
pad!(∂u∂n,ncoefficients(∂u∂n)-2)
plot(1:ncoefficients(G1),abs(coefficients(G1));yscale=:log10,label="\$ G_1(|x-y|) \$",)
plot!(1:ncoefficients(G2),abs(coefficients(G2));yscale=:log10,label="\$ G_2(|x-y|) \$",)
plot!(1:ncoefficients(uiΓ),abs(coefficients(uiΓ));yscale=:log10,label="\$ u^i|_\\Gamma \$",)
plot!(1:ncoefficients(∂u∂n),abs(coefficients(∂u∂n));yscale=:log10,label="\$ \\partial u / \\partial n \$",)

xlabel!("\$ n \$")
ylabel!("Laurent coefficients")
savefig("coefficients.pdf")
=#


text:
"""

$(TYPEDEF)

Structure that contains the detailed input options.

$(TYPEDFIELDS)

"""
@with_kw struct Options

  firstframe::Int = 1
  lastframe::Int = -1
  stride::Int = 1

  periodic::Bool = true

  irefatom::Int = -1
  n_random_samples::Int = 10

  binstep::Float64 = 0.02
  dbulk::Float64 = 10.
  cutoff::Float64 = 10.
  usecutoff::Bool = false

  # Linked cell length will be (cutoff/lcell), might be tunned for maximum
  # performance
  lcell::Int = 2

  # Sleep time between checks of threads for multple spawn
  sleep::Float64 = 0.01

  # Force garbage collection in parallel runs to avoid memory overflow,
  # whenever free memory in the system is smaller than GC_threshold
  GC::Bool = true
  GC_threshold::Float64 = 0.1

  # Random number generator options
  seed::Int = 321
  StableRNG::Bool = false

  # Manually chose on how many threads to run (-1 to use all threads)
  nthreads::Int = -1

  # Do not show any output on the screen on execution of mddf routines
  silent::Bool = false

end


text:
function mpccmodel_setup_newton_penalty(
        config::MPCCModelConfig
    )

    @unpack dimspec = config
    @unpack n, q, l, me, mi, r, s = dimspec

    # Std variables
    @variables x[1:n] pr[1:r] ps[1:s]

    # The masking variables (these will either be 0 or 1 when called)
    @variables m_ci_n[1:mi] m_ci_p[1:mi] m_F_n[1:l, 1:q] m_Fq[1:q]

    # The penalty coefficients
    @variables α_ci α_F β_ci γ

    println(" ** Doing max formulation **")

    # Construct the penalised objective function
    ϕ = config.defn.f
    for lp_ci=1:mi
        ϕ += α_ci * m_ci_n[lp_ci] * max(0, -config.defn.ci[lp_ci])
    end
    for lp_ci=1:mi
        ϕ += β_ci * m_ci_p[lp_ci] * config.defn.ci[lp_ci]
    end
    for lp_q=1:q
        for lp_l=1:l
            ϕ += α_F * m_F_n[lp_l, lp_q] * max(0, -config.defn.F[lp_l, lp_q])
        end
    end
    temp_Fq = Vector{Num}(undef, q)
    for lp_q=1:q
        temp_Fq[lp_q] = Num(1)
        for lp_l=1:l
            temp_Fq[lp_q] *= config.defn.F[lp_l, lp_q]
        end
        ϕ += γ * m_Fq[lp_q] * temp_Fq[lp_q]
    end

    ϕ_fn = build_function([ϕ], x, pr, ps, α_ci, α_F, β_ci, γ, m_ci_n, m_ci_p, m_F_n, m_Fq; expression=Val{false})

    function local_gradϕ(
            x::Vector{S}, pr::Vector{T}, ps::Vector{Int64},
            α_ci::S, α_F::S, β_ci::S, γ::S,
            m_ci_n::Vector{Float64}, m_ci_p::Vector{Float64}, m_F_n::Matrix{Float64}, m_Fq::Vector{Float64}
        ) where {S <: Real, T <: Real}  # ::Vector{S}

        return mm_fd_ntpen_gradϕ(dimspec, ϕ_fn[1], x, pr, ps, α_ci, α_F, β_ci, γ, m_ci_n, m_ci_p, m_F_n, m_Fq)
    end

    function local_hessϕ(
            x::Vector{S}, pr::Vector{T}, ps::Vector{Int64},
            α_ci::S, α_F::S, β_ci::S, γ::S,
            m_ci_n::Vector{Float64}, m_ci_p::Vector{Float64}, m_F_n::Matrix{Float64}, m_Fq::Vector{Float64}
        ) where {S <: Real, T <: Real}  # ::Vector{S}

        return mm_fd_ntpen_hessϕ(dimspec, ϕ_fn[1], x, pr, ps, α_ci, α_F, β_ci, γ, m_ci_n, m_ci_p, m_F_n, m_Fq)
    end

    return MPCCNewtonPenaltyFunctions(
            ϕ_fn[1],
            local_gradϕ,
            local_hessϕ
        )
end

# Define gradient wrt x
function mm_fd_ntpen_gradϕ(
        dimspec::MPCCDimSpec, ϕ_fn::Function,
        x::Vector{S}, pr::Vector{T}, ps::Vector{Int64},
        α_ci::S, α_F::S, β_ci::S, γ::S,
        m_ci_n::Vector{Float64}, m_ci_p::Vector{Float64}, m_F_n::Matrix{Float64}, m_Fq::Vector{Float64}
    ) where {S <: Real, T <: Real}

    local_ϕ(z) = ϕ_fn(z, pr, ps, α_ci, α_F, β_ci, γ, m_ci_n, m_ci_p, m_F_n, m_Fq)[1]
    return ForwardDiff.gradient(local_ϕ, x)
end

# Define gradient wrt x
function mm_fd_ntpen_hessϕ(
        dimspec::MPCCDimSpec, ϕ_fn::Function,
        x::Vector{S}, pr::Vector{T}, ps::Vector{Int64},
        α_ci::S, α_F::S, β_ci::S, γ::S,
        m_ci_n::Vector{Float64}, m_ci_p::Vector{Float64}, m_F_n::Matrix{Float64}, m_Fq::Vector{Float64}
    ) where {S <: Real, T <: Real}

    local_ϕ(z) = ϕ_fn(z, pr, ps, α_ci, α_F, β_ci, γ, m_ci_n, m_ci_p, m_F_n, m_Fq)[1]
    return ForwardDiff.hessian(local_ϕ, x)
end


text:
## tests todo (also see https://github.com/NREL-SIIP/WaterModelsAnalytics.jl/issues/6):
# - create graph of Epanet network parsed by WNTR with Python-only code (python/wntr_vis.py)
# - plots/pumps.jl (this should also use code in analysis/pump_bep.jl and utility.jl)

using WaterModelsAnalytics

import Memento
import JSON

#import PyCall # will likely use this to test python code

const _WM = WaterModelsAnalytics._WM
const _IM = WaterModelsAnalytics._WM._IM
const _MOI = WaterModelsAnalytics._WM._MOI

# Suppress warnings during testing.
Memento.setlevel!(Memento.getlogger(_IM), "error")
Memento.setlevel!(Memento.getlogger(_WM), "error")
WaterModelsAnalytics.logger_config!("error")

using Test

include("common.jl")

@testset "WaterModelsAnalytics" begin

    include("graph.jl")

    # empty test, so just a placeholder
    include("io.jl")

    include("simulation.jl")
    include("validation.jl");
end


text:
# Automatically generated test file
function super_exponential5()

g = [
54.598150033144236,
21.380942759123343,
9.487735836358526,
4.770733181967603,
2.718281828459045,
1.755054656960298,
1.284025416687741,
1.064494458917860,
1.000000000000000,
1.000000000000000,
0.935505541082140,
0.715974583312259,
0.244945343039702,
-0.718281828459045,
-2.770733181967603,
-7.487735836358526,
-19.380942759123343,
-52.598150033144236]

ζ_vec = [1690.5, 321.786, 61.2518, 11.6593, 2.21934, 0.42245, 0.0804132, 0.0153066, 0.00291361, 0.000554605]

I_sols = [
[1, 18],
[1, 17, 18],
[1, 3, 16, 18],
[1, 3, 15, 17, 18],
[1, 2, 4, 15, 17, 18],
Int[],
[1, 2, 3, 5, 14, 16, 17, 18],
[1, 2, 3, 4, 6, 14, 16, 17, 18],
[1, 2, 3, 4, 6, 13, 15, 16, 17, 18],
Int[],
[1, 2, 3, 4, 5, 7, 12, 14, 15, 16, 17, 18],
Int[],
[1, 2, 3, 4, 5, 6, 7, 12, 13, 14, 15, 16, 17, 18],
Int[],
[1, 2, 3, 4, 5, 6, 7, 8, 11, 12, 13, 14, 15, 16, 17, 18],
Int[],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]]

f_sols = [
2762.089239616717350,
1690.496808895726645,
164.518291814615623,
93.493059000276844,
19.818188371409633,
11.663223474130064,
2.861709058053748,
1.745479203571449,
0.472671424236978,
0.299212879232073,
0.088231472475854,
0.055714494648782,
0.020874875549453,
0.011366807450031,
0.000831907042993,
0.000693255869010,
-0.000000000002728]

return g, ζ_vec, I_sols, f_sols

end


text:
using LinearAlgebra
using SparseArrays

function global2local_transf_matrix(end_vert_u, end_vert_v)
    @assert length(end_vert_u) == length(end_vert_v)
    @assert length(end_vert_u) == 2 || length(end_vert_u) == 3
    xdim = length(end_vert_u)
    L = norm(end_vert_u-end_vert_v)
    @assert L > 1e-6
    # by convention, the new x axis is along the element's direction
    # directional cosine of the new x axis in the global world frame
    c_x = (end_vert_v[1] - end_vert_u[1])/L
    c_y = (end_vert_v[2] - end_vert_u[2])/L
    R = zeros(3,3)
    Γ = zeros(2,xdim*2)
    if 3 == xdim
        # error("Not Implemented")
        c_z = (end_vert_v[3] - end_vert_u[3]) / L
        if abs(abs(c_z) - 1.0) < eps()
            R[1, 3] = -c_z
            R[2, 2] = 1
            R[3, 1] = c_z
        else
            # local x_axis = element's vector
            new_x = [c_x, c_y, c_z]
            # local y axis = cross product with global z axis
            new_y = -cross(new_x, [0,0,1.0])
            new_y /= norm(new_y)
            new_z = cross(new_x, new_y)
            R[1, :] = new_x
            R[2, :] = new_y
            R[3, :] = new_z
        end
        Γ[1, 1:3] = R[1,:]
        Γ[2, 4:6] = R[1,:]
    elseif 2 == xdim
        R[1,:] = [c_x, c_y, 0]
        R[2,:] = [-c_y, c_x, 0]
        R[3,3] = 1.0
        Γ[1, 1:2] = R[1,1:2]
        Γ[2, 3:4] = R[1,1:2]
        # Γ = [c_x c_y 0 0
        #      0 0 c_x c_y]
    end
    return Γ
end

"""assemble element stiffness matrix into a sparse system stiffness matrix

Parameters
----------
k_list : list of element stiffness matrix
    a list of equal-sized global element stiffness matrix
nV : int
    number of vertices, used to compute total dof for sizing
    the output sparse matrix
id_map : Matrix
    (nE x (2xnode_dof)) index matrix:
        M[element index] = [nodal dof index]
    e.g. for a single beam
        M[1, :] = [1, 2, 3, 4, 5, 6]
    e.g. for a single truss
        M[1, :] = [1, 2, 3, 4]
"""
function assemble_global_stiffness_matrix(k_list, nV, id_map; exist_e_ids=undef)
    node_dof = Int(size(k_list[1].data,1)/2)
    total_dof = node_dof * Int(nV)

    if exist_e_ids==undef
        exist_e_ids = 1:length(k_list)
    else
        @assert Set(exist_e_ids) <= Set(1:length(k_list))
    end

    iT = Int
    T = Float64
    row = iT[]
    col = iT[]
    data = T[]
    for e_id in exist_e_ids
        Ke = k_list[e_id].data
        @assert size(Ke) == (node_dof*2, node_dof*2)
        for i in 1:node_dof*2
            for j in 1:node_dof*2
                if abs(Ke[i,j]) > eps()
                    push!(row,id_map[e_id, i])
                    push!(col,id_map[e_id, j])
                    push!(data, T(Ke[i,j]))
                end
            end
        end
    end
    K = sparse(row, col, data, total_dof, total_dof)
    return K
end

"""
compute a sparse permutation matrix that performs row permuation if multiply on the left
    Perm[i,j] = 1 to place j-th row of the target matrix to i-th row

Inputs
------
dof_stat: dof status vector
    dof_stat[i]  = 1 if i-th dof is fixed, 0 if free, -1 if not exists
"""
function compute_permutation(dof_stat)
    n_fixed_dof = sum(dof_stat)
    total_dof = length(dof_stat)
    n_free_dof = total_dof - n_fixed_dof + 1

    free_tail = 1
    fix_tail = n_free_dof
    id_map_RO = collect(1:total_dof)
    for i in 1:total_dof
        if 0 == dof_stat[i]
            id_map_RO[free_tail] = i
            free_tail += 1
        elseif 1 == dof_stat[i]
            id_map_RO[fix_tail] = i
            fix_tail += 1
        else
            error("Wrong value: $(dof_stat[i])")
        end
    end
    # @show id_map_RO

    iT = Int
    T = Float64
    # a row permuatation matrix (multiply left)
    perm_row = iT[]
    perm_col= iT[]
    perm_data = T[]
    for i in 1:total_dof
        push!(perm_row,i)
        push!(perm_col,id_map_RO[i])
        push!(perm_data,T(1))
    end
    Perm = sparse(perm_row, perm_col, perm_data, total_dof, total_dof)
    return Perm
end


text:
function eftAddGTE{T<:StdFloat}(a::T,b::T,c::T)
    s,t = eftAddGTE(b, c)
    x,u = eftAddGTE(a, s)
    y,z = eftAddGTE(u, t)
    x,y = eftAddGTE(x, y)
    x,y,z
end

function eftAddGTEas2{T<:StdFloat}(a::T,b::T,c::T)
    s,t = eftAddGTE(b, c)
    x,u = eftAddGTE(a, s)
    y   = u+t
    x,y = eftAddGTE(x, y)
    x,y
end

function eftAdd{T<:StdFloat}(a::T,b::T,c::T)
    s,t = eftAdd(b, c)
    x,u = eftAdd(a, s)
    y,z = eftAdd(u, t)
    x,y = eftAddGTE(x, y)
    x,y,z
end

function eftAddAs2{T<:StdFloat}(a::T,b::T,c::T)
    s,t = eftAdd(b, c)
    x,u = eftAdd(a, s)
    y   = u+t
    x,y = eftAddGTE(x, y)
    x,y
end

function eftMul{T<:StdFloat}(a::T, b::T, c::T)
    p,e = eftMul(a,b)
    x,y = eftMul(p,c)
    z,t = eftMul(e,c)
    y,z = eftAdd(y,z)
    x,y = eftAdd(x,y)
    eftAddGTE(x,y,z,t)
end

function eftMulas3{T<:StdFloat}(a::T, b::T, c::T)
    p,e = eftMul(a,b)
    x,y = eftMul(p,c)
    z,t = eftMul(e,c)
    y,z = eftAdd(y,z)
    x,y = eftAdd(x,y)
    eftAddGTEas3(x,y,z,t)
end

function eftMulAs2{T<:StdFloat}(a::T, b::T, c::T)
    p,e = eftMul(a,b)
    x,y = eftMul(p,c)
    z   = e*c
    y,z = eftAdd(y,z)
    x,y = eftAdd(x,y)
    x,y
end

function eftFMA{T<:StdFloat}(a::T, b::T, c::T)
    x = fma(a,b,c)
    u1,u2 = eftMul(a,b)
    a1,a2 = eftAdd(u2,c)
    b1,b2 = eftAdd(u1,a1)
    g = (b1-x)+b2
    y,z = eftAddGTE(g,a2)
    x,y,z
end

function eftFMAas2{T<:StdFloat}(a::T, b::T, c::T)
    x = fma(a,b,c)
    u1,u2 = eftMul(a,b)
    a1,a2 = eftAdd(u2,c)
    b1,b2 = eftAdd(u1,a1)
    g = (b1-x)+b2
    y = g+a2
    x,y
end

function eftFMS{T<:StdFloat}(a::T, b::T, c::T)
    x = fma(a,b,-c)
    u1,u2 = eftMul(a,b)
    a1,a2 = eftSub(u2,c)
    b1,b2 = eftAdd(u1,a1)
    g = (b1-x)+b2
    y,z = eftAddGTE(g,a2)
    x,y,z
end

function eftFMSas2{T<:StdFloat}(a::T, b::T, c::T)
    x = fma(a,b,-c)
    u1,u2 = eftMul(a,b)
    a1,a2 = eftSub(u2,c)
    b1,b2 = eftAdd(u1,a1)
    g = (b1-x)+b2
    y = g+a2
    x,y
end

for fn in (:eftAdd, :eftAddGTE, :eftAddAs2, :eftMul, :eftMulAs2, :eftMulAs4, :eftFMA, :eftFMS, :eftFMAas2, :eftFMSas2)
  @eval begin
    ($fn){T1<:Real,T2<:StdFloat,T3<:StdFloat}(a::T1, b::T2, c::T3) = ($fn)(convert(Float64,a), b, c)
    ($fn){T1<:StdFloat,T2<:Real,T3<:StdFloat}(a::T1, b::T2, c::T3) = ($fn)(a, convert(Float64,b), c)
    ($fn){T1<:StdFloat,T2<:StdFloat,T3<:Real}(a::T1, b::T2, c::T3) = ($fn)(a, b, convert(Float64,c))
    ($fn){T1<:StdFloat,T2<:Real,T3<:Real}(a::T1, b::T2, c::T3) = ($fn)(a, convert(Float64,b), convert(Float64,c))
    ($fn){T1<:Real,T2<:StdFloat,T3<:Real}(a::T1, b::T2, c::T3) = ($fn)(convert(Float64,a), b, convert(Float64,c))
    ($fn){T1<:Real,T2<:Real,T3<:StdFloat}(a::T1, b::T2, c::T3) = ($fn)(convert(Float64,a), convert(Float64,b), c)
    ($fn){T1<:Real,T2<:Real,T3<:Real}(a::T1, b::T2, c::T3) = ($fn)(convert(Float64,a), convert(Float64,b), convert(Float64,c))
  end
end


text:
include("sem.jl")

#msh = Nektar.nektar("meshes/ss-275D.rea")
msh = Nektar.nektar("meshes/chan2.rea")

nb, nverts, bmap = dof.basic_num_scheme(msh.elems)

tags = [b.tag for b in msh.bcs]

function get_dirichilet_nodes(bcs)

    vine = [1 2 3 4; 2 3 4 1]
    nb = length(bcs)

    dnodes = zeros(Int, 2, 3*nb)
    count = 1
    for i = 1:nb
        b = bcs[i]
        ed = b.edge
        el = b.elem
        dnodes[1,count] = el
        dnodes[2,count] = vine[1,ed]
        count += 1
        dnodes[1,count] = el
        dnodes[2,count] = ed + 4
        count += 1
        dnodes[1,count] = el
        dnodes[2,count] = vine[2,ed]
        count += 1
    end

    return dnodes
end

dnodes = get_dirichilet_nodes(msh.bcs[tags .== 'O'])

nbslv = dof.dirichilet_numbering(bmap, dnodes)

Q = 4

bmapQ, nb2, nbslv2 = dof.make_bmap(bmap, Q, nbslv)
#isvert = dof.make_bmap(bmap, Q, nbslv)

bmap2 = dof.rcm(bmap, nbslv)

bmapQ2, nbq2, nbslvq2 = dof.make_bmap(bmap2, Q, nbslv)


text:
using Plots
using FFTW

tan_taylor(x) = x+( x^3/3 )+( x^5*2/15.0 )+( x^7*17/315.0 )+( x^9*62/2835.0 )+( x^11*1382/155925.0 )+( x^13*21844/6081075.0 )+( x^15*929569/638512875.0 );

x = LinRange(-π/2+0.1, π/2-0.1, 1001);

tanx_taylor = tan_taylor.(x);
tanx_native = tan.(x);

delta = tanx_native-tanx_taylor

# Time domain comparison.
plot( 	x, tanx_native,
	 	label = "Julia native tan",
        size = (1024, 720),)
plot!( 	x, tanx_taylor,
	 	label = "8-term Taylor series",)

savefig("figures/tan/time_domain")

max_stars_repo_name: ErikBuer/Integer-Array

max_issues_repo_name: ErikBuer/Integer-Array

max_forks_repo_name: ErikBuer/Integer-Array


text:
module PlotGraphviz

using Graphs
using SimpleWeightedGraphs
using ParserCombinator
using ShowGraphviz

# Write your package code here

export
    # visualization, import, export:
    plot_graphviz, write_dot_file, read_dot_file,

    # modifier:
    set!, val, rm!,

    # utility:
    get_id,

    # data structs:
    GraphvizAttributes, Property, gvNode, gvEdge, gvEdges

include("./dots/attributes.jl")
include("./dots/graph_attributes.jl")
include("./dots/import_dots.jl")
include("./dots/export_dots.jl")
include("./dots/plot_graphviz.jl")
include("./dots/to_dot.jl")

end

max_stars_repo_name: tragisch/PlotGraphviz

max_issues_repo_name: tragisch/PlotGraphviz


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    Base.MathConstants

Module containing the mathematical constants.
See [`π`](@ref), [`ℯ`](@ref), [`γ`](@ref), [`φ`](@ref) and [`catalan`](@ref).
"""
module MathConstants

export π, pi, ℯ, e, γ, eulergamma, catalan, φ, golden

Base.@irrational π        3.14159265358979323846  pi
Base.@irrational ℯ        2.71828182845904523536  exp(big(1))
Base.@irrational γ        0.57721566490153286061  euler
Base.@irrational φ        1.61803398874989484820  (1+sqrt(big(5)))/2
Base.@irrational catalan  0.91596559417721901505  catalan

# aliases
"""
    π
    pi

The constant pi.

# Examples
```jldoctest
julia> pi
π = 3.1415926535897...
```
"""
π, const pi = π

"""
    ℯ
    e

The constant ℯ.

# Examples
```jldoctest
julia> ℯ
ℯ = 2.7182818284590...
```
"""
ℯ, const e = ℯ

"""
    γ
    eulergamma

Euler's constant.

# Examples
```jldoctest
julia> Base.MathConstants.eulergamma
γ = 0.5772156649015...
```
"""
γ, const eulergamma = γ

"""
    φ
    golden

The golden ratio.

# Examples
```jldoctest
julia> Base.MathConstants.golden
φ = 1.6180339887498...
```
"""
φ, const golden = φ

"""
    catalan

Catalan's constant.

# Examples
```jldoctest
julia> Base.MathConstants.catalan
catalan = 0.9159655941772...
```
"""
catalan

# loop over types to prevent ambiguities for ^(::Number, x)
for T in (AbstractIrrational, Rational, Integer, Number, Complex)
    Base.:^(::Irrational{:ℯ}, x::T) = exp(x)
end
Base.literal_pow(::typeof(^), ::Irrational{:ℯ}, ::Val{p}) where {p} = exp(p)

Base.log(::Irrational{:ℯ}) = 1 # use 1 to correctly promote expressions like log(x)/log(ℯ)
Base.log(::Irrational{:ℯ}, x::Number) = log(x)

end # module

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
include("../RLAlgorithms.jl")
using .RLAlgorithms
using Plots
using Flux
using Flux: Zygote.pullback, Optimise.update!
using Statistics: mean

struct Actor{S, A1, A2}
    model::S
    μ::A1
    logstd::A2
end
(m::Actor)(s) = (l = m.model(s); (m.μ(l), m.logstd(l)))
Flux.@functor Actor

struct Critic{C}
    c1::C
    c2::C
end
(m::Critic)(s, a) = (inp = vcat(s, a); (m.c1(inp), m.c2(inp)))
Flux.@functor Critic

function act(actor, s)
    μ, logstd = actor(s)
    std = exp.(logstd)
    ã = μ .+ std .* randn(Float32, size(std))
    logpã = sum(loglikelihood(ã, μ, logstd), dims=1)
    logpã -= sum(2f0 .* (log(2f0) .- ã .- softplus.(-2f0 * ã)), dims=1)
    return tanh.(ã), logpã
end

function run!(actor, aopt, critic, copt, env;
        epochs=100, steps=300, maxt=100, batchsize=128, α=2f-1, γ=99f-2, τ=5f-3)
    rewards = zeros(Float32, epochs)
    critictar = deepcopy(critic)
    memory = ReplayMemory{Float32, Float32, Float32}(
        length(env.observationspace), length(env.actionspace), 10*steps, batchsize
    )
    actorps = params(actor)
    criticps = params(critic)

    for i in 1:epochs
        episodes!(memory, env, steps, maxt) do s
            return act(actor, s)[1]
        end

        for _ in 1:steps
            s1, a1, r, s2, done = sample(memory)

            a2, logpa2 = act(actor, s2)
            Qtars2 = vec(min.(critictar(s2, a2)...) .- α .* logpa2)
            Qtar = r .+ (1 .- done) .* γ .* Qtars2
            cgs = gradient(criticps) do
                Q1s1a1, Q2s1a1 = critic(s1, a1)
                Flux.mse(vec(Q1s1a1), Qtar) + Flux.mse(vec(Q2s1a1), Qtar)
            end
            update!(copt, criticps, cgs)

            ags = gradient(actorps) do
                a, logpa = act(actor, s1)
                Qs1a = min.(critic(s1, a)...)
                return mean(α .* logpa .- Qs1a)
            end
            update!(aopt, actorps, ags)

            softupdate!(critictar, critic, τ)
        end
        rewards[i] = test(s -> tanh.(actor(s)[1]), env)
    end
    return rewards
end

env = CartPoleContinuous{Float32}()
aDim = length(env.actionspace)
oDim = length(env.observationspace)

actor = Actor(
    Chain(Dense(oDim, 32, relu), Dense(32, 16, relu)),
    Chain(Dense(16, aDim)),
    Chain(Dense(16, aDim, x -> min(max(x, typeof(x)(-20f0)), typeof(x)(2f0))))
)
critic = Critic(
    Chain(Dense(oDim+aDim, 32, relu), Dense(32, 16, relu), Dense(16, 1)),
    Chain(Dense(oDim+aDim, 32, relu), Dense(32, 16, relu), Dense(16, 1))
)

r = run!(actor, ADAM(), critic, ADAM(), env)
plt = plot(r, labels="Reward")
display(plt)
#savefig(plt, "sac")


text:
# Adapter pattern makes incompatible interfaces work with each other.

# In Julia, this can be done by implementing creating subtype from an abstract type
# or implement a wrapper interface.

module LinkedList

export Node, list, prev, next, value

mutable struct Node{T}
    prev::Union{Node,Nothing}
    next::Union{Node,Nothing}
    value::T
end

list(x) = Node(nothing, nothing, x)

# accessors
prev(n::Node) = n.prev
next(n::Node) = n.next
value(n::Node) = n.value

"Insert new elements at current position."
function Base.insert!(n::Node{T}, x::T) where T
    new_node = Node(n, n.next, x)
    if n.next !== nothing
        n.next.prev = new_node
    end
    n.next = new_node
    return n
end

function Base.show(io::IO, n::Node)
    print(io, "Node: $(n.value)")
    n.next !== nothing && (println(), show(n.next))
end

function test()
end

end #module

using .LinkedList
LL = list(1);
insert!(LL, 2);
insert!(next(LL), 3);
LL

#=
julia> AdapterExample.test()
-- constructing a new list with one element --
Node: 1
-- inserting new node --
Node: 1
Node: 2
=#

# -------------------------------------------------------

# Add indexing function as expected by AbstractArray
# Use Composition.

module AdapterExample

export MyArray

using ..LinkedList

struct MyArray{T} <: AbstractArray{T,1}
    data::Node{T}
end

function Base.getindex(n::MyArray, idx::Int)
    n = n.data
    for i in 1:(idx-1)
        next_node = next(n)
        next_node === nothing &&
            throw(BoundsError(n.data, idx))
        n = next_node
    end
    return value(n)
end

# Size is required
function Base.size(ar::MyArray)
    n = ar.data
    count = 0
    while next(n) !== nothing
        n = next(n)
        count += 1
    end
    return (1 + count, 1)
end

function test()
    println("--- create list ---")
    mylist = list(1)
    insert!(mylist, 2)
    insert!(next(mylist), 3)
    println(mylist)

    println("--- create array from list and indexing elements ---")
end

end #

using .AdapterExample

ar = MyArray(LL)
ar[1]
ar[2]
ar[end]

#=
julia> ar = MyArray(LL)
3×1 MyArray{Int64}:
 1
 2
 3

julia> ar[1]
1

julia> ar[2]
2

julia> ar[end]
3
=#


text:
using Documenter, WAdag

makedocs(modules=[WAdag],
	sitename="WAdag.jl",
	analytics="UA-132838790-1",
	pages=[
        "Introduction" => "index.md"
    ])

deploydocs(repo="github.com/queryverse/WAdag.jl.git")


text:
module MyExample

using ForwardDiff

include("extra_file.jl")

export my_f, derivative_of_my_f;

end

max_stars_repo_name: DomagojGalic/MyExample

max_issues_repo_name: DomagojGalic/MyExample

max_forks_repo_name: DomagojGalic/MyExample


text:
## Julia program for Gamma Metropolis-Hastings
## author: weiya <szcfweiya@gmail.com>
## date: 2018-08-21

## import function gamma_int
include("../../GenRV/gamma.jl")

function mh_gamma(T = 100, alpha = 1.5)
    a = Int(floor(alpha))
    b = a/alpha
    x = ones(T+1) # initial value: 1
    for t = 1:T
        yt = rgamma_int(a, b)
        rt = (yt / x[t] * exp((x[t] - yt) / alpha))^(alpha-a)
        if rt >= 1
            x[t+1] = yt
        else
            u = rand()
            if u < rt
                x[t+1] = yt
            else
                x[t+1] = x[t]
            end
        end
    end
    return(x)
end

# example
mh_gamma()

# comparison with accept-reject
res = mh_gamma(5000, 2.43)[2:end]
est = cumsum(res.^2) ./ collect(1:5000)

res2 = ones(5000)
for i = 1:5000
    res2[i] = rgamma(2.43, 1)
end
est2 = cumsum(res2.^2) ./ collect(1:5000)

using Plots
plot(est, label="Independent MH")
plot!(est2, label="Accept-Reject")
hline!([8.33], label="True value")

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
path = "docs/src/assets/"

using ISA, Plots
𝐶₀ = AMFMtriplet(t->exp(-t^2/5),t->200.0,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(ψ₀; timeaxis=0.0:0.001:1.0)
png(path*"IS_component.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->10t,t->25cos(t)+50,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(ψ₀;timeaxis=0.0:0.001:3.0)
png(path*"default1.png")

plot(ψ₀; timeaxis=0.0:0.001:10.0)
png(path*"costum0.png")

plot(ψ₀, timeaxis=0.0:0.001:3.0, camera=(20,50))
png(path*"costum1.png")

plot(ψ₀, timeaxis=0.0:0.001:3.0,
     yguide="imaginary", ymirror=false)
png(path*"costum2.png")

using ISA, Plots
ψ₀ = AMFMcomp(t->10t,t->25cos(t)+50,0.0)
plot(ψ₀, timeaxis=0.0:0.001:3.0)
png(path*"default3d.png")

plot(ψ₀,view="TR", timeaxis=0.0:0.001:3.0,
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"TRview3d.png")

plot(ψ₀, view="TI", timeaxis=0.0:0.001:3.0,
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"TIview3d.png")

plot(ψ₀, view="RI", timeaxis=0.0:0.001:3.0,
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"RIview3d.png")


text:
using DrWatson
@quickactivate "HyperspectraWithNeXL"
using Weave

using HyperspectraWithNeXL

ENV["DATADEPS_ALWAYS_ACCEPT"] = true
@warn """
    This script will download necessary data from a NIST file server the first time it is required.
    This download is over 500 MB and will require an active Internet connection.
    The data will be placed in the `joinpath(datadir(), "exp_raw") folder.  Outputs from the scripts
    are placed in `plotsdir()`, `papersdir()` and intermediate files are placed in
    `joinpath(datadir(),"exp_pro")`
    """

# To process individual files, set `doit=false` and comment `doit &&` out of the line you wish to process.
doit=true

mkpath(joinpath(papersdir(),"Figures"))
mkpath(joinpath(papersdir(),"Tables"))
mkpath(plotsdir())
mkpath(joinpath(datadir(),"exp_pro"))

doit && include("table_k2496_quant.jl")
doit && include("figure_maxpix.jl")
doit && include("figure_roimaps.jl")
doit && include("table_krmap.jl")
doit && include("fes_fes2_figure.jl")

doit && weave(joinpath(scriptsdir(), "quantify_k2496.jmd"), fig_ext=".svg")
doit && weave(joinpath(scriptsdir(), "kratiomapMnNodule.jmd"), fig_ext=".svg")
doit && weave(joinpath(scriptsdir(), "mapMnNodule.jmd"), fig_ext=".svg")
doit && weave(joinpath(scriptsdir(), "processMnNodule.jmd"), fig_ext=".svg")

notebookspath() = joinpath(projectdir(), "notebooks")

doit && weave(joinpath(notebookspath(), "K2496_linearOfit.ipynb"), fig_ext=".svg")
doit && weave(joinpath(notebookspath(), "quant.ipynb"), fig_ext=".svg")
doit && weave(joinpath(notebookspath(), "quantifyK2496_Benitoite.ipynb"), fig_ext=".svg")
doit && weave(joinpath(notebookspath(), "quantifyK2496_Benitoite_c.ipynb"), fig_ext=".svg")
doit && weave(joinpath(notebookspath(), "quantifyK2496_stds_c.ipynb"), fig_ext=".svg")
doit && weave(joinpath(notebookspath(), "quantifyK2496_stds.ipynb"), fig_ext=".svg")
doit && weave(joinpath(notebookspath(), "quantifyK2496.ipynb"), fig_ext=".svg")
doit && weave(joinpath(notebookspath(), "suitability.ipynb"), fig_ext=".svg")

# Do a little cleanup before the main show!
GC.gc(true)
doit && weave(joinpath(scriptsdir(), "quantifyMnNodule.jmd"), fig_ext=".svg")

using NeXLSpectrum: kill_weave_temporaries
kill_weave_temporaries(scriptsdir())
kill_weave_temporaries(notebookspath())


text:
using DashBootstrapComponents, DashHtmlComponents

alerts = html_div([
    dbc_alert("This is a primary alert", color="primary"),
    dbc_alert("This is a secondary alert", color="secondary"),
    dbc_alert("This is a success alert! Well done!", color="success"),
    dbc_alert("This is a warning alert... be careful...", color="warning"),
    dbc_alert("This is a danger alert. Scary!", color="danger"),
    dbc_alert("This is an info alert. Good to know!", color="info"),
    dbc_alert("This is a light alert", color="light"),
    dbc_alert("This is a dark alert", color="dark"),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
# EXAMPLE: scatterv
# here we have 3 ranks

using MPI

MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)

if rank == 0
    data = [7.0,8.0,9.0,10.0,11.0]
else
    data = Float64[]
end

counts=Int32[2,2,1]
data = MPI.Scatterv(data, counts, 0, comm)

@show rank, data

MPI.Finalize()

## output:
# (rank, data) = (2, [11.0])
# (rank, data) = (0, [7.0, 8.0])
# (rank, data) = (1, [9.0, 10.0])


text:
function vroom_cpp_unit_tests()
    c_func_vroom_input_tests()
    c_func_vroom_add_matrix_tests()
    c_func_vroom_add_vehicle_tests()
    c_func_vroom_add_job_tests()
    vroom_add_jobs_tests()
    vroom_add_vehicles_tests()
    vroom_create_input_tests()
    vroom_solve_tests()
    # rvrp_vroom_check_data_tests()
    vroom_transform_solution_tests()
    rvrp_vroom_solve_tests()
    rvrp_vroom_supported_features_tests()
end

function c_func_vroom_input_tests()
    data = Ref{Ptr{Cvoid}}()
    status = RVRP.@vroom_ccall(input, Cint, (Bool, Ref{Ptr{Cvoid}}),
                                      false, data)
    @test status == 0
end

function c_func_vroom_add_matrix_tests()
    data = Ref{Ptr{Cvoid}}()
    RVRP.@vroom_ccall input Cint (Bool, Ref{Ptr{Cvoid}},) false data
    mat = fill(0, 5, 5)
    mat_size = 5
    array::Vector{Int32} = reshape(mat, (mat_size*mat_size))
    status = RVRP.@vroom_ccall(add_matrix, Cint, (Ptr{Cvoid}, Cint,
            Ptr{Cint}), data[], Int32(mat_size), array)
    @test status == 0
end

function c_func_vroom_add_vehicle_tests()
    data = Ref{Ptr{Cvoid}}()
    RVRP.@vroom_ccall input Cint (Bool, Ref{Ptr{Cvoid}},) false data
    status = RVRP.@vroom_ccall(add_vehicle, Cint, (Ptr{Cvoid}, Cint,
            Cint, Cint, Cint, Bool, Cint, Cint), data[], UInt64(1),
            Int64(5), UInt16(1), UInt16(1), true, UInt32(5), UInt32(10))
    @test status == 0
end

function c_func_vroom_add_job_tests()
    data = Ref{Ptr{Cvoid}}()
    RVRP.@vroom_ccall input Cint (Bool, Ref{Ptr{Cvoid}},) false data
    tws::Vector{Int32} = [3, 10]
    twe::Vector{Int32} = [4, 15]
    status = RVRP.@vroom_ccall(add_job, Cint, (Ptr{Cvoid}, Cint, Cint,
            Cint, Cint, Cint, Ptr{Cint}, Ptr{Cint}), data[], UInt64(1),
            UInt16(1), UInt32(10), Int64(2), Cint(2), tws, twe)
    @test status == 0
end

# function rvrp_vroom_check_data_tests()
#     data = generate_random_vroom_data(10, 10, false)
#     @test_nowarn RVRP.check_data(data, RVRP.VroomSolver(;with_tw = false))
# end

############################ VROOM interface tests ##################

function vroom_create_input_tests()
    data = RVRP.parse_cvrplib(dirname(@__FILE__) *
                              "/../../data/cvrplib/P/P\\P-n16-k8.vrp")
    computed_data = RVRP.build_computed_data(data)
    vroom_ptr = RVRP.vroom_create_input(data, computed_data, false)

    #TODO add tests
end

function vroom_add_jobs_tests()
    data = RVRP.parse_cvrplib(dirname(@__FILE__) *
                              "/../../data/cvrplib/P/P\\P-n16-k8.vrp")
    computed_data = RVRP.build_computed_data(data)
    vroomdata = Ref{Ptr{Cvoid}}()
    RVRP.@vroom_ccall input Cint (Bool, Ref{Ptr{Cvoid}},) false vroomdata
    RVRP.vroom_add_jobs(vroomdata[], computed_data, data.requests,
            data.location_groups, data.product_specification_classes)

    #TODO add tests
end

function vroom_add_vehicles_tests()
    data = RVRP.parse_cvrplib(dirname(@__FILE__) *
                              "/../../data/cvrplib/P/P\\P-n16-k8.vrp")
    computed_data = RVRP.build_computed_data(data)
    vroomdata = Ref{Ptr{Cvoid}}()
    RVRP.@vroom_ccall input Cint (Bool, Ref{Ptr{Cvoid}},) false vroomdata
    RVRP.vroom_add_vehicles(vroomdata[], false, computed_data,
            data.vehicle_categories, data.location_groups, data.vehicle_sets)

    #TODO add tests
end

function vroom_solve_tests()
    data = RVRP.parse_cvrplib(dirname(@__FILE__)* "/../../data/cvrplib/P/P\\P-n16-k8.vrp")
    computed_data = RVRP.build_computed_data(data)
    vroom_ptr = RVRP.vroom_create_input(data, computed_data, false)
    sol = RVRP.vroom_solve(vroom_ptr, 1, 1)
    cost = RVRP.@vroom_ccall get_sol_cost Cint (Ptr{Cvoid},) sol
    service = RVRP.@vroom_ccall get_service Cint (Ptr{Cvoid},) sol
    duration = RVRP.@vroom_ccall get_duration Cint (Ptr{Cvoid},) sol
    wait = RVRP.@vroom_ccall get_waiting_time Cint (Ptr{Cvoid},) sol
    distance = RVRP.@vroom_ccall get_distance Cint (Ptr{Cvoid},) sol
    @test cost == 450
    @test service == 0
    @test duration == 450
    @test wait == 0
    @test distance == 0
end

function vroom_transform_solution_tests()
    data = RVRP.parse_cvrplib(dirname(@__FILE__) *
                              "/../../data/cvrplib/P/P\\P-n16-k8.vrp")
    computed_data = RVRP.build_computed_data(data)
    vroom_ptr = RVRP.vroom_create_input(data, computed_data, false)
    sol_ptr = RVRP.vroom_solve(vroom_ptr, 1, 1)
    solver = RVRP.VroomSolver(;with_tw = false)
    rvrp_sol = RVRP.vroom_transform_solution(data, solver, computed_data, sol_ptr)
    @test typeof(rvrp_sol) == RVRP.RvrpSolution
    @test rvrp_sol.instance_id == data.id
    @test rvrp_sol.cost == 450
    @test length(rvrp_sol.routes) == 8
    for route in rvrp_sol.routes
        @test (route.sequence[1].operation_type
               == route.sequence[end].operation_type == 0)
    end
    @test length(rvrp_sol.unassigned_request_ids) == 0
end

############################ RVRPSolver tests ################################

function rvrp_vroom_solve_tests()
    data = RVRP.parse_cvrplib(dirname(@__FILE__) *
                              "/../../data/cvrplib/P/P\\P-n16-k8.vrp")
    rvrp_sol = RVRP.solve(data, RVRP.VroomSolver(;with_tw = true))
    @test typeof(rvrp_sol) == RVRP.RvrpSolution
    @test rvrp_sol.instance_id == data.id
    @test rvrp_sol.cost == 450
    @test length(rvrp_sol.routes) == 8
    for route in rvrp_sol.routes
        @test (route.sequence[1].operation_type
               == route.sequence[end].operation_type == 0)
    end
    @test length(rvrp_sol.unassigned_request_ids) == 0

    data = RVRP.RvrpInstance()
    rvrp_sol = RVRP.solve(data, RVRP.VroomSolver(;with_tw = true))
    @test rvrp_sol.cost == 0.0
    @test rvrp_sol.routes == RVRP.Route[]
    @test rvrp_sol.unassigned_request_ids == String[]

end

function rvrp_vroom_supported_features_tests()
    features_list = RVRP.supported_features(RVRP.VroomSolver)
    @test RVRP.HAS_PICKUPONLY_REQUESTS in features_list[1]
    @test RVRP.HAS_DELIVERYONLY_REQUESTS in features_list[2]
    non_supported_features = BitSet([RVRP.HAS_PICKUPONLY_REQUESTS,
                                     RVRP.HAS_DELIVERYONLY_REQUESTS])
    @test !(issubset(non_supported_features, features_list[1]))
    @test !(issubset(non_supported_features, features_list[2]))
end


text:
module Math3D

using LinearAlgebra

export orthogonalize, orthogonalize!
export rotx, roty, rotz, rotxyz, rot, axisangle, trans, perspective, ortho

function orthogonalize(v::Vector, n::Vector)
	nUnit = normalize(n)
	v - nUnit * dot(v, nUnit)
end

function orthogonalize!(v::Vector, n::Vector)
	nUnit = normalize(n)
	v -= nUnit * dot(v, nUnit)
end

function rotx(m::Matrix, angle::Real)
	s = sin(angle)
	c = cos(angle)
	m[1,1] = 1
	m[2,1] = m[3,1] = m[1,2] = m[1,3] = 0
	m[2,2] = m[3,3] = c
	m[3,2] = s
	m[2,3] = -s
	return m
end

function roty(m::Matrix, angle::Real)
	s = sin(angle)
	c = cos(angle)
	m[1,1] = m[3,3] = c
	m[2,1] = m[1,2] = m[3,2] = m[2,3] = 0
	m[2,2] = 1
	m[3,1] = -s
	m[1,3] = s
	return m
end

function rotz(m::Matrix, angle::Real)
	s = sin(angle)
	c = cos(angle)
	m[1,1] = m[2,2] = c
	m[2,1] = s
	m[1,2] = -s
	m[3,1] = m[3,2] = m[1,3] = m[2,3] = 0
	m[3,3] = 1
	return m
end

# rotation with angles around each of the axes
function rotxyz(m::Matrix, xAngle::Real, yAngle::Real, zAngle::Real)
	# Tait-Bryan angles transform in the order X1Y2Z3
	sx = sin(xAngle)
	cx = cos(xAngle)
	sy = sin(yAngle)
	cy = cos(yAngle)
	sz = sin(zAngle)
	cz = cos(zAngle)

	m[1,1] = cy*cz
	m[2,1] = cx*sz+cz*sx*sy
	m[3,1] = sx*sz-cx*cz*sy

	m[1,2] = -cy*sz
	m[2,2] = cx*cz-sx*sy*sz
	m[3,2] = cz*sx+cx*sy*sz

	m[1,3] = sy
	m[2,3] = -cy*sx
	m[3,3] = cx*cy

	return m
end

rotx(angle::Real) = rotx(Array{typeof(angle)}(undef, 3, 3), angle)
roty(angle::Real) = roty(Array{typeof(angle)}(undef, 3, 3), angle)
rotz(angle::Real) = rotz(Array{typeof(angle)}(undef, 3, 3), angle)
rotxyz(xAngle::Real, yAngle::Real, zAngle::Real) = rotxyz(Array{typeof(xAngle)}(undef, 3, 3), xAngle, yAngle, zAngle)

function rot(m::Matrix, axis::Vector, angle::Real)
	u = normalize(axis[1:3])
	ux = u[1]
	uy = u[2]
	uz = u[3]
	s = sin(angle)
	c = cos(angle)

	m[1,1] = c+ux*ux*(1-c)
	m[2,1] = uy*ux*(1-c)+uz*s
	m[3,1] = uz*ux*(1-c)-uy*s

	m[1,2] = ux*uy*(1-c)-uz*s
	m[2,2] = c+uy*uy*(1-c)
	m[3,2] = uz*uy*(1-c)+ux*s

	m[1,3] = ux*uz*(1-c)+uy*s
	m[2,3] = uy*uz*(1-c)-ux*s
	m[3,3] = c+uz*uz*(1-c)

	return m
end

rot(axis::Vector, angle::Real) = rot(Array{typeof(angle)}(undef, 3, 3), axis, angle)

function axisangle(m::Matrix{T}) where T
	t = m[1,1] + m[2,2] + m[3,3]
	anglecos = (t-1)/2
	anglecos >= 1 - eps(T) && return T[0,0,1], zero(T)
	axis = [m[3,2] - m[2,3], m[1,3] - m[3,1], m[2,1] - m[1,2]]
	if anglecos <= -1 + eps(T)
		squares = T[(m[1,1]+1)/2, (m[2,2]+1)/2, (m[3,3]+1)/2]
		axis /= 4
		i = indmax(squares)
		axis[i] = squares[i]
		maxel = sqrt(squares[i])
		axis /= maxel
		return axis,T(pi)
	end
	axislen = norm(axis)
	anglesin = axislen/2
	angle = atan2(anglesin, anglecos)
	axis /= axislen
	axis, angle
end

function trans(m::Matrix, t::Vector)
	m[1,4] = t[1]
	m[2,4] = t[2]
	m[3,4] = t[3]
	return m
end

trans(t::Vector) = trans(Matrix{eltype(t)}(I, 4, 4), t)

function scale(m::Matrix, s::Vector)
  m[1,1] = s[1]
  m[2,2] = s[2]
  m[3,3] = s[3]
  m[1,2] = m[1,3] = m[2,1] = m[2,3] = m[3,1] = m[3,2] = 0
end

scale(s::Vector) = diagm(0 => vcat(s, 1))

function perspective(m::Matrix, left::Real, right::Real, top::Real, bottom::Real, near::Real, far::Real; leftHanded = false)
	zsign = leftHanded ? -1 : 1
	m[1:4, 1:4] = [ 2near/(right-left)				   0  (right+left)/(right-left)*zsign					  0 ;
									 0  2near/(top-bottom)  (top+bottom)/(top-bottom)*zsign					  0 ;
									 0				   0	 (-far-near)/(far-near)*zsign  -2far*near/(far-near) ;
									 0				   0						 -1*zsign					  0 ]
	return m
end

perspective(left::Real, right::Real, top::Real, near::Real, far::Real; leftHanded = false) = perspective(Array{typeof(near)}(undef, 4, 4), left, right, top, bottom, near, far, leftHanded = leftHanded)
perspective(m::Matrix, width::Real, height::Real, near::Real, far::Real; leftHanded = false) = perspective(m, -0.5width, 0.5width, -0.5height, 0.5height, near, far, leftHanded = leftHanded)
perspective(width::Real, height::Real, near::Real, far::Real; leftHanded = false) = perspective(Array{typeof(near)}(undef, 4, 4), width, height, near, far, leftHanded = leftHanded)
function persp_horizontal_fov(m::Matrix, hfov::Real, w_h_ratio::Real, near::Real, far::Real; leftHanded = false)
	width = 2near*tan(hfov/2)
	perspective(m, width, width / w_h_ratio, near, far, leftHanded = leftHanded)
end
persp_horizontal_fov(hfov::Real, w_h_ratio::Real, near::Real, far::Real; leftHanded = false) = persp_horizontal_fov(Array{typeof(near)}(undef, 4, 4), hfov, w_h_ratio, near, far, leftHanded = leftHanded)
function persp_vertical_fov(m::Matrix, vfov::Real, w_h_ratio::Real, near::Real, far::Real; leftHanded = false)
	height = 2near*tan(vfov/2)
	perspective(m, w_h_ratio * height, height, near, far, leftHanded = leftHanded)
end
persp_vertical_fov(vfov::Real, w_h_ratio::Real, near::Real, far::Real; leftHanded = false) = persp_vertical_fov(Array{typeof(near)}(undef, 4, 4), vfov, w_h_ratio, near, far, leftHanded = leftHanded)

function ortho(m::Matrix, left::Real, right::Real, top::Real, bottom::Real, near::Real, far::Real; leftHanded = false)
	zsign = leftHanded ? -1 : 1
	m[1:4, 1:4] = [ 2/(right-left)			   0				   0 (-right-left)/(right-left) ;
								 0  2/(top-bottom)				   0 (-top-bottom)/(top-bottom) ;
								 0			   0 -2/(far-near)*zsign	 (-far-near)/(far-near) ;
								 0			   0				   0						  1 ]
	return m
end

ortho(left::Real, right::Real, top::Real, bottom::Real, near::Real, far::Real; leftHanded = false) = ortho(Array{typeof(near)}(undef, 4, 4), left, right, top, bottom, near, far, leftHanded = leftHanded)
ortho(m::Matrix, width::Real, height::Real, near::Real, far::Real; leftHanded = false) = ortho(m, -0.5width, 0.5width, -0.5height, 0.5height, near, far, leftHanded = leftHanded)
ortho(width::Real, height::Real, near::Real, far::Real; leftHanded = false) = ortho(Array{typeof(near)}(undef, 4, 4), width, height, near, far, leftHanded = leftHanded)

end


text:
#=

GLSL Version | OpenGL Version |   Date     | Shader Preprocessor
1.10.59[1]	        2.0	    30 April 2004	    #version 110
1.20.8[2]	        2.1	    07 September 2006	#version 120
1.30.10[3]	        3.0	    22 November 2009	#version 130
1.40.08[4]	        3.1	    22 November 2009	#version 140
1.50.11[5]	        3.2	    04 December 2009	#version 150
3.30.6[6]	        3.3	    11 March 2010	    #version 330
4.00.9[7]	        4.0	    24 July 2010	    #version 400
4.10.6[8]	        4.1	    24 July 2010	    #version 410
4.20.11[9]	        4.2	    12 December 2011	#version 420
4.30.8[10]	        4.3	    7 February 2013	    #version 430
4.40.9[11]	        4.4	    16 June 2014	    #version 440
4.50.7[12]	        4.5	    09 May 2017         #version 450
4.60.5[13]	        4.6	    14 June 2018	    #version 460

=#

function imgui_setup(window::GLFW.Window)
    # setup Dear ImGui context
    ctx = CImGui.CreateContext()

    # setup Dear ImGui style
    CImGui.StyleColorsDark()
    # CImGui.StyleColorsClassic()
    # CImGui.StyleColorsLight()

    # load Fonts
    # - If no fonts are loaded, dear imgui will use the default font.
    fonts = CImGui.GetIO().Fonts
    default_font = CImGui.AddFontDefault(fonts)
    @assert default_font != C_NULL

    # You can also load multiple fonts and use `CImGui.PushFont/PopFont` to select them.
    # - `CImGui.AddFontFromFileTTF` will return the `Ptr{ImFont}` so you can store it if you need to select the font among multiple.
    # - If the file cannot be loaded, the function will return C_NULL. Please handle those errors in your application (e.g. use an assertion, or display an error and quit).
    # - The fonts will be rasterized at a given size (w/ oversampling) and stored into a texture when calling `CImGui.Build()`/`GetTexDataAsXXXX()``, which `ImGui_ImplXXXX_NewFrame` below will call.
    # - Read 'fonts/README.txt' for more instructions and details.

    # fonts_dir = joinpath(pathof(CImGui), "../..", "fonts")
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Cousine-Regular.ttf"), 15)
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "DroidSans.ttf"), 16)
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Karla-Regular.ttf"), 10)
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "ProggyTiny.ttf"), 10)
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Roboto-Medium.ttf"), 16)
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Recursive Mono Casual-Regular.ttf"), 16)
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Recursive Mono Linear-Regular.ttf"), 16)
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Recursive Sans Casual-Regular.ttf"), 16)
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Recursive Sans Linear-Regular.ttf"), 16)

    # setup Platform/Renderer bindings
    #ImGui_ImplGlfw_InitForOpenGL(s.window, true)
    ImGui_ImplGlfw_InitForOpenGL(window, false) # Don't install callbacks

    # Get and parse the the GLSL version string
	glsl = split(unsafe_string(glGetString(GL_SHADING_LANGUAGE_VERSION)), ['.', ' '])
    length(glsl) < 2 && error("Unexpected version number string. GLSL version string: $(glsl)")

    # Initialize CImGui with synthesized GLSL version number
    glsl_version = parse(Int, glsl[1]) * 100 + parse(Int, glsl[2])
    @info("GLSL version: $glsl_version")
    ImGui_ImplOpenGL3_Init(glsl_version)
end


text:
#=
P24 (*) Lotto: Draw N different random numbers from the set 1..M.
The selected numbers shall be returned in a list.
Example:
* (lotto-select 6 49)
(23 1 17 33 21 37)

Hint: Combine the solutions of problems P22 and P23.
=#

function lotto_select(n, limit)
    [rand(1:limit) for _ in 1:n]
end

selected = lotto_select(6, 49)
@assert length(selected) == 6
@assert all(x -> 1 <= x <= 49, selected)

println("Tests passed: JL-24.jl")


text:
module X86_32
  using ..Registers
  abstract RegisterSet <: Registers.RegisterSet
  import Base: copy

  immutable X86_32Arch <: Registers.Architecture
  end
  Registers.intptr(::X86_32Arch) = UInt32

  # See https://www.uclibc.org/docs/psABI-i386.pdf
  const dwarf_numbering = Dict{Int, Symbol}(
  0  => :eax,   1  => :ecx,   2  => :edx,
  3  => :ebx,   4  => :esp,   5  => :ebp,
  6  => :esi,   7  => :edi,   8  => :eip,
  9  => :eflags,
  # 10 is reserved
  (11+i => Symbol("st$i") for i = 0:7)...,
  # 19-20 are reserved
  (21+i => Symbol("xmm$i") for i = 0:7)...,
  (29+i => Symbol("mm$i") for i = 0:7)...,
  39 => :mxcsr,
  40 => :es, 41 => :cs, 42 => :ss,
  43 => :ds, 44 => :fs, 45 => :gs,
  # 46-47 are reserved
  48 => :tr, 49 => :ldtr
  # 50-92 are reserved
  )
  const inverse_dwarf = map(p->p[2]=>p[1], dwarf_numbering)
  const basic_regs = 0:9

  const gdb_numbering = Dict{Int, Symbol}(
    Dict(i=>dwarf_numbering[i] for i in 0:9)...,
    10 => :cs, 11 => :ss, 12 => :ds,
    13 => :es, 14 => :fs, 15 => :gs,
    (16+i => Symbol("st$i") for i = 0:7)...,
    24 => :fctrl, 25 => :fstat, 26 => :ftag,
    27 => :fiseg, 28 => :fioff, 29 => :foseg,
    30 => :fooff, 31 => :fop,
    (32+i => Symbol("xmm$i") for i = 0:7)...,
    40 => :mxcsr,
  )
  const inverse_gdb = map(p->p[2]=>p[1], gdb_numbering)

  function Registers.get_dwarf(::X86_32Arch, regs, sym)
    get_dwarf(regs, inverse_dwarf[sym])
  end

  function Registers.set_dwarf!(::X86_32Arch, RC, sym, val)
    set_dwarf!(RC, inverse_dwarf[sym], val)
  end

  function Registers.get_syscallarg(::X86_32Arch, regs, idx)
        @assert 1 <= idx <= 6
        get_dwarf(regs, [:ebx, :ecx, :edx, :esi, :edi, :ebp][idx])
  end

  # This operation can be performance critical, precompute it.
  const dwarf2gdbmap = [inverse_gdb[dwarf_numbering[regno]] for regno in basic_regs]
  dwarf2gdb(regno) = regno in basic_regs ? dwarf2gdbmap[regno+1] :
    inverse_gdb[dwarf_numbering[regno]]

end


text:
# automatically generated by the FlatBuffers compiler, do not modify

# module: Example

import FlatBuffers

FlatBuffers.@STRUCT struct Test
    a::Int16
    b::Int8
end
FlatBuffers.@ALIGN(Test, 2)

Test(buf::AbstractVector{UInt8}) = FlatBuffers.read(Test, buf)
Test(io::IO) = FlatBuffers.deserialize(io, Test)


text:
using SeisIO.SUDS, SeisIO.UW

qml_file  = path*"/SampleFiles/XML/ISC_2011-tohoku-oki.xml"
uw_file   = path*"/SampleFiles/UW/99011116541"
suds_file = path*"/SampleFiles/SUDS/eq_wvm1.sud"

printstyled("  read_quake wrapper\n", color=:light_green)

printstyled("    QML\n", color=:light_green)
H, R = read_qml(qml_file)
Ev1 = SeisEvent(hdr = H[1], source = R[1])
Ev2 = read_quake("qml", qml_file)
@test Ev1 == Ev2

printstyled("    SUDS\n", color=:light_green)
Ev1 = readsudsevt(suds_file)
Ev2 = read_quake("suds", suds_file)
@test Ev1 == Ev2

printstyled("    UW\n", color=:light_green)
Ev1 = readuwevt(uw_file)
Ev2 = read_quake("uw", uw_file)
@test Ev1 == Ev2


text:
using FactCheck

#push!(LOAD_PATH, "/Users/dsanders/Dropbox/papers/billiards/")

using BilliardModels

#using Vector2d

facts("Particle tests") do

    p = Particle([1, 2.], [3., 4.5])
    @fact p.x => Vector2D{Float64}(1,2)
    @fact p.v => Vector2D{Float64}(3.0,4.5)

end

facts("Disc tests") do
    d = Disc([1, 2], 3)
    @fact d.radius => 3.0

    d = Disc([0, 0], 1)
    p = Particle([-2, 0], [1, 0])

    @fact collision_time(p, d) => 1

end

facts("Plane tests") do
    o = Plane([0,0], [1, 1])
    p = Particle([-1, 0], [1, 0])

    @fact collision_time(p, o) => 1
end

facts("Billiard table tests") do
    table = Sinai_billiard(0.1)
    p = Particle([0.2, 0.], [1., 0.])

    next_collision = calculate_next_collision(p, table, nothing)
    @fact next_collision[2] => Vector2D(-1., 0.)  # new velocity after reflecting on vertical boundary
    @fact next_collision[3] => 0.3  # col]\

    p = Particle([0.2, 0.], [-1., 0.])
    next_collision = calculate_next_collision(p, table, nothing)
    @fact isa(next_collision[4], Disc) => true

end

facts("Finite plane") do
    f = FinitePlane([0., 0.], [1.,1.])
    p = Particle([0.5, 0.], [0., 1.])

    @fact collision_time(p, f) => 0.5

    p = Particle([1.5, 0.], [0., 1.])
    @fact collision_time(p, f) => -Inf
end

max_stars_repo_name: dpsanders/BilliardModel

max_issues_repo_name: dpsanders/BilliardModel

max_forks_repo_name: dpsanders/BilliardModel


text:
import LanguageServer: Document, get_text, get_offset, get_line_offsets, get_position_at, get_open_in_editor, set_open_in_editor, is_workspace_file, applytextdocumentchanges

s1 = """
123456
abcde
ABCDEFG
"""
d1 = Document("untitled", s1, false)
@test get_text(d1) == s1
@test get_offset(d1, 0, 4) == 4
@test get_offset(d1, 1, 2) == 9
@test get_line_offsets(d1) == [0, 7, 13, 21]
@test get_position_at(d1, 1) == (0, 1)
@test get_position_at(d1, 8) == (1, 1)
@test get_position_at(d1, 15) == (2, 2)
@test get_open_in_editor(d1) == false
set_open_in_editor(d1, true)
@test get_open_in_editor(d1) == true
set_open_in_editor(d1, false)
@test get_open_in_editor(d1) == false
@test is_workspace_file(d1) == false

s2 = """
12μ456
abηde
ABCDEFG"""
d2 = Document("untitled", s2, true)
@test get_offset(d2, 0, 4) == 5
@test get_offset(d2, 1, 2) == 10
@test get_line_offsets(d2) == [0, 8, 15]
@test get_position_at(d2, 1) == (0, 1)
@test get_position_at(d2, 9) == (1, 1)
@test get_position_at(d2, 17) == (2, 2)
@test is_workspace_file(d2) == true

applytextdocumentchanges(d2, LanguageServer.TextDocumentContentChangeEvent(Range(1), 0, "12"))
@test_broken get_line_offsets(d2) == [0, 8, 17]

s4 = "1234\r\nabcd"
d4 = Document("untitled", s4, false)
@test_broken get_line_offsets(d4) == [0, 5]

s5 = "1234\nabcd\n"
d5 = Document("untitled", s5, false)
@test get_line_offsets(d5) == [0, 5, 10]

s6 = "\n"
d6 = Document("untitled", s6, false)
@test get_line_offsets(d6) == [0,1]


text:
# Collection of utilities for special functions (Bessel/Hankel)

"Derivative of Hankel function of first kind, with respect to argument (z)"
function diffhankelh1(n,z)
    if n!=0
        0.5*(hankelh1(-1 + n, z) - hankelh1(1 + n, z))
    else
        - hankelh1(1, z)
    end
end

"Derivative of Bessel function of first kind, with respect to argument (z)"
function diffbesselj(n,z)
    if n!=0
        0.5*(besselj(-1 + n, z) - besselj(1 + n, z))
    else
        - besselj(1, z)
    end
end

fresnels(z::T) where T = T((1+im)//4*(erf((1+im)//2*sqrt(T(π))*z) - im*erf((1-im)//2*sqrt(T(π))*z) ))

fresnelc(z::T) where T = T((1-im)//4*(erf((1+im)//2*sqrt(T(π))*z) + im*erf((1-im)//2*sqrt(T(π))*z) ))


text:
module QuantumESPRESSOFormatter

include("Inputs/Inputs.jl")

end


text:
# Original ResNet from He et al. (2015): https://arxiv.org/pdf/1512.03385.pdf with optional normalization (BatchNorm, Ioffe et al. (2015), https://arxiv.org/pdf/1502.03167.pdf)
# Author: Gabrio Rizzuti, grizzuti3@gatech.edu
# Date: October 2020

export ResNet

function ResNet(n_in::Int64, n_hidden::Int64, nblocks::Int64; k::Int64=3, p::Int64=1, s::Int64=1, norm::Union{Nothing, String}="batch", n_out::Union{Nothing, Int64}=nothing)

    resnet_blocks = Array{Any, 1}(undef, nblocks)
    for i = 1:nblocks-1
        # Normalization layer
        (norm == "batch")  && (NormLayer = BatchNorm(n_hidden))
        (norm === nothing) && (NormLayer = identity)

        # Skip-connection
        resnet_blocks[i]   = SkipConnection(Chain(Conv((k, k), n_in => n_hidden; stride = s, pad = p),
                                                  NormLayer, x->relu.(x),
                                                  Conv((k, k), n_hidden => n_in; stride = s, pad = p)), +)
    end

    # Last layer
    if isnothing(n_out)
        # Normalization layer
        (norm == "batch")  && (NormLayer = BatchNorm(n_hidden))
        (norm === nothing) && (NormLayer = identity)

        # Skip-connection
        resnet_blocks[end] = SkipConnection(Chain(Conv((k, k), n_in => n_hidden; stride = s, pad = p),
                                                  NormLayer, x->relu.(x),
                                                  Conv((k, k), n_hidden => n_in; stride = s, pad = p)), +)
    else
        # Simple convolution
        resnet_blocks[end] = Conv((k, k), n_in => n_out; stride = s, pad = p)
    end
    return FluxBlock(Chain(resnet_blocks...))

end


text:
push!(LOAD_PATH, abspath("./src"))

using QMPS

# Creating Cluster-state QMPS circuit.
nBitT1 = 4
vBit1 = 1
rBit1 = 1
mpsCS = MPSC("CS", nBitT1, vBit1, rBit1)

# Creating Differentiable QMPS circuit.
nBitT2 = 4
vBit2 = 2
rBit2 = 1
depth2 = 2
mpsDC = MPSC(("DC", depth2), nBitT2, vBit2, rBit2)

# Calculating quantum differentiations.
# using Pkg; Pkg.add("Yao") # Un-comment this line if you haven't installed package Yao.
using Yao
nBitT0 = 6
depth0 = 4
# Differentiable circuit
c0 = DCbuilder(nBitT0, depth0).circuit |> markDiff
# QMPS circuit
c1 = deepcopy(mpsDC.circuit)
# MPS extended circuit
c2 = deepcopy(mpsDC.cExtend)
# Witness operators
op0 = put(nBitT0, 1=>X)
op1 = put(nqubits(c1), 1=>X)
op2 = put(nqubits(c2), 1=>X)
# Differentiable gates
dGates0 = collect_blocks(QMPS.QDiff, c0)
dGates1 = collect_blocks(QMPS.QDiff, c1)
dGates2 = collect_blocks(QMPS.QDiff, c2)
# Default registers
reg0 = zero_state(nqubits(c0))
reg1 = zero_state(nqubits(c1), nbatch=50000)
reg2 = zero_state(nqubits(c2))
# Quantum differentiations
grads = getQdiff!.( ()->(reg0 |> c0), dGates0, Ref(op0) )
grads = getQdiff!.( ()->(reg1 |> c1), dGates1, Ref(op1) )
grads = getQdiff!.( ()->(reg2 |> c2), dGates2, Ref(op2) )