text:
using PyPlot
using SpecialFunctions
Np=Int(1e4)
Nx=Int(6)
dt=0.05
tmax=1000
deg=0; #Smoothness of periodic boundary
sigmaL=4*pi;
vth=0.1;
gamma=1; # Electrons (-1), gravity 1
eps=0
pad=0

function f(x,v)
  return exp.(-0.5.*(v/vth).^2)/sqrt(2*pi)/vth.*
              exp.(-0.5.*(x/sigmaL).^2).*(1+ eps*sqrt(2)*((x/sigmaL).^2-0.5) )/sqrt(2*pi)/(1+0.5*eps)
end
function g(x,v)
  return exp.(-0.5.*(v/vth).^2)./sqrt(2*pi)/vth.*exp.(-0.5.*(x/sigmaL).^2)/sqrt(2*pi)/sigmaL
end

using Sobol
sob = SobolSeq(2)
skip(sob, 5) # Skip some entries
xk=Array(Float64,Np); vk=similar(xk)
for i=1:Np
  xk[i],vk[i]=Sobol.next(sob)
end

# Sampling
norminv=p->sqrt(2.).*erfinv.(2.*p - 1.0)
xk.=norminv(xk)*sigmaL
vk.=norminv(vk)*vth

fk=f(xk,vk)
wk=f(xk,vk)./g(xk,vk)

# Hermite functions
function HermiteH{T}(x::T,N::Int)
  Hn=Array(T, N)
  Hn[1]=1.
  Hn[2]=sqrt(2)*x
  @inbounds for idx=3:N
    n=idx-2
    Hn[idx]=x*sqrt(2/(n+1))*Hn[idx-1]- sqrt(n/(n+1))*Hn[idx-2]
  end
  return Hn*exp(x^2/2)
end

 # Hermite functions
 function HermiteHdx{T}(x::T,N::Int)
   Hdxn=Array(T, N)
   H0::T=1.0
   H1::T=sqrt(2)*x
   Hdxn[1]=-x   #-1./sqrt(2)*H1
   @inbounds for idx=2:N
     n=idx-1
     H2::T=x*sqrt(2/(n+1))*H1- sqrt(n/(n+1))*H0
     Hdxn[idx]=sqrt(n/2)*H0 -sqrt((n+1)/2)*H2
     H0=H1; H1=H2
   end
   return Hdxn.*exp(-x^2/2)
 end

function HermiteHdx_series{T}( x::Array{T,1}, u::Array{T,1} )
  E=Array(T,length(x))
  Np=length(x)
  Nx=length(u)
  @inbounds for pdx=1:Np
      E[pdx]=dot(HermiteHdx(x[pdx],Nx),u)
    end
  return E
end

function diff_Hermite{T}(H::Array{T,1})
  N=length(H)-1
  dH=zeros(T,length(H))
  dH[1]=1/sqrt(2)*H[2]
  dH[end]=-sqrt(N/2)*H[end-1]
  for idx=2:N-1
    n=idx-1
    dH[idx]=sqrt((n+1)/2)*H[idx+1] -sqrt(n/2)*H[idx-1]
  end
  return dH
end

function poissonK_Hermite(N::Int)
  K=spdiagm([-sqrt( (3:N).*((3:N)-1)*pi)/2,
             ((1:N)+0.5)*sqrt(pi)  ,
             -sqrt(((1:N-2)+2).*((1:N-2)+1)*pi)/2],
                   [2,0,-2])
  return K
end

@inbounds function  rhs_particle_Hermite{T}(xk::Array{T,1},wk::Array{T,1}, Nx::Int)
    const Np=length(xk)
    rho=zeros(T,Nx)
    # Precalculate coefficients
    alpha=Array(T,Nx-2)
    beta=Array(T,Nx-2)
    @inbounds for n=1:Nx-2
      alpha[n]=sqrt(2/(n+1))
      beta[n]=-sqrt(n/(n+1))
    end

    @inbounds for pdx=1:Np
      w::T=exp(-xk[pdx].^2/2)*wk[pdx]
      P0::T=1.0
      P1::T=sqrt(2)*P0*xk[pdx]
      rho[1]+=P0.*w
      rho[2]+=P1.*w
      for idx=1:Nx-2
        P2::T=xk[pdx].*alpha[idx].*P1+ beta[idx].*P0
        P0=P1; P1=P2;
        rho[idx+2]+=P2.*w
      end
    end
  return rho
end

@inbounds function HermiteH_series{T}(x::Array{T,1}, u::Array{T,1})
  N::Int=length(u)-1;
  alpha=Array(T,N+1)
  beta=Array(T,N+1)
  #
  for n=0:N
     alpha[n+1]=sqrt(2/(n+1))
     beta[n+1]=-sqrt((n+1)/(n+2))
  end
  E=Array(T,length(x))
  for pdx=1:length(x)
   b1::T=0
   b2::T=0
   for n=N:-1:0
     b0::T=u[n+1]+ x[pdx]*alpha[n+1]*b1 + beta[n+1]*b2
     b2=b1; b1=b0;
    end
    E[pdx]=b1*exp(-x[pdx]^2/2)
  end
  return E
end

# Coefficients for Runge Kutta
rksd=[2/3.0, -2/3.0, 1  ]
rksc=[ 7/24.0, 3/4.0, -1/24.0]

Nt=Int(ceil(tmax/dt))
fieldenergy=zeros(Float64,Nt)
kineticenergy=zeros(Float64,Nt)
momentum=zeros(Float64,Nt)
rho=zeros(Float64,Nx)
Phi=zeros(Float64,Nx+pad)
E=zeros(Float64,Nx+pad+2)
K=poissonK_Hermite(Nx+pad)

