text:
#40*sum(Z.*phi.(40*(Z-0.04)))/x.tt[end]
using Bridge
using Test
using Distributions
using Bridge: runmean, lp

#srand(1234)
Random.seed!(123)

PYPLOT = false
PYPLOT && using PyPlot
import Bridge: increment, expint, lp
import Distributions.pdf

N = 3 # number of thetas (so N+1 bins)
#b = cumsum(0.3:0.4:2.0)
if !isdefined(:simid)
    error("provide simid = {1,2,3}")
end

#sim = [:gamma, :sumgamma, :fire][simid]
sim = [:sumgamma, :sumgamma, :fire][simid]

if sim == :fire

    b = Float64[2.0]
    N = length(b)
end

T = 2000.0
n = 10000 # number of increments
m = 20 # number of augmentation points per bridge exluding the left endpoint
#m = 1
if sim == :fire
    #m = 2000
    m = 1000
end

tt = range(0.0, stop=T, length=n + 1)

iterations = [100_000, 200_000][2]

# two gamma processes
beta0a = 0.4
beta0b = beta0a/10
alpha0a = 2.0
alpha0b = alpha0a/10

transdim = [false, true, true][simid]
#fixalpha = transdim
fixalpha = false

if sim == :gamma
    simname = "gammap"
    alpha0 = alpha0a
    beta0 = beta0a
elseif sim == :sumgamma
    simname = "sumgamma"
    if transdim
        simname = "sumgammatrans"
    end
    beta0 = beta0a + beta0b
    alpha0 = (beta0a*alpha0a + beta0b*alpha0b)/(beta0a + beta0b)

elseif sim == :fire
    simname = "danishfire"

    data = readdlm("output/danish.csv", header=true)[1]

    tt_ = Float64.(Dates.value.(Date.(data[:,1])))/365
    dt = tt_[2]-tt_[1]
    tt_ -= (tt_[1] - dt)
    prepend!(tt_, 0.0)

    dxx = Float64.(data[:,2])

    n = length(tt_)-1
    T = tt_[end]
    tt = tt_

    beta0 = 1/dt*1.812 # Maximum likelihood estimate from the R script
    alpha0 = 0.588942

end

beta = beta0

P0a = GammaProcess(beta0a, alpha0a)
dxxa = rand.(increment.(diff(tt), P0a))
#X0a = sample(tt, P0a)
X0a = SamplePath(tt, Bridge.cumsum0(dxxa))

P0b = GammaProcess(beta0b, alpha0b)
#X0b = sample(tt, P0b)
dxxb = rand.(increment.(diff(tt), P0b))
X0b = SamplePath(tt, Bridge.cumsum0(dxxb))

P0 = GammaProcess(beta0, alpha0)
dxx0 = rand.(increment.(diff(tt), P0))
X0 = SamplePath(tt, Bridge.cumsum0(dxx0))

dt = mean(diff(X0a.tt))
PD = increment(dt, P0a)

@test shape(PD) ≈ beta0a*dt
@test scale(PD) ≈ 1/alpha0a
@test beta0a*dt/alpha0a ≈ mean(increment(dt, P0a))
@test beta0a*dt/alpha0a^2 ≈ var(increment(dt, P0a))

"""
    hist1(xx, r)

Count
"""
function hist1(xx, r; prob = false)
    c = zeros(Int, length(r) + 1)
    for x in xx
        c[first(searchsorted(r, x))] += 1
    end
    if prob
        c / sum(c)
    else
        c / 1
    end
end

inc(B::SamplePath) = B.tt[end] - B.tt[1], B.yy[end] - B.yy[1]

pdf(x, beta, alpha) = beta/x*exp(-alpha*x)

#lp(s, x, t, y, P::GammaProcess) = logpdf(increment(t-s, P), y-x)

function dlp(dt, dx, Po::GammaProcess, P::GammaProcess)
    assert(P.λ == Po.λ)
    alpha = dt*P.γ
    alphao = dt*Po.γ
    #log(y-x)*(alphao-1) - (y-x)*P.λ - lgamma(alphao) + log(P.λ)*alphao - (log(y-x)*(alpha-1) - (y-x)*P.λ - lgamma(alpha) + log(P.λ)*alpha)
    log(dx)*(alphao-alpha)  + lgamma(alpha) - lgamma(alphao) + log(P.λ)*(alphao - alpha)
end

function lp2(s, x, t, y, P::GammaProcess)
    alpha = (t-s)*P.γ
    log(y-x)*(alpha-1) - (y-x)*P.λ - lgamma(alpha) + log(P.λ)*alpha
end
#\frac{x^{\alpha-1} e^{-x/\theta}}{\Gamma(\alpha) \theta^\alpha},

function addprc(X, Y)
    yy = copy(X.yy)
    for i in 1:length(yy)
        yy[i] += Y.yy[last(searchsorted(Y.tt, X.tt[i]))]
    end
    SamplePath(X.tt, yy)
 end

function scalebridge!(yy, xx, u, v)
    assert(xx[1] == 0)
    yy[:] = (xx-xx[1]) .* ((v-u)/(xx[end] - xx[1])) .+ u
end

if sim == :gamma
    X = X0a
    dxx = dxxa
elseif sim == :sumgamma
    X = addprc(X0a, X0b)
    dxx = dxxa + dxxb
else
    X = SamplePath(tt, Bridge.cumsum0(dxx))
end

alphahat = mean(dxx)/var(dxx) # use estimate for alpha

# grid points
# b = quantile(increment(dt, GammaProcess(beta0,alpha0)),(1:(N))/(N+1)) # theoretical
if sim != :fire && N >= 1
    #b = cumsum(0.5:0.25:2.5)[1:N]
    b = cumsum(0:0.5:5)[2:N+1]
    b = [1.,2.,4.,8.][1:N]
    #assert(N == 5)
    #b = [0.3, 1.0, 2.1, 3.6, 5.5] #fixme
    #b = quantile(diff(X.yy), (N:2N-1)/(2N)) # first bin resembles 50% of emperical increment distributions.
elseif N == 0
    b = Float64[]
end

#println("P(Y < b1) = ", mean(diff(X.yy) .< b1))
h = hist1(diff(X.yy), b) # note that P(Y < b[1]) may be inaccurate if dt is chosen small
if any(h .< 20)
    @warn("Less than 20 observations in bin")
end

if PYPLOT
    plot(XY.tt, XY.yy)
end

## Prior

vpi = 2.0
epi = 2.0

Pi = Gamma(epi^2/vpi, vpi/epi) # alpha
Pi2 = Normal(0., 10.0) # theta
Pi3 = Normal(0., 50.0) # rho
Pi4 = Uniform(0.1, 1000.0) # beta
assert(mean(Pi) ≈ epi)
assert(var(Pi) ≈ vpi)

if sim == :fire
    epi = 0.75
    vpi = 0.6^2
    Pi = Gamma(epi^2/vpi, vpi/epi) # alpha
    Pi2 = Normal(0., 5.0) # theta
    Pi3 = Normal(0., 20.0) # rho
    epi = 90.0
    vpi = 50.0^2
    Pi4 = Gamma(epi^2/vpi, vpi/epi) # beta
    #Pi4 = Uniform(0.1, 1000.0) # beta
end

lpi(alpha, theta, rho) = logpdf(Pi, alpha) + sum(logpdf.(Pi2, theta)) + sum(logpdf.(Pi3, rho))
lpi(alpha, theta, rho, beta) = logpdf(Pi, alpha) + sum(logpdf.(Pi2, theta)) + sum(logpdf.(Pi3, rho)) + logpdf(Pi4, beta)
#initialize

#...

# alpha, beta = truth

B = Vector{typeof(X)}(n)
Z = Vector{typeof(X)}(n)
Bº = Vector{typeof(X)}(n)
Zº = Vector{typeof(X)}(n)

#Bº = Vector{typeof(X)}(n)
tt = X.tt
yy = X.yy

# Bookkeeping
mkpath(joinpath("output", simname))
try # save cp of this file as documentation
    cp(@__FILE__(), joinpath("output",simname,"$simname.jl"); force=true)
catch
end

beps = 0.0
#beps = b[1]/5
#yy = diff(vcat(B...).yy)
#var(yy[yy.< b[1]/4])

#alpha = alpha0
alpha = alphahat

beta = beta0

alpha = alpha
if transdim
    beta = 0.9*beta
end

#c = beta*(T/(n*m*alpha)) *(1-exp(-alpha*beps)) # compensator for small jumps
epsij = 1e-50

theta = zeros(N)
rho = zeros(N)

#alpha, beta = [2.83095, 162.025]
#theta = [ 1.89115, -1.9195]
#rho =  [-10.5676, 2.03732]

# prior chain step std
alphasigma = 0.025
thsigma = 0.025*ones(N)
#thsigma[1] = 0.02 #Fixme
rhosigma = 0.15*ones(N) # variance of rho1 = 0

#if N > 0 && sim != :fire
#    rhosigma[1] = 0.0
#end
betasigma = 0.01

if sim == :fire
    alphasigma = 0.05
    betasigma = 1.0
    thsigma = 0.05*ones(N)
    rhosigma = 0.1*ones(N)
end
if isdefined(:nomcmc) && nomcmc
    error("don't run: nomcmc == true")
end

open(joinpath("output", simname,"truth.txt"), "w") do f
    bn = join(["b$i" for i in 1:length(b)], " ")
    println(f, "N alpha0 beta0 T n $bn beps prior1 prior2 ")
    println(f, N, " ", join(round.([alpha0, beta0, T, n, b... ],3)," "), " $beps \"$Pi\" \"$Pi2\" ")
end

# initial augmentation

for i in 1:n
    P0 = GammaProcess(beta, alpha)
    Pº = GammaBridge(tt[i+1], yy[i+1], P0)
    tti = range(tt[i], stop=tt[i+1], length=m+1)
    Z[i] = sample(tti, P0, 0.0)
    Zº[i] = copy(Z[i])

    B[i] = copy(Z[i])
    scalebridge!(B[i].yy, Z[i].yy, yy[i], yy[i+1])
    Bº[i] = copy(B[i])
end
# workspaces
zz = copy(B[1].yy)
bb = copy(B[1].yy)
if PYPLOT
    clf()
    collect(plot(B[i].tt, B[i].yy) for i in 1:10)
    plot(X.tt[1:11], X.yy[1:11], "*")
end

# write header line for csv
open(joinpath("output", simname, "params.txt"), "w") do f
    thn = join(["theta$i" for i in 1:length(theta)], " ")
    rhon = join(["rho$i" for i in 1:length(rho)], " ")
    println(f, "n alpha beta $thn $rhon")
end

P0 = GammaProcess(beta, alpha)
P = LocalGammaProcess(P0, theta, rho, b)

mc = mcstart([alpha; beta; theta; rho])
thacc = thtotal = 0
Bacc = 0
alphaacc = alphatotal = 0
betaacc = betatotal = 0