figure()
tic()
@inbounds for tdx=1:Nt
  if (mod(tdx-1,50)==0)
    PyPlot.clf()

    scatter(xk,vk,s=2, c=fk,lw=0)
    axis([-sigmaL*5,sigmaL*5,-3.5*vth,3.5*vth])

    sleep(0.01)
  end
  @inbounds for rkdx=1:length(rksd)

    # Bin to right hand side
    rhs=rhs_particle_Hermite( xk/sigmaL,wk,Nx)/Np
    # Hermite L^2 projection
    rho=rhs./sqrt(pi)
    # Remove ion background (quick'n dirty)
    rhs[1]=0
    Phi.=\(-K, [rhs;zeros(Float64,pad)])

    E.=diff_Hermite(gamma*[Phi;0;0])

    if rkdx==1
      fieldenergy[tdx]=(Phi'*K*Phi)[1]/2*sigmaL  #sum(E.^2)*sqrt(pi) #
      kineticenergy[tdx]=0.5.*sum(vk.^2.*wk )/Np
      momentum[tdx]=sum(vk.*wk)/Np
    end

    vk.+=dt*rksc[rkdx]*HermiteH_series(xk/sigmaL,E)
    #vk.+=dt*rksc[rkdx]*HermiteHdx_series(xk/sigmaL,-Phi)
    xk.+=dt*rksd[rkdx].*vk

    # Reflecting boundary
  end

  #print(tdx)
end
toc()
ttime=(0:Nt-1)*dt

prefix="plots/hermite_sheath"
#prefix="legendre"

#
# xg=collect(linspace(minimum(xk),maximum(xk)))
# figure()
# plot(xg, HermiteH_series(xg/sigmaL,Phi))

# HermiteHdx_series(xk/sigmaL,-Phi)
# HermiteH_series(xk/sigmaL,E)
#
# figure()
# plot(abs(rhs))

# plot(abs(Phi))
# plot(abs(E))
# plot(1:Nx,abs(rho))

# figure()
# title("magnitude of coefficients")
# semilogy(abs(coefficients(rho)),label=L"density $\rho$")
# semilogy(abs(coefficients(Phi)),label=L"potential $\Phi$")
# semilogy(abs(coefficients(D1*Phi)),label=L"field $E$")
# xlabel("index")
# legend()
# grid()
# savefig("$prefix\_field_coeffs.png")

# figure()
# plot(kineticenergy)
# plot(kineticenergy+fieldenergy*sigmaL)

PyPlot.matplotlib[:rc]("font", size=18)

figure()
semilogy(ttime,fieldenergy)
omega=1.415661888604536 - 0.153359466909605im;
semilogy(ttime[ttime.<8.],
0.5*fieldenergy[1]*abs(exp(-im*omega*(ttime[ttime.<8.]-0.4))).^2,
lw=3,"r")
autoscale(tight="True");grid()
ylabel("electrostatic energy");
xlabel("time")
#tight_layout()
# savefig("$prefix\_fieldenergy.png")

# figure()
# semilogy(ttime,kineticenergy)
# grid();
# #autoscale(tight="True")
# ylabel("kinetic energy")
# xlabel("time")
# savefig("$prefix\_kineticenergy.png")

figure()
energy=kineticenergy+fieldenergy
semilogy(ttime,abs((energy-energy[1])./energy[1]))
grid()
ylabel("rel. energy error")
xlabel("time ")
# savefig("$prefix\_energy_error.png")
figure()
semilogy(ttime,abs(momentum-momentum[1]))
xlabel("time "); grid()
ylabel("absolute momentum error")
# savefig("$prefix\_momentum_error.png")

figure()
scatter(xk,vk,s=1, c=fk,lw=0)
# autoscale(tight="True")
# axis([0,L,minimum(vk),maximum(vk)])
colorbar();tight_layout()
xlabel("x");ylabel("v")
# savefig("$prefix\_particles.png")

#
# #Test OSDE HERMITE
# h=x->exp(-x.^2/2).*(x-0.5)+exp(-(x-2.2).^2)
# Dh=x->exp(-x.^2/2).*(x-0.5).*(-x) + exp(-x.^2/2)  -exp(-(x-2.2).^2).*2.*(x-2.2)
# Np=Int(1e7)
# xk=randn(Float64,Np)
#
# wk=h(xk)./(exp(-0.5.*xk.^2)/sqrt(2*pi) )
#
#
#
# hhat=rhs_particle_Hermite(xk,wk,40)/Np/sqrt(pi)
#
# x=collect(linspace(-15,15,2000))
# figure()
# plot(x,h(x))
# plot(x,HermiteH_series(x,hhat))
#
# figure()
# # plot(x,Dh(x))
# plot(x,gradient(h(x),x[2]-x[1] ))
#
# Dhhat=diff_Hermite(hhat)
#
#
# plot(x,HermiteH_series(x,Dhhat))
#  plot(x,HermiteHdx_series(x,hhat))

max_stars_repo_name: ameresj/VlasovPIF

max_issues_repo_name: ameresj/VlasovPIF

max_forks_repo_name: ameresj/VlasovPIF


text:
"""
    MMWInterval <: ConfidenceInterval

Exact confidence interval based on ordering induced by `estimator` minimizing the
expected squared width with respect to weight function `prior(p::Real)`.
"""
struct MMWInterval{TE<:Estimator,TR<:Real} <: ConfidenceInterval

    estimator::TE
    confidence::TR
    limits::Array{TR,3}
    prior::Function

    function MMWInterval{TE,TR}(
        estimator::TE,
        rho0::TR,
        prior::Function,
        solver::MathProgBase.AbstractMathProgSolver,
        confidence::TR,
        ngrid::Integer #   = 100
    ) where {TE<:Estimator,TR<:Real}
        alpha =  1 - confidence
        ngrid += 1 # number of gridpoints is number of intervals + 1
        d     = design(estimator)
        n1    = interimsamplesize(d)
        nmax::Int64  = maximum(samplesize(d))

        # construct array of possible outcomes
        x1x2_poss = zeros(Int64, 1, 2)
        for x1 in 1:n1
            for x2 in 0:convert(Int64, d.n[x1 + 1] - n1)
                x1x2_poss = [x1x2_poss; convert(Array{Int64, 2}, [x1 x2])]
            end
        end
        nposs = size(x1x2_poss, 1)
        function ispossible(x1, x2)
            sum(mapslices(x -> x == vec([x1 x2]), x1x2_poss, 2)) >= 1
        end

        # set up arrays needed for coefficients
        rho_grid = linspace(0, 1, ngrid)
        # rho_grid_refined = linspace(0, 1, (ngrid - 1)*refinement_fctr)
        probs = zeros(Float64, nposs, ngrid)
        exp_probs = zeros(Float64, nposs)
        estimates = zeros(Int64, nposs, 2)
        for o in 1:nposs
            x1 = x1x2_poss[o, 1]
            x2 = x1x2_poss[o, 2]
            # bracket estimator in grid
            i  = 1
            while rho_grid[i] < estimate(estimator, x1, x2)
                i += 1
            end
            if rho_grid[i] == estimate(estimator, x1, x2)
                estimates[o, 1] = i
                estimates[o, 2] = i
            else
                estimates[o, 1] = i
                estimates[o, 2] = i + 1
            end
            # precompute probabilities
            exp_probs[o] = QuadGK.quadgk(rho -> prior(rho)*pdf(d, x1, x2, rho), 0.0, 1.0, abstol = 10^(-3.0))[1]
            for i in 1:ngrid
                probs[o, i] = pdf(d, x1, x2, rho_grid[i])
            end
        end
        model = JuMP.Model(solver = solver)
        JuMP.@variable(model, # lower bounds
            y[o = 1:nposs, i = 1:(ngrid - 1), j = (i + 1):ngrid],
            Bin
        )
        JuMP.@variable(model,
            z_mon[o = 1:nposs, i = 1:2] >= 0
        )
        a = 0
        for i in 1:ngrid
            if rho0 > rho_grid[i]
                a = i
            end
        end
        for o in 1:nposs
            x1 = x1x2_poss[o, 1]
            x2 = x1x2_poss[o, 2]
            # monotonicity violations
            if o > 1
                x1_prev = x1x2_poss[o - 1, 1]
                x2_prev = x1x2_poss[o - 1, 2]
                if (x1 == x1_prev) & (x2 == x2_prev + 1)
                    expr1 = JuMP.AffExpr()
                    expr2 = JuMP.AffExpr()
                    for i in 1:(ngrid - 1)
                        for j in (i + 1):ngrid
                            JuMP.append!(expr1, rho_grid[i]*(y[o - 1, i, j] - y[o, i, j]))
                            JuMP.append!(expr2, rho_grid[j]*(y[o - 1, i, j] - y[o, i, j]))
                        end
                    end
                    JuMP.@constraint(model,
                        # expr1 <= 0
                        z_mon[o, 1] >= expr1
                    )
                    JuMP.@constraint(model,
                        # expr2 <= 0
                        z_mon[o, 2] >= expr2
                    )
                end
            end
            # functional constraint
            funcConstrExpr = JuMP.AffExpr()
            for i in 1:(ngrid - 1)
                for j in (i + 1):ngrid
                    JuMP.append!(funcConstrExpr, y[o, i, j])
                end
            end
            JuMP.@constraint(model,
                funcConstrExpr == 1
            )
            # estimator consistency: estimator must lie in interval!
            estConstrExpr = JuMP.AffExpr()
            for i in 1:estimates[o, 1] # lower bound must be smaller or equal to estimate
                for j in max(i + 1, estimates[o, 2]):ngrid # upper bound must be larger or equal
                    JuMP.append!(estConstrExpr, y[o, i, j]) # a possible assignment
                end
            end
            JuMP.@constraint(model,
                estConstrExpr == 1 # any of the possible assigments must be selected!
            )
            # consistency 1: x1 + x2 > c(x1) -> lower bound must be strictly larger than rho0
            if x1 + x2 > d.c[x1 + 1]
                cons1ConstrExpr = JuMP.AffExpr()
                for i in 1:min(ngrid - 1, a) # rho_grid[a] <= rho0
                    for j in (i + 1):ngrid
                        JuMP.append!(cons1ConstrExpr, y[o, i, j])
                    end
                end
                JuMP.@constraint(model,
                    cons1ConstrExpr == 0 # none of the assignments with lower bound <= rho0 can be selected
                )
            end
            # consistency 2: x1 + x2 <= c(x1) -> lower bound must be smaller or equal to rho0
            if x1 + x2 <= d.c[x1 + 1]
                cons2ConstrExpr = JuMP.AffExpr()
                for i in 1:min(ngrid - 1, a) # rho_grid[a] <= rho0
                    for j in (i + 1):ngrid
                        JuMP.append!(cons2ConstrExpr, y[o, i, j])
                    end
                end
                JuMP.@constraint(model,
                    cons2ConstrExpr == 1 # one of the assignments with lower bound <= rho0 must be selected
                )
            end
            # upper boundary bleeding, if x1 + x2 == n(x1) the upper boundary must be 1 to prevent coverage from dropping towards zero
            if x1 + x2 == d.n[x1 + 1]
                boundaryUpperConstrExpr = JuMP.AffExpr()
                for i in 1:(ngrid - 1)
                    JuMP.append!(boundaryUpperConstrExpr, y[o, i, ngrid])
                end
                JuMP.@constraint(model,
                    boundaryUpperConstrExpr == 1
                )
            end
        end
        # lower boundary bleeding, if x1 + x2 == 0 lower boundary must be 0 to prevent coverage fro mdropping too low
        boundaryLowerConstrExpr = JuMP.AffExpr()
        for j in 2:ngrid
            JuMP.append!(boundaryLowerConstrExpr, y[1, 1, j]) # x1x2_poss[1, :] = [0 0]
        end
        JuMP.@constraint(model,
            boundaryLowerConstrExpr == 1
        )
        for i in 1:ngrid
            lowerBoundCoverage = JuMP.AffExpr()
            upperBoundCoverage = JuMP.AffExpr()
            for o in 1:nposs
                for j in 1:(ngrid - 1)
                    for k in (j + 1):ngrid
                        if j <= i # lower bound <= rho_grid[i]
                            JuMP.append!(lowerBoundCoverage, probs[o, i]*y[o, j, k])
                        end
                        if k >= i # upper bound >= rho_grid[i]
                            JuMP.append!(upperBoundCoverage, probs[o, i]*y[o, j, k])
                        end
                    end
                end
            end
            JuMP.@constraint(model,
                lowerBoundCoverage >= 1 - alpha/2
            )
            JuMP.@constraint(model,
                upperBoundCoverage >= 1 - alpha/2
            )
        end
        objExpr = JuMP.AffExpr()
        for o in 1:nposs
            JuMP.append!(objExpr, (ngrid - 1)*z_mon[o, 1])
            JuMP.append!(objExpr, (ngrid - 1)*z_mon[o, 2])
            for i in 1:(ngrid - 1)
                for j in (i + 1):ngrid
                    JuMP.append!(objExpr, exp_probs[o]*(rho_grid[j] - rho_grid[i])*y[o, i, j])
                end
            end
        end
        JuMP.@objective(model, Min, objExpr)
        status = JuMP.solve(model)
        if !(status in (:Optimal, :UserLimit)) # no valid solution found!
            error("no feasible solution reached")
        end
        # create actual confidence interval object
        limits = fill(NaN, n1 + 1, nmax - n1 + 1, 2)
        for o in 1:nposs
            x1 = x1x2_poss[o, 1]
            x2 = x1x2_poss[o, 2]
            current_best_score = 0.0
            for i in 1:(ngrid - 1)
                for j in (i + 1):ngrid
                    val = JuMP.getvalue(y[o, i, j]) == 1.0
                    if val > current_best_score # in the end gives the solution where y[x1, n, c] is closest to 1
                        current_best_score = val
                        limits[x1 + 1, x2 + 1, :] = [rho_grid[i] rho_grid[j]]
                    end
                end
            end
        end
        return new(estimator, confidence, limits, prior)
    end
end

"""
    MMWInterval(
      estimator::TE,
      rho0::TR,
      prior::Function,
      solver::MathProgBase.AbstractMathProgSolver;
      confidence::TR = 0.9,
      ngrid::Integer = 100
    ) where {TE<:Estimator,TR<:Real}

Exact confidence interval based on ordering induced by `estimator` minimizing the
expected squared width with respect to weight function `prior(p::Real)`.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| estimator    | estimator object defining the sample space ordering |
| rho0         | upper boundary of null hypothesis |
| confidence   | confidence level of the interval |
| solver       | MathProgBase solver used for optimization, must support quadratic expressions |
| ngrid        | number of equally spaced grid-points on which to check coverage |
"""
function MMWInterval(
  estimator::TE,
  rho0::TR,
  prior::Function,
  solver::MathProgBase.AbstractMathProgSolver;
  confidence::TR = 0.9,
  ngrid::Integer = 100
) where {TE<:Estimator,TR<:Real}

  return MMWInterval{TE,TR}(estimator, rho0, prior, solver, confidence, ngrid)

end

Base.show(io::IO, ci::MMWInterval) = print("MMWInterval")

estimator(ci::MMWInterval) = ci.estimator

design(ci::MMWInterval) = ci |> estimator |> design

function limits{T<:Integer}(ci::MMWInterval, x1::T, x2::T)
    ispossible(design(ci), x1, x2) ? nothing : throw(InexactError())
    return convert(Vector{Float64}, view([ci.limits[x1 + 1, x2 + 1, 1] ci.limits[x1 + 1, x2 + 1, 2]], 1:2))
end

ext: jl

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
using LinearAlgebra

mutable struct NACTD{T,TP,TV,TW,OP,OV,OW} <: AbstractPolicyAgent where {T<:Real,TP<:AbstractPolicy,TV<:AbstractFuncApprox,TW<:AbstractFuncApprox,OP<:AbstractOptimizer,OV<:AbstractOptimizer,OW<:AbstractOptimizer}
    π::TP
    vf::TV
    wf::TW
    popt::OP
    vopt::OV
    wopt::OW
    γ::T
    normalizew::Bool
    gp::Array{T, 2}
    gv::Array{T, 2}
    function NACTD(π::AbstractPolicy, vf::AbstractFuncApprox, wf::AbstractFuncApprox, popt::AbstractOptimizer, vopt::AbstractOptimizer, wopt::AbstractOptimizer, γ::T, normalizew::Bool) where {T<:Real}
        tp = eltype(get_params(π))
        sp = size(get_params(π))
        sv = size(get_params(vf))

        π2 = clone(π)
        vf2 = clone(vf)
        wf2 = clone(wf)
        popt2 = clone(popt, π2)
        vopt2 = clone(vopt, vf2)
        wopt2 = clone(wopt, wf2)
        new{tp,typeof(π2),typeof(vf2),typeof(wf2),typeof(popt2),typeof(vopt2),typeof(wopt2)}(π2, vf2, wf2, popt2, vopt2, wopt2, tp(γ), normalizew, zeros(tp, sp), zeros(tp, sv))
    end
end

function act!(agent::NACTD, env::AbstractEnvironment, rng::AbstractRNG)
    v = call_gradient!(agent.gv, agent.vf, env.state, 1)
    logp = get_action_gradient_logp!(agent.gp, agent.π, env.state, rng)
    adv = call(agent.wf, vec(agent.gp), 1)
    r = step!(env, agent.π.action, rng)
    target = 0.
    vnext = 0.
    if is_terminal(env)
        target = r
    else
        vnext = call(agent.vf, env.state, 1)
        target = r + agent.γ * vnext
    end

    delta = target - v
    delta_adv = delta - adv
    update!(agent.vopt, delta, agent.gv)
    update!(agent.wopt, delta_adv, vec(agent.gp))

    w = vec(get_params(agent.wf))

    if agent.normalizew
        w = w .* (1. / norm(vec(w)))
    end

    update!(agent.popt, w)

    if is_terminal(env)
        new_episode!(agent, rng)
    end
end

function act!(agent::NACTD{T,TP,TV,TW,OP,Parl2AccumulatingTraceOptimizer{T,TV}}, env::AbstractEnvironment, rng::AbstractRNG) where {T,TP,TV,TW,OP}
    v = call_gradient!(agent.gv, agent.vf, env.state, 1)
    logp = get_action_gradient_logp!(agent.gp, agent.π, env.state, rng)
    adv = call(agent.wf, vec(agent.gp), 1)
    r = step!(env, agent.π.action, rng)
    target = 0.
    vnext = 0.
    gvnext = zeros(T, size(agent.gv))
    if is_terminal(env)
        target = r
    else
        vnext = call_gradient!(gvnext, agent.vf, env.state, 1)
        target = r + agent.γ * vnext
    end

    delta = target - v
    delta_adv = delta - adv

    update!(agent.vopt, delta, agent.gv, gvnext)
    update!(agent.wopt, delta_adv, vec(agent.gp))

    w = vec(get_params(agent.wf))

    if agent.normalizew
        w = w .* (1. / norm(vec(w)))
    end

    if is_terminal(env)
        new_episode!(agent, rng)
    end
end

function new_episode!(agent::NACTD, rng::AbstractRNG)
    new_episode!(agent.vopt)
    new_episode!(agent.wopt)
    new_episode!(agent.popt)
end

function clone(agent::NACTD)::NACTD
    a = NACTD(agent.π, agent.vf, agent.wf, agent.popt, agent.vopt, agent.wopt, agent.γ, agent.normalizew)
    return a
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
"""
    blockinfo(ssa::Array)

Infer a tree of WebAssembly blocks from ssa and return a BlockInfo struct containing
- a goto dict
- a list of which blocks each ssa index is a child of.

# Details
WebAssembly control instructions (mainly) consist of
- blocks: block, loop
- branching: br, br_if, return

# Notes concerning WebAssembly control instructions:
- We do not have Phi nodes (a block can only have one parent)
- We do not have gotos/jumps. We can only branch backwards/up the block tree.
- However, branching to a block continues at end of that block, which effectively is a forward jump (aka break)
- Special case: branching to a loop block continues at start of that block (aka continue)
- branching is specified in terms of number of levels up (br 0 goes to current block, br 1 goes to parent block and so on. return is essentialy sugar for br MAX)

So the strategy is to infer a blocktree from arbitrary gotos, insert blocks and translate gotos to branching in terms of levels up the tree.
"""
function blockinfo(ssa::Array)
    gotos = getgotos(ssa)
    if isempty(gotos)
        parents = [Int[] for _ = 1:length(ssa)]
    else
        parents = inferblocks(gotos)
    end
    return BlockInfo(gotos, parents)
end

isgoto(item) = isa(item, Expr) && (item.head == :(goto) || item.head == :(gotoif))
function getgotos(ssa)
    gotos = Dict{Int,Int}()
    for (i, line) in enumerate(ssa)
        !isa(line, Array) && continue
        for item in line
            if isgoto(item)
                gotos[i] = item.args[1]
            end
        end
    end
    return gotos
end

# inferblocks

function inferblocks(gotos::Dict{Int,Int})
    parents = initialblocks(gotos)
    parents = expandblocks(parents)
    parents = mergeblocks(parents, gotos)
    return parents
end

function initialblocks(gotos::Dict{Int,Int})
    # lists of which block(s) each ssa index belong to
    ssaparents = [Int[] for _ = 1:maximum(values(gotos)) + 1]
    #ssaparents = [Int[] for _ = 1:20]
    for (origin, target) in gotos
        a, b = minmax(origin, target)
        for j = a:b
            push!(ssaparents[j], b) # named by block end index
        end
    end
    ssaparents = unique.(sort.(ssaparents, rev=true))
    return ssaparents
end

function expandblocks(ssaparents::Array{Array{Int,1},1})
    Nblocks = length(unique(vcat(ssaparents...)))
    # make sure branching (outward) can reach the intended target
    for _ = 1:Nblocks
        for (i, parents) in enumerate(ssaparents)
            length(parents) < 1 && continue
            for (level, parent) in enumerate(parents)
                length(ssaparents[parent]) < level && continue
                # a parent is named according to where it ends,
                requiredparent = ssaparents[parent][level] # so parent in this line refers to where the parent ends
                if parent != requiredparent
                    # expand required block (upward) to dominate current block
                    insert!(ssaparents[i], level, requiredparent)
                end
            end
        end
    end

    return ssaparents
end

function mergeblocks(ssaparents::Array{Array{Int,1},1}, gotos::Dict{Int,Int})
    # gotos are named according to max(origin, target) for sorting and block expanding purposes,
    # this means backward jumps (unlike forward jumps) will be named
    # according to "origin" instead of "target" and will thus produce extra blocks.
    # So conceptually; this function renames loop blocks to be named according to start
    # index and removes extra blocks produced by a "continue" jump.
    backjumps = filter(kv -> kv[2] < kv[1], gotos)
    for (i, parents) in enumerate(ssaparents)
        for (j, name) in enumerate(parents)
            if haskey(backjumps, name)
                ssaparents[i][j] = backjumps[name] # rename
            end
        end
    end
    return unique.(ssaparents) # merge
end

# translate

function goto2brlevel(bi::BlockInfo, i::Int, target::Int)
    parents = bi.parents[i]
    !(target in parents) && error("Target block not accessible: ssa index $i is not a child of target block $target")

    for (level, parent) in enumerate(reverse(parents))
        target == parent && return level - 1
    end
end

translate_gotos(bi::BlockInfo, i::Int, items::Array) = translate_gotos.((bi,), (i,), items)
function translate_gotos(bi::BlockInfo, i::Int, item)
    if isa(item, Expr) && item.head == :(goto)
        brlevel = goto2brlevel(bi, i, item.args[1])
        return "(br $brlevel)"
    elseif isa(item, Expr) && item.head == :(gotoif)
        brlevel = goto2brlevel(bi, i, item.args[1])
        return "br_if $brlevel"
    else
        return item
    end
end

function addblocks(bi::BlockInfo, wat::Array)
    # where and what blocks to insert into wat
    blocks = [[] for _ = 1:length(wat)]
    backtargets = values(filter(kv -> kv[2] < kv[1], bi.gotos)) # v<k means backjump (target less than origin)
    for (i, targets) in enumerate(bi.parents)
        for (level, target) in enumerate(targets)
            prevmaxlevel = length(bi.parents[i - 1])
            nextmaxlevel = length(bi.parents[i + 1])
            if prevmaxlevel < level || bi.parents[i - 1][level] != target
                target in backtargets ? push!(blocks[i], "(loop") : push!(blocks[i], "(block")
            elseif nextmaxlevel < level || bi.parents[i + 1][level] != target
                target in backtargets ? push!(blocks[i + 1], ")") : push!(blocks[i], ")") # loop end at end of line i aka first at i+1
            end
        end
    end

    # add them
    for (i, v) in enumerate(blocks)
        length(v) < 1 && continue
        if isnothing(wat[i])
            wat[i] = [spacedjoin(v)]
        elseif isa(wat[i], TypedSlot)
            # wat[i] = [v,wat[i]]
            wat[i] = spacedjoin(v)
        else
            pushfirst!(wat[i], spacedjoin(v))
        end
    end
    return wat
end


text:
abstract type Regularizer end

# regularizers are separable
function reg(R::Regularizer, regweights::Array{Float64,1}, theta::Array{Float64,1})
    s = 0.0
    for i=1:length(theta)
        if regweights[i] > 0 # test otherwise 0*Inf = NaN. We would want 0
            s += regweights[i]*reg(R, theta[i])
        end
    end
    return s
end

# construct prox of vector by applying prox elementwise
# include regweights
function prox(R::Regularizer, gamma::Float64,
              regweights::Array{Float64, 1},
              v::Array{Float64,1})
    n = length(v)
    p = zeros(v)
    for i=1:n
        if regweights[i] > 0
            p[i] = prox(R,  gamma/regweights[i], v[i])
        else
            p[i] = v[i]
        end
    end
    return p
end

##############################################################################

#########################################
# L1 Regularizer
#########################################
"""
`L1Reg()` constructs the L1 regularizer. Use with `Model()`
"""
struct L1Reg <: Regularizer end
reg(R::L1Reg, a::Float64) = abs(a)
cvxreg(R::L1Reg, a) = abs(a)

# returns arg min ( R(x) + gamma*\norm(x-v)^2 )
function prox(R::L1Reg, gamma::Float64, v::Float64)
    s = 0.5/gamma
    if v > s
        return v-s
    elseif v < -s
        return v+s
    end
    return 0
end

#########################################
# L2 Regularizer
#########################################
"`L2Reg()` constructs the L2 regularizer. Use with `Model()`"
struct L2Reg <: Regularizer end
reg(R::L2Reg, a::Float64) = a*a
prox(R::L2Reg, gamma::Float64, v::Float64) = v*gamma/(1+gamma)
cvxreg(R::L2Reg, a) = a*a

#########################################
# Sqrt Regularizer
#########################################
"`SqrtReg()` computes the square root regularizer `sqrt(abs(x_i))`. Use with `Model()`"
struct SqrtReg <: Regularizer end
reg(R::SqrtReg, u::Float64) = sqrt(abs(u))

# horrific approach to computing prox
function prox(R::SqrtReg, ga::Float64, vv::Float64)
    if vv<0
        return -prox(R, ga,-vv)
    end

    ff(x) = sqrt(abs(x)) + ga*(x-vv)*(x-vv)

    g = complex(ga)
    v = complex(vv)
    # roots of  4 g x^3 - 4  g v  x  + 1 = 0
    z1 = (2*g*v)/(3^(1/3)*(sqrt(3)*sqrt(27*g^4 - 64*g^6*v^3) - 9*g^2)^(1/3)) + (sqrt(3)*sqrt(27*g^4 - 64*g^6*v^3) - 9*g^2)^(1/3)/(2*3^(2/3)*g)
    z2 = -((1 + im*sqrt(3))*g*v)/(3^(1/3)*(sqrt(3)*sqrt(27*g^4 - 64*g^6*v^3) - 9*g^2)^(1/3)) - ((1 - im*sqrt(3))*(sqrt(3)*sqrt(27*g^4 - 64*g^6*v^3) - 9*g^2)^(1/3))/(4*3^(2/3)*g)
    z3 = -((1 - im*sqrt(3))*g*v)/(3^(1/3)*(sqrt(3)*sqrt(27*g^4 - 64*g^6*v^3) - 9*g^2)^(1/3)) - ((1 + im*sqrt(3))*(sqrt(3)*sqrt(27*g^4 - 64*g^6*v^3) - 9*g^2)^(1/3))/(4*3^(2/3)*g)

    x1 = z1*z1
    x2 = z2*z2
    x3 = z3*z3

    xmin = 0.0
    fmin = ff(xmin)
    if abs(imag(x1)) < 1e-10
        t = ff(real(x1))
        if t < fmin
            fmin = t
            xmin = real(x1)
        end
    end
    if abs(imag(x2)) < 1e-10
        t = ff(real(x2))
        if t < fmin
            fmin = t
            xmin = real(x2)
        end
    end
    if abs(imag(x3)) < 1e-10
        t = ff(real(x3))
        if t < fmin
            fmin = t
            xmin = real(x3)
        end
    end
    return xmin
end

#########################################
# Nonneg Regularizer
#########################################
"`NonnegReg()` is the nonnegative regularizer. Use with `Model()`"
struct NonnegReg <: Regularizer end
function reg(R::NonnegReg, a)
    if a >= 0
        return 0
    end
    return Inf
end
function prox(R::NonnegReg, gamma::Float64, v::Float64)
    if v>0
        return v
    end
    return 0
end


text:
function hashjoin(A::Array, ja::Int, B::Array, jb::Int)
    M = Dict(t[jb] => filter(l -> l[jb] == t[jb], B) for t in B)
    return collect([a, b] for a in A for b in get(M, a[ja], ()))
end

table1 = [(27, "Jonah"),
          (18, "Alan"),
          (28, "Glory"),
          (18, "Popeye"),
          (28, "Alan")]
table2 = [("Jonah", "Whales"),
          ("Jonah", "Spiders"),
          ("Alan", "Ghosts"),
          ("Alan", "Zombies"),
          ("Glory", "Buffy")]

for r in hashjoin(table1, 2, table2, 1)
    println(r)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module Predix

end


text:
using ParameterizedFunctions, Test
using OrdinaryDiffEq, Calculus, ForwardDiff

f = @ode_def begin
    dx = -a
    dy = b
end a b

cb = ContinuousCallback((u,t,i) -> u[1], (integrator)->(println("Stopped.");integrator.p[2]=zero(integrator.p[2])))
function test_f(p)
  prob = ODEProblem(f,eltype(p).([1.0,0.0]),eltype(p).((0.0,1.0)),copy(p))
  integrator = init(prob,Tsit5(),abstol=1e-14,reltol=1e-14,callback=cb)
  step!(integrator)
  solve!(integrator).u[end]
end
p = [2.0, 1.0]
findiff = Calculus.finite_difference_jacobian(test_f,p)
fordiff = ForwardDiff.jacobian(test_f,p)
@test findiff ≈ fordiff

f2 = @ode_def begin
    dx = -x
    dy = b
end a b

function test_f2(p)
  prob = ODEProblem(f2,eltype(p).([1.0,0.0]),eltype(p).((0.0,1.0)),copy(p))
  integrator = init(prob,Tsit5(),abstol=1e-14,reltol=1e-14,callback=cb)
  step!(integrator)
  solve!(integrator).u[end]
end
p = [2.0, 1.0]
findiff = Calculus.finite_difference_jacobian(test_f2,p)
fordiff = ForwardDiff.jacobian(test_f2,p)
@test findiff ≈ fordiff


text:
"Autumn SubExpressions"
module SubExpressions

using ..AExpressions
using ..Util

export SubExpr,
       subexpr,
       resolve,
       update,
       subexprdfs,
       subexprs

"Subexpression of `parent::AE` indicated by pointer `p::P`"
struct SubExpr{AExpr, P}
  parent::AExpr
  pointer::P
end

subexpr(aexpr, id) = SubExpr(aexpr, id)
subexpr(aexpr, id::Integer) = SubExpr(aexpr, [id])
subexpr(::SubExpr) = error("Cannot take subexpression of subexpression, yet")

"Returns subexpressions that are children"
function AExpressions.args(subexpr::SubExpr)
  q = g(resolve(subexpr))
  [SubExpr(subexpr.parent, append(subexpr.pointer, i)) for i = 1:length(q)]
end
g(aex::AExpr) = aex.args
g(x) = []
append(xs::AbstractVector, x) = [xs; x]

"Resolve Value pointed to by `subexpr`"
function resolve(subexpr::SubExpr)
  ex = subexpr.parent
  for id in subexpr.pointer
    ex = arg(ex, id)
  end
  ex
end

"""Update subexpr.parent such that `subexpr` is `newexpr`

```
prog = au\"\"\"
(program
  (: x Int)
  (= x 3)
  (= y (initnext (+ 1 2) (/ 3 this)))
  (= z (f 1 2))
)\"\"\"
.push(w)

subexpr_ = subexpr(prog, [2])

prog2 = au\"\"\"(= x 4000)\"\"\"
update(subexpr_, prog2)
```
\"\"\"

```
"""
function update(subexpr::SubExpr, newexpr)
  function subchild(expr, pos)
    # @show pos, subexpr.pointer
    # @show pos == subexpr.pointer
    pos == subexpr.pointer ? newexpr : expr
  end
  postwalkpos(subchild, subexpr.parent)
end

"Depth first search of subexpressions"
function subexprdfs(subexpr::SubExpr)
  s = SubExpr[]
  res = SubExpr[]
  push!(s, subexpr)
  discovered = Set{SubExpr}()
  while !isempty(s)
    v = pop!(s)
    push!(res, v)
    if v ∉ discovered
      push!(discovered, v)
      for subexpr in args(v)
        push!(s, subexpr)
      end
    end
  end
  res

  # What's wrong with this
  # 1. we're not using the result
  # 2. push pop needs better data structure
  # 3. not iterator friedly
  # 4. will be hashing these sexpr, which is slow
  # 5. args(v) undefiined on most things ERR
  # 6. args(v) doesn't returb subexpression ERR
end

# Use DFS by default
subexprs(x) = subexprdfs(x)

# Base.iterate(S::SubExpr, state)

subexprdfs(aexpr::AExpr) =
  subexprdfs(subexpr(aexpr, Int[]))

Base.show(io::IO, subexpr::SubExpr) =
  print(io, "Subexpression @ ", subexpr.pointer, ":\n", subexpr.parent, " => \n", resolve(subexpr), "\n")

end


text:
mutable struct BSpline
	degree   :: Int
	numKnots :: Int
	lower    :: Float64
	upper    :: Float64
	knots    :: Array{Float64,1}

	# constructor with equally spaced interior knots
	function BSpline(nKnots,deg,lb,ub)
		knots = zeros(nKnots+2*deg)
		if (nKnots < 2*(deg+1)-1)
			throw(ArgumentError("need at least 2*(deg+1) - 1 = $(2*(deg+1) - 1) knots"))
		end
		if deg < 0
			throw(ArgumentError("degree must be non-negative integer"))
		end
		if lb > ub
			throw(ArgumentError("require lb < ub"))
		end

		h = (ub - lb) / (nKnots - 1)
		for i=1:length(knots)
			if i < deg + 1
				knots[i] = lb
			elseif i > nKnots + deg
				knots[i] = ub
			else
				knots[i] = lb + (i-deg-1)*h
			end
		end
		# sp = sparsevec(zeros(num_coefs))
		new(deg,nKnots,lb,ub,knots)
	end

	# constructor with user supplied interior knots
	function BSpline(knots::Array,deg::Integer)
		if !issorted(knots)
			throw(ArgumentError("knots must be sorted"))
		end
		lb = knots[1]
		ub = knots[end]
		# extend knot vector by degree on both sides
		knots_new = [ [lb for i=1:deg]; knots; [ub for i=1:deg] ]
		numKnots = length(knots)	# number of INTERIOR knots
		new(deg,numKnots,lb,ub,knots_new)
	end

end

function show(io::IO, b::BSpline)
	print(io,"BSpline object with\n")
	print(io,"degree: $(b.degree)\n")
	print(io,"number of knots: $(b.numKnots)\n")
	print(io,"[lower,upper]: [$(b.lower),$(b.upper)]\n")
	print(io,"knot vector: $(b.knots)\n")
end

function getNumKnots(b::BSpline) return b.numKnots end
function getNumCoefs(b::BSpline) return b.numKnots + b.degree -1 end
function getCoefs(b::BSpline) return b.knots[ (b.degree+1):(length(b.knots)-b.degree)] end
function getDegree(b::BSpline) return b.degree end

# """
# 	compute value of the ``d + 1`` B-splines at ``x``

# returns

# ```math
# b = \right(B_{\mu-d,d}(x),\dots,B_{\mu,d}(x)\left)^T
# ```
# """
# function getBasis2(x::Float64,b::BSpline)

# 	num_nodes = getNumCoefs(b)
# 	deg       = b.degree

# 	B = zeros(num_nodes)  # Bspline

# end

#
function getBasis(x::Float64,b::BSpline)

	num_nodes = getNumCoefs(b)
	deg       = b.degree

	# tmp
	d = 0.0
	e = 0.0

	# create a basis function
	# return a colvec because currntly only CSC format
	# a 1-row sparse matrix is dense
	# bs = spzeros(num_nodes,1)
	bs = spzeros(num_nodes,1)

	# check x
	if x < b.lower
		@warn("x < lb: x=$x, lb=$(b.lower). be careful!")
	elseif x > b.upper
		@warn("x > ub: x=$x, ub=$(b.upper). be careful!")
	end

	# get mu s.t. knot_mu < knot_{mu+1} and x in [knot_mu, knot_{mu+1})
	# i.e. get the index of the lower knot in the active knot span

	# fix bound behaviour
	if x <= b.lower
		mu = deg+1
	elseif x>=b.upper
		mu = num_nodes
	else
		mu = searchsortedlast(b.knots,x)
	end

	# set 0-degree basis function
	# 0-deg basis is an indicator function
	# that is 1.0 in the active knot span and 0.0 else.
	bs[mu] = 1.0

	# loop over degrees
	for k=1:deg

		# loop over basis functions
		for j in mu-k:mu

			# take care of "division by zero" issue
			# dividing by zero must return 0.0
			if j+k <= deg +1
				d = 0.0
			elseif b.knots[j+k]-b.knots[j] == 0
				d = 0.0
			else
				d = bs[j] * (x - b.knots[j]) / (b.knots[j+k]-b.knots[j])
			end

			if j+1 >= num_nodes+1
				e = 0.0
			elseif b.knots[j+k+1]-b.knots[j+1] == 0
				e = 0.0
			else
				e = bs[j+1] * (b.knots[j+k+1] - x)/(b.knots[j+k+1]-b.knots[j+1])
			end

			bs[j] = d + e

		end

	end
	return bs
end

# vector of points
function getBasis(x::Vector{Float64},b::BSpline)

	n = length(x)
	num_nodes = getNumCoefs(b)
	deg       = b.degree

	# if you evaluate a degree 1 basis at its interior knots ("coefs"), you get identity
	if (deg == 1) && (n == num_nodes) && (sum(abs.(x .- getCoefs(b))) < 1e-8)
		return sparse(I*1.0,n,n)
	else
	# sort x?

		# tmp
		d = 0.0
		e = 0.0

		# create a basis function
		# return a colvec because currntly only CSC format
		# a 1-row sparse matrix is dense
		bs = spzeros(n,num_nodes)

		for xi in 1:n

			# check x
			if x[xi] < b.lower
				@warn("x < lb: x=$(x[xi]), lb=$(b.lower). be careful!")
			elseif x[xi] > b.upper
				@warn("x > ub: x=$(x[xi]), ub=$(b.upper). be careful!")
			end

			# get mu s.t. knot_mu < knot_{mu+1} and x in [knot_mu, knot_{mu+1})
			# i.e. get the index of the lower knot in the active knot span

			# fix bound behaviour
			if x[xi] <= b.lower
				mu = deg+1
			elseif x[xi]>=b.upper
				mu = num_nodes
			else
				mu = searchsortedlast(b.knots,x[xi])
			end

			# set 0-degree basis function
			# 0-deg basis is an indicator function
			# that is 1.0 in the active knot span and 0.0 else.
			bs[xi,mu] = 1.0

			# loop over degrees
			for k=1:deg

				# loop over basis functions
				for j in mu-k:mu

					# take care of "division by zero" issue
					# dividing by zero must return 0.0
					if j+k <= deg +1
						d = 0.0
					elseif b.knots[j+k]-b.knots[j] == 0
						d = 0.0
					else
						d = bs[xi,j] * (x[xi] - b.knots[j]) / (b.knots[j+k]-b.knots[j])
					end

					if j+1 >= num_nodes+1
						e = 0.0
					elseif b.knots[j+k+1]-b.knots[j+1] == 0
						e = 0.0
					else
						e = bs[xi,j+1] * (b.knots[j+k+1] - x[xi])/(b.knots[j+k+1]-b.knots[j+1])
					end
					bs[xi,j] = d + e
				end

			end
		end
		return bs

	end

end


text:
module TestConvert
    using Base.Test
    using CategoricalArrays

    pool = CategoricalPool([1, 2, 3])
    @test convert(CategoricalPool{Int, CategoricalArrays.DefaultRefType}, pool) === pool
    @test convert(CategoricalPool{Int}, pool) === pool
    @test convert(CategoricalPool, pool) === pool
    convert(CategoricalPool{Float64, UInt8}, pool)
    convert(CategoricalPool{Float64}, pool)
    convert(CategoricalPool, pool)

    v1 = CategoricalValue(1, pool)
    v2 = CategoricalValue(2, pool)
    v3 = CategoricalValue(3, pool)

    convert(Int32, v1)
    convert(Int32, v2)
    convert(Int32, v3)

    convert(UInt8, v1)
    convert(UInt8, v2)
    convert(UInt8, v3)

    @test get(v1) === 1
    @test get(v2) === 2
    @test get(v3) === 3

    @test promote(1, v1) === (1, 1)
    @test promote(1.0, v1) === (1.0, 1.0)
    @test promote(0x1, v1) === (1, 1)
end


text:
function predict_new_rbf_deriv_deriv(a_all,a_layer1,x_new,obs,neurons,n_per_part,inds_all,x_old,vars,cc,var_seq)

    layer1=zeros(obs,neurons)
    for i=1:neurons
        x1=copy(x_old[inds_all[n_per_part[i]+1:n_per_part[i+1]],:])
        obs1=size(x1,1)
        phi1 = calc_phi_deriv_deriv(x1,x_new,obs1,obs,cc,vars,var_seq)
        # if sum(abs.(a_all[i,1:obs1]))>0
        #     layer1[:,i].=phi1*a_all[i,1:obs1]
        # end
        layer1[:,i].=phi1*a_all[i]
    end

    predl1=[layer1 ones(obs)]*a_layer1

    return predl1

end


text:
"""
     SparseInSto

Storage type for in-place sparse calculations
"""
type SparseInSto
     Xh
     Xj
     Yh
     Yj
     Zh
     Zj
     nx
end
SparseInSto() = SparseInSto([],[],[],[],[],[],[])

"""
     Preconditioner(h,X,P;jac="User")

Directly applies inverse preconditioning matrix.
"""
function Preconditioner(h,X,P;jac="User")
  J = h(X,P)
  #println("J:   ",J)
  if (length(X)>1)
    #println("mid.(J)", mid.(J))
    Y = inv(mid.(J))
  else
    Y = 1.0/(mid(J[1]))
  end
  return Y
end

"""
     Sparse_Forward_Elimination!(x::Vector{Tq},L::SparseMatrixCSC{Tv,Ti},
                                 b::Vector{Tq},nx::Int64) where {Tv,Tq,Ti}

Solves for `Lx=b` via forward elimination. A must be a lower triangular sparse
matrix of CSC format.
"""
function Sparse_Forward_Elimination!(x::Vector{Tq},L::SparseMatrixCSC{Tv,Ti},
                                     b::Vector{Tq},nx::Int64) where {Tv,Tq,Ti}
     # converts to CSR (expensive if dense ... much cheaper calc if sparse)
     A::SparseMatrixCSC{Tv,Ti} = L'
     # standard row-oriented forward elimination
     x[1,:] = b[1,:]/A.nzval[1]
     for i=2:(nx)
          for k=(A.colptr[i]):(A.colptr[i+1]-2)
               b[i,:] = b[i,:] - A.nzval[k]*x[A.rowval[k],:]
          end
          x[i,:] = b[i,:]/A.nzval[A.colptr[i+1]-1]
     end
end

"""
     Sparse_Forward_Elimination!(x::SparseMatrixCSC{Tq,Ti},U::SparseMatrixCSC{Tv,Ti},
                                  b::SparseMatrixCSC{Tq,Ti},nx::Int64)

Solves for `Lx=b` via forward elimination. A must be a lower triangular sparse
matrix of CSC format.
"""
# TO DO: Improve row access for b
function Sparse_Forward_Elimination!(x::SparseMatrixCSC{Tq,Ti},L::SparseMatrixCSC{Tv,Ti},
                                     b::SparseMatrixCSC{Tq,Ti},nx::Int64) where {Tv,Tq,Ti}
     # converts to CSR (expensive if dense ... much cheaper calc if sparse)
     A::SparseMatrixCSC{Tv,Ti} = L'
     # standard row-oriented forward elimination
     x[1,:] = b[1,:]/A.nzval[1]
     for i=2:(nx)
          for k=(A.colptr[i]):(A.colptr[i+1]-2)
               b[i,:] = b[i,:] - A.nzval[k]*x[A.rowval[k],:]
          end
          x[i,:] = b[i,:]/A.nzval[A.colptr[i+1]-1]
     end
end

"""
     Sparse_Back_Substitution!(x::Vector{Tq},U::SparseMatrixCSC{Tv,Ti},
                              b::Vector{Tq},nx::Int64)

Solves for `Ux=b` via backsubstitution. A must be a upper triangular sparse
matrix of CSC format.
"""
function Sparse_Back_Substitution!(x::Vector{Tq},U::SparseMatrixCSC{Tv,Ti},
                                   b::Vector{Tq},nx::Int64) where {Tv,Tq,Ti}
     # converts to CSR (expensive if dense ... much cheaper calc if sparse)
     A::SparseMatrixCSC{Tv,Ti} = U'
     # standard row-oriented back substituion
     x[end,:] = b[end,:]/A.nzval[end]
     for i=(nx-1):-1:1
          for k=(A.colptr[i+1]-1):-1:(A.colptr[i]+1)
               b[i,:] = b[i,:] - A.nzval[k]*x[A.rowval[k],:]
          end
          x[i,:] = b[i,:]/A.nzval[A.colptr[i]]
     end
end

"""
     Sparse_Back_Substitution!(x::SparseMatrixCSC{Tq,Ti},U::SparseMatrixCSC{Tv,Ti},
                               b::SparseMatrixCSC{Tq,Ti},nx::Int64)

Solves for `Ux=b` via backsubstitution. A must be a upper triangular sparse
matrix of CSC format.
"""
# TO DO: Improve row access for b
function Sparse_Back_Substitution!(x::SparseMatrixCSC{Tq,Ti},U::SparseMatrixCSC{Tv,Ti},
                                   b::SparseMatrixCSC{Tq,Ti},nx::Int64) where {Tv,Tq,Ti}
     # converts to CSR (expensive if dense ... much cheaper calc if sparse)
     A::SparseMatrixCSC{Tv,Ti} = U'
     # standard row-oriented back substituion
     x[end,:] = b[end,:]/A.nzval[end]
     for i=(nx-1):-1:1
          for k=(A.colptr[i+1]-1):-1:(A.colptr[i]+1)
               b[i,:] = b[i,:] - A.nzval[k]*x[A.rowval[k],:]
          end
          x[i,:] = b[i,:]/A.nzval[A.colptr[i]]
     end
end

"""
     Sparse_Precondition!(H::Vector{Ta},J::SparseMatrixCSC{Ta,Ti},
                          P::SparseMatrixCSC{Tp,Ti},st::SparseInSto)

Preconditions the H & J to inv(P)H and inv(P)J using a sparse LU factorization
method with full pivoting. J and P must be of size nx-by-nx and H must be of
size nx. st is the inplace storage type.
"""
function Sparse_Precondition!(H::Vector{Ta},J::SparseMatrixCSC{Ta,Ti},
                              P::SparseMatrixCSC{Tp,Ti},st::SparseInSto) where {Ta,Tp,Ti}

     # generate LU-PDQ factorization
     lu = lufact(P)

     # solves Lz = PDH for z
     Sparse_Forward_Elimination!(st.Zh,lu[:L],(lu[:Rs].*H)[lu[:p]],st.nx)
     Sparse_Forward_Elimination!(st.Zj,lu[:L],(sparse(convert.(Ta,lu[:Rs])).*J)[lu[:p],:],st.nx)

     # solves Uy = z for y
     Sparse_Back_Substitution!(st.Yh,lu[:U],st.Zh,st.nx)
     Sparse_Back_Substitution!(st.Yj,lu[:U],st.Zj,st.nx)

     # solves x = Qy
     st.Xh = st.Yh[lu[:q]]
     st.Xj = st.Yj[lu[:q],:]

     # stores the preconditioned matrices back in place
     H[:] = st.Xh[:,1]
     J[:] = st.Xj[:,1:(st.nx)]
end

"""
     Dense_Precondition!(H::Vector{Ta},J::Array{Ta,2},P::Array{Tp,2})

Preconditions the H & J to inv(P)H and inv(P)J using a sparse LU factorization
method with full pivoting. J and P must be of size nx-by-nx and H must be of
size nx. st is the inplace storage type.
"""
function Dense_Precondition!(H::Vector{Ta},J::Array{Ta,2},P::Array{Tp,2},nx) where {Ta,Tp}

     # generate PLU factorization
     lu = lufact(P)
     JHmerge = [J H]
     yH = lu[:L]\JHmerge[lu[:p],:]
     xH = lu[:U]\yH
     # stores the preconditioned matrices back in place
     H[:] = xH[:,nx+1]
     J[:] = xH[:,1:nx]
end


text:
"""
Various helpers for data representation.
"""
module Data
  export INT_SIZE, CHAR_SIZE, BLOCK_SIZE
  export string_size

  const BLOCK_SIZE = 400
  const INT_SIZE = sizeof(Int)
  const CHAR_SIZE = 1
  # sizeof(Char) == 4, but most alphabet is 1. This will
  # also cause issue when reading from penis buffer as it will read 4 bytes
  # (hence getstring() hack).
  #
  # TODO#charsize: Implement either
  # 1. (4 - charsize)-byte 0x00 padding (naive).
  # 2. proper utf handling.

  """
  Size of a string, as formatted and stored on disk by db, in bytes.
  Each string is represented as an array of chars, preceded by an Int
  representing the number of chars:

      |3|A|B|C|

  REF#stringsize
  """
  function string_size(str_length)
    (CHAR_SIZE * str_length) + INT_SIZE
  end
end


text:
function  producecode(app,mesh)

# search compilers and set options
app = Gencode.setcompilers(app);

# generate input files and store them in datain folder
app, mesh, master, dmd = preprocessing(app,mesh);

# generate source codes and store them in app folder
Gencode.gencode(app);

# compile source codes to build an executable file and store it in app folder
compilerstr = Gencode.compilecode(app);

return compilerstr,app,mesh,master,dmd

end

max_stars_repo_name: rloekvh/Exasim

max_issues_repo_name: rloekvh/Exasim

max_forks_repo_name: rloekvh/Exasim


text:
@testset "DeepONet Training Accuracy" begin
    ϵ = Burgers.train_don()

    @test ϵ < 0.4
end


text:
__precompile__()

module AIMACore
include("errors.jl")
include("utils.jl")
include("agents.jl")
include("search.jl")
end # module


text:
# # von Mises plasticity
#
# ![Shows the von Mises stress distribution in a cantilever beam.](plasticity.png)
# *Figure 1.* A coarse mesh solution of a cantilever beam subjected to a load
# causing plastic deformations. The initial yield limit is 200 MPa but due to
# hardening it increases up to approximately 240 MPa.

#
# ## Introduction
#
# This example illustrates the use of a nonlinear material model in Ferrite.
# The particular model is von Mises plasticity (also know as J₂-plasticity) with
# isotropic hardening. The model is fully 3D, meaning that no assumptions like *plane stress*
# or *plane strain* are introduced.
#
# Also note that the theory of the model is not described here, instead one is
# referred to standard textbooks on material modeling.
#
# To illustrate the use of the plasticity model, we setup and solve a FE-problem
# consisting of a cantilever beam loaded at its free end. But first, we shortly
# describe the parts of the implementation deadling with the material modeling.

#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`plasticity.ipynb`](@__NBVIEWER_ROOT_URL__/examples/plasticity.ipynb)
#-

# ## Material modeling
# This section describes the `struct`s and methods used to implement the material
# model

# ### Material parameters and state variables
#
# Start by loading some necessary packages
using Ferrite, Tensors, SparseArrays, LinearAlgebra, Printf

# We define a J₂-plasticity-material, containing material parameters and the elastic
# stiffness Dᵉ (since it is constant)
struct J2Plasticity{T, S <: SymmetricTensor{4, 3, T}}
    G::T  # Shear modulus
    K::T  # Bulk modulus
    σ₀::T # Initial yield limit
    H::T  # Hardening modulus
    Dᵉ::S # Elastic stiffness tensor
end;

# Next, we define a constructor for the material instance.
function J2Plasticity(E, ν, σ₀, H)
    δ(i,j) = i == j ? 1.0 : 0.0 # helper function
    G = E / 2(1 + ν)
    K = E / 3(1 - 2ν)

    Isymdev(i,j,k,l) = 0.5*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k)) - 1.0/3.0*δ(i,j)*δ(k,l)
    temp(i,j,k,l) = 2.0G *( 0.5*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k)) + ν/(1.0-2.0ν)*δ(i,j)*δ(k,l))
    Dᵉ = SymmetricTensor{4, 3}(temp)
    return J2Plasticity(G, K, σ₀, H, Dᵉ)
end;

#md # !!! note
#md #     Above, we defined a constructor `J2Plasticity(E, ν, σ₀, H)` in terms of the more common
#md #     material parameters ``E`` and ``ν`` - simply as a convenience for the user.
#md #

# Define a `struct` to store the material state for a Gauss point.
struct MaterialState{T, S <: SecondOrderTensor{3, T}}
    ## Store "converged" values
    ϵᵖ::S # plastic strain
    σ::S # stress
    k::T # hardening variable
end

# Constructor for initializing a material state. Every quantity is set to zero.
function MaterialState()
    return MaterialState(
                zero(SymmetricTensor{2, 3}),
                zero(SymmetricTensor{2, 3}),
                0.0)
end

# For later use, during the post-processing step, we define a function to
# compute the von Mises effective stress.
function vonMises(σ)
    s = dev(σ)
    return sqrt(3.0/2.0 * s ⊡ s)
end;

# ## Constitutive driver
#
# This is the actual method which computes the stress and material tangent
# stiffness in a given integration point.
# Input is the current strain and the material state from the previous timestep.
function compute_stress_tangent(ϵ::SymmetricTensor{2, 3}, material::J2Plasticity, state::MaterialState)
    ## unpack some material parameters
    G = material.G
    H = material.H

    ## We use (•)ᵗ to denote *trial*-values
    σᵗ = material.Dᵉ ⊡ (ϵ - state.ϵᵖ) # trial-stress
    sᵗ = dev(σᵗ)         # deviatoric part of trial-stress
    J₂ = 0.5 * sᵗ ⊡ sᵗ   # second invariant of sᵗ
    σᵗₑ = sqrt(3.0*J₂)   # effective trial-stress (von Mises stress)
    σʸ = material.σ₀ + H * state.k # Previous yield limit

    φᵗ  = σᵗₑ - σʸ # Trial-value of the yield surface

    if φᵗ < 0.0 # elastic loading
        return σᵗ, material.Dᵉ, MaterialState(state.ϵᵖ, σᵗ, state.k)
    else # plastic loading
        h = H + 3G
        μ =  φᵗ / h   # plastic multiplier

        c1 = 1 - 3G * μ / σᵗₑ
        s = c1 * sᵗ           # updated deviatoric stress
        σ = s + vol(σᵗ)       # updated stress

        ## Compute algorithmic tangent stiffness ``D = \frac{\Delta \sigma }{\Delta \epsilon}``
        κ = H * (state.k + μ) # drag stress
        σₑ = material.σ₀ + κ  # updated yield surface

        δ(i,j) = i == j ? 1.0 : 0.0
        Isymdev(i,j,k,l)  = 0.5*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k)) - 1.0/3.0*δ(i,j)*δ(k,l)
        Q(i,j,k,l) = Isymdev(i,j,k,l) - 3.0 / (2.0*σₑ^2) * s[i,j]*s[k,l]
        b = (3G*μ/σₑ) / (1.0 + 3G*μ/σₑ)

        Dtemp(i,j,k,l) = -2G*b * Q(i,j,k,l) - 9G^2 / (h*σₑ^2) * s[i,j]*s[k,l]
        D = material.Dᵉ + SymmetricTensor{4, 3}(Dtemp)

        ## Return new state
        Δϵᵖ = 3/2 * μ / σₑ * s # plastic strain
        ϵᵖ = state.ϵᵖ + Δϵᵖ    # plastic strain
        k = state.k + μ        # hardening variable
        return σ, D, MaterialState(ϵᵖ, σ, k)
    end