for iter in 1:iterations
    # logging
    mc = mcnext(mc, [alpha; beta; theta; rho])
    open(joinpath("output", simname, "params.txt"), "a") do f
        println(f, iter, " ", join(round.([alpha; beta; theta; rho], 8), " "))
    end

    # compensator for small jumps
    #c = beta*(T/(n*m*alpha)) * (1 - exp(-alpha*beps))
    c = 0.0

    P0 = GammaProcess(beta, alpha)
    P = LocalGammaProcess(P0, theta, rho, b)
    # sample bridges
    for i in 1:n
        #Pº = GammaBridge(tt[i+1], yy[i+1], P0)
        #Bº = SamplePath(B[i].tt, zz)
        #sample!(Bº, Pº, yy[i])
        sample!(Zº[i], P0, 0.0)
        scalebridge!(Bº[i].yy, Zº[i].yy, yy[i], yy[i+1])

        ll = llikelihood(Bº[i], P, c) - llikelihood(B[i], P, c)
       # print(i, " ", ll, " ")
        if rand() < exp(ll)
            B[i], Bº[i] = Bº[i], B[i]
            Z[i], Zº[i] = Zº[i], Z[i]
            Bacc += 1
        end
    end

    # sample parameters
    # update theta and rho
    if !transdim || iter % 5 < 2 # remember to update formula for acceptance rates
        alphaº = alpha + alphasigma*randn()

        P0 = GammaProcess(beta, alpha)
        P0º = GammaProcess(beta, alphaº)

        if N > 0
            thetaº = theta + thsigma.*randn(length(theta)) - (alphaº - alpha)
            rhoº = rho + rhosigma.*randn(length(rho))
        else
            rhoº = rho
            thetaº = theta
        end
        thtotal += 1

        if N == 0 && alphaº < eps()
        elseif N > 0 && (thetaº[end] + alphaº < eps() || alphaº < eps())
            # reject
        elseif N == 0
            P = LocalGammaProcess(P0, theta, rho, b)
            Pº = LocalGammaProcess(P0º, thetaº, rhoº, b)

            ll = 0.0
            for i in 1:n
                ll += llikelihood(B[i], Pº, P, c)
            end
            print("$iter \t\t\t\t paramsº: ", round(ll, digits=5), " ", round.([alphaº; thetaº; rhoº], 3))
            if rand() < exp(ll + lpi(alphaº, thetaº, rhoº, beta) - lpi(alpha, theta, rho, beta))
                print("✓")
                theta = thetaº
                rho = rhoº
                alpha = alphaº
                thacc += 1
            end
            println()
        else
            P = LocalGammaProcess(P0, theta, rho, b)
            Pº = LocalGammaProcess(P0º, thetaº, rhoº, b)

            ll = 0.0
            for i in 1:n
                ll += llikelihood(B[i], Pº, P, c)
            end
            print("$iter \t\t\t\t paramsº: ", round(ll, digits=5), " ", round.([alphaº; thetaº; rhoº], 3))
            if rand() < exp(ll + lpi(alphaº, thetaº, rhoº, beta) - lpi(alpha, theta, rho, beta))
                print("✓")
                theta = thetaº
                rho = rhoº
                alpha = alphaº
                thacc += 1
            end
            println()
        end
    end

    if transdim && iter % 5 >= 2
        betaº = beta + betasigma*randn()
        betatotal += 1
        if betaº <= 0
            #reject
        else
            ll = 0.0
            P0 = GammaProcess(beta, alpha)
            P0º = GammaProcess(betaº, alpha)
            P = LocalGammaProcess(P0, theta, rho, b)
            Pº = LocalGammaProcess(P0º, theta, rho, b)

            if betaº >= beta
                PD = GammaProcess(betaº - beta, alpha)

                for i in 1:n
                    sample!(Zº[i], PD, 0.0)
                    Zº[i].yy[:] += Z[i].yy

                    scalebridge!(Bº[i].yy, Zº[i].yy, yy[i], yy[i+1])

                    ll += llikelihood(Bº[i], Pº, c) - llikelihood(B[i], P, c) + dlp(tt[i+1] - tt[i], dxx[i], P0º, P0)
                end
            else # thinning
                for i in 1:n

                    Zº[i].yy[1] = 0.0
                    dtti = step(range(tt[i], stop=tt[i+1], length=m+1))

                    Be = Beta(dtti*betaº, dtti*(beta-betaº))
                    for j in 2:length(Zº[i].tt)
                        #dt = Z[i].tt[j] - Z[i].tt[j-1]
                        #Be = Beta(dt*betaº, dt*(beta-betaº))
                        Zº[i].yy[j] = Zº[i].yy[j-1] + (Z[i].yy[j] - Z[i].yy[j-1])*rand(Be)
                    end

                    scalebridge!(Bº[i].yy, Zº[i].yy, yy[i], yy[i+1])
                    ll += llikelihood(Bº[i], Pº, c) - llikelihood(B[i], P, c) + dlp(tt[i+1] - tt[i], dxx[i], P0º, P0)
                end

            end
            print("$iter \t ")
            printstyled("betaº: ", round(ll, digits=5), " [", round(betaº, digits=3), "]", color=:green)
            #println(llº, " ", ll, " ", lpi(alpha, theta, rho, betaº), " ", lpi(alpha, theta, rho, beta))
            if rand() < exp(ll + lpi(alpha, theta, rho, betaº) - lpi(alpha, theta, rho, beta))
                for i in 1:n
                    B[i], Bº[i] = Bº[i], B[i]
                    Z[i], Zº[i] = Zº[i], Z[i]
                end
                beta = betaº
                betaacc += 1
                printstyled("✓", color=:green)
            end
            println()
        end
    end
    #=
    if  !fixalpha && iter % 5 == 2 # remember to update formula for acceptane rates
        alphaº = alpha + alphasigma*randn()

        if alphaº < 0 || (N > 0 && theta[end] + alphaº < eps()) # according to Wilkinson
            # reject
        else
            P0º = GammaProcess(beta, alphaº)
            P0 = GammaProcess(beta, alpha)

            Pº = LocalGammaProcess(P0º, theta, rho, b)
            P = LocalGammaProcess(P0, theta, rho, b)

            ll = 0.0
            for i in 1:n
                ll += llikelihood(B[i], Pº, P)
            end
            print("$iter \t alphaº: ", round(ll, 5), " [", round(alphaº, 3), "]")
            if rand() < exp(ll + lpi(alphaº, theta, rho, beta) - lpi(alpha, theta, rho, beta))
                alpha = alphaº
                alphaacc += 1
                print("✓")
            end
            alphatotal += 1
            println()
        end
    end
    =#

end

println("alpha acc ", alphaacc/max(1,alphatotal))
println("beta acc ", betaacc/betatotal)
println("theta acc ", thacc/thetatotal)
println("posterior band")
display(hcat(mcband(mc)...))
println("posterior band for mean")
display(hcat(mcbandmean(mc)...))

params = readdlm(joinpath("output", simname, "params.txt"), Float64; skipstart=1)[:,2:end];

function plotparams(simname = "gammap")
    params = readdlm(joinpath("output", simname, "params.txt"), Float64; skipstart=1)[:,2:end];
    clf();
    plot(params[:,1], color=:blue , lw = 0.7)
    plot(params[:,2], color=:blue , lw = 0.7)

    plot(params[:,3:end], color=:lightblue , lw = 0.4)
    plot(runmean(params), color=:darkgrey, lw = 0.4)

    params
end

Ph = GammaProcess( mean(params,1)[2:-1:1]... )
dxxh = rand.(increment.(diff(tt), Ph))
Xh = SamplePath(tt, Bridge.cumsum0(dxxh))

#assert(c==0)

figure();
z = sort(vcat((diff(cumsum(shuffle(diff(X2.yy)))[1:m:end]) for i in 1:100)...))

plot(z[1:100:end], range(0, stop=1, length=length(z[1:100:end])), label="augmented")
z = sort(diff(X.yy))
plot(z, range(0, stop=1, length=length(z)), label="obs")
#plot(z, cdf.(Gamma( var(z)/mean(z), mean(z)^2/var(z)), z))
#plot(z, cdf.(Gamma(beta0*dt, 1/alpha0), z), label="Gamma ML")
legend()


text:
# --------- SOLUTION --------------------------------------------
# ----------------------------- LOGISIC ------------------------
using JuMP,Ipopt,DataFrames, GLM, ForwardDiff
m=nothing
m = Model(solver=IpoptSolver())  # sum( -1 * y .* (x*beta) + log(1+exp((x*beta))) )
dfd = readtable("/home/iriadmin/g/JLib/src/StatOpt/logreg.csv")
x = Array(dfd[:,2:end])
y = Array(dfd[:,1])
(rows, cols) = size(x)
g = glm(y ~ -1 + x0 + x1 + x2 + x3 + x4 + x6 + x7, dfd,  Bernoulli(), LogitLink())
@variable(m, B[z=1:cols], start=0)  # getvalue(B)  @defVar( m, 0.0 <= b[1:(K1-1)] <= 100.0 )
@NLobjective( m, Min, sum{ ( -1*y[r]*sum{x[r,c]*B[c],c in 1:cols} ) + log(1+exp(sum{x[r,cc]*B[cc],cc in 1:cols})) ,r in 1:rows} )

statusM = solve( m )
df = DataFrame( beta = getvalue( B ) )

function mynllf(beta::Vector)
	nll = 0.0
	for r in 1:rows
		dbl_part1 = 0.0
		for c in 1:cols
			dbl_part1 = dbl_part1 + x[r,c] * beta[c]
		end
		nll = nll -1.0 * y[r] * dbl_part1 + log(1+exp(dbl_part1))
	end
	return nll
end

sqrt(diag(inv(  ForwardDiff.hessian(mynllf,getvalue(B))  )))
getvalue(B)

# ----------------------- POISON --------------------------------
using JuMP,Ipopt,DataFrames, GLM
m=nothing
m = Model(solver=IpoptSolver())
dfd = readtable("/home/iriadmin/g/JLib/src/StatOpt/ZD/poison_gamma/y.csv")
names!(dfd,[:y,:x1, :x2, :x3,:x4,:x5])
x = Array(dfd[:,2:end])
y = Array(dfd[:,1])
(rows, cols) = size(x)
g = glm(y ~ -1 + x1 + x2 + x3 + x4 +x5 , dfd,  Poisson(), LogLink())
@variable(m, B[z=1:cols], start=0)
@NLobjective( m, Min, sum{ exp(sum{x[r,cc]*B[cc],cc in 1:cols}) - (y[r]*sum{x[r,c]*B[c],c in 1:cols}) ,r in 1:rows} )
statusM = solve( m )
df = DataFrame(beta=getvalue(B))
mynllf(beta::Vector) = sum(exp(x*beta) .- ((x*beta).*y))
sqrt(diag(inv(ForwardDiff.hessian(mynllf,getvalue(B)))))

# ----------------------- GAMMA - WORKING On THIS --------------------------------
using JuMP,Ipopt,DataFrames, GLM
m=nothing
m = Model(solver=IpoptSolver())
dfd = readtable("/home/iriadmin/g/JLib/src/StatOpt/ZD/poison_gamma/g.csv")
names!(dfd,[:y,:x1, :x2, :x3,:x4,:x5])
x = Array(dfd[:,2:end])
y = Array(dfd[:,1])
(rows, cols) = size(x)
g = glm(y ~ -1 + x1 + x2 + x3 + x4 +x5 , dfd,  Gamma(), LogLink())
@variable(m, B[z=1:cols], start=0)
@NLobjective( m, Min, sum{ exp(-1.0*sum{x[r,cc]*B[cc],cc in 1:cols}) + (sum{x[r,c]*B[c],c in 1:cols}) ,r in 1:rows} )
statusM = solve( m )
df = DataFrame(beta=getvalue(B))
mynllf(beta::Vector) = sum(exp(-1.0*(x*beta)) .+ ((x*beta)))
sqrt(diag(inv(ForwardDiff.hessian(mynllf,getvalue(B)))))

str_path <- "C:\\Users\\mdzhd\\Documents\\Bayesian_Approach_with_Lampros\\"

myd <- read.csv( paste0( str_path, "gamma_regression_ex2.csv" ) )

myglm <- glm( y ~ x1 + x2 + x3, data = myd, family = Gamma( link = log ) )

x <- as.matrix( myd[, -c(1)] )
y <- as.matrix( myd[, 1] )

NegLogLike <- function( beta )
{
	xbeta <- x %*% beta
	Li <- y * exp( -1.0 * xbeta ) +  xbeta
	sum( Li )
}

beta_initial <- as.matrix( rep( 0, ncol( x ) ), ncol(x), 1 )

myoptim <- optim( beta_initial, f = NegLogLike, g = NULL, method = "BFGS", hessian = TRUE)

u <- exp( x %*% as.matrix( myoptim$par , ncol( x ), 1 ) )

v <- ( y - u )/u

v2 <- v * v

phi <- sum( v2 ) / ( nrow( x ) - ncol(x ) )

myhessian <- myoptim$hessian / phi

var_cov_m <- solve( myhessian )

beta_se <- sqrt( diag( var_cov_m ) )

zvalue <- myoptim$par / beta_se

pvalue <- 2 * ( 1 - pnorm(abs( zvalue ) ) )

coeff_outputs <- cbind( myoptim$par, beta_se, zvalue, round( pvalue, 6 ) )

colnames( coeff_outputs ) <- c( "Estimates", "Std. Error", "z value", "Pr(>|z|)" )
rownames( coeff_outputs ) <- c()

print( coeff_outputs )
summary( myglm )$coefficients

# --------- SOLUTION END --------------------------------------------

# ----------------------------------------------------------------------
# -------- OLD WORKING.... but wrong -----------------------------------
# ----------------------------------------------------------------------

using JuMP;
using Ipopt;
using DataFrames;
my_table = readtable("/home/iriadmin/g/jlib/src/StatOpt/newdata.csv")
( N, K1 ) = size( my_table )
m = Model(solver=IpoptSolver())
@defVar( m, 0.0 <= b[1:(K1-1)] <= 100.0 )
@NLobjective( m, Min, sum{ ( my_table[i,1] - sum{ my_table[i,j] * b[ j - 1] , j in 2:K1} )^2 , i in 1:N} )
statusM = solve( m )
df = DataFrame( beta = getvalue( b ) )

# ----------------------- GT LOGISTIC --------------------------------
using JuMP,Ipopt,DataFrames, GLM
m=nothing
m = Model(solver=IpoptSolver())  # sum( -1 * y .* (x*beta) + log(1+exp((x*beta))) )
dfd = readtable("/home/iriadmin/g/jlib/src/StatOpt/logreg.csv")
x = Array(dfd[:,2:end])
y = Array(dfd[:,1])
(rows, cols) = size(x)
g = glm(y ~ -1 + x0 + x1 + x2 + x3 + x4 + x6 + x7, dfd,  Bernoulli(), LogitLink()) # X0 s the replacement for intercept - otherwise glm will generate a [1...] arr
@variable(m, B[z=1:cols], start=0)  # getvalue(B)  @defVar( m, 0.0 <= b[1:(K1-1)] <= 100.0 )
@NLobjective( m, Min, sum{ ( -1*y[r]*sum{x[r,c]*B[c],c in 1:cols} ) + log(1+exp(sum{x[r,cc]*B[cc],cc in 1:cols})) ,r in 1:rows} )
#..... sum( -1 * y .* (x*B) + log(1+exp((x*B))) )

d = JuMP.NLPEvaluator(m)
MathProgBase.initialize(d, [:Grad,:HessVec,:Hess])
statusM = solve( m )  #statusM = solve( m, IpoptOptions=[("tol",1e-6)] )
df = DataFrame( beta = getvalue( B ) )

d = JuMP.NLPEvaluator(m)
MathProgBase.initialize(d, [:Grad,:HessVec,:Hess])
hess = MathProgBase.hesslag_structure(d)
diag(hess)
beta_se <- sqrt( diag( var_cov_m ) )

ext: jl

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct ImagesApi <: SwaggerApi
    client::Swagger.Client
end