end

# ## FE-problem
# What follows are methods for assembling and and solving the FE-problem.
function create_values(interpolation)
    ## setup quadrature rules
    qr      = QuadratureRule{3,RefTetrahedron}(2)
    face_qr = QuadratureRule{2,RefTetrahedron}(3)

    ## create geometric interpolation (use the same as for u)
    interpolation_geom = Lagrange{3,RefTetrahedron,1}()

    ## cell and facevalues for u
    cellvalues_u = CellVectorValues(qr, interpolation, interpolation_geom)
    facevalues_u = FaceVectorValues(face_qr, interpolation, interpolation_geom)

    return cellvalues_u, facevalues_u
end;

# ### Add degrees of freedom
function create_dofhandler(grid, interpolation)
    dh = DofHandler(grid)
    dim = 3
    push!(dh, :u, dim, interpolation) # add a displacement field with 3 components
    close!(dh)
    return dh
end

# ### Boundary conditions
function create_bc(dh, grid)
    dbcs = ConstraintHandler(dh)
    ## Clamped on the left side
    dofs = [1, 2, 3]
    dbc = Dirichlet(:u, getfaceset(grid, "left"), (x,t) -> [0.0, 0.0, 0.0], dofs)
    add!(dbcs, dbc)
    close!(dbcs)
    return dbcs
end;

# ### Assembling of element contributions
#
# * Residual vector `r`
# * Tangent stiffness `K`
function doassemble(cellvalues::CellVectorValues{dim},
                    facevalues::FaceVectorValues{dim}, K::SparseMatrixCSC, grid::Grid,
                    dh::DofHandler, material::J2Plasticity, u, states, states_old, t) where {dim}
    r = zeros(ndofs(dh))
    assembler = start_assemble(K, r)
    nu = getnbasefunctions(cellvalues)
    re = zeros(nu)     # element residual vector
    ke = zeros(nu, nu) # element tangent matrix

    for (i, cell) in enumerate(CellIterator(dh))
        fill!(ke, 0)
        fill!(re, 0)
        eldofs = celldofs(cell)
        ue = u[eldofs]
        state = @view states[:, i]
        state_old = @view states_old[:, i]
        assemble_cell!(ke, re, cell, cellvalues, facevalues, grid, material,
                       ue, state, state_old, t)
        assemble!(assembler, eldofs, re, ke)
    end
    return K, r
end

# Compute element contribution to the residual and the tangent.
#md # !!! note
#md #     Due to symmetry, we only compute the lower half of the tangent
#md #     and then symmetrize it.
#md #
function assemble_cell!(Ke, re, cell, cellvalues, facevalues, grid, material,
                        ue, state, state_old, t)
    n_basefuncs = getnbasefunctions(cellvalues)
    reinit!(cellvalues, cell)

    for q_point in 1:getnquadpoints(cellvalues)
        ## For each integration point, compute stress and material stiffness
        ϵ = function_symmetric_gradient(cellvalues, q_point, ue) # Total strain
        σ, D, state[q_point] = compute_stress_tangent(ϵ, material, state_old[q_point])

        dΩ = getdetJdV(cellvalues, q_point)
        for i in 1:n_basefuncs
            δϵ = shape_symmetric_gradient(cellvalues, q_point, i)
            re[i] += (δϵ ⊡ σ) * dΩ # add internal force to residual
            for j in 1:i # loop only over lower half
                Δϵ = shape_symmetric_gradient(cellvalues, q_point, j)
                Ke[i, j] += δϵ ⊡ D ⊡ Δϵ * dΩ
            end
        end
    end
    symmetrize_lower!(Ke)

    ## Add traction as a negative contribution to the element residual `re`:
    for face in 1:nfaces(cell)
        if onboundary(cell, face) && (cellid(cell), face) ∈ getfaceset(grid, "right")
            reinit!(facevalues, cell, face)
            for q_point in 1:getnquadpoints(facevalues)
                dΓ = getdetJdV(facevalues, q_point)
                for i in 1:n_basefuncs
                    δu = shape_value(facevalues, q_point, i)
                    re[i] -= (δu ⋅ t) * dΓ
                end
            end
        end
    end

end

# Helper function to symmetrize the material tangent
function symmetrize_lower!(K)
    for i in 1:size(K,1)
        for j in i+1:size(K,1)
            K[i,j] = K[j,i]
        end
    end
end;

# Define a function which solves the FE-problem.
function solve()
    ## Define material parameters
    E = 200.0e9 # [Pa]
    H = E/20   # [Pa]
    ν = 0.3     # [-]
    σ₀ = 200e6  # [Pa]
    material = J2Plasticity(E, ν, σ₀, H)

    L = 10.0 # beam length [m]
    w = 1.0  # beam width [m]
    h = 1.0  # beam height[m]
    n_timesteps = 10
    u_max = zeros(n_timesteps)
    traction_magnitude = 1.e7 * range(0.5, 1.0, length=n_timesteps)

    ## Create geometry, dofs and boundary conditions
    n = 2
    nels = (10n, n, 2n) # number of elements in each spatial direction
    P1 = Vec((0.0, 0.0, 0.0))  # start point for geometry
    P2 = Vec((L, w, h))        # end point for geometry
    grid = generate_grid(Tetrahedron, nels, P1, P2)
    interpolation = Lagrange{3, RefTetrahedron, 1}() # Linear tet with 3 unknowns/node

    dh = create_dofhandler(grid, interpolation) # JuaFEM helper function
    dbcs = create_bc(dh, grid) # create Dirichlet boundary-conditions

    cellvalues, facevalues = create_values(interpolation)

    ## Pre-allocate solution vectors, etc.
    n_dofs = ndofs(dh)  # total number of dofs
    u  = zeros(n_dofs)  # solution vector
    Δu = zeros(n_dofs)  # displacement correction
    r = zeros(n_dofs)   # residual
    K = create_sparsity_pattern(dh); # tangent stiffness matrix

    ## Create material states. One array for each cell, where each element is an array of material-
    ## states - one for each integration point
    nqp = getnquadpoints(cellvalues)
    states = [MaterialState() for _ in 1:nqp, _ in 1:getncells(grid)]
    states_old = [MaterialState() for _ in 1:nqp, _ in 1:getncells(grid)]

    ## Newton-Raphson loop
    NEWTON_TOL = 1 # 1 N
    print("\n Starting Netwon iterations:\n")

    for timestep in 1:n_timesteps
        t = timestep # actual time (used for evaluating d-bndc)
        traction = Vec((0.0, 0.0, traction_magnitude[timestep]))
        newton_itr = -1
        print("\n Time step @time = $timestep:\n")
        update!(dbcs, t) # evaluates the D-bndc at time t
        apply!(u, dbcs)  # set the prescribed values in the solution vector

        while true; newton_itr += 1

            if newton_itr > 8
                error("Reached maximum Newton iterations, aborting")
                break
            end
            K, r = doassemble(cellvalues, facevalues, K, grid, dh, material, u,
                             states, states_old, traction);
            norm_r = norm(r[Ferrite.free_dofs(dbcs)])

            print("Iteration: $newton_itr \tresidual: $(@sprintf("%.8f", norm_r))\n")
            if norm_r < NEWTON_TOL
                break
            end

            apply_zero!(K, r, dbcs)
            Δu = Symmetric(K) \ r
            u -= Δu
        end

        ## Update the old states with the converged values for next timestep
        states_old .= states

        u_max[timestep] = max(abs.(u)...) # maximum displacement in current timestep
    end

    ## ## Postprocessing
    ## Only a vtu-file corrsponding to the last time-step is exported.
    ##
    ## The following is a quick (and dirty) way of extracting average cell data for export.
    mises_values = zeros(getncells(grid))
    κ_values = zeros(getncells(grid))
    for (el, cell_states) in enumerate(eachcol(states))
        for state in cell_states
            mises_values[el] += vonMises(state.σ)
            κ_values[el] += state.k*material.H
        end
        mises_values[el] /= length(cell_states) # average von Mises stress
        κ_values[el] /= length(cell_states)     # average drag stress
    end
    vtk_grid("plasticity", dh) do vtkfile
        vtk_point_data(vtkfile, dh, u) # displacement field
        vtk_cell_data(vtkfile, mises_values, "von Mises [Pa]")
        vtk_cell_data(vtkfile, κ_values, "Drag stress [Pa]")
    end

    return u_max, traction_magnitude