"""

Create or update an image.
Param: resourceGroupName::String (required)
Param: imageName::String (required)
Param: parameters::Image (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Image
"""
function _swaggerinternal_imagesCreateOrUpdate(_api::ImagesApi, resourceGroupName::String, imageName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", Image, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/images/{imageName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "imageName", imageName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function imagesCreateOrUpdate(_api::ImagesApi, resourceGroupName::String, imageName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_imagesCreateOrUpdate(_api, resourceGroupName, imageName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function imagesCreateOrUpdate(_api::ImagesApi, response_stream::Channel, resourceGroupName::String, imageName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_imagesCreateOrUpdate(_api, resourceGroupName, imageName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Deletes an Image.
Param: resourceGroupName::String (required)
Param: imageName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_imagesDelete(_api::ImagesApi, resourceGroupName::String, imageName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/images/{imageName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "imageName", imageName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function imagesDelete(_api::ImagesApi, resourceGroupName::String, imageName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_imagesDelete(_api, resourceGroupName, imageName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function imagesDelete(_api::ImagesApi, response_stream::Channel, resourceGroupName::String, imageName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_imagesDelete(_api, resourceGroupName, imageName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets an image.
Param: resourceGroupName::String (required)
Param: imageName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: expand::String
Return: Image
"""
function _swaggerinternal_imagesGet(_api::ImagesApi, resourceGroupName::String, imageName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Image, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/images/{imageName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "imageName", imageName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function imagesGet(_api::ImagesApi, resourceGroupName::String, imageName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_imagesGet(_api, resourceGroupName, imageName, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function imagesGet(_api::ImagesApi, response_stream::Channel, resourceGroupName::String, imageName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_imagesGet(_api, resourceGroupName, imageName, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets the list of Images in the subscription. Use nextLink property in the response to get the next page of Images. Do this till nextLink is null to fetch all the Images.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ImageListResult
"""
function _swaggerinternal_imagesList(_api::ImagesApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ImageListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/images", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function imagesList(_api::ImagesApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_imagesList(_api, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function imagesList(_api::ImagesApi, response_stream::Channel, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_imagesList(_api, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets the list of images under a resource group.
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ImageListResult
"""
function _swaggerinternal_imagesListByResourceGroup(_api::ImagesApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ImageListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/images", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function imagesListByResourceGroup(_api::ImagesApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_imagesListByResourceGroup(_api, resourceGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function imagesListByResourceGroup(_api::ImagesApi, response_stream::Channel, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_imagesListByResourceGroup(_api, resourceGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Update an image.
Param: resourceGroupName::String (required)
Param: imageName::String (required)
Param: parameters::ImageUpdate (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Image
"""
function _swaggerinternal_imagesUpdate(_api::ImagesApi, resourceGroupName::String, imageName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PATCH", Image, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/images/{imageName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "imageName", imageName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function imagesUpdate(_api::ImagesApi, resourceGroupName::String, imageName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_imagesUpdate(_api, resourceGroupName, imageName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function imagesUpdate(_api::ImagesApi, response_stream::Channel, resourceGroupName::String, imageName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_imagesUpdate(_api, resourceGroupName, imageName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export imagesCreateOrUpdate, imagesDelete, imagesGet, imagesList, imagesListByResourceGroup, imagesUpdate


text:
using Test
@isdefined(TestHelpers) || include("TestHelpers.jl")

TestHelpers.@timed_testset "net_components/" begin
    include("net_components/core_ops.jl")
    include("net_components/layers.jl")
    include("net_components/nets.jl")
end


text:
# Index
# =====
#
# Index types for genomic intervals.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

module Indexes

import BGZFStreams
import BioCore
import BufferedStreams
import GenomicFeatures: Interval

include("chunk.jl")
include("bgzfindex.jl")
include("tabix.jl")
include("overlap.jl")

end # module


text:
module geom

using LinearAlgebra
using Quaternions
using ArgCheck
using PyCall

pysys = pyimport("sys")

# add python directory to python path if available
pushfirst!(PyVector(pysys."path"), joinpath("."));
isdir(joinpath(pwd(), "pyfiles/")) && pushfirst!(PyVector(pysys."path"), joinpath(pwd(), "pyfiles"));
isdir(joinpath(pwd(), "viz/pyfiles/")) && pushfirst!(PyVector(pysys."path"), joinpath(pwd(), "viz/pyfiles"));

# check we can find the relevant py files.
while !any([isfile(joinpath(p, "Quaternions.py")) for p in pysys."path"])
    @warn "Dan Holden python processing scripts not found."
    println("Please run this julia file from a directory which contains these files,\n" *
             "or contains a 'pyfiles' directory with them.")
    println("\nAlternatively please provide a path on which python can find them:")
    userinput = chomp(readline())
    if isdir(userinput)
        pushfirst!(PyVector(pysys."path"), userinput);
    else
        throw(ErrorException("Cannot find file Quaternions.py"))
    end
end
Quatpy = pyimport("Quaternions");

_rowmaj_reshape_3d(x, ix, iy, iz) = (x=reshape(x, ix, iz, iy); permutedims(x, [1,3,2]);)

# trajectory is represented as (x,z), but most quaternion ops need y (=0) too.
cat_zero_y(x::T, z::T) where T <: Number = [x, 0, z]
cat_zero_y(x::Vector{T}, z::Vector{T}) where T = hcat(x, zeros(T, length(x)), z)
cat_zero_y(X::Matrix{T}) where T = begin; @argcheck size(X,2)==2;
    hcat(X[:,1], zeros(T, size(X,1)), X[:,2]); end

# Quaternion related utils
_qimag = Quaternions.imag
_quat_list(x) = [quat(x[i,:]) for i in 1:size(x,1)]
_quat_list_to_mat(x) = reduce(vcat, [_qimag(xx)' for xx in x])
_quaterion_angle_axis_w_y(θ) = quat(cos(θ/2), 0, sin(θ/2), 0)
_xz_plane_angle(q::Quaternion) = begin; @assert q.v1 ≈ 0 && q.v3 ≈ 0; atan(q.v2, q.s)*2; end  # inverse of above
_apply_rotation(x, qrot) = qrot * x * conj(qrot)
_vec4_to_quat(x::Vector) = begin; @argcheck length(x) == 4; Quaternion(x[1], x[2], x[3], x[4]); end  # quat constructor is WEIRD o.w.
_vec2_to_quat(x::Vector) = begin; @argcheck length(x) == 2;
    Quatpy.Quaternions.between([0,0,1], cat_zero_y(x[1], x[2])).qs[1,:] |> _vec4_to_quat; end

#= ---------------------------------------------------------------------------
          FORWARD KINEMATICS (for Lagrangian skeleton representation)
       These functions used some of Holden's view.py as a starting point.
------------------------------------------------------------------------------=#

reconstruct_modelled(Y::Matrix) = reconstruct(Y, :modelled)
reconstruct_modelled64(Y::Matrix) = reconstruct(Y, :modelled64)
reconstruct_raw(Y::Matrix) = reconstruct(Y, :raw)
reconstruct_root(Y::Matrix) = reconstruct(Y, :root)[:,1,:]

"""
    reconstruct(Y, input_type)

Reconstruct the absolute positions of joints in a contiguous set of frames.
This proceeds by applying forward kinematics using the root rotation from the
Lagrangian representation, and the x-z root velocities of the first few dims of
the processed matrix.

See also shortcuts:
    reconstruct_modelled(Y)
    reconstruct_raw(Y)
    reconstruct_root(Y)
"""
function reconstruct(Y::Matrix{T}, input_type::Symbol) where T <: Number
    Y = convert(Matrix{Float64}, Y)   # reduce error propagation from iterative scheme
    if input_type == :raw
        (size(Y, 2) < 70) && @warn "expecting matrix with >= 70 columns"
        root_r, root_x, root_z, joints = Y[:,1], Y[:,2], Y[:,3], Y[:,8:(63+7)]

    elseif input_type == :modelled64
        (size(Y, 2) != 64) && @warn "expecting matrix with exactly 64 columns"
        N = size(Y, 1)
        root_r, root_x, root_z, joints = Y[:,1], Y[:,2], Y[:,3], Y[:,5:end]
        rootjoint = reduce(hcat,  (zeros(T, N, 1), Y[:,4:4], zeros(T, N, 1)))
        joints = hcat(rootjoint, joints)

    elseif input_type == :modelled
        @assert size(Y, 2) == 67 "expecting matrix with exactly 67 columns. See `reconstruct_modelled64`"
        N = size(Y, 1)
        root_r, root_x, root_z = (Y[:,i] + Y[:,i+3] for i in 1:3)
        rootjoint = reduce(hcat,  (zeros(T, N, 1), Y[:,7:7], zeros(T, N, 1)))
        joints = hcat(rootjoint, Y[:,8:end])

    elseif input_type == :root
        (size(Y, 2) != 3) && @warn "expecting matrix with exactly 3 columns"
        root_r, root_x, root_z = Y[:,1], Y[:,2], Y[:,3]
        joints = zeros(size(Y,1), 3)   # (implicitly starting at the origin.)
    end

    return _joints_fk(joints, root_x, root_z, root_r)
end

function get_traj(Y::Matrix{T}, ixs::AbstractVector; offset_xz=zeros(T, 2),
        offset_r::T=T(0)) where T
    @assert size(Y, 2) == 67 "expecting matrix with exactly 67 columns."
    base_r, base_x, base_z = Y[:,1], Y[:,2], Y[:,3]

    if ixs[1] > 1
        root_r, root_x, root_z = (Y[:,i] + Y[:,i+3] for i in 1:3)

        # Calculate difference between base coods and root over i < ixs[1]
        begin_ix = 1:ixs[1]-1
        start_angle_base = -sum(base_r[begin_ix])
        start_angle_root = -sum(root_r[begin_ix])
        start_angle_δ = let θ=start_angle_base - start_angle_root + offset_r; [sin(θ), cos(θ)]; end

        start_base = _traj_fk(base_x[begin_ix], base_z[begin_ix], base_r[begin_ix],
            start=cat_zero_y(offset_xz...), start_angle=[sin(offset_r), cos(offset_r)])
        start_base = reduce(vcat, v[end] for v in start_base)
        start_root = _traj_fk(root_x[begin_ix], root_z[begin_ix], root_r[begin_ix]) |> x -> reduce(vcat, v[end] for v in x)
        start_δ = cat_zero_y((start_base - start_root + offset_xz)...);

        traj = _traj_fk(base_x[ixs], base_z[ixs], base_r[ixs], start=start_δ ,
                start_angle=start_angle_δ)
    else
        traj = _traj_fk(base_x[ixs], base_z[ixs], base_r[ixs], start=cat_zero_y(offset_xz...),
                start_angle=[sin(offset_r), cos(offset_r)])
    end
    return traj
end

function _traj_fk(root_x::Vector{T}, root_z::Vector{T}, root_r::Vector{T};
        start::Vector=zeros(T, 3), start_angle::Vector=[0,T(1)]) where T <: Number

    n = length(root_x)
    @argcheck norm(start_angle) ≈ 1
    rotation = _vec2_to_quat(start_angle)
    translation = start
    # translation = _qimag(_apply_rotation(quat(start), rotation))  # No. o.w. start pos != start.

    traj = Matrix{T}(undef, 3, n+1);
    traj[1, 1] = translation[1]
    traj[3, 1] = translation[3]

    for i = 1:n
        rotation = _quaterion_angle_axis_w_y(-root_r[i]) * rotation
        next_step = quat(cat_zero_y(root_x[i], root_z[i]))
        translation = translation + _qimag(_apply_rotation(next_step, rotation))
        traj[1, i+1] = translation[1]
        traj[3, i+1] = translation[3]
    end

    return traj[1,:], traj[3,:]
end

function _joints_fk(joints::Matrix{T}, root_x::Vector{T}, root_z::Vector{T},
        root_r::Vector{T}; start::Vector{T}=zeros(T, 3),
        start_angle::Vector{T}=[0,T(1)]) where T <: Number

    n = size(joints, 1)
    njoints = size(joints, 2) ÷ 3
    @assert (njoints * 3 == size(joints, 2)) "number of columns must be div. 3"
    joints = _rowmaj_reshape_3d(joints, n, njoints, 3)

    rotation = _vec2_to_quat(start_angle)
    translation = start

    for i = 1:n
        # Apply the rotation to the skeleton about the y axis
        # so that the azimuth at the **END** of the step is correct.
        joints[i,:,:] = _apply_rotation(_quat_list(joints[i,:,:]), rotation) |> _quat_list_to_mat

        # Move from origin --> calculated next step (see below for def.)
        # Note that this does not happen wrt current azimuth, and so is
        # **INDEPENDENT OF THE ABOVE ROTATION**.
        joints[i,:,1] = joints[i,:,1] .+ translation[1]
        joints[i,:,3] = joints[i,:,3] .+ translation[3]

        # Increment the cumulative y-axis rotation by next rotational difference.
        rotation = _quaterion_angle_axis_w_y(-root_r[i]) * rotation

        # Rotate the next delta (step; x,z diff) and express relative to origin.
        next_step = quat(cat_zero_y(root_x[i], root_z[i]))
        rotated_step = _qimag(_apply_rotation(next_step, rotation))
        translation += rotated_step
    end

    # note that vs the trajectory FK, we don't obtain an (n+1) output. This is
    # because while we know the position of the root at time (n+1), we don't
    # have the relative position of the joints, and so discard this final transformation.
    return joints
end

#= ---------------------------------------------------------------------------
        User facing reconstruction of ground path FK from data matrix
------------------------------------------------------------------------------=#
"""
    fk_path(Y, path_ixs, offset_xzr; for_viz=true, approx_to=nothing)

calculate the path trajectory from a given sequence `Y` (corresponding to an
**entire file**), the the time indices to extract from this sequence `path_ixs`.
The offset_xzr must be specified in order to obtain the initial differences in
position and rotation between the true and smoothed path trajectories. This function
outputs an N x 12 matrix for use in `mocapviz.create_animation` with appropriate
step-ahead entries, but `for_viz=false` will output just the N x 2 path.

To correct the rotation relative to a known trajectory, supply the N x 2 matrix in
`approx_to`, and the `correct_rotation` function will be used to align them.
"""
function fk_path(Y::AbstractMatrix, path_ixs::AbstractVector, offset_xzr::AbstractVector;
        for_viz::Bool=true, approx_to::Union{Nothing, AbstractMatrix}=nothing)
    N = length(path_ixs)
    T = promote_type(eltype(Y), eltype(offset_xzr))
    Y, offset_xzr = T.(Y), T.(offset_xzr)
    _traj = hcat(get_traj(Y, path_ixs, offset_xz=offset_xzr[1:2], offset_r=offset_xzr[3])...)

    if !(approx_to === nothing)
        _traj = _optimise_start(_traj, approx_to, vcat(zeros(T, 2), [randn(T, 2)*10 for i in 1:50]...))
    end
    !for_viz && return _traj

    path = randn(N, 12)*0.01
    for i in 1:6, j in 0:1
        path[:,i + j*6] = vcat(_traj[1+(5*i):N,1+j], _traj[N,1+j] .+ randn(5*i)*0.01)
    end
    return path
end

"""
    correct_rotation(A, B)
    correct_rotation(approx, target)

rotate path matrix A (n x 2) to closest match with target path
matrix B (n x 2) using orthogonal procrustes. This function
returns the corrected path matrix Â.
"""
function correct_rotation(A, B)
    @argcheck size(A, 2) == size(B, 2) == 2
    prod_svd = svd(B' * A)
    return A * prod_svd.V * prod_svd.U'
end

function _procrustes_matrix(A, B)
    prod_svd = svd(B' * A)
    return prod_svd.V * prod_svd.U'
end

function _optimise_start(A, B, xz::Vector)
    err = [sum(x->x^2, A[1:end-1,:] * _procrustes_matrix(A[1:end-1,:] .+ pos', B) - B) for pos in xz]
    return A * (_procrustes_matrix(A[1:end-1,:] .+ xz[argmin(err)]', B))
end
end

max_stars_repo_name: ornithos/pytorch-mtds-mocap

max_issues_repo_name: ornithos/pytorch-mtds-mocap

max_forks_repo_name: ornithos/pytorch-mtds-mocap


text:
using ModelingToolkit
using Base.Test

@IVar t
@DVar x(t) y(t) z(t)
@test_throws AssertionError @DVar w

@Deriv D'~t # Default of first derivative, Derivative(t,1)
@Param σ ρ β
@Const c=0

# Default values
p = Parameter(:p, 1)
u = DependentVariable(:u, [1], dependents = [t])

s = JumpVariable(:s,3,dependents=[t])
n = NoiseVariable(:n,dependents=[t])

σ*(y-x)
D*x
D*x ~ -σ*(y-x)
D*y ~ x*(ρ-z)-sin(y)

@test D*t == Constant(1)


text:
abstract type Side end
struct Left <: Side end
struct Right <: Side end

abstract type SourceDestination end
struct Src <: SourceDestination end
struct Dst <: SourceDestination end

keep(x::VecOrMat, y) = keep.(x, Ref(y))
keep(x::VecOrMat, y, z) = keep.(x, Ref(y), Ref(z))

keep(x::AbstractNode, ::Type{T}) where {T <: Node} = keep(x.node, T)
keep(x::AbstractNode, ::Type{T}) where {T <: AbstractNode} = x
keep(x::Node, ::Type{T}) where {T <: Node} = x
keep(x::Node, ::Type{T}) where {T <: AbstractNode} = x

keep(x::Symbol, y, z) = keep(convert(AbstractNode, x), y, z)
keep(x::Symbol, y) = keep(convert(AbstractNode, x), y)

keep(x::AbstractEdge, ::Type{T}) where {T <: Edge} = keep(x.edge, T)
keep(x::AbstractEdge, ::Type{T}) where {T <: AbstractEdge} = keep(x.edge, T)
keep(x::Edge, ::Type{T}) where {T <: Edge} = x
keep(x::Edge, ::Type{T}) where {T <: AbstractEdge} = x

keep(x::AbstractEdge, ::Type{T}) where {T <: Src} = x.src
keep(x::AbstractEdge, ::Type{T}) where {T <: Dst} = x.dst

keep(x::Arrow, ::Type{T}) where {T <: AbstractEdge} = x.edges
keep(x::Arrow, ::Type{T}) where {T <: Edge} = keep(x.edges, T)

keep(x::Arrow, ::Type{T}) where {T <: Left} = x.lhs
keep(x::Arrow, ::Type{T}) where {T <: Right} = x.rhs

function keep(x::AbstractNode, ::Type{T1}, ::Type{T2}) where {T1 <: Union{Side, SourceDestination}, T2 <: AbstractNode}
    keep(x, T2)
end

function keep(x::Arrow, ::Type{T1}, ::Type{T2}) where {T1 <: Union{Side, SourceDestination}, T2 <: AbstractNode}
    keep(keep(x, T1), T2)
end

keep(x::Arrow, ::Type{T}) where {T <: Src} = keep(x.edges, T)
keep(x::Arrow, ::Type{T}) where {T <: Dst} = keep(x.edges, T)

keep(x::Arrow, _) = x

unmeta(x::AbstractEdge) = keep(x, Edge)
unmeta(x::AbstractNode) = keep(x, Node)

unarrow(x) = keep(x, AbstractEdge)


text:
function __init__()
    # ensure that operations executed by the REPL back-end finish before returning,
    # because displaying values happens on a different task
    if isdefined(Base, :active_repl_backend)
        push!(Base.active_repl_backend.ast_transforms, synchronize_metal_tasks)
    end
end

function synchronize_metal_tasks(ex)
    quote
        try
            $(ex)
        finally
            if haskey($task_local_storage(), :MtlDevice)
                $device_synchronize()
            end
        end
    end
end


text:
#portfolio composition maps

#Plots
using Random
Random.seed!(111)
tickers = ["IBM", "Google", "Apple", "Intel"]
N = 10
D = length(tickers)
weights = rand(N, D)
weights ./= sum(weights, dims = 2)
returns = sort!((1:N) + D * randn(N))
portfoliocomposition(weights, returns, labels = permutedims(tickers))

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
# n-queens (nqueens) sowva, fow nsquawesx-by-nsquawesy boawd

stwuct Queen
    x::Int
    y::Int
end
hitshowz(queena, queenb) = queena.x == queenb.x
hitsvewt(queena, queenb) = queena.y == queenb.y
hitsdiag(queena, queenb) = abs(queena.x - queenb.x) == abs(queena.y - queenb.y)
hitshvd(qa, qb) = hitshowz(qa, qb) || hitsvewt(qa, qb) || hitsdiag(qa, qb)
hitsany(testqueen, queens) = any(q -> hitshvd(testqueen, q), queens)

function twysowve(nsquawesx, nsquawesy, nqueens, pwesqueens = ())
    nqueens == 0 && wetuwn pwesqueens
    fow xsquawe in 1:nsquawesx
        fow ysquawe in 1:nsquawesy
            testqueen = Queen(xsquawe, ysquawe)
            if !hitsany(testqueen, pwesqueens)
                twyqueens = (pwesqueens..., testqueen)
                maybesow = twysowve(nsquawesx, nsquawesy, nqueens - 1, twyqueens)
                maybesow !== nothing && wetuwn maybesow
            end
        end
    end
    wetuwn nothing
end

max_stars_repo_licenses: OML RSA-MD

max_issues_repo_licenses: OML RSA-MD

max_forks_repo_licenses: OML RSA-MD


text:
begin
    dir1_name = "test_readdirpath"
    mkdir(dir1_name)
    file1_path = joinpath(dir1_name, "file1.txt")
    touch(file1_path)
    dir2_name = joinpath(dir1_name, "dir2")
    mkdir(dir2_name)
    file2_path = joinpath(dir2_name, "file2.txt")
    touch(file2_path)

    @test readdirpath(dir2_name) == String[file2_path]
    @test readdirpath(dir1_name) == String[dir2_name, file1_path]

    rm(dir1_name,recursive=true)
end


text:
#This file will have 0 dependencies...
"""
    rbinomial( p, size... )

Makes an N-dimensional array of size(s) `size` with a probability of being a 1 over a 0 of 1 `p`.

Suggested by Baggepinnen on Discourse!
"""
rbinomial( p, size... ) = rand( size... ) .< p

"""
    EmpiricalQuantiles(X, quantiles)

Finds the column-wise `quantiles` of 2-Array `X` and returns them in a 2-Array of size `quantiles` by `variables`.
*Note: This copies the array... Use a subset if memory is the concern. *
"""
function EmpiricalQuantiles(X, quantiles)
    @assert all((quantiles .>= 0.0) .& (quantiles .< 1.0))
    ( obs, vars ) = size( X )
    quantilevalues = zeros( length( quantiles ), vars )
    for v in 1 : vars
        Potentials = unique(X[:, v])
        if length(Potentials) > 2
            quantilevalues[:, v] = Statistics.quantile( X[:,v], quantiles )
        else
            if length(Potentials) == 2
                quantilevalues[1:2, v] = Potentials[1 : 2]
            else
                quantilevalues[:, v] .= Potentials[1]
            end
        end
    end
    return quantilevalues
end

mutable struct RunningMean
    mu::Float64
    p::Int
end

"""
    RunningMean(x)

Constructs a running mean object with an initial scalar value of `x`.
"""
RunningMean(x) = RunningMean( x, 1 )
"""
    Update!(RM::RunningMean, x)

Adds new observation(`x`) to a RunningMean object(`RM`) in place.
"""
function Update!(RM::RunningMean, x)
    RM.p += 1
    RM.mu += (x - RM.mu) / RM.p
end

"""
    Update(RM::RunningMean, x)

Adds new observation(`x`) to a RunningMean object(`RM`) and returns the new object.
"""
Update(RM::RunningMean, x) = RunningMean( RM.mu + (x - RM.mu) / ( RM.p + 1 ), RM.p + 1 )

"""
    Remove!(RM::RunningMean, x)

Removes an observation(`x`) from a RunningMean object(`RM`) and reculates the mean in place.
"""
function Remove!(RM::RunningMean, x)
    RM.mu = (RM.p * RM.mu - x) / (RM.p - 1)
    RM.p -= 1
end

"""
    Remove!(RM::RunningMean, x)

Removes an observation(`x`) from a RunningMean object(`RM`) and recuturns the new RunningMean object.
"""
Remove(RM::RunningMean, x) = RunningMean( (RM.p * RM.mu - x) / (RM.p - 1), RM.p - 1 )

mutable struct RunningVar
    m::RunningMean
    v::Float64
end

"""
    RunningVar(x)

Constructs a RunningVar object with an initial scalar value of `x`.
*Note: RunningVar objects implicitly calculate the running mean.*
"""
RunningVar(x) = RunningVar( RunningMean( x, 1 ), 0.0 )

"""
    Update!(RV::RunningVar, x)

Adds new observation(`x`) to a RunningVar object(`RV`) and updates it in place.
"""
function Update!(RV::RunningVar, x)
    OldMean = copy(RV.m.mu)
    Update!(RV.m, x)
    RV.v = ( (RV.v * (RV.m.p - 2)) + ( (x - OldMean) * ( x - RV.m.mu ) ) ) / (RV.m.p - 1.0)
end

"""
    Variance(rv::RunningVar)

Returns the current variance inside of a RunningVar object.
"""
Variance(rv::RunningVar) = rv.v

"""
    Mean(rv::RunningVar)

Returns the current mean inside of a RunningVar object.
"""
Mean(rv::RunningVar) = rv.m.mu

"""
    Mean(rv::RunningMean)

Returns the current mean inside of a RunningMean object.
"""
Mean(rm::RunningMean) = rm.mu

"""
    Skewness(X)

returns a measure of skewness for a population vector `X`.

Joanes, D. N., and C. A. Gill. 1998. “Comparing Measures of Sample Skewness and Kurtosis”. The Statistician 47(1): 183–189.
"""
Skewness(X) = (sum( (X .- Statistics.mean(X)) .^ 3) / length(X)) / (Statistics.var(X) ^ (1.5))

"""
    SampleSkewness(X)

returns a measure of skewness for vector `X` that is corrected for a sample of the population.

Joanes, D. N., and C. A. Gill. 1998. “Comparing Measures of Sample Skewness and Kurtosis”. The Statistician 47(1): 183–189.
"""
function SampleSkewness(X)
    N = length(X)
    @assert N > 2
    return ( sqrt( N * (N - 1) ) / (N - 2) ) * Skewness( X )
end

"""
    CorrelationMatrix(X; DOF_used = 0)

Returns the Pearson correlation matrix from a centered covariance matrix.

This is only included because finding a legible implementation was hard for me
to find some years ago (for the reader). But, also I don't like assumptions on
whether or not we should use all N, N-1, etc for scaling (hence `DOF_used`).
"""
function CorrelationMatrix(X; DOF_used = 0)
    obs, _ = size(X)
    obs -= DOF_used
    C = LinearAlgebra.I - ( ( 1 / obs ) .* ones(obs) * ones(obs)' )
    Xs = C * X #Same thing as subtracting the column means
    XstXs = (1 / obs) * (Xs' * Xs) #get covariance matrix
    D = LinearAlgebra.Diagonal( 1 ./ sqrt.( LinearAlgebra.diag( XstXs ) ) )#define the scaling matrix
    return D * XstXs * D
end

"""
    CorrelationVectors( A, B )

Returns the Pearson correlation of 2 vectors.

This is only included because finding a legible implementation was hard for me
to find some years ago (for the reader).
"""
function CorrelationVectors( A, B )
    obs = length( A )
    @assert( obs == length( B ), "Vectors must have the same length." )
    A = ( A .- mean(A) )
    B = ( B .- mean(B) )
    return ( A' * B ) * ( 1 / ( (obs - 1) * std( A ) * std( B )) )
end

#This was written for an algorithm and didn't fit in anywhere so for now it's kept
#but it may not have use...
struct PermutedVectorPair{A,B,C}
    vec1::A
    vec2::B
    operation::C
    i::Int
    length::Int
end

"""
    PermutedVectorPair(vec1, vec2; op = +)

Returns an iterator which applies each element in vec2 to vec1 via the user selected operator(op)
"""
function PermutedVectorPair(vec1, vec2; op = +)
    return PermutedVectorPair(vec1, vec2, op, 1, length(vec2))
end

function Base.iterate(iterator::PermutedVectorPair, state = 1)
    if state > iterator.length
        return nothing
    else
        return ( broadcast(iterator.operation, iterator.vec1, iterator.vec2[state]) , state + 1)
    end
end


text:
function PMM_settings(D_ul::Int, D_ll::Int, fnum::Int)

    settings = zeros(Int32, 5)
    ccall((:PMM_config, bilevelBenchmark),
          Cvoid,
        (Int32,Int32,Ptr{Int32}, Int32),
        D_ul, D_ll,settings, fnum)

    _, _,m,lenG, leng = settings

    return m,lenG, leng
end

function PMM_Ψ(x::Array{Float64}, D_ll::Int, fnum::Int)

    D_ul = length(x)

    m,lenG, leng = PMM_settings(D_ul,D_ll,fnum)

    y = zeros(D_ll)

    ccall((:PMM_Psi, bilevelBenchmark),
          Cvoid,
        (Int32, Int32,Int32, Ptr{Cdouble}, Ptr{Cdouble}, Int32),
        D_ul, D_ll, m,       x,           y,        fnum)

    return y
end

function PMM_leader(x::Array{Float64}, y::Array{Float64}, fnum::Int)
    if fnum > 6
        @error "PMM1 to PMM6 are alowed."
    end
    F = [0.0]

    D_ul = length(x)
    D_ll = length(y)

    m,lenG, leng = PMM_settings(D_ul,D_ll,fnum)
    G = zeros(lenG)

    ccall((:PMM_leader, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Int32),
        D_ul, D_ll, x, y, F, G, fnum)

    if lenG == 0
        return F[1]
    end

    return F[1], G

end

function PMM_follower(x::Array{Float64}, y::Array{Float64}, fnum::Int)
    if fnum > 6
        @error "PMM1 to PMM6 are alowed."
    end
    f = [0.0]

    D_ul = length(x)
    D_ll = length(y)

    m,lenG, leng = PMM_settings(D_ul,D_ll,fnum)
    g = zeros(leng)

    ccall((:PMM_follower, bilevelBenchmark),
          Cvoid,
        (Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Int32),
        D_ul, D_ll, x, y, f, g, fnum)

    if leng == 0
        return f[1]
    end

    return f[1], g
end

function PMM_ranges(D_ul, D_ll, fnum=0)
    bounds_ul = Array([-10ones(D_ul) 10ones(D_ul)]')
    bounds_ll = Array([-10ones(D_ll) 10ones(D_ll)]')
    bounds_ul, bounds_ll
end

function PMM_test(D_ul, D_ll, fnum)
    x = -10 .+ 20*rand(D_ul)
    y = PMM_Ψ(x, D_ll, fnum)#rand(10)
    if fnum <= 6
        F = PMM_leader(x, y, fnum)
        f = PMM_follower(x, y, fnum)
    else
        F, _ = PMM_leader(x, y, fnum)
        f, _ = PMM_follower(x, y, fnum)
    end

    return F, f
end

function PMM_optimum(fnum)
    return 0.0, 0.0

end


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

import Base.StackTraces: lookup

# Test location information for inlined code (ref issues #1334 #12544)
module test_inline_bt
using Test
import ..lookup

function get_bt_frames(functionname, bt)
    for i = 1:length(bt)
        lkup = lookup(bt[i])
        lkup[end].func == functionname && return lkup
    end
    return StackTraces.StackFrame[]
end

# same-file inline
eval(Expr(:function, Expr(:call, :test_inline_1),
          Expr(:block, Expr(:line, 99, Symbol("backtrace.jl")),
                       Expr(:block, Expr(:line, 42),
                                    Expr(:meta, :push_loc, Symbol("backtrace.jl"), :inlfunc),
                                    Expr(:line, 37),
                                    Expr(:call, :throw, "foo"),
                                    Expr(:meta, :pop_loc),
                                    Expr(:line, 99)))))

@test functionloc(test_inline_1) == ("backtrace.jl", 99)
try
    test_inline_1()
    error("unexpected")
catch err
    lkup = get_bt_frames(:test_inline_1, catch_backtrace())
    @test length(lkup) == 2
    @test endswith(string(lkup[2].file), "backtrace.jl")
    @test lkup[2].line == 42
    # TODO: we don't support surface AST locations with inlined function names
    @test_broken lkup[1].func == :inlfunc
    @test endswith(string(lkup[1].file), "backtrace.jl")
    @test lkup[1].line == 37
end

# different-file inline
const absfilepath = Sys.iswindows() ? "C:\\foo\\bar\\baz.jl" : "/foo/bar/baz.jl"
eval(Expr(:function, Expr(:call, :test_inline_2),
          Expr(:block, Expr(:line, 81, Symbol("backtrace.jl")),
                       Expr(:block, Expr(:meta, :push_loc, Symbol(absfilepath)),
                                    Expr(:line, 111),
                                    Expr(:call, :throw, "foo"),
                                    Expr(:meta, :pop_loc),
                                    Expr(:line, 99)))))

@test functionloc(test_inline_2) == ("backtrace.jl", 81)
try
    test_inline_2()
    error("unexpected")
catch err
    lkup = get_bt_frames(:test_inline_2, catch_backtrace())
    @test length(lkup) == 2
    @test endswith(string(lkup[2].file), "backtrace.jl")
    @test lkup[2].line == 81
    @test string(lkup[1].file) == absfilepath
    @test lkup[1].line == 111
end

end # module

#issue 12977: line numbers for kwarg methods.
linenum = @__LINE__; f12977(; args...) = ()
loc = functionloc(f12977)
@test endswith(loc[1], "backtrace.jl")
@test loc[2] == linenum

@noinline function test_throw_commoning(x)
    if x==1; throw(AssertionError()); end
    if x==2; throw(AssertionError()); end
end

let
    local b1, b2
    try
        test_throw_commoning(1)
    catch
        b1 = stacktrace(catch_backtrace())
    end
    try
        test_throw_commoning(2)
    catch
        b2 = stacktrace(catch_backtrace())
    end
    i1 = findfirst(frame -> frame.func === :test_throw_commoning, b1)
    i2 = findfirst(frame -> frame.func === :test_throw_commoning, b2)
    @test i1 > 0 && i2 > 0
    @test b1[i1].line != b2[i2].line
end

module BackTraceTesting
using Test
import ..lookup

@inline bt2() = backtrace()
@inline bt1() = bt2()
bt() = bt1()

lkup = map(lookup, bt())
hasbt = hasbt2 = false
for sfs in lkup
    for sf in sfs
        if sf.func == :bt
            global hasbt = true
        end
        if sf.func == :bt2
            global hasbt2 = true
        end
    end
end
@test hasbt
@test hasbt2

function btmacro()
    ret = @timed backtrace()
    ret[1]
end
lkup = map(lookup, btmacro())
hasme = hasbtmacro = false
for sfs in lkup
    for sf in sfs
        if sf.func == Symbol("macro expansion")
            global hasme = true
        end
        if sf.func == :btmacro
            global hasbtmacro = true
        end
    end
end
@test hasme
@test hasbtmacro

end

# Interpreter backtraces
bt = eval(quote
    try
        error()
    catch
        catch_backtrace()
    end
end)
lkup = map(lookup, bt)
hastoplevel = false
for sfs in lkup
    for sf in sfs
        if sf.linfo isa Core.CodeInfo
            global hastoplevel = true
        end
    end
end
@test hastoplevel

# issue #23971
let
    for i = 1:1
        global bt23971 = backtrace()
    end
end
let st = stacktrace(bt23971)
    @test StackTraces.is_top_level_frame(st[1])
    @test string(st[1].file) == @__FILE__
    @test !occursin("missing", string(st[2].file))
end

# issue #27959
let bt, found = false
    @testset begin
        bt = backtrace()
    end
    for frame in map(lookup, bt)
        if frame[1].line == @__LINE__() - 3 && frame[1].file == Symbol(@__FILE__)
            found = true; break
        end
    end
    @test found
end

# issue 28618
let bt, found = false
    @info ""
    bt = backtrace()
    for frame in map(lookup, bt)
        if frame[1].line == @__LINE__() - 2 && frame[1].file == Symbol(@__FILE__)
            found = true; break
        end
    end
    @test found
end

# Syntax error locations appear in backtraces
let trace = try
        include_string(@__MODULE__,
            """

            )

            """, "a_filename")
    catch
        stacktrace(catch_backtrace())
    end
    @test trace[1].func == Symbol("top-level scope")
    @test trace[1].file == :a_filename
    @test trace[1].line == 2
end
let trace = try
        include_string(@__MODULE__,
            """

            incomplete_syntax(

            """, "a_filename")
    catch
        stacktrace(catch_backtrace())
    end
    @test trace[1].func == Symbol("top-level scope")
    @test trace[1].file == :a_filename
    @test trace[1].line == 2
end

# issue #29695 (see also test for #28442)
let code = """
    f29695(c) = g29695(c)
    g29695(c) = c >= 1000 ? (return backtrace()) : f29695(c + 1)
    bt = f29695(1)
    meth_names = [ip.code.def.name for ip in bt
                  if ip isa Base.InterpreterIP && ip.code isa Core.MethodInstance]
    num_fs = sum(meth_names .== :f29695)
    num_gs = sum(meth_names .== :g29695)
    print(num_fs, ' ', num_gs)
    """

    @test read(`$(Base.julia_cmd()) --startup-file=no --compile=min -e $code`, String) == "1000 1000"
end

# Test that modules make it into InterpreterIP for top-level code
let code = """
    module A
    foo() = error("Expected")
    try
        foo()
    catch
        global bt = catch_backtrace()
    end
    end

    foreach(println, A.bt)
    """

    bt_str = read(`$(Base.julia_cmd()) --startup-file=no --compile=min -e $code`, String)
    @test occursin("InterpreterIP in MethodInstance for foo", bt_str)
    @test occursin("InterpreterIP in top-level CodeInfo for Main.A", bt_str)
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
using JFVM, PyPlot, JPhreeqc

# Physical properties
n_comp = 6 # number of components
Di = [1.0e-19, 2.0e-9, 3.0e-9] # m^2/s
u_inj = 1.0/(24*3600) # m/s
c_init = [0.0, 0.5, 0.1] # mol/m^3
phi_val = 0.4 # [-] porosity
k_val = 0.001e-12 # m^2 permeability
c_bound = [1.0, 0.0, 2.0] # mol/m^3 left boundary

# Define the domain
Lx = 1.0 # m
Nx = 100 # [-]
nxyz = Nx
m = createMesh1D(Nx, Lx)
v_cell = cellVolume(m)
sw = createCellVariable(m, 1.0)
phi = createCellVariable(m, phi_val)

# m = createMesh2D(Nx, 10, Lx, Lx/10)

nthreads=1 # number of threads
id=RM_Create(nxyz, nthreads)
setDefaultPhreeqcProperties(id) # set properties
setDefaultPhreeqcUnits(id)
status = RM_SetUnitsSurface(id, 0)       # 0, mol/L cell 1, mol/L water 2 mol/L rock
status = RM_SetUnitsSolution(id, 2)
status = RM_UseSolutionDensityVolume(id, 2)
status = RM_SetSelectedOutputOn(id, 1)       # enable selected output
grid2chem = collect(Int, 0:nxyz-1)
status = RM_CreateMapping(id, grid2chem)
status = RM_SetPrintChemistryOn(id, 0, 1, 0) # workers, initial_phreeqc, utility

RM_LoadDatabase(id, "PHREEQC.DAT")
vol=internalCells(v_cell)[:]*1000 # cel volume in liter
RM_SetRepresentativeVolume(id, vol)
por=internalCells(phi)[:]
RM_SetPorosity(id, por)
# Set initial saturation
sat = internalCells(sw)[:]
status = RM_SetSaturation(id, sat)
RM_SetDensity(id, ones(nxyz))
RM_RunFile(id, 1,1,1, "phreeqc_input.pqi")

RM_FindComponents(id)

ic1 = zeros(Int, nxyz, 7)
ic2 = zeros(Int, nxyz, 7)
f11 = zeros(Float64, nxyz, 7)
ic1[:,1] = 1:nxyz       # Solution 1
ic1[:,2] = -1      # Equilibrium phases none
ic1[:,3] = -1       # Exchange none
ic1[:,4] = -1      # Surface 1
ic1[:,5] = -1      # Gas phase none
ic1[:,6] = -1      # Solid solutions none
ic1[:,7] = -1      # Kinetics none
ic2[:,1] = -1      # Solution none
ic2[:,2] = -1      # Equilibrium phases none
ic2[:,3] = -1      # Exchange none
ic2[:,4] = -1      # Surface none
ic2[:,5] = -1      # Gas phase none
ic2[:,6] = -1      # Solid solutions none
ic2[:,7] = -1      # Kinetics none
f11[:,1]  = 1.0      # Mixing fraction ic1 Solution
f11[:,2] = 1.0      # Mixing fraction ic1 Equilibrium phases
f11[:,3] = 1.0      # Mixing fraction ic1 Exchange 1
f11[:,4] = 1.0      # Mixing fraction ic1 Surface
f11[:,5] = 1.0      # Mixing fraction ic1 Gas phase
f11[:,6] = 1.0      # Mixing fraction ic1 Solid solutions
f11[:,7] = 1.0      # Mixing fraction ic1 Kinetics
status = RM_InitialPhreeqc2Module(id, ic1[:], ic2[:], f11[:])

time1 = 0.0 # changed time to time1 to avoid conflicts with julia time function
time_step = 0.0
status = RM_SetTime(id, time1)
status = RM_SetTimeStep(id, time_step)
# get elements
# status = RM_SetSpeciesSaveOn(id, 0)
status = RM_RunCells(id)
element_list= getComponentList(id) # include charge
n_element= RM_GetComponentCount(id) # include charge
c = zeros(Float64, n_element*nxyz)
status = RM_GetConcentrations(id, c) # get the aqueous concentrations of elements
c=reshape(c, nxyz, n_element)

c_left=zeros(Float64, n_element)
n_bound=1
status= RM_InitialPhreeqc2Concentrations(id, c_left, n_bound, zeros(Int, n_bound),
  zeros(Int, n_bound), ones(Float64, n_bound)) # solution 0

# boundary conditions
bc = Array{BoundaryCondition}(n_comp)
for i in 1:n_comp
  bc[i] = createBC(m)
  bc[i].left.a[:] = 0.0
  bc[i].left.b[:] = 1.0
  bc[i].left.c[:] = c_left[i]
end

# Initial conditions
c_old = Array{CellValue}(n_comp)
c_new = Array{CellValue}(n_comp)
for i in 1:n_comp
  c_old[i] = createCellVariable(m, c[:,i][:], bc[i])
end

# velocity and diffusion coefficients
# D = Array{FaceValue}(n_comp)
u = createFaceVariable(m, u_inj)
D = createFaceVariable(m, Di[1])

# solver options
t_final = 0.5*Lx/(u_inj/phi_val)
dt = Lx/(u_inj/phi_val)/Nx

# Discretization
# M_diff = Array{Any}(3)
M_conv = convectionTerm(u)
M_diff = diffusionTerm(D)
M_bc = Array{Any}(n_comp)
RHS_bc = Array{Any}(n_comp)
for i in 1:n_comp
  M_bc[i], RHS_bc[i] = boundaryConditionTerm(bc[i])
end

M_t = Array{Any}(n_comp)
RHS_t = Array{Any}(n_comp)
c_reac= zeros(Float64, nxyz, n_element)

for t in dt:dt:t_final
  for i in 1:n_comp
    M_t[i], RHS_t[i] = transientTerm(c_old[i], dt, phi_val)
    c_new[i] = solveLinearPDE(m, M_t[i]+M_bc[i]+M_conv-phi_val*M_diff, RHS_bc[i]+RHS_t[i])
    c_reac[:,i]=internalCells(c_new[i])[:]
  end

  status= RM_SetConcentrations(id, c_reac[:])
  sat = internalCells(sw)[:]
  status = RM_SetSaturation(id, sat)
  status = RM_SetTime(id, t)
  status = RM_SetTimeStep(id, dt)
  status = RM_RunCells(id)
  c = zeros(Float64, n_element*nxyz)
  status = RM_GetConcentrations(id, c)

  c=reshape(c, nxyz, n_element)
  # update the concentrations
  # conc_init=Array{Any}(n_element+1)
  for i in 1:n_element
    c_old[i]=createCellVariable(m, c[:,i][:], bc[i])
  end
  println("progress: $(t/t_final*100) %")
end
figure(4)
visualizeCells(c_old[4])
visualizeCells(c_old[5])
visualizeCells(c_old[6])
savefig("profiles1.png")

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using Logging
using Gadfly
using JLSO
using Revise
includet("./utils.jl")
using .PlottingUtils: get_plot_data, sr_plots,
                        box_error_plots, oe_plots

"""
Load serialized results and save the
most relevant plots
"""
function plot_main_results(results::Dict,
                           path::String,
                           methods::Vector=[:SS,:DSS])::Nothing

    first_key = collect(keys(results))[1]
    first_result = results[first_key]
    samples = collect(keys(first_result))
    noise_levels = Array{Symbol}(collect(keys(first_result[samples[1]])))
    #methods = Array{Symbol}(collect(keys(first_result[samples[1]][noise_levels[1]])))

    # cur_colors = [colorant"#e6612f",colorant"#42b3d5",colorant"seagreen"]
    # [Flame, Maximum Blue, Seagreen]
    cur_colors = [colorant"#ff595e",colorant"#ffca3a",colorant"#8ac926",colorant"#1982c4"]
    # [Sizzling Red, Sunglow, Yellow Green, Green Blue Crayola]

    method_label = Dict()
    for m in methods
        method_label[m] = String(m)
    end

    method_color = Dict()
    for (m,c) in zip([:SS,:SSR,:DS,:DSS],cur_colors)
        method_color[m] = c
    end

    for (problem,result) in results
        nrmse_log = "$(problem)\n"
        sr_log = "$(problem)\n"
        oe_log = "$(problem)\n"

        dir_path = joinpath(path,string(problem))
        mkdir(dir_path)

        for sam in samples
            nrmse_log *= "Sample size:\t$(sam)\n"
            sr_log *= "Sample size:\t$(sam)\n"
            oe_log *= "Sample size:\t$(sam)\n"

            res = result[sam]
            plot_data = get_plot_data(res,noise_levels,methods)

            sr_log *= sr_plots(plot_data,noise_levels,methods,
                    method_label,method_color,sam,dir_path)

            oe_log *= oe_plots(plot_data,noise_levels,methods,
                    method_label,method_color,sam,dir_path)

            if :DS in methods
                oe_log *= oe_plots(plot_data,
                                   noise_levels,filter(x -> x != :DS, methods),
                                   method_label,method_color,sam,dir_path)
            end

            for noise in noise_levels
                plot_data = get_plot_data(res,[noise],methods)

                nrmse_log *= box_error_plots(plot_data,noise,methods,
                                  method_label,method_color,sam,dir_path)

                #parameter_plots(plot_data,noise,methods,
                #                   method_label,method_color,sam,dir_path)

                sr_log *= sr_plots(plot_data,noise,methods,
                        method_label,method_color,sam,dir_path)

                oe_log *= oe_plots(plot_data,[noise],methods,
                        method_label,method_color,sam,dir_path)

            end
        end
        @info nrmse_log
        @info sr_log
        @info oe_log
    end
    nothing
end

function main(args::Array{<:String})::Nothing
    jlso_file_or_path = string(args[1])
    save_path = string(args[2])

    println(jlso_file_or_path)
    io = open(joinpath(save_path, "plots.log"), "a+")
    logger = SimpleLogger(io)
    global_logger(logger)

    if isdir(jlso_file_or_path)
        filter_jlso = files -> filter(f -> occursin(".jlso", f), files)

        prepend_full_path = file -> joinpath(jlso_file_or_path, file)
        generate_jlso_paths = files -> map(prepend_full_path, files)

        get_jlso_files(path) =
            path |> readdir |> filter_jlso |> generate_jlso_paths

        jlso_files = get_jlso_files(jlso_file_or_path)
    else
        jlso_files = [jlso_file_or_path]
    end

    # base_dir = "/home/andrew/git/DEParamEst/"
    # save_path = joinpath(base_dir,"results/remote/results_low/")
    # jlso_file = joinpath(save_path,"experiment_results.jlso")
    results = map(JLSO.load, jlso_files)

    for r in results
        plot_main_results(r, save_path, [:SS, :SSR, :DS, :DSS])
    end
    close(io)
    nothing
end

main(ARGS)

max_stars_repo_name: atgmello/DEParamEst

max_issues_repo_name: atgmello/DEParamEst

max_forks_repo_name: atgmello/DEParamEst


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

include("../perfutil.jl")

include("indexing.jl")

briefname(A) = typeof(A).name.name

# Small array tests
sz = (3,5)
Alist = makearrays(Int, sz)
for Ar in Alist
    @timeit sumelt(Ar, 10^5) string("sumeltIs ", briefname(Ar)) string("for a in A indexing, ", briefname(Ar)) sz
    @timeit sumeach(Ar, 10^5) string("sumeachIs ", briefname(Ar)) string("for I in eachindex(A), ", briefname(Ar)) sz
    @timeit sumlinear(Ar, 10^5) string("sumlinearIs ", briefname(Ar)) string("for I in 1:length(A), ", briefname(Ar)) sz
    @timeit sumcartesian(Ar, 10^5) string("sumcartesianIs ", briefname(Ar)) string("for I in CartesianIndices(size(A)), ", briefname(Ar)) sz
    @timeit sumcolon(Ar, 10^5) string("sumcolonIs ", briefname(Ar)) string("colon indexing, ", briefname(Ar)) sz
    @timeit sumrange(Ar, 10^5) string("sumrangeIs ", briefname(Ar)) string("range indexing, ", briefname(Ar)) sz
    @timeit sumlogical(Ar, 10^5) string("sumlogicalIs ", briefname(Ar)) string("logical indexing, ", briefname(Ar)) sz
    @timeit sumvector(Ar, 10^5) string("sumvectorIs ", briefname(Ar)) string("vector indexing, ", briefname(Ar)) sz
end

Alist = makearrays(Float32, sz)   # SIMD-able
for Ar in Alist
    @timeit sumelt(Ar, 10^5) string("sumeltFs ", briefname(Ar)) string("for a in A indexing, ", briefname(Ar)) sz
    @timeit sumeach(Ar, 10^5) string("sumeachFs ", briefname(Ar)) string("for I in eachindex(A), ", briefname(Ar)) sz
    @timeit sumlinear(Ar, 10^5) string("sumlinearFs ", briefname(Ar)) string("for I in 1:length(A), ", briefname(Ar)) sz
    @timeit sumcartesian(Ar, 10^5) string("sumcartesianFs ", briefname(Ar)) string("for I in CartesianIndices(size(A)), ", briefname(Ar)) sz
    @timeit sumcolon(Ar, 10^5) string("sumcolonFs ", briefname(Ar)) string("colon indexing, ", briefname(Ar)) sz
    @timeit sumrange(Ar, 10^5) string("sumrangeFs ", briefname(Ar)) string("range indexing, ", briefname(Ar)) sz
    @timeit sumlogical(Ar, 10^5) string("sumlogicalFs ", briefname(Ar)) string("logical indexing, ", briefname(Ar)) sz
    @timeit sumvector(Ar, 10^5) string("sumvectorFs ", briefname(Ar)) string("vector indexing, ", briefname(Ar)) sz
end

# Big array tests
sz = (300,500)
Alist = makearrays(Int, sz)
for Ar in Alist
    @timeit sumelt(Ar, 100) string("sumeltIb ", briefname(Ar)) string("for a in A indexing, ", briefname(Ar)) sz
    @timeit sumeach(Ar, 100) string("sumeachIb ", briefname(Ar)) string("for I in eachindex(A), ", briefname(Ar)) sz
    @timeit sumlinear(Ar, 100) string("sumlinearIb ", briefname(Ar)) string("for I in 1:length(A), ", briefname(Ar)) sz
    @timeit sumcartesian(Ar, 100) string("sumcartesianIb ", briefname(Ar)) string("for I in CartesianIndices(size(A)), ", briefname(Ar)) sz
    @timeit sumcolon(Ar, 100) string("sumcolonIb ", briefname(Ar)) string("colon indexing, ", briefname(Ar)) sz
    @timeit sumrange(Ar, 100) string("sumrangeIb ", briefname(Ar)) string("range indexing, ", briefname(Ar)) sz
    @timeit sumlogical(Ar, 100) string("sumlogicalIb ", briefname(Ar)) string("logical indexing, ", briefname(Ar)) sz
    @timeit sumvector(Ar, 100) string("sumvectorIb ", briefname(Ar)) string("vector indexing, ", briefname(Ar)) sz
end

Alist = makearrays(Float32, sz)   # SIMD-able
for Ar in Alist
    @timeit sumelt(Ar, 100) string("sumeltFb ", briefname(Ar)) string("for a in A indexing, ", briefname(Ar)) sz
    @timeit sumeach(Ar, 100) string("sumeachFb ", briefname(Ar)) string("for I in eachindex(A), ", briefname(Ar)) sz
    @timeit sumlinear(Ar, 100) string("sumlinearFb ", briefname(Ar)) string("for I in 1:length(A), ", briefname(Ar)) sz
    @timeit sumcartesian(Ar, 100) string("sumcartesianFb ", briefname(Ar)) string("for I in CartesianIndices(size(A)), ", briefname(Ar)) sz
    @timeit sumcolon(Ar, 100) string("sumcolonFb ", briefname(Ar)) string("colon indexing, ", briefname(Ar)) sz
    @timeit sumrange(Ar, 100) string("sumrangeFb ", briefname(Ar)) string("range indexing, ", briefname(Ar)) sz
    @timeit sumlogical(Ar, 100) string("sumlogicalFb ", briefname(Ar)) string("logical indexing, ", briefname(Ar)) sz
    @timeit sumvector(Ar, 100) string("sumvectorFb ", briefname(Ar)) string("vector indexing, ", briefname(Ar)) sz
end

max_stars_repo_name: JuliaTagBot/Benchmarks

max_issues_repo_name: JuliaArchive/Benchmarks

max_forks_repo_name: JuliaArchive/Benchmarks


text:
module SpringCollab2020CustomRespawnTimeRefill

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/CustomRespawnTimeRefill" CustomRespawnTimeRefill(x::Integer, y::Integer, twoDash::Bool=false, respawnTime::Number=2.5)

const placements = Ahorn.PlacementDict(
    "Refill (Custom Respawn Time) (Spring Collab 2020)" => Ahorn.EntityPlacement(
        CustomRespawnTimeRefill
    ),

    "Refill (Two Dashes, Custom Respawn Time) (Spring Collab 2020)" => Ahorn.EntityPlacement(
        CustomRespawnTimeRefill,
        "point",
        Dict{String, Any}(
            "twoDash" => true
        )
    )
)

spriteOneDash = "objects/refill/idle00"
spriteTwoDash = "objects/refillTwo/idle00"

function getSprite(entity::CustomRespawnTimeRefill)
    twoDash = get(entity.data, "twoDash", false)

    return twoDash ? spriteTwoDash : spriteOneDash
end

function Ahorn.selection(entity::CustomRespawnTimeRefill)
    x, y = Ahorn.position(entity)
    sprite = getSprite(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CustomRespawnTimeRefill, room::Maple.Room)
    sprite = getSprite(entity)
    Ahorn.drawSprite(ctx, sprite, 0, 0)
end

end


text:
struct OptimizerStats
    simulation_step::Int
    execution_index::Int
    objective_value::Float64
    termination_status::Int
    primal_status::Int
    dual_status::Int
    solver_solve_time::Float64
    timed_solve_time::Float64
    solve_bytes_alloc::Union{Nothing, Float64}
    sec_in_gc::Union{Nothing, Float64}
end

_SIMULATION_FIELDS = Set((:simulation_step, :execution_index))
_BASE_FIELDS = setdiff(fieldnames(OptimizerStats), _SIMULATION_FIELDS)

"""
Construct OptimizerStats when the OperationsProblem is part of a simulation.
"""
function OptimizerStats(problem, simulation_step)
    timed_log = get_solve_timed_log(problem)
    execution_index = get_execution_count(problem)
    model = get_jump_model(problem)

    return OptimizerStats(
        simulation_step,
        execution_index,
        JuMP.objective_value(model),
        Int(JuMP.termination_status(model)),
        Int(JuMP.primal_status(model)),
        Int(JuMP.dual_status(model)),
        _get_solver_time(model),
        timed_log[:timed_solve_time],
        timed_log[:solve_bytes_alloc],
        timed_log[:sec_in_gc],
    )
end

"""
Construct OptimizerStats when the OperationsProblem is not part of a simulation.
"""
function OptimizerStats(problem)
    timed_log = get_solve_timed_log(problem)
    model = get_jump_model(problem)

    return OptimizerStats(
        0,
        0,
        JuMP.objective_value(model),
        Int(JuMP.termination_status(model)),
        Int(JuMP.primal_status(model)),
        Int(JuMP.dual_status(model)),
        _get_solver_time(model),
        timed_log[:timed_solve_time],
        timed_log[:solve_bytes_alloc],
        timed_log[:sec_in_gc],
    )
end

"""
Construct OptimizerStats from a vector that was serialized to HDF5.
"""
function OptimizerStats(data::Vector{Float64})
    return OptimizerStats(data...)
end

"""
Convert OptimizerStats to an array of floats that can be serialized to HDF5.
"""
function to_array(stats::OptimizerStats)
    return [Float64(getfield(stats, x)) for x in fieldnames(OptimizerStats)]
end

function to_dataframe(stats::OptimizerStats)
    df = DataFrames.DataFrame([to_namedtuple(stats)])
    if !_part_of_simulation(stats)
        DataFrames.select!(df, _BASE_FIELDS)
    end

    return df
end

function to_dict(stats::OptimizerStats)
    data = Dict()
    for field in fieldnames(typeof(stats))
        !_part_of_simulation(stats) && field in _SIMULATION_FIELDS && continue
        data[String(field)] = getfield(stats, field)
    end

    return data
end

function _get_solver_time(model)
    solver_solve_time = NaN
    try
        solver_solve_time = MOI.get(model, MOI.SolveTime())
    catch
        @warn "SolveTime() property not supported by the Solver"
    end

    return solver_solve_time
end

function get_column_names(::Type{OptimizerStats})
    return collect(string.(fieldnames(OptimizerStats)))
end

_part_of_simulation(stats::OptimizerStats) = stats.simulation_step != 0


text:
module GeometryPredicates
export orient_fast, orient_exact, orient, incircle_fast, incircle_exact, incircle

include("adaptive_arithmetic.jl")

const EPS_ORIENT_A_FLOAT64 = (3.0 + 16.0 * eps(Float64)) * eps(Float64)
const EPS_ORIENT_B_FLOAT64 = (2.0 + 12.0 * eps(Float64)) * eps(Float64)
const EPS_ORIENT_C1_FLOAT64 = (3.0 + 8.0 * eps(Float64)) * eps(Float64)
const EPS_ORIENT_C2_FLOAT64 = (9.0 + 64.0 * eps(Float64)) * eps(Float64) * eps(Float64)
const EPS_ORIENT_A_FLOAT32 = (3.0 + 16.0 * eps(Float32)) * eps(Float32)
const EPS_ORIENT_B_FLOAT32 = (2.0 + 12.0 * eps(Float32)) * eps(Float32)
const EPS_ORIENT_C1_FLOAT32 = (3.0 + 8.0 * eps(Float32)) * eps(Float32)
const EPS_ORIENT_C2_FLOAT32 = (9.0 + 64.0 * eps(Float32)) * eps(Float32) * eps(Float32)
const EPS_ORIENT_A_FLOAT16 = (3.0 + 16.0 * eps(Float16)) * eps(Float16)
const EPS_ORIENT_B_FLOAT16 = (2.0 + 12.0 * eps(Float16)) * eps(Float16)
const EPS_ORIENT_C1_FLOAT16 = (3.0 + 8.0 * eps(Float16)) * eps(Float16)
const EPS_ORIENT_C2_FLOAT16 = (9.0 + 64.0 * eps(Float16)) * eps(Float16) * eps(Float16)
const EPS_INCIRCLE_A_FLOAT64 = (10.0 + 96.0 * eps(Float64)) * eps(Float64)
const EPS_INCIRCLE_A_FLOAT32 = (10.0 + 96.0 * eps(Float32)) * eps(Float32)
const EPS_INCIRCLE_A_FLOAT16 = (10.0 + 96.0 * eps(Float16)) * eps(Float16)

function orient_fast(a::NTuple{2,T}, b::NTuple{2,T}, c::NTuple{2,T}) where {T}
    ac1, ac2 = a[1] - c[1], a[2] - c[2]
    bc1, bc2 = b[1] - c[1], b[2] - c[2]
    return ac1 * bc2 - ac2 * bc1
end

function orient_exact(a::NTuple{2,T}, b::NTuple{2,T}, c::NTuple{2,T}) where {T}
    ab = two_prod(a[1], b[2])
    ac = two_prod(a[1], c[2])
    bc = two_prod(b[1], c[2])
    ba = two_prod(b[1], a[2])
    ca = two_prod(c[1], a[2])
    cb = two_prod(c[1], b[2])
    u = expansion_sum(expansion_diff(ab, ac), expansion_diff(bc, ba))
    v = expansion_sum(u, expansion_diff(ca, cb))
    return v[end]
end

@inline orient_test_A(A::Float64, detsum::Float64) = abs(A) >= EPS_ORIENT_A_FLOAT64 * detsum
@inline orient_test_A(A::Float32, detsum::Float32) = abs(A) >= EPS_ORIENT_A_FLOAT32 * detsum
@inline orient_test_A(A::Float16, detsum::Float16) = abs(A) >= EPS_ORIENT_A_FLOAT16 * detsum
@inline orient_test_B(B::Float64, detsum::Float64) = abs(B) >= EPS_ORIENT_B_FLOAT64 * detsum
@inline orient_test_B(B::Float32, detsum::Float32) = abs(B) >= EPS_ORIENT_B_FLOAT32 * detsum
@inline orient_test_B(B::Float16, detsum::Float16) = abs(B) >= EPS_ORIENT_B_FLOAT16 * detsum
@inline function orient_test_C(C::Float64, B::Float64, detsum::Float64)
    return abs(C) >= EPS_ORIENT_C1_FLOAT64 * abs(B) + EPS_ORIENT_C2_FLOAT64 * detsum
end
@inline function orient_test_C(C::Float32, B::Float32, detsum::Float32)
    return abs(C) >= EPS_ORIENT_C1_FLOAT32 * abs(B) + EPS_ORIENT_C2_FLOAT32 * detsum
end
@inline function orient_test_C(C::Float16, B::Float16, detsum::Float16)
    return abs(C) >= EPS_ORIENT_C1_FLOAT16 * abs(B) + EPS_ORIENT_C2_FLOAT16 * detsum
end

function orient(a::NTuple{2,T}, b::NTuple{2,T}, c::NTuple{2,T}) where {T}
    ac1, ac2 = a[1] - c[1], a[2] - c[2]
    bc1, bc2 = b[1] - c[1], b[2] - c[2]
    det1 = ac1 * bc2
    det2 = ac2 * bc1
    det = det1 - det2
    detsum = det1 + det2
    # The result is only questionable when det is O(epsilon)
    if det1 > 0.0
        if det2 <= 0.0
            return det
        end
    elseif det1 < 0.0
        if det2 >= 0.01
            return det
        else
            detsum *= -1.0
        end
    else
        return det
    end
    if orient_test_A(det, detsum)
        return det
    else
        return orient_adapt(a, b, c, ac1, ac2, bc1, bc2, detsum)
    end
end

function orient_adapt(a::NTuple{2,T}, b::NTuple{2,T}, c::NTuple{2,T},
                      ac1::T, ac2::T, bc1::T, bc2::T, detsum::T) where {T}
    det1 = two_prod(ac1, bc2)
    det2 = two_prod(ac2, bc1)
    B = expansion_diff(det1, det2)
    B_prime = approximate(B)
    if orient_test_B(B_prime, detsum)
        return B_prime
    end

    ac1_error = two_diff_tail(a[1], c[1], ac1)
    ac2_error = two_diff_tail(a[2], c[2], ac2)
    bc1_error = two_diff_tail(b[1], c[1], bc1)
    bc2_error = two_diff_tail(b[2], c[2], bc2)
    if (ac1_error == 0.0) && (ac2_error == 0.0) && (bc1_error == 0.0) && (bc2_error == 0.0)
        return B_prime
    end

    C = B_prime + (ac1 * bc2_error + bc2 * ac1_error) - (ac2 * bc1_error + bc1 * ac2_error)
    if orient_test_C(C, B_prime, detsum)
        return C
    end

    s = two_prod(ac1_error, bc2)
    t = two_prod(ac2_error, bc1)
    u = expansion_diff(s, t)
    C1 = expansion_sum(B, u)

    s = two_prod(ac1, bc2_error)
    t = two_prod(ac2, bc1_error)
    u = expansion_diff(s, t)
    C2 = expansion_sum(C1, u)

    s = two_prod(ac1_error, bc2_error)
    t = two_prod(ac2_error, bc1_error)
    u = expansion_diff(s, t)
    D = expansion_sum(C2, u)

    return D[end]
end

function incircle_fast(a::NTuple{2,T}, b::NTuple{2,T}, c::NTuple{2,T}, d::NTuple{2,T}) where {T}
    ad1, ad2 = a[1] - d[1], a[2] - d[2]
    bd1, bd2 = b[1] - d[1], b[2] - d[2]
    cd1, cd2 = c[1] - d[1], c[2] - d[2]
    ab_det = ad1 * bd2 - bd1 * ad2
    bc_det = bd1 * cd2 - cd1 * bd2
    ca_det = cd1 * ad2 - ad1 * cd2
    a_lift = ad1 * ad1 + ad2 * ad2
    b_lift = bd1 * bd1 + bd2 * bd2
    c_lift = cd1 * cd1 + cd2 * cd2
    return a_lift * bc_det + b_lift * ca_det + c_lift * ab_det
end

function incircle_exact(a::NTuple{2,T}, b::NTuple{2,T}, c::NTuple{2,T}, d::NTuple{2,T}) where {T}
    ab_1, ab_2 = two_prod(a[1], b[2]), two_prod(b[1], a[2])
    bc_1, bc_2 = two_prod(b[1], c[2]), two_prod(c[1], b[2])
    cd_1, cd_2 = two_prod(c[1], d[2]), two_prod(d[1], c[2])
    da_1, da_2 = two_prod(d[1], a[2]), two_prod(a[1], d[2])
    ac_1, ac_2 = two_prod(a[1], c[2]), two_prod(c[1], a[2])
    bd_1, bd_2 = two_prod(b[1], d[2]), two_prod(d[1], b[2])

    ab = expansion_diff(ab_1, ab_2)
    bc = expansion_diff(bc_1, bc_2)
    cd = expansion_diff(cd_1, cd_2)
    da = expansion_diff(da_1, da_2)
    ac = expansion_diff(ac_1, ac_2)
    bd = expansion_diff(bd_1, bd_2)

    ca = (-ac[1], -ac[2], -ac[3], -ac[4])
    db = (-bd[1], -bd[2], -bd[3], -bd[4])

    cda = expansion_sum(expansion_sum(cd, da), ac)
    dab = expansion_sum(expansion_sum(da, ab), bd)
    abc = expansion_sum(expansion_sum(ab, bc), ca)
    bcd = expansion_sum(expansion_sum(bc, cd), db)

    adet_1 = expansion_scale(expansion_scale(bcd, a[1]), a[1])
    adet_2 = expansion_scale(expansion_scale(bcd, a[2]), a[2])
    bdet_1 = expansion_scale(expansion_scale(cda, b[1]), -b[1])
    bdet_2 = expansion_scale(expansion_scale(cda, b[2]), -b[2])
    cdet_1 = expansion_scale(expansion_scale(dab, c[1]), c[1])
    cdet_2 = expansion_scale(expansion_scale(dab, c[2]), c[2])
    ddet_1 = expansion_scale(expansion_scale(abc, d[1]), -d[1])
    ddet_2 = expansion_scale(expansion_scale(abc, d[2]), -d[2])

    adet = expansion_sum(adet_1, adet_2)
    bdet = expansion_sum(bdet_1, bdet_2)
    cdet = expansion_sum(cdet_1, cdet_2)
    ddet = expansion_sum(ddet_1, ddet_2)

    det = expansion_sum(expansion_sum(adet, bdet), expansion_sum(cdet, ddet))
    return det[end]
end

@inline incircle_test_A(A::Float64, alpha::Float64) = abs(A) >= EPS_INCIRCLE_A_FLOAT64 * alpha
@inline incircle_test_A(A::Float32, alpha::Float32) = abs(A) >= EPS_INCIRCLE_A_FLOAT32 * alpha
@inline incircle_test_A(A::Float16, alpha::Float16) = abs(A) >= EPS_INCIRCLE_A_FLOAT16 * alpha

# TODO: Add adaptive method instead of falling back to exact
function incircle(a::NTuple{2,T}, b::NTuple{2,T}, c::NTuple{2,T}, d::NTuple{2,T}) where {T}
    ad1, ad2 = a[1] - d[1], a[2] - d[2]
    bd1, bd2 = b[1] - d[1], b[2] - d[2]
    cd1, cd2 = c[1] - d[1], c[2] - d[2]
    adbd, bdad = ad1 * bd2, bd1 * ad2
    bdcd, cdbd = bd1 * cd2, cd1 * bd2
    cdad, adcd = cd1 * ad2, ad1 * cd2
    a_lift = ad1 * ad1 + ad2 * ad2
    b_lift = bd1 * bd1 + bd2 * bd2
    c_lift = cd1 * cd1 + cd2 * cd2
    A = (bdcd - cdbd) * a_lift + (cdad - adcd) * b_lift + (adbd - bdad) * c_lift
    alpha = (
        (abs(bdcd) + abs(cdbd)) * a_lift +
        (abs(cdad) + abs(adcd)) * b_lift +
        (abs(adbd) + abs(bdad)) * c_lift
    )
    if incircle_test_A(A, alpha)
        return A
    else
        return incircle_exact(a, b, c, d)
    end
end

# Test if a is in the outer half plane of bc
function incircle(a::NTuple{2,T}, b::NTuple{2,T}, c::NTuple{2,T}) where {T}
    # Check open half plane
    in_open_half_plane = orient(a, b, c)
    if in_open_half_plane == zero(T)
        # Check if a lies along the segment bc
        if (b[1] == c[1]) && (a[2] > min(b[2], c[2])) && (a[2] < max(b[2], c[2]))
            in_open_half_plane = one(T)
        elseif (a[1] > min(b[1], c[1])) && (a[1] < max(b[1], c[1]))
            in_open_half_plane = one(T)
        end
    end
    return in_open_half_plane
end

end # module

max_stars_repo_name: PaulMatlashewski/DelaunayTriangulations

max_issues_repo_name: PaulMatlashewski/DelaunayTriangulations

max_forks_repo_name: PaulMatlashewski/DelaunayTriangulations


text:
module FHIRClient

import Base64
import Dates
import HTTP
import JSON3
import SaferIntegers
import StructTypes
import TimeZones

include("types.jl")

include("r4.jl")

include("credentials.jl")
include("fhir-to-julia.jl")
include("headers.jl")
include("other-fhir-versions.jl")
include("requests.jl")

end # end module FHIRClient


text:
##############
# maximize_f #
##############

"""
Optimizes f using Newton's method and exact Hessians.

Returns:
  - iter_count: The number of iterations taken
  - max_f: The maximum function value achieved
  - max_x: The optimal function input
  - ret: The return code of optimize()
"""
function maximize_f{F}(f::F, ea::ElboArgs, transform::DataTransform;
                       omitted_ids=Int[],
                       use_default_optim_params=false,
                       xtol_rel=1e-7,
                       ftol_abs=1e-6,
                       verbose=false,
                       max_iters=50,
                       fast_hessian=true)
    # Make sure the model parameters are within the transform bounds
    enforce_bounds!(ea.vp, ea.active_sources, transform)
    @assert ea.active_sources == transform.active_sources

    f_evals = 0
    n_active_sources::Int = length(transform.active_sources)
    kept_ids::Vector{Int} = setdiff(1:length(UnconstrainedParams), omitted_ids)

    all_kept_ids::Vector{Int} = Int[]
    for i in eachindex(transform.active_sources)
        append!(all_kept_ids, kept_ids + (i - 1) * length(kept_ids))
    end

    x0 = vec(Transform.vp_to_array(transform, ea.vp, omitted_ids))
    last_sf::SensitiveFloat{Float64} = SensitiveFloat{Float64}(
                    length(UnconstrainedParams), n_active_sources, true, true)
    last_x::Vector{Float64} = fill(NaN, size(x0))

    f_wrapped_nocache! = (x::Vector) -> begin
        # Evaluate in the constrained space and then unconstrain again.
        reshaped_x = reshape(x, length(kept_ids), n_active_sources)
        Transform.array_to_vp!(transform, reshaped_x, ea.vp, kept_ids)
        f_res = f(ea)
        f_evals += 1
        verbose && log_f_eval(f_evals, ea, f_res)
        Transform.transform_sensitive_float!(transform, last_sf, f_res, ea.vp, ea.active_sources)
    end

    f_wrapped_cached! = (x::Vector) -> begin
        if x != last_x
            copy!(last_x, x)
            f_wrapped_nocache!(x)
        end
        return nothing
    end

    neg_f_value = (x::Vector) -> begin
        f_wrapped_cached!(x)
        return -(last_sf.v[])
    end

    neg_f_grad! = (x::Vector, grad::Vector) -> begin
        f_wrapped_cached!(x)
        for i in kept_ids
            for j in 1:n_active_sources
                grad[i, j] = -last_sf.d[i, j]
            end
        end
        return grad
    end

    neg_f_hessian! = (x::Vector, hess::Matrix) -> begin
        f_wrapped_cached!(x)
        for i in all_kept_ids
            for j in all_kept_ids
                hess[i, j] = last_sf.h[i, j]
            end
        end
        Transform.symmetrize!(hess, -0.5)
        return hess
    end

    if use_default_optim_params
        tr_method = Optim.NewtonTrustRegion()
    else
        tr_method = Optim.NewtonTrustRegion(
            initial_delta=10.0,
            delta_hat=1e9,
            eta=0.1,
            rho_lower=0.25,
            rho_upper=0.75)
    end

    options = Optim.OptimizationOptions(;
        x_tol = xtol_rel, f_tol = ftol_abs, g_tol = 1e-8,
        iterations = max_iters, store_trace = verbose,
        show_trace = false, extended_trace = verbose)

    nm_result = Optim.optimize(neg_f_value,
                               neg_f_grad!,
                               neg_f_hessian!,
                               x0,
                               tr_method,
                               options)

    reshaped_min = reshape(Optim.minimizer(nm_result), size(x0))
    Transform.array_to_vp!(transform, reshaped_min, ea.vp, kept_ids)
    max_f = -1.0 * Optim.minimum(nm_result)
    max_x = Optim.minimizer(nm_result)

    #Log.debug("elbo is $max_f after $(nm_result.iterations) Newton steps")
    return f_evals, max_f, max_x, nm_result, transform
end

function maximize_f{F}(f::F, ea::ElboArgs;
                       loc_width=1.5e-3,
                       loc_scale=1.0,
                       omitted_ids=Int[],
                       xtol_rel=1e-7,
                       ftol_abs=1e-6,
                       verbose=false,
                       max_iters=50,
                       fast_hessian=true,
                       use_default_optim_params=false)
    transform = get_mp_transform(ea.vp, ea.active_sources,
                                 loc_width=loc_width, loc_scale=loc_scale)

    maximize_f(f, ea, transform;
                omitted_ids=omitted_ids,
                xtol_rel=xtol_rel,
                ftol_abs=ftol_abs,
                verbose=verbose,
                max_iters=max_iters,
                fast_hessian=fast_hessian,
                use_default_optim_params=use_default_optim_params)
end

function log_f_eval(f_evals, ea::ElboArgs, f_res)
    # TODO: Add an option to print either the transformed or
    # free parameterizations.
    Log.debug("f_evals=$(f_evals) | value=$(f_res.v[])")

    iter_vp = ea.vp[ea.active_sources]

    if length(iter_vp[1]) == length(ids_names)
        state_df = DataFrames.DataFrame(names=ids_names)
    elseif length(iter_vp[1]) == length(ids_free_names)
        state_df = DataFrames.DataFrame(names=ids_free_names)
    else
        state_df = DataFrames.DataFrame(names=[ "x$i" for i=1:length(iter_vp[1, :])])
    end

    for s in eachindex(iter_vp)
        state_df[Symbol(string("val", s))] = iter_vp[s]
    end

    for s in eachindex(iter_vp)
        state_df[Symbol(string("grad", s))] = f_res.d[:, s]
    end

    Log.debug(repr(state_df))
    Log.debug("=======================================\n")
end


text:
@testset "PTBLM" begin
    n_obs = 5574
    n_targets = 1
    n_features = ()
    Tx=String
    Ty=String

    d = SMSSpamCollection()

    test_supervised_array_dataset(d;
        n_obs, n_targets, n_features,
        Tx, Ty)
end

@testset "SMSSpamCollection" begin
    n_obs = 5574
    n_targets = 1
    n_features = ()
    Tx=String
    Ty=String

    d = SMSSpamCollection()

    test_supervised_array_dataset(d;
        n_obs, n_targets, n_features,
        Tx, Ty)
end

@testset "UD_English" begin
    n_features = ()
    Tx = Vector{Vector{String}}
    for (n_obs, split) in [(12543, :train), (2077, :test), (2001, :dev)]
        d = UD_English(split)

        test_unsupervised_array_dataset(d;
            n_obs, n_features, Tx)
    end
end


text:
using Base.Test

@testset "Lagrange" begin include("test_lagrange.jl") end
@testset "Benders" begin include("test_benders.jl") end


text:
using RefractiveIndex
using Base.Test

# write your own tests here
materials = RefractiveIndex.Material[]
push!(materials, getdata("main", "au", "babar")) # tabulated nk
push!(materials, getdata("glass", "ami", "amtir-2")) # formula 1 / tabulated k
push!(materials, getdata("glass", "schott", "k7")) # formula 2 / tabulated k
push!(materials, getdata("glass", "hoya", "bac4")) # formula 3 / tabulated k
push!(materials, getdata("main", "Ag3AsS3", "Hulme-e.yml")) # formula 4
push!(materials, getdata("glass", "misc", "soda-lime/Rubin-bronze.yml")) # formula 5
push!(materials, getdata("main", "Ar", "Bideau-Mehu.yml")) # formula 6
push!(materials, getdata("main", "Si", "Edwards.yml")) # formula 7
push!(materials, getdata("main", "AgBr", "Schroter.yml")) # formula 8
push!(materials, getdata("organic", "CH4N2O - urea", "Rosker-e.yml")) # formula 9

λ =     [1.937,   1.937,     1.937,        1.0,       1.3,     1.3,         0.5677,     5.3,     0.5,     0.5]
nvals = [0.25991, 2.8104,    1.4909,       1.5578,    2.5608,  1.5100,      1.00028201, 3.4254,  2.3095,  1.6167]
kvals = [13.099,  0.00000023, 0.00000137, 0.00000002, nothing, 0.00001430, nothing,    nothing, nothing, nothing]
num_dec_k = [3, 8, 8, 8, nothing, 8, nothing, nothing, nothing, nothing]

function getalldata()
    print("|")
    for (i, shelf) in enumerate(getshelves())
        for book in getbooks(shelf), page in getpages(shelf,book)
            # println("\"$shelf\", \"$book\", \"$page\"")
            @test typeof(getdata(shelf, book, page)) == RefractiveIndex.Material
        end
    end
end

function check_refractiveindices(materials, λ, nvals, kvals)
    for i = 1:length(materials)
        num_dec_n = length(split(string(nvals[i]), ".")[2])
        @test round(materials[i].nfun(λ[i]), num_dec_n) ≈ nvals[i]
        kvals[i] == nothing || @test round(materials[i].kfun(λ[i]), num_dec_k[i]) ≈ kvals[i]
    end
end

@testset "Data Reachability Test" begin getalldata() end
@testset "Refractive Index Tests" begin check_refractiveindices(materials, λ, nvals, kvals) end


text:
# Autogenerated wrapper script for Aria2_jll for x86_64-linux-gnu-cxx03
export aria2c, libaria2

using Cares_jll
using LibSSH2_jll
using OpenSSL_jll
using XML2_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("Aria2")
JLLWrappers.@declare_executable_product(aria2c)
JLLWrappers.@declare_library_product(libaria2, "libaria2.so.0")
function __init__()
    JLLWrappers.@generate_init_header(Cares_jll, LibSSH2_jll, OpenSSL_jll, XML2_jll, Zlib_jll)
    JLLWrappers.@init_executable_product(
        aria2c,
        "bin/aria2c",
    )

    JLLWrappers.@init_library_product(
        libaria2,
        "lib/libaria2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()