end

# Solve the FE-problem and for each time-step extract maximum displacement and
# the corresponding traction load. Also compute the limit-traction-load
u_max, traction_magnitude = solve();

# Finally we plot the load-displacement curve.
using Plots
plot(
    vcat(0.0, u_max),                # add the origin as a point
    vcat(0.0, traction_magnitude),
    linewidth=2,
    title="Traction-displacement",
    label=nothing,
    markershape=:auto
    )
ylabel!("Traction [Pa]")
xlabel!("Maximum deflection [m]")

# *Figure 2.* Load-displacement-curve for the beam, showing a clear decrease
# in stiffness as more material starts to yield.

## test the result                       #src
using Test                               #src
@test norm(u_max[end]) ≈ 0.254452645     #src

#md # ## [Raw source](@id plasticity-raw-code)
#md #
#md # Below follows a version of the program without any comments.
#md # The file is also available here: [plasticity.jl](plasticity.jl)
#md #
#md # ```julia
#md # @__CODE__
#md # ```


text:
using Test
include("../required.jl")

# Checkout 41a9e9 & load data
verify_jld2_data("data/jld2"; extract_recent=true)
@load joinpath("data", "jld2", "raw-distributions.jld2") distributions spenames specommon speindex
Y = Ymatrix(distributions)
inds_obs = _indsobs(Y)
Yobs = _Yobs(Y, inds_obs)
env_df = CSV.read(
    joinpath("data", "proc", "distributions_env_full.csv"),
    DataFrame;
    header=true,
    delim="\t",
)
# Save new data
new = (
    distributions=distributions,
    spenames=spenames,
    Y=Y,
    inds_obs=inds_obs,
    Yobs=Yobs,
    env_df=env_df,
)
# Checkout a1cbac, re-run, & save old data
old = (
    distributions=distributions,
    spenames=spenames,
    Y=Y,
    inds_obs=inds_obs,
    Yobs=Yobs,
    env_df=env_df,
)

# Investigate sites with NA for land cover
allowmissing!(old.env_df)
for col in eachcol(old.env_df)
    replace!(col, NaN => missing)
end
@test nrow(old.env_df[old.inds_obs, :]) == nrow(dropmissing(old.env_df[old.inds_obs, :]))

# Investigate number of occurrence differece
old.distributions
new.distributions
comp = DataFrame(; old=length.(old.distributions), new=length.(new.distributions))
insertcols!(comp, :diff => comp.new .- comp.old)
insertcols!(comp, :reldiff => comp.diff ./ comp.old)
show(comp; allrows=true)
# Verify order of species
@test old.spenames == new.spenames
insertcols!(comp, :oldnames => old.spenames, :newnames => new.spenames)
insertcols!(comp, :equal => isequal.(comp.oldnames, comp.newnames))
show(comp; allrows=true)

# Compare by species
olddf = DataFrame(; sp=old.spenames, n_old=length.(old.distributions))
newdf = DataFrame(; sp=new.spenames, n_new=length.(new.distributions))
comp2 = outerjoin(olddf, newdf; on=:sp)
insertcols!(comp2, :diff => comp2.n_new .- comp2.n_old)
insertcols!(comp2, :reldiff => comp2.diff ./ comp2.n_old)
show(comp2; allrows=true)
show(sort(comp2, :diff); allrows=true)
show(sort(comp2, :reldiff); allrows=true)

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
module Layer1Tests

using Test
using HDF5

using PicoQuant
include("test_utils.jl")

@testset "Layer 1 tests" begin
    @testset "Test executing dsl commands from file" begin
        qasm_str = """OPENQASM 2.0;
                    include "qelib1.inc";
                    qreg q[3];
                    h q[0];
                    cx q[0],q[1];
                    cx q[1],q[2];"""

        circ = load_qasm_as_circuit(qasm_str)
        tng = convert_qiskit_circ_to_network(circ, DSLBackend("contract_network.tl",
                                                              "tensor_data.h5",
                                                              "",
                                                              true))
        add_input!(tng, "000")
        add_output!(tng, "000")
        plan = random_contraction_plan(tng)

        try
            # Try executing the dsl commands created by contract_network!().
            contract_network!(tng, plan)
            execute_dsl_file(Array{ComplexF64}, "contract_network.tl", "tensor_data.h5")

            # Check if the correct amplitube was computed
            @test begin
                result = h5open("tensor_data.h5", "r") do file
                    read(file, "result")
                end
                result ≈ 1/sqrt(2)
            end

        finally
            # Clean up any files created.
            rm("contract_network.tl", force=true)
            rm("tensor_data.h5", force=true)
        end
    end

    @testset "Test tensor functions" begin

        # Create a tensor A to manipulate.
        A = [1+1im 1im; -1im 2.0]
        tensors = Dict{Symbol, Array{<:Number}}(:A => A)

        # Use the transpose and conjugate functions to get the adjoint of A.
        tensors[:A] = transpose_tensor(tensors[:A], [2,1])
        tensors[:A] = conjugate_tensor(tensors[:A])

        # Test if the correct tensor was created
        @test tensors[:A] == A'

        # Test reshaping the tensor.
        tensors[:A] = reshape_tensor(tensors[:A], 4)
        @test tensors[:A] == [1-1im, -1im, 1im, 2.0]
    end
end
end


text:
using Libdl

const depsfile = joinpath(@__DIR__, "deps.jl")

libpath = get(ENV, "JULIA_HDF5_PATH",
                get(ENV, "JULIA_HDF5_LIBRARY_PATH", nothing)) # legacy env variable for compatibility

new_contents = if libpath === nothing
    # By default, use HDF5_jll
    """
    # This file is automatically generated
    # Do not edit
    using HDF5_jll
    check_deps() = nothing
    """
else
    @info "using system HDF5"

    libhdf5 = find_library("libhdf5", [libpath, joinpath(libpath, "lib"), joinpath(libpath, "lib64")])
    libhdf5_hl = find_library("libhdf5_hl", [libpath, joinpath(libpath, "lib"), joinpath(libpath, "lib64")])

    isempty(libhdf5) && error("libhdf5 could not be found")
    isempty(libhdf5_hl) && error("libhdf5_hl could not be found")

    libhdf5_size = filesize(dlpath(libhdf5))

    """
    # This file is automatically generated
    # Do not edit

    function check_deps()
        if libhdf5_size != filesize(Libdl.dlpath(libhdf5))
            error("HDF5 library has changed, re-run Pkg.build(\\\"HDF5\\\")")
        end
        if h5_get_libversion() < v"1.10.4"
            error("HDF5.jl requires ≥ v1.10.4 of the HDF5 library.")
        end
    end
    $(:(const libhdf5 = $libhdf5))
    $(:(const libhdf5_hl = $libhdf5_hl))
    $(:(const libhdf5_size = $libhdf5_size))
    """
end

if !isfile(depsfile) || new_contents != read(depsfile, String)
    # only write file if contents have changed to avoid triggering re-precompilation each build
    open(depsfile, "w") do io
        print(io, new_contents)
    end
end


text:
# Sample postprocessing script for customized figure.
#
# Plotting of the same style from a series of input files. Combined with ffmpeg,
# we can easily make animations from data.
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, Glob, PyPlot, Printf
using Vlasiator: set_args, prep2d, set_colorbar, set_plot

files = glob("out/bulk*.vlsv")
nfile = length(files)

meta = load(files[1])
vardict = Dict("v"=>"proton/vg_v", "rho"=>"proton/vg_rho", "b"=>"vg_b_vol", "b2"=>"fg_b")
varname = "rho"

fig, ax = plt.subplots()

op = :z
axisunit = RE
colorscale = Log
addcolorbar = true
cmap = matplotlib.cm.turbo
vmin = 7.0e4
vmax = 2.5e6

pArgs = set_args(meta, vardict[varname], axisunit; normal=:none)

cnorm, cticks = set_colorbar(colorscale, vmin, vmax)

for (i, file) in enumerate(files)
   @info "$i out of $nfile"
   local meta = load(file)

   var = meta[vardict[varname]]
   t = readparameter(meta, "time")

   if i == 1
      x, y = Vlasiator.get_axis(pArgs)
      data = prep2d(meta, vardict[varname], op)'

      c = ax.pcolormesh(x, y, data, norm=cnorm, cmap=cmap, shading="nearest")

      set_plot(c, ax, pArgs, cticks, addcolorbar)
   else
      data = prep2d(meta, vardict[varname], op)'

      c = ax.pcolormesh(x, y, data, norm=cnorm, cmap=cmap, shading="nearest")
   end

   xlabel("X")
   ylabel("Y")
   str_title = @sprintf "t= %4.1fs" t
   title(str_title)

   savefig("out/"*lpad(i, 4, '0')*".png", bbox_inches="tight")

   ax.cla()
end

close()


text:
module NonconvexAugLagLab

export AugLag2, AugLag2Options

import Optim
using Reexport, Parameters, ChainRulesCore
using Setfield, NonconvexMMA
@reexport using NonconvexCore
using NonconvexCore: @params, AbstractOptimizer, AbstractModel
using NonconvexCore: ConvergenceState, ConvergenceCriteria
using NonconvexCore: VecModel, Tolerance, AbstractFunction
using NonconvexCore: Solution, GenericCriteria, Trace, NoCallback
using NonconvexCore: getnineqconstraints, debugging, value_gradient
using NonconvexCore: assess_convergence!, hasconverged, GenericResult

import NonconvexCore: optimize!, Workspace, getmin, getmax, getdim
import NonconvexCore: getobjective, getresiduals
import NonconvexCore: getineqconstraints, geteqconstraints
import NonconvexCore: getobjectiveconstraints

include("aggregations.jl")
include("model.jl")
include("algorithm.jl")

end


text:
# Rawdata Import
function LoadData(mpath, filename)
    # data import
    if contains(filename,".mat")
        ratdata = matread(joinpath(mpath,filename))
    else
        ratdata = matread(joinpath(mpath,*(filename,".mat")))
    end
    println("rawdata from ", filename, " imported" )

    # number of trials
    ntrials = size(ratdata["rawdata"]["leftbups"],2)

    return ratdata, ntrials
end

# Get Trial Data
function TrialData(rawdata, trial::Int)
    if rawdata["pokedR"][trial] > 0
        rat_choice = 1;  # "R"
    else
        rat_choice = -1; # "L"
    end;

    if typeof(rawdata["rightbups"][trial]) <: Array
        rvec = vec(rawdata["rightbups"][trial])::Array{Float64,1};
    else
        rvec = Float64[rawdata["rightbups"][trial]]
    end
    if typeof(rawdata["leftbups"][trial]) <: Array
        lvec = vec(rawdata["leftbups"][trial])::Array{Float64,1};
    else
        lvec = Float64[rawdata["leftbups"][trial]]
    end

    return rvec, lvec,
    rawdata["T"][trial]::Float64, rat_choice
end

# Write File
function WriteFile(mpath, filename, D)
   saveto_filename = joinpath(mpath,*("julia_out_",filename))
   matwrite(saveto_filename, D)
end


text:
# Pose3Pose3 evaluation functions

mutable struct PriorPose3 <: IncrementalInference.FunctorSingleton
    Zi::Distribution
    PriorPose3() = new()
    PriorPose3(st::FloatInt, sr::Float64) = new( MvNormal(zeros(6), [[st*eye(3);zeros(3,3)];[zeros(3);sr*eye(3)]] )  )
    PriorPose3(s::Distribution) = new(s)
end
function getSample(p3::PriorPose3, N::Int=1)
  # mv = Distributions.MvNormal(veeEuler(p3.Zi), p3.Cov)
  return (rand(p3.Zi, N),)
end
mutable struct PackedPriorPose3  <: IncrementalInference.PackedInferenceType
    vecZi::Array{Float64,1} # 0rotations, 1translation in each column
    vecCov::Array{Float64,1}
    dimc::Int
    PackedPriorPose3() = new()
    PackedPriorPose3(x...) = new(x[1], x[2], x[3])
end
function convert(::Type{PriorPose3}, d::PackedPriorPose3)
  Zi = SE3(d.vecZi[1:3], Quaternion(d.vecZi[4],d.vecZi[5:7]))
  Cov = reshapeVec2Mat(d.vecCov, d.dimc)
  return PriorPose3( MvNormal( veeEuler(Zi), Cov) )
end
function convert(::Type{PackedPriorPose3}, d::PriorPose3)
  tf = SE3(d.Zi.μ[1:3], Euler(d.Zi.μ[4:6]...) )
  v1 = veeQuaternion(tf)
  v2 = d.Zi.Σ.mat[:];
  return PackedPriorPose3(v1, v2, size(d.Zi.Σ.mat,1))
end

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

mutable struct PP3REUSE
  wTi::SE3
  wTj::SE3
  iTi::SE3
  PP3REUSE() = new(SE3(0),SE3(0),SE3(0))
end

function fastpose3pose3residual!(reusethrid::PP3REUSE,
      res::Array{Float64},
      idx::Int,
      meas::Tuple,
      wXi::Array{Float64,2},
      wXj::Array{Float64,2}  )
  #
  reusethrid.wTi.t[1:3] = wXi[1:3,idx]
  TransformUtils.convert!(reusethrid.wTi.R, Euler(wXi[4,idx],wXi[5,idx],wXi[6,idx]))
  reusethrid.wTj.t[1:3] = wXj[1:3,idx]
  TransformUtils.convert!(reusethrid.wTj.R, Euler(wXj[4,idx],wXj[5,idx],wXj[6,idx]))

  # TODO -- convert to in place convert! functions, many speed-ups possible here
  jTi = SE3( matrix(reusethrid.wTj)\matrix(reusethrid.wTi) )
  # also wasted memory here, should operate directly on iTi and not be assigning new memory
  reusethrid.iTi = (SE3(meas[1][1:3,idx],Euler(meas[1][4:6,idx]...)) * jTi)
  res[:] = veeEuler(reusethrid.iTi)
  nothing
end

mutable struct Pose3Pose3 <: RoME.BetweenPoses
    Zij::Distribution
    reuse::Vector{PP3REUSE}
    Pose3Pose3() = new()
    Pose3Pose3(s) = new(s, fill(PP3REUSE(), Threads.nthreads() )  )
end
function getSample(pp3::Pose3Pose3, N::Int=1)
  return (rand(pp3.Zij, N), )
end
function (pp3::Pose3Pose3)(res::Array{Float64},
      idx::Int,
      meas::Tuple,
      wXi::Array{Float64,2},
      wXj::Array{Float64,2}  )
  #
  reusethrid = pp3.reuse[Threads.threadid()]
  fastpose3pose3residual!(reusethrid, res, idx, meas, wXi, wXj)
  nothing
end

mutable struct PackedPose3Pose3 <: IncrementalInference.PackedInferenceType
  vecZij::Array{Float64,1} # 3translations, 3rotation
  vecCov::Array{Float64,1}
  dimc::Int
  PackedPose3Pose3() = new()
  PackedPose3Pose3(x...) = new(x[1], x[2], x[3])
end
function convert(::Type{Pose3Pose3}, d::PackedPose3Pose3)
  qu = Quaternion(d.vecZij[4], d.vecZij[5:7])
  se3val = SE3(d.vecZij[1:3], qu)
  cov = reshapeVec2Mat(d.vecCov, d.dimc)
  return Pose3Pose3( MvNormal(veeEuler(se3val), cov) )
end
function convert(::Type{PackedPose3Pose3}, d::Pose3Pose3)
  val = d.Zij.μ
  se3val = SE3(val[1:3], Euler(val[4:6]...))
  v1 = veeQuaternion(se3val)
  v2 = d.Zij.Σ.mat
  return PackedPose3Pose3(v1[:], v2[:], size(v2,1) )
end

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

mutable struct Pose3Pose3NH <: IncrementalInference.FunctorPairwiseNH
    Zij::Distribution
    nullhypothesis::Distributions.Categorical
    reuse::Vector{PP3REUSE}
    Pose3Pose3NH() = new()
    Pose3Pose3NH(s::Distribution, vh::Vector{Float64}) = new(s, Distributions.Categorical(vh), fill(PP3REUSE(), Threads.nthreads() )  )
    # Pose3Pose3NH(s::SE3, c::Array{Float64,2}, vh::Float64) = new(s,c, Distributions.Categorical([(1.0-vh);vh]),SE3(0),SE3(0),SE3(0))
    # Pose3Pose3NH(st::FloatInt, sr::Float64;vh::Float64=1.0) = new(SE3(0), [[st*eye(3);zeros(3,3)];[zeros(3);sr*eye(3)]], Distributions.Categorical([(1.0-vh);vh]),SE3(0),SE3(0),SE3(0))
end
function getSample(pp3::Pose3Pose3NH, N::Int=1)
  return (rand(pp3.Zij, N), )
end
function (pp3::Pose3Pose3NH)(res::Array{Float64},
      idx::Int,
      meas::Tuple,
      wXi::Array{Float64,2},
      wXj::Array{Float64,2}  )
  #
  reusethrid = pp3.reuse[Threads.threadid()]
  fastpose3pose3residual!(reusethrid, res, idx, meas, wXi, wXj)
  nothing
end

mutable struct PackedPose3Pose3NH <: IncrementalInference.PackedInferenceType
  vecZij::Vector{Float64} # 3translations, 3rotation
  vecCov::Vector{Float64}
  dimc::Int
  nullhypothesis::Vector{Float64}
  PackedPose3Pose3NH() = new()
  PackedPose3Pose3NH(x1::Vector{Float64},x2::Vector{Float64},x3::Int,x4::Vector{Float64}) = new(x1, x2, x3, x4)
end

function convert(::Type{Pose3Pose3NH}, d::PackedPose3Pose3NH)
  qu = Quaternion(d.vecZij[4], d.vecZij[5:7])
  se3val = SE3(d.vecZij[1:3], qu)
  cov = reshapeVec2Mat(d.vecCov, d.dimc)
  return Pose3Pose3NH( MvNormal(veeEuler(se3val), cov), d.nullhypothesis )
end
function convert(::Type{PackedPose3Pose3NH}, d::Pose3Pose3NH)
  val = d.Zij.μ
  se3val = SE3(val[1:3], Euler(val[4:6]...))
  v1 = veeQuaternion(se3val)
  v2 = d.Zij.Σ.mat
  return PackedPose3Pose3NH(v1[:], v2[:], size(v2,1), d.nullhypothesis.p )
end

# TODO -- stronger type safety required here
# Project all particles (columns) Xval with Z, that is for all  SE3(Xval[:,i])*Z
function projectParticles(Xval::Array{Float64,2}, Z::Distribution)
  # TODO optimize for more speed with threads and better memory management
  r,c = size(Xval)
  RES = zeros(r,c) #*cz

  ent, x = SE3(0), SE3(0)
  ENT = rand( Z, c )
  # ENT = rand( MvNormal(zeros(6), Cov), c )
  j=1
  # for j in 1:cz
  for i in 1:c
    x.R = TransformUtils.convert(SO3,Euler(Xval[4,i],Xval[5,i],Xval[6,i]))
    x.t = Xval[1:3,i]
    ent.R =  TransformUtils.convert(SO3, Euler(ENT[4:6,i]...)) # so3
    ent.t = ENT[1:3,i]
    # newval = Z*ent
    # res = x*newval
    res = x*ent
    RES[1:r,i*j] = veeEuler(res)
  end
  # end
  #
  return RES
end

⊕(Xpts::Array{Float64,2}, z::Pose3Pose3) = projectParticles(Xpts, z.Zij)
⊕(Xvert::Graphs.ExVertex, z::Pose3Pose3) = ⊕(getVal(Xvert), z)

# 3 dimensionsal evaluation functions.
# Transforms code is scattered, and will be moved to better location with stable 3D Examples
# also worth it to consolidate if other affine transforms package is available

# using Euler angles for linear sampling in product of potentials
# Gaussian model for prior
# function evalPotential(obs::PriorPose3, Xi::Array{Graphs.ExVertex,1}; N::Int=200)
#   mu = veeEuler(obs.Zi)
#   return rand( MvNormal(mu, obs.Cov), N)
# end
#
#
# # Still limited to linear sampler, then reprojected onto ball -- TODO upgrade manifold sampler
# function evalPotential(odom::Pose3Pose3, Xi::Array{Graphs.ExVertex,1}, Xid::Int; N::Int=200)
#     # rz,cz = size(odom.Zij)
#     Xval = Array{Float64,2}()
#     # implicit equation portion -- bi-directional pairwise function made explicit here
#     if Xid == Xi[1].index #odom.
#         # reverse direction
#         Z = inverse(odom.Zij)
#         Xval = getVal(Xi[2])
#     elseif Xid == Xi[2].index
#         # forward direction
#         Z = odom.Zij
#         Xval = getVal(Xi[1])
#     else
#         error("Bad evalPairwise Pose3Pose3")
#     end
#
#     return projectParticles(Xval, Z, odom.Cov)
# end
#
# # Still limited to linear sampler, then reprojected onto ball -- TODO upgrade manifold sampler
# function evalPotential(odom::Pose3Pose3NH,
#       Xi::Array{Graphs.ExVertex,1},
#       Xid::Int;
#       N::Int=200)
#    #
#   # rz,cz = size(odom.Zij)
#   Xval = Array{Float64,2}()
#   # implicit equation portion -- bi-directional pairwise function made explicit here
#   if Xid == Xi[1].index #odom.
#       # reverse direction
#       Z = inverse(odom.Zij)
#       Xval = getVal(Xi[2])
#       oldval = getVal(Xi[1])
#   elseif Xid == Xi[2].index
#       # forward direction
#       Z = odom.Zij
#       Xval = getVal(Xi[1])
#       oldval = getVal(Xi[2])
#   else
#       error("Bad evalPairwise Pose3Pose3")
#   end
#
#   skipdos = rand(odom.ValidHypot, N)
#   dos, donts = skipdos .== 2, skipdos .== 1
#
#   projted = Array{Float64}(6, N)
#   projted[:,dos] = projectParticles(Xval[:,dos], Z, odom.Cov)
#
#   projted[:,donts] = oldval[:,donts]
#
#   return projted
# end

#


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

mutable struct IPConfigurationBgpPeeringAddress <: SwaggerModel
    ipconfigurationId::Any # spec type: Union{ Nothing, String } # spec name: ipconfigurationId
    defaultBgpIpAddresses::Any # spec type: Union{ Nothing, Vector{String} } # spec name: defaultBgpIpAddresses
    customBgpIpAddresses::Any # spec type: Union{ Nothing, Vector{String} } # spec name: customBgpIpAddresses
    tunnelIpAddresses::Any # spec type: Union{ Nothing, Vector{String} } # spec name: tunnelIpAddresses

    function IPConfigurationBgpPeeringAddress(;ipconfigurationId=nothing, defaultBgpIpAddresses=nothing, customBgpIpAddresses=nothing, tunnelIpAddresses=nothing)
        o = new()
        validate_property(IPConfigurationBgpPeeringAddress, Symbol("ipconfigurationId"), ipconfigurationId)
        setfield!(o, Symbol("ipconfigurationId"), ipconfigurationId)
        validate_property(IPConfigurationBgpPeeringAddress, Symbol("defaultBgpIpAddresses"), defaultBgpIpAddresses)
        setfield!(o, Symbol("defaultBgpIpAddresses"), defaultBgpIpAddresses)
        validate_property(IPConfigurationBgpPeeringAddress, Symbol("customBgpIpAddresses"), customBgpIpAddresses)
        setfield!(o, Symbol("customBgpIpAddresses"), customBgpIpAddresses)
        validate_property(IPConfigurationBgpPeeringAddress, Symbol("tunnelIpAddresses"), tunnelIpAddresses)
        setfield!(o, Symbol("tunnelIpAddresses"), tunnelIpAddresses)
        o
    end
end # type IPConfigurationBgpPeeringAddress

const _property_map_IPConfigurationBgpPeeringAddress = Dict{Symbol,Symbol}(Symbol("ipconfigurationId")=>Symbol("ipconfigurationId"), Symbol("defaultBgpIpAddresses")=>Symbol("defaultBgpIpAddresses"), Symbol("customBgpIpAddresses")=>Symbol("customBgpIpAddresses"), Symbol("tunnelIpAddresses")=>Symbol("tunnelIpAddresses"))
const _property_types_IPConfigurationBgpPeeringAddress = Dict{Symbol,String}(Symbol("ipconfigurationId")=>"String", Symbol("defaultBgpIpAddresses")=>"Vector{String}", Symbol("customBgpIpAddresses")=>"Vector{String}", Symbol("tunnelIpAddresses")=>"Vector{String}")
Base.propertynames(::Type{ IPConfigurationBgpPeeringAddress }) = collect(keys(_property_map_IPConfigurationBgpPeeringAddress))
Swagger.property_type(::Type{ IPConfigurationBgpPeeringAddress }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IPConfigurationBgpPeeringAddress[name]))}
Swagger.field_name(::Type{ IPConfigurationBgpPeeringAddress }, property_name::Symbol) =  _property_map_IPConfigurationBgpPeeringAddress[property_name]

function check_required(o::IPConfigurationBgpPeeringAddress)
    true
end

function validate_property(::Type{ IPConfigurationBgpPeeringAddress }, name::Symbol, val)
end


text:
using DomainSets, GridArrays, Base.Broadcast, Test

@testset "broadcast" begin

    broadcast(in, PeriodicEquispacedGrid(10,-1,1) , UnitInterval()) == [i> 5 for i in 1:10]

    uniondomain = UnionDomain(0.0..0.1,0.5..1.0)
    @test broadcast(in, PeriodicEquispacedGrid(10,-1,1) , uniondomain) == [6<=i<=6 || 9<=i<=10 for i in 1:10]

    intersectdomain = IntersectionDomain(-1.0..0.1,0.0..1.0)
    @test broadcast(in, PeriodicEquispacedGrid(10,-1,1) , intersectdomain) == [i==6 for i in 1:10]

    diffdomain = DifferenceDomain(-1.0..0.1,0.0..1.0)
    @test broadcast(in, PeriodicEquispacedGrid(10,-1,1) , diffdomain) == [i<6 for i in 1:10]
end


text:
# # Ice Flow over a bed with variable slipperiness and viscosity.

# In this example, we consider flow in two horizontal dimensions, with variable ice viscosity and basal slipperiness.
# This example is similar in spirit to the MISMIP 3D Intercomparison exercise (http://homepages.ulb.ac.be/~fpattyn/mismip3d/)
# This example demonstrates how to apply spatially heterogenous basal slipperiness and ice viscosity fields.

# ## Install dependencies
#
# First let's make sure we have all required packages installed.

#using Pkg
#Pkg.add("https://github.com/RJArthern/WAVI.jl"), Pkg.add(Plots)
using WAVI, Plots

# ## Model Setup
# The grid has length 800 km in the x-direction and 50km in the y-direction. The grid spacing is 2.5km
dx = 2500., dy = 2500.,
nx = round(Int, 800*1e3 / dx)
ny = round(Int, 50*1e3 / dy)
grid = Grid(nx=nx,ny=ny,dx=2000.,dy=2000.);

# Since we're running to steady state, we'll only do one velocity solve iteration per timestep:
solver_params = SolverParams(maxiter_picard = 1);

# and we'll take a uniform accumulation rate
params = Params( accumulation_rate = 0.3);

# ## Basal slipperiness
# We set the basal slipperiness via the model's initial conditions. In the MISMIP 3D experiment, the sliding coefficient is a Gaussian bump:
# $C = C_0 \left[1 - a  \exp (-(x- x_b)^2 / 2x_c^2  - (y - y_b)^2 / 2y_c^2) \right] where a = 0.75, x_b = 400km, ...


text:
struct EarthMoonBarycenter <: abstractBarycenter end

naif_id(::Type{EarthMoonBarycenter}) = 3
parent(::Type{EarthMoonBarycenter}) = SSB


text:
using Distributions, SpecialFunctions, Roots, LinearAlgebra

"""
GVarSmile: Variational SMiLe for the Gaussian task (section 2.3.1)
"""
struct GVarSmile <: GLearner
    m::Float64
    χ_0::Float64
    ν_0::Float64
    sigma::Float64 # std of observations
    χ_n::Array{Float64, 1}
    ν_n::Array{Float64, 1}
    mu_n::Array{Float64, 1} #parameter we are interested in
    Sgm_n::Array{Float64, 1} # Bayes Factor Surprise
    var_n::Array{Float64, 1}
    γ_n::Array{Float64, 1}
    Ssh_n::Array{Float64, 1} # Shannon Surprise
    changeprobability::Float64 # Used only for Shannon calculation
end
function GVarSmile(;m = .1, mu_0 = 0., sigma_0 = 1., sigma = 0.1,
                        changeprobability = 0.05)
    ν_0 = calcν(sigma, sigma_0)
    χ_0 = calcχ(mu_0, ν_0, sigma)
    var_n_0 = sigma^2 / ν_0
    mu_n_0 = sigma * χ_0 / ν_0
    GVarSmile(m, χ_0, ν_0, sigma, [χ_0], [ν_0], [mu_n_0], [1.], [var_n_0],
                    [m/(1. + m)], [1.], changeprobability)
end
export GVarSmile
function update!(learnerG::GVarSmile, y)
    calcSgm!(learnerG, y)
    calcγ!(learnerG)
    calcShannon!(learnerG, y)
    update_χ_ν!(learnerG, y)
    computemu_n!(learnerG)
    computevar_n!(learnerG)

end
export update!
function update_χ_ν!(learnerG::GVarSmile, y)
    learnerG.χ_n[1] = (1. - learnerG.γ_n[1]) * learnerG.χ_n[1] + learnerG.γ_n[1] * learnerG.χ_0 + calcφ(learnerG, y)
    learnerG.ν_n[1] = (1. - learnerG.γ_n[1]) * learnerG.ν_n[1] + learnerG.γ_n[1] * learnerG.ν_0 + 1.
end


text:
export rref, rref!, solve, solve_augmented

"""
`swap_rows!(A,i,j)` swaps rows `i` and `j` in the matrix `A`.
"""
function swap_rows!(A::Array{T,2}, i::Int, j::Int) where T
  if i==j
    return nothing
  end
  A[ [i,j], :] = A[ [j,i], :]
  return nothing
end

"""
`add_row_to_row!(A,i,j)` adds row `i` to row `j` in the matrix `A`.
"""
function add_row_to_row!(A::Array{T,2},i::Int,j::Int) where T<:Number
  A[j,:] += A[i,:]
  return nothing
end

# The following functions developed by Tara Abrishami

"""
`rref!(A)` overwrites `A` with its row reduced echelon form.
"""
function rref!(A::Array{GF2,2})
  r, c = size(A)
  s = 0
  for x in 1:r
    b = false
    while !b && x + s <= c
      if A[x, x+s] == 1
        break
      elseif A[x, x + s] == 0
        for y in x:r
          if A[y, x + s] == 1
            swap_rows!(A, y, x)
            b = true
            break
          end
        end
      end
      if !b
        s = s + 1
      end
    end
    for m in 1:r
      if x + s <= c && m != x && A[m, x+s] == 1
        add_row_to_row!(A, x, m)
      end
    end
  end
end

"""
`rref(A)` returns the row reduced echelon form of `A`.
"""
function rref(A::Array{GF2,2})
  AA = copy(A)
  rref!(AA)
  return AA
end

"""
`solve(A,b)` returns a solution `x` to the linear system
`A*x == b` or throws an error if no solution can be found.
"""
function solve(A::Array{GF2, 2}, b::Array{GF2, 1})
  r, c = size(A)
  if r != size(b)[1]
    error("Dimensionally incorrect input")
  end
  C = [A b]
  return solve_augmented(C)
end

# returns a single solution to the system with matrix C1
function solve_augmented(C1::Array{GF2, 2})
  r, c = size(C1)
  D = copy(C1)
  rref!(D)
  x = 0
  for a in 1:r
   in = true
   for b in 1:c-1
     if D[a, b] != 0
       in = false
     end
   end
   if in && D[a, c] != 0
     error("Inconsistent system")
   end
  end
  ret = zeros(GF2, c-1)
  for p in 1:r
    if D[p, c] == 1
      for n in 1:c-1
        if D[p, n] == 1
          ret[n] = 1
          break
        end
      end
    end
  end
  return ret
end

import Base.inv

function inv(A::Array{GF2,2})
  n,m = size(A)
  if n!= m
    error("Cannot invert a matrix that isn't square.")
  end
  if det(A)==0
    error("Cannot invert a singular matrix.")
  end
  In = Matrix{GF2}(I,n,n)

  AB = [A  In]
  rref!(AB)

  B = AB[:,n+1:end]
  return B
end

function LinearAlgebra.nullspace(A::Array{GF2, 2})
  r, c = size(A)
  M = rref(A)
  ret = zeros(GF2, c)
  s = 0
  x = 1
  left = false
  while x <= c
    if x > r
      left = true
      break
    end
    if x + s > c
      break
    end
    if M[x, x+s] == 1
      x = x + 1
      continue
    else
      p = zeros(GF2, c)
      p[x + s] = 1
      for t in 1:r
        if M[t, x + s] == 1
          for q in 1:c
            if M[t, q] == 1
              p[q] = 1
              break
            end
          end
        end
      end
      s = s + 1
      ret = hcat(ret, p)
    end
  end
  if left
    while x + s <= c
      p = zeros(GF2, c)
      p[x + s] = 1
      for t in 1:r
        if M[t, x + s] == 1
          for q in 1:c
            if M[t, q] == 1
              p[q] = 1
              break
            end
          end
        end
      end
      s = s + 1
      ret = hcat(ret, p)
    end
  end
  ret = ret[:, 1:size(ret,2) .!= 1]
  return ret
end

export solve_all
"""
`solve_all(A,b)` returns a solution to `A*x==b` together with
a basis for the nullspace of `A`.
"""
function solve_all(A::Array{GF2, 2}, b::Array{GF2, 1})
  return solve(A, b), nullspace(A)
end

export nullity
"""
`nullity(A)` returns the dimension of the nullspace of `A`.
"""
function nullity(A::Array{GF2,2})
  NS = nullspace(A)
  (x,n) = size(NS)
  return n
end

function LinearAlgebra.rank(A::Array{GF2,2})
  r,c = size(A)
  n = nullity(A)
  return c-n
end


text:
#test file

x = 3
y = 4

x + y

using DataFrames
df = DataFrame(a = 1:20,b = rand([2,5,10],20), x = rand(20))
using Plots
plot(rand(10))

max_stars_repo_name: snowdj/CoursePack

max_issues_repo_name: snowdj/CoursePack

max_forks_repo_name: snowdj/CoursePack


text:
@forwardfun autocov
@forwardfun autocorr
@forwardfun make_ufunc
@forwardfun wrap_xarray_ufunc


text:
export ishexcolor

function ishexcolor(str::AbstractString)::Bool
    hexcolorReg = r"^#?([0-9A-F]{3}|[0-9A-F]{4}|[0-9A-F]{6}|[0-9A-F]{8})$"i

    return contains(str, hexcolorReg)
end


text:
"""
FOLDER_PATH

Container to keep track of where Franklin is being run.
"""
const FOLDER_PATH = Ref{String}()

"""
IGNORE_FILES

Collection of file names that will be ignored at compile time.
"""
const IGNORE_FILES = [".DS_Store", ".gitignore", "LICENSE.md", "README.md",
                      "franklin", "franklin.pub", "node_modules/"]

"""
INFRA_EXT

Collection of file extensions considered as potential infrastructure files.
"""
const INFRA_FILES = [".html", ".css"]

"""
PATHS

Dictionary for the paths of the input folders and the output folders. The simpler case only
requires the main input folder to be defined i.e. `PATHS[:src]` and infers the others via the
`set_paths!()` function.
"""
const PATHS = LittleDict{Symbol,String}()

"""
$(SIGNATURES)

This assigns all the paths where files will be read and written with root the `FOLDER_PATH`
which is assigned at runtime.
"""
function set_paths!()::LittleDict{Symbol,String}
    @assert isassigned(FOLDER_PATH) "FOLDER_PATH undefined"
    @assert isdir(FOLDER_PATH[])    "FOLDER_PATH is not a valid path"

    # NOTE it is recommended not to change the names of those paths.
    # Particularly for the output dir. If you do, check for example that
    # functions such as Franklin.publish point to the right dirs/files.

    if FD_ENV[:STRUCTURE] < v"0.2"
        PATHS[:folder]    = normpath(FOLDER_PATH[])
        PATHS[:src]       = joinpath(PATHS[:folder], "src")
        PATHS[:src_pages] = joinpath(PATHS[:src],    "pages")
        PATHS[:src_css]   = joinpath(PATHS[:src],    "_css")
        PATHS[:src_html]  = joinpath(PATHS[:src],    "_html_parts")
        PATHS[:pub]       = joinpath(PATHS[:folder], "pub")
        PATHS[:css]       = joinpath(PATHS[:folder], "css")
        PATHS[:libs]      = joinpath(PATHS[:folder], "libs")
        PATHS[:assets]    = joinpath(PATHS[:folder], "assets")
        PATHS[:literate]  = joinpath(PATHS[:folder], "scripts")
        PATHS[:tag]       = joinpath(PATHS[:pub],    "tag")
    else
        PATHS[:folder]   = normpath(FOLDER_PATH[])
        PATHS[:site]     = joinpath(PATHS[:folder], "__site")    # mandatory
        PATHS[:assets]   = joinpath(PATHS[:folder], "_assets")   # mandatory
        PATHS[:css]      = joinpath(PATHS[:folder], "_css")      # mandatory
        PATHS[:layout]   = joinpath(PATHS[:folder], "_layout")   # mandatory
        PATHS[:libs]     = joinpath(PATHS[:folder], "_libs")     # mandatory
        PATHS[:literate] = joinpath(PATHS[:folder], "_literate") # optional
        PATHS[:tag]      = joinpath(PATHS[:site],   "tag")
    end

    return PATHS
end

"""
    path(s)

Return the paths corresponding to `s` e.g. `path(:folder)`.
"""
path(s) = PATHS[Symbol(s)]

"""
Pointer to the `/output/` folder associated with an eval block (see also
[`@OUTPUT`](@ref)).
"""
const OUT_PATH = Ref("")

"""
This macro points to the `/output/` folder associated with an eval block.
So for instance, if an eval block generates a plot, you could save the plot
with something like `savefig(joinpath(@OUTPUT, "ex1.png"))`.
"""
macro OUTPUT()
    return OUT_PATH[]
end


text:
export TransposeOperator

mutable struct TransposeOperator{T<:Number,B<:Operator} <: Operator{T}
    op::B
end

TransposeOperator(B::Operator{T}) where {T<:Number}=TransposeOperator{T,typeof(B)}(B)

convert(::Type{Operator{T}},A::TransposeOperator) where {T}=TransposeOperator(convert(Operator{T},A.op))

domainspace(P::TransposeOperator)=rangespace(P.op)
rangespace(P::TransposeOperator)=domainspace(P.op)

domain(P::TransposeOperator)=domain(P.op)

bandwidths(P::TransposeOperator) = reverse(bandwidths(P.op))

getindex(P::TransposeOperator,k::Integer,j::Integer) =
    P.op[j,k]

function BandedMatrix(S::SubOperator{T,TO}) where {T,TO<:TransposeOperator}
    kr,jr=parentindices(S)
    transpose(BandedMatrix(view(parent(S).op,jr,kr)))
end

Base.transpose(A::Operator)=TransposeOperator(A)
adjoint(A::Operator{T}) where {T<:Real}=TransposeOperator(A)


text:
module FluxTikz

using Flux: Chain
import TikzPictures: TikzPicture
export TikzPicture

include("layer.jl")
include("network.jl")
include("tikzpicture.jl")

end