text:
# ANNEX 2C – THICKNESS, MAWP AND STRESS EQUATIONS FOR A FFS ASSESSMENT

# 2C.2.5 Treatment of Weld and Riveted Joint Efficiency, and Ligament Efficiency
#+++ condiser weld joint efficient for damage over welds +++ to add

# 2C.2.6 Treatment of Damage in Formed Heads
#+++ If damage (e.g. corrosion/erosion, pitting, etc.) occurs in the center section of an elliptical or torispherical head +++ to add

#2C.3 Pressure Vessels and Boiler Components – Internal Pressure
# 2C.3.3 Cylindrical Shells
# 2C.3.4 Spherical Shell or Hemispherical Head
# 2C.3.5 Elliptical Head
# 2C.3.6 Torispherical Head
# 2C.3.7 Conical Shell
# 2C.3.8 Toriconical Head
# 2C.3.9 Conical Transition
# 2C.3.10 Nozzles Connections in Shells
#2C.3.12 Other Components
# 2C.3.11 Junction Reinforcement Requirements at Conical Transitions
# Tubesheets: VIII-1 Part UHX, VIII-2, Part 4, paragraph 4.18 or TEMA.
# Flat head to cylinder connections: VIII-1 paragraph UG-34 or VIII-2, Part 4, paragraph 4.3.
# Bolted Flanges: VIII-1, Appendix 2 or VIII-2, Part 4, paragraph 4.16.
#2C.4 Pressure Vessels and Boiler Components – External Pressure
# 2C.5.3 Required Thickness and MAWP – Straight Pipes Subject To Internal Pressure
@doc """
    Pipingtcmin(P::Float64; Do::Float64=0.0, S::Float64=0.0, E::Float64=0.0, Yb31::Float64=0.0, MA::Float64=0.0)

2C.5 Piping Components and Boiler Tubes - equation (2C.146)\n
Circumferential stress (Longitudinal Joints)\n
Minimum required thickness based on the circumferential membrane stress for a cylinder or cone, as applicable\n
Do # Outside Diameter\n
P # internal design pressure\n
S # allowable stress\n
E # weld joint efficiency or quality factor from the original construction code, if unknown use 0.7\n
MA # mechanical allowances (thread or groove depth); for threaded components, the nominal thread depth (dimension h of ASME B.1.20.1) shall apply\n
Yb31 # coefficient from ASME B31 Piping codes used for determining the pipe wall thickness, the coefficient can be determined from the following table that is valid for tmin < Do / 6 Annex 2C\n
t = trd # thickness of the shell or pipe adjusted for mill tolerance, LOSS and FCA , or cylinder thickness at a conical transition for a junction reinforcement calculation adjusted for mill tolerance, LOSS and FCA , as applicable\n
""" ->
function Pipingtcmin(P::Float64; Do::Float64=0.0, S::Float64=0.0, E::Float64=0.0, Yb31::Float64=0.0, MA::Float64=0.0)
    tcmin = ((P*Do) / (2*((S*E) + (P*Yb31)))) + MA # (2C.146)
    return tcmin
end

@doc """
    PipingMAWPc(S::Float64; E::Float64=0.0, t::Float64=0.0, MA::Float64=0.0, Do::Float64=0.0, Yb31::Float64=0.0)

2C.5 Piping Components and Boiler Tubes - equation (2C.147)\n
Circumferential stress (Longitudinal Joints)\n
Maximum allowable working pressure based on circumferential stress\n
Do # Outside Diameter\n
S # allowable stress\n
E # weld joint efficiency or quality factor from the original construction code, if unknown use 0.7\n
MA # mechanical allowances (thread or groove depth); for threaded components, the nominal thread depth (dimension h of ASME B.1.20.1) shall apply\n
Yb31 # coefficient from ASME B31 Piping codes used for determining the pipe wall thickness, the coefficient can be determined from the following table that is valid for tmin < Do / 6 Annex 2C\n
t = trd # thickness of the shell or pipe adjusted for mill tolerance, LOSS and FCA , or cylinder thickness at a conical transition for a junction reinforcement calculation adjusted for mill tolerance, LOSS and FCA , as applicable\n
""" ->
function PipingMAWPc(S::Float64; E::Float64=0.0, t::Float64=0.0, MA::Float64=0.0, Do::Float64=0.0, Yb31::Float64=0.0)
    MAWPc = (2*S*E*(t-MA))/(Do-(2*Yb31*(t-MA))) # eq (2C.147)
    return MAWPc
end

@doc """
    PipingOcm(P::Float64; E::Float64=0.0, Do::Float64=0.0, t::Float64=0.0, MA::Float64=0.0, Yb31::Float64=0.0)

2C.5 Piping Components and Boiler Tubes - equation (2C.148)\n
Circumferential stress (Longitudinal Joints)\n
Nominal circumferential membrane stress for a cylinder or cone, as applicable\n
Do # Outside Diameter\n
P # internal design pressure\n
S # allowable stress\n
E # weld joint efficiency or quality factor from the original construction code, if unknown use 0.7\n
MA # mechanical allowances (thread or groove depth); for threaded components, the nominal thread depth (dimension h of ASME B.1.20.1) shall apply\n
Yb31 # coefficient from ASME B31 Piping codes used for determining the pipe wall thickness, the coefficient can be determined from the following table that is valid for tmin < Do / 6 Annex 2C\n
t = trd # thickness of the shell or pipe adjusted for mill tolerance, LOSS and FCA , or cylinder thickness at a conical transition for a junction reinforcement calculation adjusted for mill tolerance, LOSS and FCA , as applicable\n

""" ->
function PipingOcm(P::Float64; E::Float64=0.0, Do::Float64=0.0, t::Float64=0.0, MA::Float64=0.0, Yb31::Float64=0.0)
    Ocm = (P/E)*((Do/(2*(t-MA)))-Yb31) # eq (2C.148)
    return Ocm
end

@doc """
    Pipingtlmin(P::Float64; Do::Float64=0.0, S::Float64=0.0, E::Float64=0.0, Yb31::Float64=0.0, tsl::Float64=0.0, MA::Float64=0.0)

2C.5 Piping Components and Boiler Tubes - equation (2C.149)\n
Longitudinal stress (Circumferential Joints)\n
Minimum required thickness based on the longitudinal membrane stress for a cylinder or cone, as applicable.e\n
Do # Outside Diameter\n
P # internal design pressure\n
S # allowable stress\n
E # weld joint efficiency or quality factor from the original construction code, if unknown use 0.7\n
MA # mechanical allowances (thread or groove depth); for threaded components, the nominal thread depth (dimension h of ASME B.1.20.1) shall apply\n
tsl # supplemental thickness for mechanical loads other than pressure that result in longitudinal stress; this thickness is usually obtained from the results of a weight case in a stress analysis of the piping system (see paragraph 2C.2.7)\n
Yb31 # coefficient from ASME B31 Piping codes used for determining the pipe wall thickness, the coefficient can be determined from the following table that is valid for tmin < Do / 6 Annex 2C\n
t = trd # thickness of the shell or pipe adjusted for mill tolerance, LOSS and FCA , or cylinder thickness at a conical transition for a junction reinforcement calculation adjusted for mill tolerance, LOSS and FCA , as applicable\n
""" ->
function Pipingtlmin(P::Float64; Do::Float64=0.0, S::Float64=0.0, E::Float64=0.0, Yb31::Float64=0.0, tsl::Float64=0.0, MA::Float64=0.0)
    tlmin =(P*Do)/(4*((S*E)+(P*Yb31)))+tsl+MA # (2C.149)
    return tlmin
end

@doc """
    PipingMAWPl(S::Float64; E::Float64=0.0, t::Float64=0.0, tsl::Float64=0.0, MA::Float64=0.0, Do::Float64=0.0, Yb31::Float64=0.0)

2C.5 Piping Components and Boiler Tubes - equation (2C.150)\n
Longitudinal stress (Circumferential Joints)\n
Maximum allowable working pressure based on circumferential stress\n
Do # Outside Diameter\n
S # allowable stress\n
E # weld joint efficiency or quality factor from the original construction code, if unknown use 0.7\n
MA # mechanical allowances (thread or groove depth); for threaded components, the nominal thread depth (dimension h of ASME B.1.20.1) shall apply\n
tsl # supplemental thickness for mechanical loads other than pressure that result in longitudinal stress; this thickness is usually obtained from the results of a weight case in a stress analysis of the piping system (see paragraph 2C.2.7)\n
Yb31 # coefficient from ASME B31 Piping codes used for determining the pipe wall thickness, the coefficient can be determined from the following table that is valid for tmin < Do / 6 Annex 2C\n
t = trd # thickness of the shell or pipe adjusted for mill tolerance, LOSS and FCA , or cylinder thickness at a conical transition for a junction reinforcement calculation adjusted for mill tolerance, LOSS and FCA , as applicable\n
""" ->
function PipingMAWPl(S::Float64; E::Float64=0.0, t::Float64=0.0, tsl::Float64=0.0, MA::Float64=0.0, Do::Float64=0.0, Yb31::Float64=0.0)
    MAWPl = (4*S*E*(t-tsl-MA))/(Do-(4*Yb31*(t-tsl-MA))) # eq (2C.150)
    return MAWPl
end

@doc """
    PipingOlm(P::Float64; E::Float64=0.0, Do::Float64=0.0, t::Float64=0.0, tsl::Float64=0.0, MA::Float64=0.0, Yb31::Float64=0.0)

2C.5 Piping Components and Boiler Tubes - equation (2C.151)\n
Longitudinal stress (Circumferential Joints)\n
Nominal circumferential membrane stress for a cylinder or cone, as applicable\n
Do # Outside Diameter\n
P # internal design pressure\n
S # allowable stress\n
E # weld joint efficiency or quality factor from the original construction code, if unknown use 0.7\n
MA # mechanical allowances (thread or groove depth); for threaded components, the nominal thread depth (dimension h of ASME B.1.20.1) shall apply\n
Yb31 # coefficient from ASME B31 Piping codes used for determining the pipe wall thickness, the coefficient can be determined from the following table that is valid for tmin < Do / 6 Annex 2C\n
t = trd # thickness of the shell or pipe adjusted for mill tolerance, LOSS and FCA , or cylinder thickness at a conical transition for a junction reinforcement calculation adjusted for mill tolerance, LOSS and FCA , as applicable\n

""" ->
function PipingOlm(P::Float64; E::Float64=0.0, Do::Float64=0.0, t::Float64=0.0, tsl::Float64=0.0, MA::Float64=0.0, Yb31::Float64=0.0)
    Olm = (P/E)*((Do/(4*(t-tsl-MA)))-Yb31) # eq (2C.151)
    return Olm
end

# ASME B31.8 P
@doc """
    B31_8Piping_P(S::Float64; F::Float64=0.0, E::Float64=0.0, T::Float64=0.0, Do::Float64=0.0, t::Float64=0.0)

    ASME B31.8\n
    841.1 Steel Piping Systems Design Requirements\n
    ð18Þ 841.1.1 Steel Pipe Design Formula\n
Do # nominal outside diameter of pipe\n
S # specified minimum yield strength\n
F # design factor obtained from Table 841.1.6-1. In setting the values of the design factor, F, due consideration has been given and allowance has been made for the various underthickness \ntolerances provided for in the pipe specifications listed and approved for usage in this Code (B31.8)\n
E # longitudinal joint factor obtained from Table 841.1.7-1 [see also para. 817.1.3(d)]\n
T # temperature derating factor obtained from Table 841.1.8-1 (B31.8)\n
t = nominal wall thickness\n
P = design pressure\n
""" ->
function B31_8Piping_P(S::Float64; F::Float64=0.0, E::Float64=0.0, T::Float64=0.0, Do::Float64=0.0, t::Float64=0.0)
    P = (2*S*t / Do) * F * E * T # eq (2C.147)
    return P
end

# 2C.5.4 Required Thickness and MAWP – Boiler Tubes
# 2C.5.5 Required Thickness and MAWP – Pipe Bends Subject To Internal Pressure
#=
# API 579 - 2C.5.5 Required Thickness and MAWP – Pipe Bends Subject To Internal Pressure
P = 1900 # design pressure
tnom = .719 # nominal WT
LOSS = 0.0 # uniform metal loss away from the damage area at the time of the inspection.
FCA = 0.05 # future corrosion allowance
MA = 0 # mechanical allowances (thread or groove depth); for threaded components, the nominal thread depth (dimension h of ASME B.1.20.1) shall apply.
t = tnom - LOSS - FCA # adjuted for LOSS and FCA
NPS = 10 # Nominal pipe size
Do = 10.75 # Outside Diameter
Di = Do-(2*tnom)+2*(LOSS+FCA) # inside diameter
S = 20000 # code allowable material stress
E = 1.0 # weld joint efficiency
W = 1.0 # wed joint reduction factor
Yb31 = 0.4 # B31.3 temperature coefficient
Rb = 1.5 * NPS   # Centerline bend radius (see Figure 2C.1).
Rm = (Do+Di)/4 # Mean radius of the component; use the large end radius for a conical shell .

## Check applicability
if (Rm/t <= 10)
    print("Considered thin wall piping and 2C.5.5 is applicable")
else
    print("Considered thick wall piping - consult DIN 2413, Parts 1 and 2.")
end

##

θo = 20 # The circumferential angle measured from the crown of the elbow (+90° at the extrados, -90° at the intrados)
Lf = ((Rb/Rm)+(sind(θo)/2))/((Rb/Rm)+sind(θo)) # The Lorenz factor is computed using Equation (2C.161).
tCmin=(P*Do)/(2*(S*E/Lf)+(P*Yb31))+MA # Minimum required thickness based on the circumferential membrane stress for a cylinder or cone, as applicable
MAWPc =(2*((S*E)/Lf)*(t-MA))/(Do-(2*Yb31*(t-MA))) # Maximum allowable working pressure based on circumferential stress

print("\nMinimum required thickness based on the circumferential membrane stress: ", round(tCmin,digits=3))
print("\nMaximum allowable working pressure based on circumferential stress: ", round(MAWPc,digits=3),"psig")
print("\nMinimum required thickness based on the circumferential membrane stress @ PSV Set pressure 1900psig: ", round(tCmin,digits=3))
print("\nMaximum allowable working pressure based on circumferential stress @ PSV Set pressure 1900psig: ", round(MAWPc,digits=3),"psig")
=#
# 2C.5.6 Required Thickness and MAWP for External Pressure
# 2C.5.7 Branch Connections

#2C.6 API 650 Storage Tanks
# 2C.6.3 Required Thickness and MFH for Liquid Hydrostatic Loading
#


text:
# plane rotator
import SDELab2
include("common.jl")
########################## problem parameters
T=1
K=1;
tau=1/8; sig=sqrt(2*tau)
init_mu_sig_sq=[pi/2 3*pi/4]
quad0=gauss_quadrature.Gauss_Hermite_per(init_mu_sig_sq,40,2*pi)
#########################################
function drift_ff(u,K)
  x=u[1]; mfs=u[2]; mfc=u[3]
  return K*mfs*cos(x)-(K*mfc+1)*sin(x)
end
#########################################
function diff_gg(u,sig)
  return sig
end
#########################################
function mf_r(x)
  return [sin(x); cos(x)]
end
#########################################
fcn=SDELab2.set_fcn_mf(u->drift_ff(u,K),
          u->diff_gg(u,sig),
          x->mf_r(x))
#########################################
######################################### method parameters
solver="MFEM2"
dt=1e-2
params=SDELab2.set_opt(dt,solver)
params["const"]=1           # for no_Gauss_points
params["length_scale"]=2*pi # lam=1/length^2
params["skiphalf"]=true
params["periodic"]=true
params["no_partition"]=5
#########################################
###################### computing reference value
print("Computing a reference value...")
phi=x->[sin(x)^2;sin(x)]
tic()
params["extrap"]=false
dtpts=0.1./(2.0).^(1:6)
params["MaxStepSize"]=dtpts[end]
exact,cpu_time,stats,q1,q0=get_soln(fcn,quad0,params,phi,false)
println(stats["no_pts"])
println("exact=", exact)
#########################################
initial_final_cdf_pdf(q0,q1,[100,40],[0.08,0.01])

#########################################
dtpts=dtpts[1:end-1] # omit last value, used for reference
errors=zeros(length(dtpts),2); errors_no=deepcopy(errors); errors_nd=deepcopy(errors)
cpu_times=zeros(length(dtpts)); cpu_times_no=deepcopy(cpu_times);cpu_times_nd=deepcopy(cpu_times)
params["MaxStepSize"]=dtpts[1]
######################
params["Solver"]="MFEMe"
dtpts_no=dtpts/4
for i=1:(length(dtpts))
    params["MaxStepSize"]=dtpts_no[i]
    qoi,cpu_time,stats,q=get_soln(fcn,quad0,params,phi,false)
    errors_no[i,:],cpu_times_no[i]=helper(qoi, exact, cpu_time, "normal")
end
########################
params["Solver"]="MFEM0"
for i=1:(length(dtpts))
    params["MaxStepSize"]=dtpts[i]
    qoi,cpu_time,stats,q=get_soln(fcn,quad0,params,phi,true)
    errors[i,:],cpu_times[i]=helper(qoi, exact, cpu_time, "normal")
end
################
params["Solver"]="MFEM2"
for i=1:(length(dtpts))
    params["MaxStepSize"]=dtpts[i]
    qoi,cpu_time,stats,q=get_soln(fcn,quad0,params,phi,false)
    errors_nd[i,:],cpu_times_nd[i]=helper(qoi, exact, cpu_time, "normal")
end
#########################################
str1="--" # extrap
str2="k-." # 2nd-order
figure(1);
# first functional
i=1
time_step_err(dtpts,errors[:,i],
  dtpts_no,errors_no[:,i],
  dtpts,errors_nd[:,i],[1,1],utf8("mom1.pdf"))
err_cpu(errors[:,i],cpu_times[:],
  errors_no[:,i],cpu_times_no[:],
  errors_nd[:,i],cpu_times_nd[:],[1,1], utf8("mom1_.pdf"))
# second functional
i=2
time_step_err(dtpts,errors[:,i],
      dtpts_no,errors_no[:,i],
      dtpts,errors_nd[:,i],[1,1],utf8("mom2.pdf"))
err_cpu(errors[:,i], cpu_times[:],
      errors_no[:,i],cpu_times_no[:],
      errors_nd[:,i],cpu_times_nd[:],[1,1], utf8("mom2_.pdf"))
#########################################
println("finished preparing plots")

ext: jl


text:
include("./aa_properties.jl")
include("./which_natural_aminoacid.jl")
include("./selmacros.jl")
include("./select.jl")

export isacidic, isaliphatic, isaromatic, isbackbone, isbasic,
       ischarged, ishydrophobic, isneutral, isnonpolar, ispolar,
       isprotein, issidechain, iswater

export select, selindex


text:
function deep_nnm(tmp_weights::Array{Float64,1},tmp_xx::Array{Float64,2},nodes::Array{Int64,2},obs::Int64,vars::Int64,vars_depe::Int64)

    n_bias=sum(nodes)+vars_depe
    w_bias=(tmp_weights[end-n_bias+1:end])
    w1=(tmp_weights[1:vars*nodes[1]])
    w11=reshape(w1,vars,nodes[1])
    a_sig=sigm1.(tmp_xx*w11+ones(obs,nodes[1]).*w_bias[1:nodes[1]]')

    ind=length(w1)
    for ii=2:length(nodes)
        w1=(tmp_weights[ind+1:ind+nodes[ii-1]*nodes[ii]])
        ind+=length(w1)
        w11=reshape(w1,nodes[ii-1],nodes[ii])
        a_sig=sigm1.(a_sig*w11+ones(obs,nodes[ii]).*w_bias[sum(nodes[1:ii-1])+1:sum(nodes[1:ii])]')
    end

    w1=(tmp_weights[ind+1:end-n_bias])
    w11=reshape(w1,nodes[end],vars_depe)
    out=(a_sig*w11+ones(obs,vars_depe).*w_bias[end-vars_depe+1:end]')

    return out

end

#


text:
__precompile__()

module LCMPolynomials

using LCMCore
using DrakeLCMTypes: polynomial_t, piecewise_polynomial_t, polynomial_matrix_t
using MultivariatePolynomials
import MultivariatePolynomials: variables, nvariables, polynomial, differentiate
using TypedPolynomials

export polynomial_t,
       polynomial_matrix_t,
       piecewise_polynomial_t,
       PiecewiseFunction,
       PiecewisePolynomial,
       cubic_spline

include("piecewise.jl")
using .PiecewiseFunctions

const PiecewisePolynomial = PiecewiseFunction{<:AbstractPolynomial}

include("lcm_conversions.jl")
include("splines.jl")

end # module


text:
using Pkg
Pkg.add("SMM")
using Optim
using HTTP
using GLM
using LinearAlgebra
using Random
using Statistics
using DataFrames
using DataFramesMeta
using CSV

#Question1#
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS1-julia-intro/nlsw88.csv"
df = CSV.read(HTTP.get(url).body, DataFrame)
X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.married.==1

function ols_gmm(beta, X, y)
    ssr = (y.-X*beta)'*(y.-X*beta)
    return ssr
end
g=Y-X*beta
end
beta_optim = optimize(a -> ols_gmm(a, X, y), rand(size(X,2)), LBFGS(), Optim.Options(g_tol=1e-8, iterations=100_000))
println(beta_optim.minimizer)

#Question2#part a#
X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.occupation

function mlogit(alpha, X, y)

    K = size(X,2)
    J = length(unique(y))
    N = length(y)
    bigY = zeros(N,J)
    for j=1:J
        bigY[:,j] = y.==j
    end
    bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

    num = zeros(N,J)
    dem = zeros(N)
    for j=1:J
        num[:,j] = exp.(X*bigAlpha[:,j])
        dem .+= num[:,j]
    end

    P = num./repeat(dem,1,J)

    loglike = -sum( bigY.*log.(P) )

    return loglike
end

alpha_zero = zeros(6*size(X,2))
alpha_rand = rand(6*size(X,2))
alpha_true = [.1910213,-.0335262,.5963968,.4165052,-.1698368,-.0359784,1.30684,-.430997,.6894727,-.0104578,.5231634,-1.492475,-2.26748,-.0053001,1.391402,-.9849661,-1.398468,-.0142969,-.0176531,-1.495123,.2454891,-.0067267,-.5382892,-3.78975]
alpha_start = alpha_true.*rand(size(alpha_true))
println(size(alpha_true))
alpha_hat_optim = optimize(a -> mlogit(a, X, y), alpha_start, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
alpha_hat_mle = alpha_hat_optim.minimizer
println(alpha_hat_mle)

#partb#
unction mlogit_gmm(alpha, X, y)

        K = size(X,2)
        J = length(unique(y))
        N = length(x)
        g = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
        bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

        num = zeros(N,J)
        dem = zeros(N)
        for j=1:J
            num[:,j] = exp.(X*bigAlpha[:,j])
            dem .+= num[:,j]
        end

        P = num./repeat(dem,1,J)

        loglike = -sum( bigY.*log.(P) )

        return loglike
    end
    #part c#

    #we would ideally differentiate to see if the function is globally concave#

    #Question 3#
    randn()

    #Question4#
    using SMM, DataFrames
MA = SMM.parallelNormal() # Note: this line may take up to 5 minutes to execute
dc = SMM.history(MA.chains[1])
dc = dc[dc[:accepted].==true, :]
println(describe(dc))

#Question5#
function ols_smm(θ, X, y, D)
    K = size(X,2)
    N = size(y,1)
    β = θ[1:end-1]
    σ = θ[end]
    if length(β)==1
        β = β[1]
    end
    # N+1 moments in both model and data
    gmodel = zeros(N+1,D)
    # data moments are just the y vector itself
    # and the variance of the y vector
    gdata  = vcat(y,var(y))
    #### !!!!!!!!!!!!!!!!!!!!!!!!!!!!! ####
    # This is critical!                   #
    Random.seed!(1234)                    #
    # You must always use the same ε draw #
    # for every guess of θ!               #
    #### !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ###
    # simulated model moments
    for d=1:D
        ε = σ*randn(N)
        ỹ = X*β .+ ε
        gmodel[1:end-1,d] = ỹ
        gmodel[  end  ,d] = var(ỹ)
    end
    # criterion function
    err = vec(gdata .- mean(gmodel; dims=2))
    # weighting matrix is the identity matrix
    # minimize weighted difference between data and moments
    J = err'*I*err
    return J
end


text:
#using JuMP, Clp

function idm(d :: Array{Int64,1}, s=2.0 :: Float64, logspace=true :: Bool)
    n = length(d)
    p = zeros(2,n)
    p[1,:] = d
    p[2,:] = d .+ s
    p = p ./ (sum(d) .+ s)
    return logspace ? [log.(x) for x in p] : p
end

function lp(c :: Array{Float64,1}, p :: Array{Float64,2}, maximize=true :: Bool, logspace=true :: Bool)
    x = p[:,1]

    i = sortperm(c,rev=maximize)[1]
    x[i] = getindex(p[i,:], 2)
    return logspace ? c'*exp.(x) : c'*x
end

function add_missingness(data :: PlainXData{Bool,BitArray{2}} , ratio :: Float64)
    data2 = convert.(Int8,data.x)
    data2[rand(Float16,size(data.x)) .< ratio] .= -1
    return data2
end
#=
function minimi(coeff :: Array{Float64,1}, l_bounds :: Array{Float64,1}, u_bounds :: Array{Float64,1})
    my_model = Model(Clp.Optimizer)
    set_optimizer_attribute(my_model, "LogLevel", 0)
    set_optimizer_attribute(my_model, "Algorithm", 4)
    @variable(my_model, l_bounds[i] <= x[i = 1:length(l_bounds)]<= u_bounds[i])
    @objective(my_model,Min,coeff' *x) # perché anche * ?
    @constraint(my_model, normalization, ones(Float64, length(l_bounds))'x == 1)
    optimize!(my_model)
    return JuMP.objective_value(my_model)
end
=#

max_stars_repo_name: IDSIA/cjuice

max_issues_repo_name: IDSIA/cjuice

max_forks_repo_name: IDSIA/cjuice


text:
@testset "Checking io-equations: single output" begin
    test_cases = []

    # 2-compartiment model
    push!(test_cases, Dict(
        "ode" => @ODEmodel(
            x0'(t) = -(a01 + a21) * x0(t) + a12 * x1(t) + u(t),
            x1'(t) = a21 * x0(t) - a12 * x1(t),
            y(t) = x0(t)
        ),
        "correct" => "y_2 + y_1*a01 + y_1*a21 + y_1*a12 + y_0*a01*a12 - u_0*a12 - u_1"
    ))

    #---------------------------------------
    # Chen-Lee model
    correct = "9*y1_3^2*y1_0^2 - 18*y1_3*y1_2*y1_1*y1_0 - 18*y1_3*y1_2*y1_0^2*a - 36*y1_3*y1_2*y1_0^2*b - 36*y1_3*y1_2*y1_0^2*c + 18*y1_3*y1_1^2*y1_0*a + 18*y1_3*y1_1^2*y1_0*b + 18*y1_3*y1_1^2*y1_0*c - 24*y1_3*y1_1*y1_0^4 + 18*y1_3*y1_1*y1_0^2*a*b + 18*y1_3*y1_1*y1_0^2*a*c + 18*y1_3*y1_1*y1_0^2*b^2 + 36*y1_3*y1_1*y1_0^2*b*c + 18*y1_3*y1_1*y1_0^2*c^2 + 24*y1_3*y1_0^5*a - 18*y1_3*y1_0^3*a*b^2 - 36*y1_3*y1_0^3*a*b*c - 18*y1_3*y1_0^3*a*c^2 + 9*y1_2^2*y1_1^2 + 18*y1_2^2*y1_1*y1_0*a + 36*y1_2^2*y1_1*y1_0*b + 36*y1_2^2*y1_1*y1_0*c + 9*y1_2^2*y1_0^2*a^2 + 36*y1_2^2*y1_0^2*a*b + 36*y1_2^2*y1_0^2*a*c + 27*y1_2^2*y1_0^2*b^2 + 90*y1_2^2*y1_0^2*b*c + 27*y1_2^2*y1_0^2*c^2 - 18*y1_2*y1_1^3*a - 18*y1_2*y1_1^3*b - 18*y1_2*y1_1^3*c + 24*y1_2*y1_1^2*y1_0^3 - 18*y1_2*y1_1^2*y1_0*a^2 - 72*y1_2*y1_1^2*y1_0*a*b - 72*y1_2*y1_1^2*y1_0*a*c - 54*y1_2*y1_1^2*y1_0*b^2 - 108*y1_2*y1_1^2*y1_0*b*c - 54*y1_2*y1_1^2*y1_0*c^2 + 48*y1_2*y1_1*y1_0^4*b + 48*y1_2*y1_1*y1_0^4*c - 18*y1_2*y1_1*y1_0^2*a^2*b - 18*y1_2*y1_1*y1_0^2*a^2*c - 18*y1_2*y1_1*y1_0^2*a*b^2 - 108*y1_2*y1_1*y1_0^2*a*b*c - 18*y1_2*y1_1*y1_0^2*a*c^2 - 18*y1_2*y1_1*y1_0^2*b^3 - 126*y1_2*y1_1*y1_0^2*b^2*c - 126*y1_2*y1_1*y1_0^2*b*c^2 - 18*y1_2*y1_1*y1_0^2*c^3 - 24*y1_2*y1_0^5*a^2 - 48*y1_2*y1_0^5*a*b - 48*y1_2*y1_0^5*a*c + 18*y1_2*y1_0^3*a^2*b^2 + 36*y1_2*y1_0^3*a^2*b*c + 18*y1_2*y1_0^3*a^2*c^2 + 18*y1_2*y1_0^3*a*b^3 + 126*y1_2*y1_0^3*a*b^2*c + 126*y1_2*y1_0^3*a*b*c^2 + 18*y1_2*y1_0^3*a*c^3 + 9*y1_1^4*a^2 + 18*y1_1^4*a*b + 18*y1_1^4*a*c + 9*y1_1^4*b^2 + 18*y1_1^4*b*c + 9*y1_1^4*c^2 - 24*y1_1^3*y1_0^3*a - 24*y1_1^3*y1_0^3*b - 24*y1_1^3*y1_0^3*c + 18*y1_1^3*y1_0*a^2*b + 18*y1_1^3*y1_0*a^2*c + 36*y1_1^3*y1_0*a*b^2 + 72*y1_1^3*y1_0*a*b*c + 36*y1_1^3*y1_0*a*c^2 + 18*y1_1^3*y1_0*b^3 + 54*y1_1^3*y1_0*b^2*c + 54*y1_1^3*y1_0*b*c^2 + 18*y1_1^3*y1_0*c^3 + 16*y1_1^2*y1_0^6 + 24*y1_1^2*y1_0^4*a^2 - 12*y1_1^2*y1_0^4*b^2 - 72*y1_1^2*y1_0^4*b*c - 12*y1_1^2*y1_0^4*c^2 - 18*y1_1^2*y1_0^2*a^2*b^2 - 18*y1_1^2*y1_0^2*a^2*c^2 - 18*y1_1^2*y1_0^2*a*b^3 + 18*y1_1^2*y1_0^2*a*b^2*c + 18*y1_1^2*y1_0^2*a*b*c^2 - 18*y1_1^2*y1_0^2*a*c^3 + 36*y1_1^2*y1_0^2*b^3*c + 72*y1_1^2*y1_0^2*b^2*c^2 + 36*y1_1^2*y1_0^2*b*c^3 - 32*y1_1*y1_0^7*a + 24*y1_1*y1_0^5*a^2*b + 24*y1_1*y1_0^5*a^2*c + 24*y1_1*y1_0^5*a*b^2 + 144*y1_1*y1_0^5*a*b*c + 24*y1_1*y1_0^5*a*c^2 - 72*y1_1*y1_0^3*a^2*b^2*c - 72*y1_1*y1_0^3*a^2*b*c^2 - 72*y1_1*y1_0^3*a*b^3*c - 144*y1_1*y1_0^3*a*b^2*c^2 - 72*y1_1*y1_0^3*a*b*c^3 + 16*y1_0^8*a^2 - 12*y1_0^6*a^2*b^2 - 72*y1_0^6*a^2*b*c - 12*y1_0^6*a^2*c^2 + 36*y1_0^4*a^2*b^3*c + 72*y1_0^4*a^2*b^2*c^2 + 36*y1_0^4*a^2*b*c^3"

    push!(test_cases, Dict(
        "ode" => @ODEmodel(
            x0'(t) = a * x0(t) - x1(t) * x2(t),
            x1'(t) = b * x1(t) + x0(t) * x2(t),
            x2'(t) = c * x2(t) + 1//3 * x0(t) * x1(t),
            y1(t) = x0(t)
        ),
        "correct" => correct
    ))

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

    # predator-prey model
    correct = "y1_2*y1_0 - y1_1^2 - y1_1*y1_0^2*d + y1_1*y1_0*c + y1_0^3*a*d - y1_0^2*a*c"

    push!(test_cases, Dict(
        "ode" => @ODEmodel(
            x0'(t) = a * x0(t) - b * x0(t) * x1(t),
            x1'(t) = -c * x1(t) + d * x0(t) * x1(t),
            y1(t) = x0(t)
        ),
        "correct" => correct
    ))

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

    for case in test_cases
        ode = case["ode"]
        io_eq = collect(values(find_ioequations(ode)))[1]
        gen_symbols = symbols(parent(io_eq))
        correct_expr = Meta.parse(case["correct"])
        for (i, s) in enumerate(gen_symbols)
            correct_expr = StructuralIdentifiability.MacroTools.postwalk(ex -> ex == s ? :(polyvars[$i]) : ex, correct_expr)
        end
        # dirty hack to eval within the local scope
        correct_function = @eval function(eq)
            polyvars = gens(parent(eq))
            $correct_expr
        end
        correct = correct_function(io_eq)
        divisibility, remainder = divides(io_eq, correct)
        @test divisibility
        @test total_degree(remainder) == 0
    end
end


text:
"""
    codegen(fun, argtypes; optimize_lowering = true, triple = nothing, datalayout = nothing)

Return the bitcode for `fun` with argument types `argtypes`.
`argtypes` is a tuple type (e.g. `Tuple{Float64, Int}`).

The optional argument `optimize_lowering` determines whether
`code_typed` uses optimization when infering types.

The optional arguments `triple` and `datalayout` are strings that
set the LLVM triple and data layout for the module being generated.

`codegen` creates a `CodeCtx` and runs `codegen!(cg)`.
Direct calls to that can be done for more flexibility.

The return value is an LLVM module. This can be written to a bitcode
files with `write(mod, filepath)`. It can be optimized with
`optimize!(mod)`.
"""
function codegen(@nospecialize(fun), @nospecialize(argtypes); optimize_lowering = true, triple = nothing, datalayout = nothing)
    ci, dt = code_typed(fun, argtypes, optimize = optimize_lowering)[1]
    sig = first(methods(fun, argtypes)).sig
    funname = getfunname(fun, argtypes)
    cg = CodeCtx(funname, ci, dt, argtypes, fun, sig)
    return codegen!(cg)
end
getfunname(fun, argtypes) = string(basename(fun), "_", join(collect(argtypes.parameters), "_"))
basename(f::Function) = Base.function_name(f)
basename(f::Core.IntrinsicFunction) = Symbol(unsafe_string(ccall(:jl_intrinsic_name, Cstring, (Core.IntrinsicFunction,), f)))
basename(x::GlobalRef) = x.name
basename(m::Core.MethodInstance) = basename(m.def)
basename(m::Method) = m.name == :Type ? m.sig.parameters[1].parameters[1].name.name : m.name

#
# Main entry
#

"""
    codegen!(cg::CodeCtx)
    codegen!(cg::CodeCtx, ...)

This is the main function for dispatching various types of code generation.
"""
function codegen!(cg::CodeCtx)
    @info "## $(cg.name)"
    ci = cg.code_info
    lastsig(x) = last(x.parameters)
    lastsig(x::UnionAll) = lastsig(x.body)
    siglen(x) = length(x.parameters)
    siglen(x::UnionAll) = siglen(x.body)
    sigend = lastsig(cg.sig)
    hasvararg = isa(sigend, UnionAll) && sigend.body.name.name == :Vararg
    if hasvararg
        siglength = siglen(cg.sig) - 1
        argtypes = LLVMType[llvmtype(p) for p in cg.argtypes.parameters[1:siglength-1]]
        push!(argtypes, llvmtype(Tuple{cg.argtypes.parameters[siglength:end]...}))
    else
        argtypes = LLVMType[llvmtype(p) for p in cg.argtypes.parameters]
    end
    for i in 1:length(argtypes)
        if isa(argtypes[i], LLVM.VoidType)
            argtypes[i] = llvmtype(Int)
        end
    end
    func_type = LLVM.FunctionType(llvmtype(cg.result_type), argtypes)
    cg.func = LLVM.Function(cg.mod, cg.name, func_type)
    cg.nargs = length(argtypes)
    LLVM.linkage!(cg.func, LLVM.API.LLVMExternalLinkage)
    entry = LLVM.BasicBlock(cg.func, "entry", ctx)
    LLVM.position!(cg.builder, entry)
    for (i, param) in enumerate(LLVM.parameters(cg.func))
        argname = string(ci.slotnames[i + 1], "_as", i+1)
        if !isa(LLVM.llvmtype(param), LLVM.VoidType)
            alloc = LLVM.alloca!(cg.builder, LLVM.llvmtype(param), argname)
            store!(cg, param, alloc)
            cg.slotlocs[i + 1] = alloc
        end
    end
    for i in cg.nargs+2:length(ci.slotnames)
        varname = string(ci.slotnames[i], "_s", i)
        vartype = llvmtype(ci.slottypes[i])
        if !isa(vartype, LLVM.VoidType)
            alloc = LLVM.alloca!(cg.builder, vartype, varname)
            cg.slotlocs[i] = alloc
        end
    end
    for (i, node) in enumerate(ci.code)
        @debug "$(cg.name): ***** node $i/$(length(ci.code)) ***** " node
        codegen!(cg, node)
    end
    if !has_terminator(entry)
        LLVM.unreachable!(cg.builder)
    end
    # LLVM.verify(cg.func)
    LLVM.dispose(cg.builder)
    return cg.mod
end

function codegen!(cg::CodeCtx, @nospecialize(fun), @nospecialize(argtypes); optimize_lowering = true)
    ci, dt = code_typed(fun, argtypes, optimize = optimize_lowering)[1]
    funname = getfunname(fun, argtypes)
    sig = first(methods(fun, argtypes)).sig
    return codegen!(CodeCtx(cg, funname, ci, dt, argtypes, fun, sig))
end

#
# Expressions
#
function codegen!(cg::CodeCtx, e::Expr)
    # Slow dispatches here but easy to write and to customize
    codegen!(cg, Val(e.head), e.args, e.typ)
end

#
# Variable assignment and reading
#

function codegen!(cg::CodeCtx, x::SlotNumber)
    varname = string(cg.code_info.slotnames[x.id], "_s", x.id)
    # cg.slots[v.id]
    # p = get(current_scope(cg), varname, nothing)
    p = cg.slotlocs[x.id]
    p == nothing && error("did not find variable $(varname)")
    return LLVM.load!(cg.builder, p, varname*"_")
end

codegen!(cg::CodeCtx, x::SSAValue) = LLVM.load!(cg.builder, cg.ssas[x.id])

function codegen!(cg::CodeCtx, ::Val{:(=)}, args, typ)
    result = codegen!(cg, args[2])
    if isa(args[1], SlotNumber)
        varname = string(cg.code_info.slotnames[args[1].id], "_s", args[1].id)
        @debug "$(cg.name): Assigning slot $(args[1].id)" varname args[2] _typeof(cg, args[1]) _typeof(cg, args[2])
        p = cg.slotlocs[args[1].id]
        unboxed_result = emit_unbox!(cg, result, _typeof(cg, args[1]))
        if _typeof(cg, args[2]) == Union{}
            return
        end
        store!(cg, unboxed_result, p)
        cg.slots[args[1].id] = unboxed_result
        return unboxed_result
    end
    if isa(args[1], SSAValue)
        t = _typeof(cg, args[2])
        # if !isconcrete(t) # || sizeof(t) == 0
        #     return
        # end
        @debug "$(cg.name): Assigning SSA $(args[1].id)" typ llvmtype(_typeof(cg, args[1]))
        # @info "$(cg.name): Assigning SSA $(args[1].id) $typ $result $(llvmtype(_typeof(cg, args[1])))"
        # dump(args)
        unboxed_result = emit_unbox!(cg, result, _typeof(cg, args[1]))
        p = alloca!(cg.builder, llvmtype(_typeof(cg, args[1])))
        store!(cg, unboxed_result, p)
        cg.ssas[args[1].id] = p
        return unboxed_result
    end
end

function codegen!(cg::CodeCtx, f::GlobalRef)
    if haskey(LLVM.functions(cg.mod), string(f.name))
        return LLVM.functions(cg.mod)[string(f.name)]
    end
    evf = eval(f)
    if isa(evf, Type)
        return codegen!(cg, evf)
    else
        # If we refer to something global that we don't have, need to add it as a global variable.
        # return emit_box!(cg, Int32(999)) # KLUDGE - WRONG
        return codegen!(cg, 999) # KLUDGE - WRONG
    end
end

codegen!(cg::CodeCtx, x::QuoteNode) = codegen!(cg, x.value)
codegen!(cg::CodeCtx, x::Symbol) = emit_symbol!(cg, x)
function codegen!(cg::CodeCtx, x::T) where {T}
    if isbits(T) && sizeof(T) == 0
        return codegen!(cg, 0)   # I'm not sure how to handle singleton types
    end
    error("Unsupported type: $T")
end
function codegen!(cg::CodeCtx, x::T) where {T<:Union{LLVM.LoadInst,LLVM.ConstantInt}}
    return x
end

#
# Function calls
#
function codegen!(cg::CodeCtx, ::Val{:call}, args, typ)
    llvmargs = LLVM.Value[]
    llvmargs = Any[]
    fun = eval(args[1])
    name = string(args[1])
    if isa(fun, Core.IntrinsicFunction)
        @debug "$(cg.name): calling intrinsic: $name"
        # @info "$(cg.name): calling intrinsic: $name"
        # dump(args)
        # @show basename(args[1])
        return emit_intrinsic!(cg, basename(args[1]), args[2:end])
    end
    if isa(fun, Core.Builtin)
        @debug "$(cg.name): calling builtin: $name"
        return emit_unbox!(cg, emit_builtin!(cg, basename(args[1]), args[2:end], typ), typ)
    end
    @debug "$(cg.name): calling other method: $name" args
    ## NOTE: everything past here may be wrong!
    argstypetuple = Tuple{(gettypes(cg, a) for a in args[2:end])...}
    @debug "$(cg.name): more info" argstypetuple
    # argstypetuple = Tuple{(Any for a in args[2:end])...}
    global method = which(fun, argstypetuple)
    # codegen!(cg, Val(:invoke), Any[method.specializations.func, args[2:end]...], typ)
    codegen!(cg, Val(:invoke), Any[method, args[2:end]...], typ)
    # error("Function $name not supported.")
end
gettypes(cg::CodeCtx, x::SlotNumber) = cg.code_info.slottypes[x.id]
# gettypes(cg::CodeCtx, x::GlobalRef) = eval(x)
gettypes(cg::CodeCtx, x::GlobalRef) = Type{Any}
gettypes(cg::CodeCtx, x::Expr) = Type{Any}
gettypes(cg::CodeCtx, x) = typeof(x)

function codegen!(cg::CodeCtx, ::Val{:invoke}, args, typ)
    # name = string(Base.function_name(args[1]))
    argtypes = getargtypes(args[1])
    name = getfunname(args[1], argtypes)
    @info "$(cg.name): invoking $name"
    @debug "$(cg.name): invoking $name" args typ argtypes name
    # dump(args, maxdepth=2)
    if haskey(LLVM.functions(cg.mod), name)
        func = LLVM.functions(cg.mod)[name]
    elseif name == "throw_overflowerr_binaryop_Symbol_Int32_Int32"
        return codegen!(cg, Int32(888)) # KLUDGE - WRONG
    else
        fun = getfun(args[1])
        global MI = args[1]
        @debug "$(cg.name): invoke argtypes" args argtypes
        if isa(args[1], Core.MethodInstance)  && isdefined(args[1], :inferred) && args[1].inferred != nothing
            MI = args[1]
            ci = Base.uncompressed_ast(MI.def, MI.inferred)
            dt = MI.rettype
            sig = MI.def.sig
        else
            methodtable = code_typed(fun, argtypes, optimize = true)
            if length(methodtable) < 1
                @show fun
                @show argtypes
                dump(args, maxdepth=2)
                error("Method not found")
            end
            ci, dt = methodtable[1]
            sig = first(methods(fun, argtypes)).sig
        end
        newcg = CodeCtx(cg, name, ci, dt, argtypes, fun, sig)
        codegen!(newcg)
        func = newcg.func
    end
    llvmargs = LLVM.Value[]
    startpos = isa(args[1], Core.MethodInstance) ? 3 : 2
    for (i,v) in enumerate(args[startpos:end])
        a = codegen!(cg, v)
        if i <= length(argtypes.parameters) && llvmtype(argtypes.parameters[i]) != LLVM.llvmtype(a)  # kludgey?
            a = emit_box!(cg, v)
        end
        push!(llvmargs, a)
    end
    return LLVM.call!(cg.builder, func, llvmargs)
end

# getname(x::Core.MethodInstance) = getfield(x.def.module, x.def.name)
function getname(x::Core.MethodInstance)
    nm = string(getfield(x.def.module, x.def.name))
    if nm != "Type"
        return nm
    end
    return getname(first(x.def.sig.parameters))
end

function getname(x::Method)
    nm = string(getfield(x.module, x.name))
    if nm != "Type"
        return nm
    end
    first(x.sig.parameters)
    getname(first(x.sig.parameters))
    return getname(first(x.sig.parameters))
end
# getname(x::Method) = getfield(x.module, getname(x.sig.parameters[1]))
getname(::Type{T}) where T = string(T.parameters[1].name.name)
getname(::T) where T = string(T.name.name)
getargtypes(x::Core.MethodInstance) = Tuple{x.specTypes.parameters[2:end]...}
getargtypes(x::Method) = Tuple{x.sig.parameters[2:end]...}
getfun(x::Core.MethodInstance) = getfun(x.def)
getfun(x::Method) = getfield(x.module, basename(x))

is_nothing(x) = false
is_nothing(::Nothing) = true
is_nothing(x::GlobalRef) = x.name == :nothing

function codegen!(cg::CodeCtx, ::Val{:return}, args, typ)
    if length(args) == 1 && !is_nothing(args[1])
        res = codegen!(cg, args[1])
        if LLVM.llvmtype(res) == int1_t
            res = LLVM.zext!(cg.builder, res, int8_t)
        elseif LLVM.llvmtype(res) != llvmtype(cg.result_type)
            res = emit_unbox!(cg, res, cg.result_type)
        end
        LLVM.ret!(cg.builder, res)
    else
        LLVM.ret!(cg.builder)
    end
end

#
# Constants / values
#

codegen!(cg::CodeCtx, x::T) where T <: Base.IEEEFloat =
    LLVM.ConstantFP(llvmtype(T), x)

codegen!(cg::CodeCtx, x::T) where T <: Base.BitInteger =
    LLVM.ConstantInt(llvmtype(T), x)

codegen!(cg::CodeCtx, x::Bool) =
    LLVM.ConstantInt(bool_t, Int8(x))

function codegen!(cg::CodeCtx, s::String)
    # strinst = globalstring!(builder, s)
    gs = globalstring_ptr!(cg.builder, s)
    return LLVM.call!(cg.builder, cg.extern[:jl_pchar_to_string], LLVM.Value[gs, codegen!(cg, Int32(length(s)+1))])
end

function codegen!(cg::CodeCtx, x::Tuple)
    if isbits(x)
        loc = alloca!(cg.builder, llvmtype(typeof(x)))
        for i in 1:length(x)
            p = LLVM.struct_gep!(cg.builder, loc, i-1)
            store!(cg, codegen!(cg, x[i]), p)
        end
        return LLVM.load!(cg.builder, loc)
    end
end

function codegen!(cg::CodeCtx, x::Array{T,N}) where {T,N}
    # typ = llvmtype(T)
    # array_type = LLVM.call!(cg.builder, cg.extern[:jl_apply_array_type], LLVM.Value[typ, codegen!(cg, UInt32(N))])
    # tuple_type = LLVM.call!(cg.builder, cg.extern[:jl_apply_tuple_type], LLVM.Value[])
    # dims_tuple = jl_new_struct(tuple_type, jl_box...(arg1), jl_box...(arg2), ...);
    # array = LLVM.call!(cg.builder, cg.extern[:jl_new_array], LLVM.Value[array_type, codegen!(cg, UInt32(N))])
    # array = jl_new_array_1d(array_type, length(x))
    # typ = LLVM.load!(cg.builder, cg.datatype[T])
    # return
end

#
# Miscellaneous
#
codegen!(cg::CodeCtx, ::LineNumberNode) = nothing

codegen!(cg::CodeCtx, ::NewvarNode) = nothing

codegen!(cg::CodeCtx, ::Nothing) = cg.datatype[Nothing]

codegen!(cg::CodeCtx, ::Tuple{}) = cg.datatype[Tuple{}]

codegen!(cg::CodeCtx, ::Val{:meta}, args, typ) = nothing

codegen!(cg::CodeCtx, ::Val{:simdloop}, args, typ) = nothing
codegen!(cg::CodeCtx, ::Val{:gc_preserve_begin}, args, typ) = codegen!(cg, args[1])
codegen!(cg::CodeCtx, ::Val{:gc_preserve_end}, args, typ) = codegen!(cg, args[1])

function codegen!(cg::CodeCtx, ::Val{:static_parameter}, args, typ)
    @debug "$(cg.name) emit static_parameter found: $args"
    # Adapted from https://github.com/SimonDanisch/Sugar.jl/blob/afa3da222b3b3eb1df3fe1be76b4d892e91c9d98/src/methods.jl#L350-L357
    world = typemax(UInt)
    y = Base._methods(cg.jfun, cg.argtypes, -1, world)
    isempty(y) && return ()
    x = first(y)
    val = tuple(x[2]...)[args[1]]
    codegen!(cg, val)
end

function codegen!(cg::CodeCtx, ::Val{:boundscheck}, args, typ)
    return codegen!(cg, 999) # KLUDGE - WRONG
end

#
# ccall
#
function codegen!(cg::CodeCtx, ::Val{:foreigncall}, args, typ)
    name = isa(args[1], Tuple) ? nameof(args[1][1]) : nameof(args[1])
    if haskey(cg.extern, name)
        @debug "$(cg.name) ccall found: $name"
        func = cg.extern[name]
    else
        @debug "$(cg.name) ccall creating: $name"
        func = extern!(cg.mod, name, llvmtype(args[2]), LLVMType[llvmtype(x) for x in collect(args[3])])
        cg.extern[name] = func
    end
    llvmargs = LLVM.Value[]
    for v in args[6:5+args[5]]
        push!(llvmargs, codegen!(cg, v))
    end

    return LLVM.call!(cg.builder, func, llvmargs)
end
nameof(x::QuoteNode) = nameof(x.value)
nameof(x::String) = Symbol(x)
nameof(x) = x

function codegen!(cg::CodeCtx, ::Val{:extcall}, name, rettype, args)
    llvmargs = LLVM.Value[]
    llvmargtypes = LLVM.Value[]
    for v in args
        push!(llvmargs, llvmtype(typeof(v)))
        push!(llvmargs, codegen!(cg, v))
    end

    func_type = LLVM.FunctionType(llvmtype(rettype), llvmargtypes)
    func = LLVM.Function(cg.mod, name, func_type)
    LLVM.linkage!(func, LLVM.API.LLVMExternalLinkage)

    return LLVM.call!(cg.builder, func, llvmargs)
end

#
# Control flow
#
function codegen!(cg::CodeCtx, ln::LabelNode)
    if !haskey(cg.labels, ln.label)
        func = LLVM.parent(LLVM.position(cg.builder))
        cg.labels[ln.label] = LLVM.BasicBlock(func, "L", ctx)
    end
    if !has_terminator(position(cg.builder))
        br!(cg.builder, cg.labels[ln.label])
    end
    position!(cg.builder, cg.labels[ln.label])
end

function codegen!(cg::CodeCtx, gn::GotoNode)
    if !haskey(cg.labels, gn.label)
        func = LLVM.parent(LLVM.position(cg.builder))
        cg.labels[gn.label] = LLVM.BasicBlock(func, "L", ctx)
    end
    br!(cg.builder, cg.labels[gn.label])
end

function codegen!(cg::CodeCtx, ::Val{:gotoifnot}, args, typ)
    condv = emit_condition!(cg, codegen!(cg, args[1]))
    func = LLVM.parent(LLVM.position(cg.builder))
    ifso = LLVM.BasicBlock(func, "if", ctx)
    ifnot = LLVM.BasicBlock(func, "L", ctx)
    cg.labels[args[2]] = ifnot
    LLVM.br!(cg.builder, condv, ifso, ifnot)
    position!(cg.builder, ifso)
end

emit_condition!(cg, condv) = LLVM.width(LLVM.llvmtype(condv)) > 1 ?
        LLVM.trunc!(cg.builder, condv, int1_t) : condv

#
# New - structure creation
#
function codegen!(cg::CodeCtx, ::Val{:new}, args, typ)
    typ = eval(args[1])
    if isbits(typ)
        @debug "$(cg.name): creating bitstype" args typ llvmtype(typ)
        loc = alloca!(cg.builder, llvmtype(typ))
        for i in 2:length(args)
            p = LLVM.struct_gep!(cg.builder, loc, i-2)
            store!(cg, codegen!(cg, args[i]), p)
        end
        return LLVM.load!(cg.builder, loc)
    else
        dt = load_and_emit_datatype!(cg, args[1])
        @debug "$(cg.name): creating composite" args typ
        llvmargs = LLVM.Value[dt]
        for a in args[2:end]
            push!(llvmargs, emit_box!(cg, a))
        end
        return LLVM.call!(cg.builder, cg.extern[:jl_new_struct], llvmargs)
    end
end

#
# Types - emit and return a stored type or create a new type and store it
#

# The following creates the appropriate structures in memory
function codegen!(cg::CodeCtx, ::Type{Array{T,N}}) where {T,N}
    typ = LLVM.load!(cg.builder, cg.datatype[T])
    return LLVM.call!(cg.builder, cg.extern[:jl_apply_array_type], LLVM.Value[typ, codegen!(cg, UInt32(N))])
end

codegen!(cg::CodeCtx, x::Type{T}) where T = load_and_emit_datatype!(cg, x)

function load_and_emit_datatype!(cg, ::Type{JT}) where JT
    if haskey(cg.datatype, JT)
        return LLVM.load!(cg.builder, cg.datatype[JT])
        # return cg.datatype[JT]
    end
    name = string(JT)
    # @info "$(cg.name): emitting new type: $name"
    typeof(JT) != DataType && error("Not supported, yet")
        # JL_DLLEXPORT jl_value_t *jl_type_union(jl_value_t **ts, size_t n);

    lname = emit_symbol!(cg, JT.name)
    mod = LLVM.load!(cg.builder, cg.extern[:jl_main_module_g])
    super = load_and_emit_datatype!(cg, JT.super)
    params = LLVM.call!(cg.builder, cg.extern[:jl_svec],
        LLVM.Value[codegen!(cg, length(JT.parameters)), [load_and_emit_datatype!(cg, t) for t in JT.parameters]...])
    if isconcrete(JT)
        fnames = LLVM.call!(cg.builder, cg.extern[:jl_svec],
            LLVM.Value[codegen!(cg, length(fieldnames(JT))), [emit_symbol!(cg, s) for s in fieldnames(JT)]...])
        ftypes = LLVM.call!(cg.builder, cg.extern[:jl_svec],
            LLVM.Value[codegen!(cg, length(JT.types)), [load_and_emit_datatype!(cg, t) for t in JT.types]...])
    else
        fnames = ftypes = LLVM.call!(cg.builder, cg.extern[:jl_svec], LLVM.Value[codegen!(cg, 0)])
    end
    abstrct = codegen!(cg, UInt32(JT.abstract))
    mutabl = JT.mutable ? codegen!(cg, UInt32(1)) : codegen!(cg, UInt32(0))
    ninitialized = codegen!(cg, UInt32(JT.ninitialized))
    dt = LLVM.call!(cg.builder, cg.extern[:jl_new_datatype],
        LLVM.Value[lname, mod, super, params, fnames, ftypes, abstrct, mutabl, ninitialized])
        # LLVMType[jl_sym_t_ptr, jl_module_t_ptr, jl_datatype_t_ptr, jl_svec_t_ptr, jl_svec_t_ptr, jl_svec_t_ptr, int32_t, int32_t, int32_t])
    loc = LLVM.GlobalVariable(cg.mod, jl_datatype_t_ptr, string(JT))
    LLVM.linkage!(loc, LLVM.API.LLVMCommonLinkage)
    # LLVM.initializer!(loc, null(jl_value_t_ptr))
    LLVM.API.LLVMSetInitializer(LLVM.ref(loc), LLVM.ref(null(jl_value_t_ptr)))
    store!(cg, dt, loc)
    result = LLVM.load!(cg.builder, loc)
    # cg.datatype[JT] = result
    cg.datatype[JT] = loc
    return result
end

load_and_emit_datatype!(cg, x::GlobalRef) = load_and_emit_datatype!(cg, eval(x))

load_and_emit_datatype!(cg, x) = codegen!(cg, x)

emit_symbol!(cg, x) =
    LLVM.call!(cg.builder, cg.extern[:jl_symbol], [globalstring_ptr!(cg.builder, string(x))])


text:
# Starting example file for working with chars and strings in Julia

# Julia has a specific character type
myChar = 'x' # use single quotation mark for character
println(Int(myChar))
println(Char(120))

# Strings are defined using double quotes or triple quotes
myStr = "This is a sample string in Julia"
myOtherStr = """
    Hello There
    This is a string
    """

# Get the length of a string
println(length(myStr))

# Access individual characters - note that they are 1-indexed
println(myStr[1]) # obtains the first character of myStr
println(myStr[end]) # obtains the last character of myStr
println(myStr[end - 3]) # obtains the fourth last character of myStr

# Slicing strings is used with the : notation
sliceMyStr = myStr[2:end - 1] # slices myStr from the second to the second to last character
println(typeof(sliceMyStr)) # sliced type is still S tring

# Iterate over characters
for c in myStr
    print(c)
end
println()

# String concatenation using *
w1 = "Hello"
w2 = "world"
println(w1 * ", " * w2 * "! ")

# String interpolation
a = 5
b = 10
println("The result of $a + $b is $(a+b)")

max_stars_repo_name: scaotravis/Julia-Course

max_issues_repo_name: scaotravis/Julia-Course

max_forks_repo_name: scaotravis/Julia-Course


text:
using Documenter, ImageQuilting

makedocs()

deploydocs(
  deps  = Deps.pip("mkdocs"),
  repo  = "github.com/juliohm/ImageQuilting.jl.git",
  julia = "release"
)


text:
#!/usr/bin/env julia

x = [3.0, 4.0]
n = Cint(length(x))
out = Ref{Cdouble}(0)

# TODO: Follow this
# http://docs.julialang.org/en/stable/manual/calling-c-and-fortran-code/#passing-pointers-for-modifying-inputs
# Not obvious how to do this with arrays rather than scalars.
ccall((:c_l2norm, "l2norm"), Void
      , (Ref{Cdouble}, Ref{Cint}, Ref{Cdouble})
      , Ref(x), Ref(n), Ref(out)
      )
      #, Ref(x), Ref(n), Ref(out))

print(out, "\n")

max_stars_repo_name: clarkfitzg/templates

max_issues_repo_name: clarkfitzg/templates

max_forks_repo_name: clarkfitzg/templates


text:
using Test,KnetLayers, Knet
gpu() < 0 ?  KnetLayers.settype!(Array{Float64}) : KnetLayers.settype!(KnetArray{Float64})
println("Testing with $(KnetLayers.arrtype)")


text:
"""
    GenerativeBeliefMDP(pomdp, updater)

Create a generative model of the belief MDP corresponding to POMDP `pomdp` with belief updates performed by `updater`.
"""
struct GenerativeBeliefMDP{P<:POMDP, U<:Updater, B, A} <: MDP{B, A}
    pomdp::P
    updater::U
end

function GenerativeBeliefMDP(pomdp::P, up::U) where {P<:POMDP, U<:Updater}
    # XXX hack to determine belief type
    b0 = initialize_belief(up, initialstate(pomdp))
    GenerativeBeliefMDP{P, U, typeof(b0), actiontype(pomdp)}(pomdp, up)
end

function POMDPs.gen(bmdp::GenerativeBeliefMDP, b, a, rng::AbstractRNG)
    s = rand(rng, b)
    if isterminal(bmdp.pomdp, s)
        bp = gbmdp_handle_terminal(bmdp.pomdp, bmdp.updater, b, s, a, rng::AbstractRNG)::typeof(b)
        return (sp=bp, r=0.0)
    end
    sp, o, r = @gen(:sp, :o, :r)(bmdp.pomdp, s, a, rng) # maybe this should have been generate_or?
    bp = update(bmdp.updater, b, a, o)
    return (sp=bp, r=r)
end

actions(bmdp::GenerativeBeliefMDP{P,U,B,A}, b::B) where {P,U,B,A} = actions(bmdp.pomdp, b)
actions(bmdp::GenerativeBeliefMDP) = actions(bmdp.pomdp)

isterminal(bmdp::GenerativeBeliefMDP, b) = all(isterminal(bmdp.pomdp, s) for s in support(b))

discount(bmdp::GenerativeBeliefMDP) = discount(bmdp.pomdp)

const warned_about_gbmdp_terminal=false

# override this if you want to handle it in a special way
function gbmdp_handle_terminal(pomdp::POMDP, updater::Updater, b, s, a, rng)
    global warned_about_gbmdp_terminal
    if !warned_about_gbmdp_terminal
        @warn("""
             Sampled a terminal state for a GenerativeBeliefMDP transition - not sure how to proceed, but will try.

             See $(@__FILE__) and implement a new method of POMDPToolbox.gbmdp_handle_terminal if you want special behavior in this case.

             """)
        warned_about_gbmdp_terminal = true
    end
    sp, o, r = @gen(:sp, :o, :r)(pomdp, s, a, rng)
    bp = update(updater, b, a, o)
    return bp
end

function initialstate(bmdp::GenerativeBeliefMDP)
    return Deterministic(initialize_belief(bmdp.updater, initialstate_distribution(bmdp.pomdp)))
end

# deprecated in POMDPs v0.9
function initialstate(bmdp::GenerativeBeliefMDP, rng::AbstractRNG)
    return initialize_belief(bmdp.updater, initialstate_distribution(bmdp.pomdp))
end


text:
@testset "Dictionary keys" begin
    s = """
Dict(:a=>1, :b=>2, :a=>3)
"""
    msgs = lintstr(s)
    @test msgs[1].code == :E334

    s = """
Dict{Symbol,Int}(:a=>1, :b=>"")
"""
    msgs = lintstr(s)
    @test msgs[1].code == :E532

    s = """
Dict{Symbol,Int}(:a=>1, "b"=>2)
"""
    msgs = lintstr(s)
    @test msgs[1].code == :E531
end


text:
export HDPHyperparam

"Hierarchical Dirichlet Process Mixture Model Hyperparameters"
struct HDPHyperparam <: AbstractHyperparam

  γ_a::Float64
  γ_b::Float64

  α_a::Float64
  α_b::Float64

  "default values"
  HDPHyperparam(;γ_a = 1.0, γ_b = 1.0, α_a = 1.0, α_b = 1.0) = new(γ_a, γ_b, α_a, α_b)

end

"Hierarchical Dirichlet Process Mixture Model Data Object"
mutable struct HDPData <: AbstractModelData

  # Energy
  energy::Float64

  # Distributions
  distributions::Vector{ConjugatePostDistribution}

  # Assignments
  assignments::Vector{Vector{Int}}

	# Weights
  weights::Vector{Vector{Float64}}

end

mutable struct HDPBuffer{T <: Real} <: AbstractModelBuffer

  # samples
  X::Vector{Vector{T}}

  # number of groups
  N0::Int

  # indecies of groups
  N0idx::Vector{Int}

  # number of samples per group
  Nj::Vector{Int}

  # indecies of samples per group
  Njidx::Array{Vector}

  # assignments
  Z::Vector{Vector{Int}}

  # number of active cluster
  K::Int

  # number of samples per group per cluster
  C::Array{Int, 2}

  # total number of tables
  totalnt::Array{Int}

  # number of clusters per table
  classnt::Array{Int}

  # distributions
  G::Vector{ConjugatePostDistribution}

  # base distribution
  G0::ConjugatePostDistribution

  # concentration parameter
  α::Float64

  # concentration parameters of clusters
  β::Array{Float64}

  # gamma
  γ::Float64
end

"HDP initialization using random assignments."
function init(X::Vector{Vector{T}}, model::HDP, init::RandomInitialisation) where T <: Real

    Z = Vector{Vector{Int}}(length(X))
    G = ConjugatePostDistribution[]

    K = init.k

    for i in 1:length(X)
      N = length(X[i])
      Z[i] = rand(1:K, N)

      for c in 1:K
        idx = findall(Z[i] .== c)
        if length(idx) > 0
          if size(G, 1) < c
            push!(G, add(model.H, X[i][idx]))
          else
            add!(G[c], X[i][idx])
          end
        end
      end
    end

    N0 = length(X)
    N0idx = collect(1:N0)

    Nj = [length(x) for x in X]
    Njidx = [collect(1:N) for N in Nj]

    C = reduce(hcat, [counts(z, 1:K) for z in Z])

    # init step
    β = ones(K) ./ K
    classnt = randnumtable(model.α * β[:,ones(Int, N0)]', C')
    totalnt = sum(classnt, 1)

    # set alpha vector of Dirichlet Distribution to sample β
    a = zeros(K + 1)
    a[1:end-1] = totalnt
    a[end] = model.γ

    # update beta
    β = rand(Dirichlet(a));

    return HDPBuffer{T}(
      X,
      N0,
      N0idx,
      Nj,
      Njidx,
      Z,
      K,
      C,
      totalnt,
      classnt,
      G,
      model.H,
      model.α,
      β,
      model.γ)
end

"Single Gibbs sweep of HDP training using beta variables."
function gibbs!(B::HDPBuffer)

    shuffle!(B.N0idx)
    prob = zeros(B.K + 1) * -Inf

    for j in B.N0idx

        # get samples of group
        data = B.X[j]
        z = B.Z[j]

        shuffle!(B.Njidx[j])

        for i in B.Njidx[j]

            cluster = z[i]

            # remove data item from model
            remove!(B.G[cluster], data[i])
            B.C[cluster, j] -= 1

            if isdistempty(B.G[cluster])

                # remove cluster
                B.G = B.G[[1:cluster-1; cluster+1:end]]

                for jj in 1:B.N0
                    B.Z[jj][B.Z[jj] .> cluster] -= 1
                end

                B.β = B.β[[1:cluster-1; cluster+1:end]]

                B.classnt = B.classnt[:,[1:cluster-1; cluster+1:end]]
                B.totalnt = B.totalnt[[1:cluster-1; cluster+1:end]]'

                B.K -= 1

                B.C = B.C[[1:cluster-1; cluster+1:end],:]
                prob = zeros(B.K + 1) * -Inf
            end

            # compute log likelihood
            for k in 1:B.K
              llh = logpostpred(B.G[k], data[i])[1]
              prior = B.C[k, j] + B.β[k] * B.α
              prob[k] = llh + log( prior )
            end

            prob[B.K + 1] = logpostpred(B.G0, data[i])[1] + log( B.β[B.K+1] * B.α )
            prob = exp(prob - maximum(prob))

            c = randomindex(prob)

            # add data to model
            B.Z[j][i] = c

            if c > B.K
                # create new cluster
                B.K += 1
                B.G = vcat(B.G, deepcopy(B.G0))
                b = rand(Dirichlet([1, B.γ]))
                b = b * B.β[end]
                B.β = vcat(B.β, 1)
                B.β[end-1:end] = b
                B.C = vcat(B.C, zeros(Int, 1, B.N0))
                prob = zeros(B.K + 1) * -Inf
            end

            B.C[c, j] += 1
            add!(B.G[c], [data[i]])

        end

    end

    # sample number of tables
    kk = maximum([0, B.K - length(B.totalnt)])
    B.totalnt = hcat(B.totalnt - sum(B.classnt, 1), zeros(Int, 1, kk))
    B.classnt = randnumtable(B.α .* B.β[:,ones(Int, B.N0)]', B.C')
    B.totalnt = B.totalnt + sum(B.classnt, 1)

    # update beta weights
    a = zeros(B.K + 1)
    a[1:end-1] = B.totalnt
    a[end] = B.γ

    B.β = rand(Dirichlet(a))

    B
end

function updateenergy!(B::HDPData, X::Vector{Vector{T}}) where T <: Real

  E = 0
  numl = 0

  for j in length(X)

      # get samples of group
      data = X[j]

      pp = 0.0
      for i in length(data)
        for k = 1:length(B.distributions)
          pp += logpostpred( B.distributions[k], data[i] )[1] + log(B.weights[j][k])
        end
      end
      E += pp
      numl += 1

  end

  B.energy = E - log(numl)

end

function sampleparameters!(B::HDPBuffer, P::HDPHyperparam)
  totalnt = sum(B.classnt, 1)
  B.γ = resampleα(B.γ, sum(B.totalnt), B.K, k = P.γ_a, Θ = P.γ_b, maxiter = 10)
  B.α = resampleα(B.α, B.Nj, B.totalnt, k = P.α_a, Θ = P.α_b, maxiter = 10)
end

function extractpointestimate(B::HDPBuffer)
  W = map(j -> map(k -> B.C[k, j] + B.β[k] * B.α, 1:B.K), B.N0idx)
  model = HDPData(0.0, deepcopy(B.G), deepcopy(B.Z), deepcopy(W) )
  updateenergy!(model, B.X)
  model
end


text:
"""
    MetaAbstractPurse{T}

Supertype for abstract purses.  This type should be used when defining fallback methods of
functions for purses that unwrap the value and does not make use of the cache.
"""
abstract type MetaAbstractPurse{T} end

"""
    AbstractPurse{T,F<:Tuple}

Supertype for purses (or purse-like types) wrapping a value of type `T` and caching
functions of types `F`.  To create a subtype which conforms to the AbstractPurse interface,
[`value`](@ref) and [`cache`](@ref) must be implemented.
"""
abstract type AbstractPurse{T,F<:Tuple} <: MetaAbstractPurse{T} end

Base.convert(::Type{Any}, x::T) where {T<:AbstractPurse} = x
Base.convert(::Type{T}, x::T) where {T<:AbstractPurse} = x
Base.convert(::Type{T}, x::AbstractPurse) where {T} = convert(T, value(x))
Base.convert(::Type{T}, x) where {T<:AbstractPurse} = T(x)
Base.convert(::Type{T}, x::AbstractPurse) where {T<:AbstractPurse} = T(value(x))

cache_signature(::Type{<:AbstractPurse{<:Any,F}}) where {F} = F

"""
    Purse{T,F<:Tuple,S<:Tuple}

Purse with a value and cache.

# Fields
- `value::T`: the value stored in the purse
- `cache::S`: the cached values stored as a tuple in the same order as `F`
"""
struct Purse{T,F<:Tuple,S<:Tuple} <: AbstractPurse{T,F}
    value::T
    cache::S
end

"""
    Purse(value, fs...)

Return `value` wrapped in a `Purse`.  Each `f` in `fs` can be either a callable or a pair a
callable and a value.  If `f` is a callable, the cache will automatically be calculated by
applying each `f` to `value`.  If `f` is a pair, the cache for the first element of the pair
is set to the value of the second element.

# Examples
```jldoctest
julia> Purse(1.0)
Purse{Float64,Tuple{},Tuple{}}(1.0, ())

julia> Purse(2, inv)
Purse{Int64,Tuple{typeof(inv)},Tuple{Float64}}(2, (0.5,))

julia> Purse(2, inv => 0.2)
Purse{Int64,Tuple{typeof(inv)},Tuple{Float64}}(2, (0.2,))
```
"""
function Purse(value, fs...)
    Fs = map(f -> _register!(f isa Pair ? first(f) : f), fs)
    cache = map(f -> f isa Pair ? last(f) : f(value), fs)
    return Purse(value, Fs, cache)
end

Purse{T,F}(value) where {T,F<:Tuple} = Purse{T}(value, map(f -> _register!(f.instance), F.parameters)...)
Purse{T}(value, fs...) where {T} = Purse(convert(T, value), fs...)
Purse(value::T, ::F, cache::S) where {T,F<:Tuple,S<:Tuple} = Purse{T,F,S}(value, cache)

function Purse{T,F,S}(value) where {T,F<:Tuple,S<:Tuple}
    self = Purse{T}(value, map(f -> _register!(f.instance), F.parameters)...)
    return Purse{T,F,S}(self.value, convert(S, self.cache))
end

"""
    value(x)

If `x` is a purse, return its stored value, otherwise return x itself.

# Examples
```jldoctest
julia> Purses.value(Purse(2.0))
2.0
```
"""
@inline value(x) = x
@inline value(x::Purse) = x.value

"""
    cache(x[, idx])

Return the cache associated with the purse `x`.  If the optional argument `idx` is supplied,
return the cache item stored at index `idx`.

# Examples
```jldoctest
julia> purse = Purse(1.0, -);

julia> Purses.cache(purse, 1)
-1.0
```
"""
@inline cache(x::Purse) = x.cache
@inline cache(x, idx) = cache(x)[idx]


text:
#__precompile__()

module MuJoCo

using Reexport
using Libdl

const depsfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(depsfile)
   include(depsfile)
else
   error("MuJoCo was not installed correctly.")
end

include("./mj.jl")
@reexport using .mj

include("./export_all.jl") # a list of all function and struct names

function teardown()
   ACTIVATED[] = false
   mj_deactivate()
end

const ACTIVATED = Base.RefValue(false)
isactivated() = ACTIVATED[]
function __init__()
   if !isactivated()
      if Sys.islinux()
         mjglew = Libdl.dlopen_e(libglew, Libdl.RTLD_LAZY | Libdl.RTLD_DEEPBIND | Libdl.RTLD_GLOBAL)
         mjGL   = Libdl.dlopen_e(libGL, Libdl.RTLD_LAZY | Libdl.RTLD_DEEPBIND | Libdl.RTLD_GLOBAL)
      end

      key = ""
      try
         key = ENV["MUJOCO_KEY_PATH"]

         mjlib = Libdl.dlopen_e(libmujoco, Libdl.RTLD_LAZY | Libdl.RTLD_DEEPBIND | Libdl.RTLD_GLOBAL)
         if mjlib == C_NULL
            @error("Please build MuJoCo.jl")
         end
         ACTIVATED[] = Bool( ccall(Libdl.dlsym(mjlib, :mj_activate),Cint,(Cstring,),key) )
         if !isactivated()
            println("MuJoCo not activated; call mj_activate with path to license key.")
         end
      catch e
         println("Set MUJOCO_KEY_PATH environment variable, please.")
      end

      atexit(teardown)
   end
end

end


text:
function optimize!(mod::LLVM.Module)
    # everying except unroll, slpvec, loop-vec
    # then finish Julia GC
    ModulePassManager() do pm
        add_library_info!(pm, triple(mod))
        add_transform_info!(pm, tm[])

        propagate_julia_addrsp!(pm)
        scoped_no_alias_aa!(pm)
        type_based_alias_analysis!(pm)
        basic_alias_analysis!(pm)
        cfgsimplification!(pm)
        # TODO: DCE (doesn't exist in llvm-c)
        scalar_repl_aggregates!(pm) # SSA variant?
        mem_cpy_opt!(pm)
        always_inliner!(pm)
        alloc_opt!(pm)
        instruction_combining!(pm)
        cfgsimplification!(pm)
        scalar_repl_aggregates!(pm) # SSA variant?
        instruction_combining!(pm)
        jump_threading!(pm)
        instruction_combining!(pm)
        reassociate!(pm)
        early_cse!(pm)
        alloc_opt!(pm)
        loop_idiom!(pm)
        loop_rotate!(pm)
        lower_simdloop!(pm)
        licm!(pm)
        loop_unswitch!(pm)
        instruction_combining!(pm)
        ind_var_simplify!(pm)
        loop_deletion!(pm)
        # SimpleLoopUnroll -- not for Enzyme
        alloc_opt!(pm)
        scalar_repl_aggregates!(pm) # SSA variant?
        instruction_combining!(pm)
        gvn!(pm)
        mem_cpy_opt!(pm)
        sccp!(pm)
        # TODO: Sinking Pass
        # TODO: LLVM <7 InstructionSimplifier
        instruction_combining!(pm)
        jump_threading!(pm)
        dead_store_elimination!(pm)
        alloc_opt!(pm)
        cfgsimplification!(pm)
        loop_idiom!(pm)
        loop_deletion!(pm)
        jump_threading!(pm)
        # SLP_Vectorizer -- not for Enzyme
        aggressive_dce!(pm)
        instruction_combining!(pm)
        # Loop Vectorize -- not for Enzyme
        # InstCombine

        # GC passes
        barrier_noop!(pm)
        lower_exc_handlers!(pm)
        gc_invariant_verifier!(pm, false)
        late_lower_gc_frame!(pm)
        final_lower_gc!(pm)
        # TODO: DCE doesn't exist in llvm-c
        lower_ptls!(pm, #=dump_native=# false)

        # FIXME: Currently crashes printing
        # remove_julia_addrspaces!(pm)
        cfgsimplification!(pm)
        instruction_combining!(pm) # Extra for Enzyme

        run!(pm, mod)
    end
end

function post_optimze!(mod)
    # run second set of optimizations post enzyme
    ModulePassManager() do pm
        add_library_info!(pm, triple(mod))
        add_transform_info!(pm, tm[])

        scoped_no_alias_aa!(pm)
        type_based_alias_analysis!(pm)
        basic_alias_analysis!(pm)
        cfgsimplification!(pm)
        # TODO: DCE (doesn't exist in llvm-c)
        scalar_repl_aggregates!(pm) # SSA variant?
        mem_cpy_opt!(pm)
        always_inliner!(pm)
        instruction_combining!(pm)
        cfgsimplification!(pm)
        scalar_repl_aggregates!(pm) # SSA variant?
        instruction_combining!(pm)
        jump_threading!(pm)
        instruction_combining!(pm)
        reassociate!(pm)
        early_cse!(pm)
        loop_idiom!(pm)
        loop_rotate!(pm)
        lower_simdloop!(pm)
        licm!(pm)
        loop_unswitch!(pm)
        instruction_combining!(pm)
        ind_var_simplify!(pm)
        loop_deletion!(pm)
        # SimpleLoopUnroll -- not for Enzyme
        scalar_repl_aggregates!(pm) # SSA variant?
        instruction_combining!(pm)
        gvn!(pm)
        mem_cpy_opt!(pm)
        sccp!(pm)
        # TODO: Sinking Pass
        # TODO: LLVM <7 InstructionSimplifier
        instruction_combining!(pm)
        jump_threading!(pm)
        dead_store_elimination!(pm)
        cfgsimplification!(pm)
        loop_idiom!(pm)
        loop_deletion!(pm)
        jump_threading!(pm)
        # SLP_Vectorizer -- not for Enzyme
        aggressive_dce!(pm)
        instruction_combining!(pm)
        # Loop Vectorize -- not for Enzyme
        # InstCombine

        cfgsimplification!(pm)
        instruction_combining!(pm)
        # CombineMulAddPass will run on second pass

        run!(pm, mod)
    end
end


text:
# quick reference to commonly used modules
const Proto = typedb.protocol

# Convenient types

const Optional{T} = Union{T,Nothing}
const Transaction_Res_All = Union{Proto.Transaction_Res, Proto.Transaction_ResPart}

# protobuf types
const EnumType = Int32
const Bytes = Array{UInt8,1}


text:
################################################################################
#
#  Factor base over (Euclidean) Rings
#
################################################################################

# Note that T must admit gcd's and base must consist of elements x for which
# valuation(_, x) is definied.
# works (at least) for fmpz and nmod_poly, so it can be used for the
# smoothness test

function _compose(a::node{T}, b::node{T}, check = false) where T
  if check && !isone(gcd(a.content, b.content))
    error("input not coprime")
  end
  return node{T}(a.content * b.content, a, b)
end

# assume that the set or array consists of pairwise coprime elements
function FactorBase(x::Union{Set{T}, AbstractVector{T}}; check::Bool = true) where T
  if length(x)==0
    z = FactorBase{T}(T(1), x)
    return z
  end
  ax = [ node{T}(p) for p in x]
  while length(ax) > 1
    bx = [ _compose(ax[2*i-1], ax[2*i], check) for i=1:div(length(ax), 2)]
    if isodd(length(ax))
      push!(bx, ax[end])
    end
    ax = bx
  end
  z = FactorBase{T}(ax[1].content, x)
  z.ptree = ax[1]
  return z
end

function show(io::IO, B::FactorBase{T}) where T
  print(io, "Factor base with \n$(B.base) and type $T")
end

function ring(B::FactorBase)
  return parent(B.prod)
end

function issmooth(c::FactorBase{T}, a::T) where T
  @assert a != 0
  g = gcd(c.prod, a)
  while g != 1
    a = div(a, g)
    g = gcd(g, a)
  end
  return a == 1 || a==-1
end

function issmooth!(c::FactorBase{fmpz}, a::fmpz)
  @assert a != 0
  g = gcd(c.prod, a)
  if g==1
    return a==1 || a==-1, a
  end
  b = copy(a)
  while g != 1
    divexact!(b, b, g)
    gcd!(g, g, b)
  end
  return b == 1 || b==-1, b
end

function isleaf(a::node{T}) where T
  return !(isdefined(a, :left) || isdefined(a, :right))
end

function _split(c::node{T}, a::T) where T
  if isleaf(c)
    return [gcd(a, c.content)]
  end
  if isdefined(c, :left)
    l = gcd(a, c.left.content)
    if l != 1
      ls = _split(c.left, l)
    else
      ls = Vector{T}()
    end
  else
    ls = Vector{T}()
  end
  if isdefined(c, :right)
    r = gcd(a, c.right.content)
    if r != 1
      rs = _split(c.right, r)
    else
      rs = Vector{T}()
    end
  else
    rs = Vector{T}()
  end
  return vcat(ls, rs)
end

function factor(c::FactorBase{T}, a::T, do_error::Bool = true) where T
  @assert a != 0
  f = Dict{T, Int}()
  isunit(a) && return f
  lp = _split(c.ptree, a)
  for i in lp
    if mod(a, i)==0  ## combine: use divmod and do val of rest
                     ## to save a division
      v = remove(a, i)
      f[i] = v[1]
      a = v[2]
      if isunit(a)
        break
      end
    end
  end
  @assert (!do_error || a==1 || a==-1)
  return f
end

function factor(c::FactorBase{fmpz}, a::fmpq)
  @assert a != 0
  a = deepcopy(a)
  f = Dict{fmpz, Int}()
  n = abs(numerator(a))
  d = denominator(a)
  lp = _split(c.ptree, n*d)
  for i in lp
    v = remove!(a, i)
    if isdefined(f, :i)
      f[i] += v
    else
      f[i] = v
    end
    if isone(a)
      break
    end
  end
  @hassert :ClassGroup 1 isone(a)
  return f
end


text:
@safetestset "isafter" begin
    using Validator

    # date
    @test Validator.isafter("2022-02-28", "2022-02-27") === true
    @test Validator.isafter("2022-02-28", "2022-02-28") === false
    @test Validator.isafter("2022-02-28", "2022-03-01") === false

    # datetime
    @test Validator.isafter("2022-02-28T12:00:00", "2022-02-27T12:00:00") === true
    @test Validator.isafter("2022-02-28T12:00:00", "2022-02-28T09:00:00") === true
    @test Validator.isafter("2022-02-28T12:00:00", "2022-02-28T12:00:00") === false
    @test Validator.isafter("2022-02-28T12:00:00", "2022-03-01T12:00:00") === false

    # default
    @test Validator.isafter("2022-02-27") === false
    @test Validator.isafter("2050-02-28") === true
    @test Validator.isafter("2022-02-28T12:00:00") === false
    @test Validator.isafter("2050-02-28T09:00:00") === true
end


text:
module SpringCollab2020CancelLightningRemoveRoutineTrigger

using ..Ahorn, Maple

@mapdef Trigger "SpringCollab2020/CancelLightningRemoveRoutineTrigger" CancelLightningRemoveRoutineTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight)

const placements = Ahorn.PlacementDict(
    "Cancel Lightning Remove Routine Trigger (Spring Collab 2020)" => Ahorn.EntityPlacement(
        CancelLightningRemoveRoutineTrigger,
        "rectangle",
    ),
)

end


text:
#Zero Momentum

#without a maxiter

function dkF1Tilde3(k, m, T)
    -1 / 48 * (
        k^4 * (
            6 * T * coth(Epi(k, m) / (2 * T)) +
            csch(Epi(k, m) / (2 * T))^2 * Epi(k, m)
        )
    ) / (pi^2 * T * Epi(k, m)^5)
end

function dkF2Tilde3(k, m, T)
    -1 / 48 * (
        k^4 *
        csch(Epi(k, m) / (2 * T))^2 *
        (2 * T + coth(Epi(k, m) / (2 * T)) * Epi(k, m))
    ) / (pi^2 * T^2 * Epi(k, m)^4)
end

function PvdkF1Tildeps(p0, qsmax, k, m, T; kwargs...)
    if k > qsmax / 2
        return quadgk_PV(
            x -> ppfunps(x, qsmax, k, m, T),
            2 * Epi(k, m),
            p0,
            Epi(k + qsmax, m) + Epi(k, m);
            kwargs...,
        )
    elseif k <= qsmax / 2
        return quadgk_PV(
            x -> ppfunps(x, qsmax, k, m, T),
            Epi(k, m) + Epi(qsmax - k, m),
            p0,
            Epi(k + qsmax, m) + Epi(k, m);
            kwargs...,
        )
    end
end

function dkF1Tildeps_delta1(p0, qsmax, k, m, T)
    if k > qsmax / 2
        deltasumps(p0, qsmax, k, m, T)
    elseif k <= qsmax / 2
        0.0
    end
end

function dkF1Tildeps_delta2(p0, qsmax, k, m, T)
    if k > qsmax / 2
        ppfunps_delta(p0, qsmax, k, m, T)
    elseif k <= qsmax / 2
        0.0
    end
end

function PvdkF2Tildeps(p0, qsmax, k, m, T; kwargs...)
    return quadgk_PV(
        x -> pmfunps(x, qsmax, k, m, T),
        0.0,
        p0,
        Epi(k + qsmax, m) - Epi(k, m);
        kwargs...,
    )[1] + pmfunps_zerofix(p0, qsmax, k, m, T, 9)
end

#dkF1Tilde_intcostheqs has a delta function contribution in ps < 2 * k + qsmax brunch
#but we have seprated it

@doc raw"""
    PvdkF1Tildeps(p0, ps, qsmax, k, m, T; kwargs...)

Piecewise hints:

We only integrate the `p0'` domain which the  `flowpp_intcostheqs` is non zero.
"""
function PvdkF1Tildeps(p0, ps, qsmax, k, m, T; kwargs...)
    if ps >= 2 * k + qsmax
        return quadgk_PV2(
            x -> flowpp_intcostheqs(x, ps, qsmax, k, m, T),
            Epi(ps - k - qsmax, m) + Epi(k, m),
            p0,
            Epi(k + ps + qsmax, m) + Epi(k, m);
            kwargs...,
        )
    elseif ps < 2 * k + qsmax
        return quadgk_PV3(
            x -> flowpp_intcostheqs(x, ps, qsmax, k, m, T),
            2 * Epi(k, m),
            p0,
            Epi(k + ps + qsmax, m) + Epi(k, m);
            kwargs...,
        )
    end
end

function PvdkF1Tildeps_Compensate(p0, ps, qsmax, k, m, T)
    if ps >= 2 * k + qsmax
        return 0.0
    elseif ps < 2 * k + qsmax
        if Epi(k + ps + qsmax, m) + Epi(k, m) > p0 > 2 * Epi(k, m)
            return flowpp_intcostheqs(p0, ps, qsmax, k, m, T) * log(
                abs(
                    ((Epi(k + ps + qsmax, m) + Epi(k, m))^2 - p0^2) /
                    ((2 * Epi(k, m))^2 - p0^2),
                ),
            )
        else
            return flowpp_intcostheqs(nextfloat(2 * Epi(k, m)), ps, qsmax, k, m, T) * log(
                abs(
                    ((Epi(k + ps + qsmax, m) + Epi(k, m))^2 - p0^2) /
                    ((2 * Epi(k, m))^2 - p0^2),
                ),
            )
        end
    end
end

function dkF1Tildeps_delta1(p0, ps, qsmax, k, m, T)
    if ps >= 2 * k + qsmax
        0.0
    elseif ps < 2 * k + qsmax
        delta1_intcosthqs(p0, ps, qsmax, k, m, T,0.5)
    end
end

function dkF1Tildeps_delta2(p0, ps, qsmax, k, m, T)
    if ps >= 2 * k + qsmax
        0.0
    elseif ps < 2 * k + qsmax
        delta2Smooth_intcosthqs(p0, ps, qsmax, k, m, T)
    end
end

function PvdkF2Tildeps(p0, ps, qsmax, k, m, T; kwargs...)
    if ps > 2 * k + qsmax
        return quadgk_PV2(
            x -> flowpm_intcostheqs(x, ps, qsmax, k, m, T),
            Epi(ps - k - qsmax, m) - Epi(k, m),
            p0,
            Epi(k + ps + qsmax, m) - Epi(k, m);
            kwargs...,
        )
    elseif ps < 2 * k + qsmax
        return quadgk_PV2(
            x -> flowpm_intcostheqs(x, ps, qsmax, k, m, T),
            0.0,
            p0,
            Epi(k + ps + qsmax, m) - Epi(k, m);
            kwargs...,
        )
    end
end

@doc raw"""
    dkF1TildeintqsAll(p0, qsmax, k, m, T; kwargs...)

`costh` is not integrated we need an extra `2` at somewhere
"""
function dkF1TildeintqsAll(p0, qsmax, k, m, T; kwargs...)
    return PvdkF1Tildeps(abs(p0), qsmax, k, m, T; kwargs...) +
           dkF1Tildeps_delta1(abs(p0), qsmax, k, m, T) +
           dkF1Tildeps_delta2(abs(p0), qsmax, k, m, T)
end
@doc raw"""
    dkF2TildeintqsAll(p0, qsmax, k, m, T; kwargs...)

`costh` is not integrated we need an extra `2` at somewhere
"""
function dkF2TildeintqsAll(p0, qsmax, k, m, T; kwargs...)
    PvdkF2Tildeps(abs(p0), qsmax, k, m, T; kwargs...)
end

function dkF1TildeintqsAll(p0, ps, qsmax, k, m, T; kwargs...)
    return PvdkF1Tildeps(abs(p0), ps, qsmax, k, m, T; kwargs...)
    # +dkF1Tildeps_delta1(abs(p0), ps, qsmax, k, m, T)+dkF1Tildeps_delta2(abs(p0), ps, qsmax, k, m, T)
    # +PvdkF1Tildeps(abs(p0), ps, qsmax, k, m, T; kwargs...)
    # +
end

function dkF1TildeintqsAll_delta1(p0, ps, qsmax, k, m, T)
    return dkF1Tildeps_delta1(abs(p0), ps, qsmax, k, m, T)
end
function dkF1TildeintqsAll_delta2(p0, ps, qsmax, k, m, T)
    return dkF1Tildeps_delta2(abs(p0), ps, qsmax, k, m, T)
end
function dkF1TildeintqsAll_Compensate(p0, ps, qsmax, k, m, T)
    return PvdkF1Tildeps_Compensate(abs(p0), ps, qsmax, k, m, T)
end

function dkF2TildeintqsAll(p0, ps, qsmax, k, m, T; kwargs...)
    PvdkF2Tildeps(abs(p0), ps, qsmax, k, m, T; kwargs...)
end

function dkF1TildeAll(k, m, T)
    dkF1Tilde3(k, m, T)
end
function dkF2TildeAll(k, m, T)
    dkF2Tilde3(k, m, T)
end


text:
install_term_repr()

@testset "REPR renderable repr" begin
    p = string(Panel())
    @test sprint(show, Panel()) ==
          "\e[38;5;117mPanel <: AbstractRenderable\e[0m \e[2m(w:88, h:2)\e[0m"
end

@testset "REPR @with_repr" begin
    # just check that creating a @with_repr is ok
    @with_repr struct Rocket
        width::Int
        height::Int
        mass::Float64

        manufacturer::String
    end

    obj = Rocket(10, 50, 5000, "NASA")

    io = IOBuffer()
    show(IOContext(io), "text/plain", obj)
    s = String(take!(io))

    correct_s = "\e[38;2;155;179;224m╭──────────────────────────────────────╮\e[39m\n\e[0m\e[38;2;155;179;224m│\e[39m\e[0m             \e[1m\e[38;2;224;219;121mwidth\e[22m\e[39m\e[38;2;187;134;219m::Int64\e[39m\e[2m\e[38;2;126;157;217m│\e[22m\e[39m\e[0m \e[38;2;179;212;255m10\e[39m         \e[0m\e[38;2;155;179;224m│\e[39m\e[0m\n\e[0m\e[38;2;155;179;224m│\e[39m\e[0m            \e[1m\e[38;2;224;219;121mheight\e[22m\e[39m\e[38;2;187;134;219m::Int64\e[39m\e[2m\e[38;2;126;157;217m│\e[22m\e[39m\e[0m \e[38;2;179;212;255m50\e[39m         \e[0m\e[38;2;155;179;224m│\e[39m\e[0m\n\e[0m\e[38;2;155;179;224m│\e[39m\e[0m            \e[1m\e[38;2;224;219;121mmass\e[22m\e[39m\e[38;2;187;134;219m::Float64\e[39m\e[2m\e[38;2;126;157;217m│\e[22m\e[39m\e[0m \e[38;2;179;212;255m5000.0\e[39m     \e[0m\e[38;2;155;179;224m│\e[39m\e[0m\n\e[0m\e[38;2;155;179;224m│\e[39m\e[0m     \e[1m\e[38;2;224;219;121mmanufacturer\e[22m\e[39m\e[38;2;187;134;219m::String\e[39m\e[2m\e[38;2;126;157;217m│\e[22m\e[39m\e[0m \e[38;2;179;212;255mNASA\e[39m       \e[0m\e[38;2;155;179;224m│\e[39m\e[0m\n\e[38;2;155;179;224m╰─────────────────────────── \e[38;2;227;172;141mRocket\e[39m\e[38;2;155;179;224m\e[38;2;155;179;224m ───╯\e[39m\e[0m\e[39m\e[38;2;155;179;224m\e[0m\n"
    @test s == correct_s
    @test sprint(termshow, obj) == correct_s

    termshow(devnull, Rocket)  # coverage

    @with_repr struct T end

    @test sprint(termshow, T()) ==
          "\e[38;2;155;179;224m╭──────────────────────────────────────╮\e[39m\n\e[0m\e[38;2;155;179;224m│\e[39m\e[0m                T()\e[38;2;187;134;219m::T\e[39m                \e[0m\e[38;2;155;179;224m│\e[39m\e[0m\n\e[38;2;155;179;224m╰──────────────────────────────── \e[38;2;227;172;141mT\e[39m\e[38;2;155;179;224m\e[38;2;155;179;224m ───╯\e[39m\e[0m\e[39m\e[38;2;155;179;224m\e[0m\n"
end

objs = if VERSION >= v"1.7.1"
    (
        (1, [1, 2, 3]),
        (2, Dict(:x => [1, 2, 3], "a" => Dict(:z => "a"))),
        (3, Dict(i => i for i in 1:100)),
        (4, zeros(120, 300)),
        (5, zeros(200)),
        (6, zeros(3, 3, 3)),
        (7, termshow),
        (8, :(x / y + √9)),
    )
else
    (
        (1, [1, 2, 3]),
        (2, Dict(:x => [1, 2, 3], "a" => Dict(:z => "a"))),
        (3, Dict(i => i for i in 1:100)),
        (7, termshow),
    )
end

for (i, t) in objs
    t = sprint(termshow, t)
    tofile(string(t), "./txtfiles/termshow_$i.txt")
end

# @testset "TERMSHOW for types" begin
#     for (i, t) in objs
#         t = sprint(termshow, t)
#         @test fromfile("./txtfiles/termshow_$i.txt") == t
#     end
# end

@testset "Term automatic repr" begin
    repr_show(io, x) = show(io, MIME("text/plain"), x)
    @test sprint(repr_show, 1) == "\e[38;2;144;202;249m1\e[39m"

    @test sprint(repr_show, Dict(1 => :x)) ==
          "\e[38;2;155;179;224m╭──── \e[38;2;227;172;141mDict {Int64, Symbol} \e[39m\e[38;2;155;179;224m\e[38;2;155;179;224m ───────────╮\e[39m\e[0m\e[39m\e[38;2;155;179;224m\n\e[0m\e[38;2;155;179;224m│\e[39m\e[0m                                      \e[0m\e[38;2;155;179;224m│\e[39m\e[0m\n\e[0m\e[38;2;155;179;224m│\e[39m\e[0m    \e[2m\e[38;2;187;134;219m {Int64} \e[22m\e[39m\e[2m\e[38;2;126;157;217m│\e[22m\e[39m\e[0m \e[1m\e[38;2;224;219;121m1\e[22m\e[39m \e[1m\e[31m=>\e[22m\e[39m \e[1m\e[38;2;179;212;255mx\e[22m\e[39m \e[2m\e[38;2;126;157;217m│\e[22m\e[39m\e[0m\e[2m\e[38;2;187;134;219m {Symbol} \e[22m\e[39m     \e[0m\e[38;2;155;179;224m│\e[39m\e[0m\n\e[0m\e[38;2;155;179;224m│\e[39m\e[0m                                      \e[0m\e[38;2;155;179;224m│\e[39m\e[0m\n\e[38;2;155;179;224m╰────────────────────────── \e[1m\e[37m1\e[22m\e[39m\e[22m items\e[22m\e[39m\e[38;2;155;179;224m ───╯\e[39m\e[0m\e[39m\e[0m\n"

    @test sprint(repr_show, :(x + y)) ==
          "\e[38;2;155;179;224m╭──────────────────────────────────────╮\e[39m\n\e[0m\e[38;2;155;179;224m│\e[39m\e[0m                \e[32mx \e[38;2;239;83;80m+\e[39m y\e[39m                 \e[0m\e[38;2;155;179;224m│\e[39m\e[0m\n\e[0m\e[38;2;155;179;224m│\e[39m\e[0m    \e[22m─────────────────────────────\e[0m\e[22m     \e[0m\e[38;2;155;179;224m│\e[39m\e[0m\n\e[0m\e[38;2;155;179;224m│\e[39m\e[0m    \e[1m\e[38;2;224;219;121mhead\e[22m\e[39m\e[38;2;187;134;219m::Symbol\e[39m\e[2m\e[38;2;126;157;217m│\e[22m\e[39m\e[0m \e[38;2;179;212;255mcall\e[39m                \e[0m\e[38;2;155;179;224m│\e[39m\e[0m\n\e[0m\e[38;2;155;179;224m│\e[39m\e[0m    \e[1m\e[38;2;224;219;121margs\e[22m\e[39m\e[38;2;187;134;219m::Vector\e[39m\e[38;2;187;134;219m\e[2m\e[38;2;126;157;217m│\e[22m\e[39m\e[0m \e[38;2;179;212;255mAny[:+, :x, :y]\e[39m     \e[0m\e[38;2;155;179;224m│\e[39m\e[0m\n\e[38;2;155;179;224m╰───────────────────────────── \e[38;2;227;172;141mExpr\e[39m\e[38;2;155;179;224m\e[38;2;155;179;224m ───╯\e[39m\e[0m\e[39m\e[38;2;155;179;224m\e[0m\n"
end


text:
@testset "Testing lae.jl" begin

    # copied from the lae lecture
    s = 0.2
    δ = 0.1
    a_σ = 0.4       # A = exp(B) where B ~ N(0, a_sigma)
    α = 0.4         # We set f(k) = k**alpha
    ϕ = LogNormal(a_σ)

    function p(x, y)
        d = s * x.^α

        pdf_arg = clamp((y .- (1-δ) .* x) ./ d, eps(), Inf)
        return pdf(ϕ, pdf_arg) ./ d
    end

    # other data
    n_a, n_b, n_y = 50, (5, 5), 20
    a = rand(n_a) + 0.01
    b = rand(n_b...) + 0.01

    y = linspace(0, 10, 20)

    lae_a = LAE(p, a)
    lae_b = LAE(p, b)

    laes = [lae_a, lae_b]

    # test stuff
    for l in laes
        # make sure X is made a column vector by constructor
        @test size(l.X, 2) == 1

        # make sure X is indeed a column vector (2 dims)
        @test ndims(l.X) == 2

        # Make sure we get 1d estimate out
        @test size(lae_est(l, y)) == (n_y, )
    end
end  # @testset


text:
include("learn_Plotly.jl")

app = dash(external_stylesheets=[
    dbc_themes.SPACELAB,
    "https://cdnjs.cloudflare.com/ajax/libs/font-awesome/4.7.0/css/font-awesome.min.css",
    #"/myassets/spacelab.css",
    #"myassets/my.css"
    ])

# Uncomment one of the following to run
#app = chapter_configuration_options(app=app)
#app = chapter_scatter_plots(app=app)
#app = chapter_line_charts(app=app)
#app = chapter_bar_charts(app=app)
#app = chapter_pie_charts(app=app)
#app = chapter_bubble_charts(app=app)
app = chapter_histograms(app=app)
#app = chapter_contour_plots(app=app)
#app = chapter_heatmaps(app=app)
#app = chapter_3d_scatter_plots(app=app)
#app = chapter_3d_surface_plots(app=app)

try # Enable Ctrl-C
    run_server(
        app,
        "0.0.0.0",
        8050,
        debug=true, # enables hot reload and more
    )
catch e
    @show e
end


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 6df3beed-24a7-4b26-a315-0520f4863190
develop=false

# ╔═╡ 9701cbe0-d048-11eb-151b-67dda7b72b71
begin
   using Pkg
   Pkg.activate(mktempdir())
   Pkg.add("Revise")
   using Revise
   Pkg.add(["PyPlot","ExtendableGrids","PlutoUI","Plots"])
	if develop
	   Pkg.develop(["PlutoVista","GridVisualize"])
	else
	   Pkg.add(["PlutoVista","GridVisualize"])
	end
   using PyPlot,PlutoVista,GridVisualize,ExtendableGrids,PlutoUI,Plots
end

# ╔═╡ 68e2c958-b417-4ba1-9577-697304fe140a
TableOfContents()

# ╔═╡ b35d982d-1fa9-413d-b008-892b4f241097
md"""
# Test notebook for PlutoVista backend of GridVisualize
"""

# ╔═╡ 00b04f6b-34a6-4e30-9864-d273305281d4
md"""
## 1D scalar plot
"""

# ╔═╡ a20d74c9-16da-408a-b247-0c17321888f9
function testplot1(Plotter)
	grid=simplexgrid(0:0.01:10)
	scalarplot(grid,map(sin,grid),Plotter=Plotter,resolution=(600,200),markershape=:star5,markevery=20, xlabel="x",ylabel="z",legend=:rt,label="sin")
end

# ╔═╡ cc17187f-404c-4c31-8625-fa067eea7273
testplot1(PyPlot)

# ╔═╡ ad1c4dd6-7c6b-4433-a1ac-b2f817ba5d81
testplot1(Plots)

# ╔═╡ 33482af8-3542-4723-ae43-770a789b69b3
testplot1(PlutoVista)

# ╔═╡ c4eeb06f-932e-4acc-8e5b-f2a7f9242a42
function testplot2(Plotter;t=0)
	p=GridVisualizer(Plotter=Plotter, resolution=(500,200),legend=:rt,xlabel="x")
	grid=simplexgrid(0:0.01:10)
	scalarplot!(p,grid,map(x->sin(x-t),grid),Plotter=Plotter,color=:red,label="sin(x-$(t))",linestyle=:dash)
	scalarplot!(p,grid,map(cos,grid),Plotter=Plotter,color=:green,clear=false,label="cos",linestyle=:dashdot,linewidth=3)
	reveal(p)
end

# ╔═╡ 29ca4775-6ba5-474c-bd2c-8f770b09addd
testplot2(PyPlot)

# ╔═╡ c0ab77e8-01ea-436d-85f2-34e253944f11
testplot2(Plots)

# ╔═╡ 84192945-d4b6-4949-8f06-d94e04a7a56d
testplot2(PlutoVista)

# ╔═╡ 63fe3259-7d79-40ec-98be-e0592e40ee6b
@bind t2 PlutoUI.Slider(0:0.1:5,show_value=true)

# ╔═╡ 4de6b5c9-4d2d-4bcb-bc88-c6f50a23f9b6
testplot2(PlutoVista,t=t2)

# ╔═╡ 2061e7fd-c740-4d4b-af5b-7a3a9444aafd
md"""
### Changeable data (experimental)

This just updates the data. In the 1D case the difference seems to be not critical.
"""

# ╔═╡ f84beb4f-4136-4e5a-ba43-279b703fc75f
begin
	X2=0:0.001:10
	grid2=simplexgrid(collect(X2))
	f2(t)=map( (x)->sin(x^2-t),grid2)
end

# ╔═╡ c1278fb2-3e75-445f-893a-b8b8a7e931d3
p=PlutoVista.PlutoVistaPlot(resolution=(600,200))

# ╔═╡ 29fa4467-65ee-4dad-a660-5197864ddbdc
md"""
t4: $(@bind t4 PlutoUI.Slider(-10:0.1:10, default=0, show_value=true))
"""

# ╔═╡ 661531f7-f740-4dd4-9a59-89ddff06ba5c
PlutoVista.plot!(p,X2,f2(t4))

# ╔═╡ 9bb243cc-c69a-405b-bb35-6cddfde8fd30
begin
	myplot(t)=scalarplot!(vis2,grid2,f2(t),show=true)
	vis2=GridVisualizer(resolution=(600,200),Plotter=PlutoVista,datadim=1)
end

# ╔═╡ f4c78c61-19b1-4889-aa20-c6a3b157d435
md"""
t3: $(@bind t3 PlutoUI.Slider(-10:0.1:10, default=0, show_value=true))
"""

# ╔═╡ be369a01-a0c2-4a6b-831b-f716cc807240
myplot(t3)

# ╔═╡ ed9b80e5-9678-4ba6-bb36-c2e0674ed9ba
md"""
## 1D grid plot
"""

# ╔═╡ 9ce4f63d-cd96-48d7-a637-07cb84fa88ab
function testgridplot(Plotter)
	grid=simplexgrid(0:1:10)
	cellmask!(grid,[0.0],[5],2)
	bfacemask!(grid,[5],[5],3)
	gridplot(grid, Plotter=Plotter,resolution=(600,200),legend=:rt)
end

# ╔═╡ 77eeefc7-e416-426b-8f87-1bc8439dae6d
testgridplot(PyPlot)

# ╔═╡ d503ee1e-1e1f-4235-b286-dc3137a2c96a
testgridplot(PlutoVista)

# ╔═╡ ae1fe1ab-4a0e-4c80-bd6f-912201fb4bb4
md"""
## 2D Scalar plot
"""

# ╔═╡ d5258595-60e4-406f-a71e-69111cdad8b9
function testplot3(Plotter)
	X=0:0.1:10
	grid=simplexgrid(X,X)
	f=map( (x,y)->sin(x)*atan(y),grid)
	scalarplot(grid,f,Plotter=Plotter,resolution=(300,300),flimits=(-π/2,π/2),backend=:vtk)
end

# ╔═╡ c98a90bf-1a3e-4681-a3b0-663c6844df6b
testplot3(PyPlot)

# ╔═╡ a0c3067b-3aa5-493e-b132-89746483b5ce
testplot3(Plots)

# ╔═╡ 0998a9a7-d57a-476e-aacd-bee9396e9b8f
testplot3(PlutoVista)

# ╔═╡ cefb38c1-159e-42db-8088-294573fcece2
md"""
### Changeable data (experimental)

Here we observe a more profound advantage, and vtk.js also has a rather understandable way how data can be updated.

Generally, with plutovista we need two cells - one with the graph shown, and a second one which triggers the modification.
"""

# ╔═╡ a9f4f98f-ec2f-42d6-88da-4a8a6f727e93
begin
	X=0:0.1:10
	grid=simplexgrid(X,X)
	f(t)=map( (x,y)->sin(x-t)*atan(y)*cos((y-t)),grid)
end

# ╔═╡ 412c905f-050c-4b78-a66f-0d03978e7edf
begin
	vis=GridVisualizer(resolution=(300,300),Plotter=PlutoVista)
	myplot2(t)=scalarplot!(vis,grid,f(t),show=true,clear=true,flimits=(-π/2,π/2))
end

# ╔═╡ 6f1707ed-79ab-42dc-8ad8-d66a9e1a65b3
md"""
t= $(@bind t PlutoUI.Slider(-10:0.1:10, default=0, show_value=true))
"""

# ╔═╡ 461481ef-f88b-4e4e-b57d-ce003abbfdf1
myplot2(t)

# ╔═╡ Cell order:
# ╠═6df3beed-24a7-4b26-a315-0520f4863190
# ╠═9701cbe0-d048-11eb-151b-67dda7b72b71
# ╠═68e2c958-b417-4ba1-9577-697304fe140a
# ╟─b35d982d-1fa9-413d-b008-892b4f241097
# ╟─00b04f6b-34a6-4e30-9864-d273305281d4
# ╠═a20d74c9-16da-408a-b247-0c17321888f9
# ╠═cc17187f-404c-4c31-8625-fa067eea7273
# ╠═ad1c4dd6-7c6b-4433-a1ac-b2f817ba5d81
# ╠═33482af8-3542-4723-ae43-770a789b69b3
# ╠═c4eeb06f-932e-4acc-8e5b-f2a7f9242a42
# ╠═29ca4775-6ba5-474c-bd2c-8f770b09addd
# ╠═c0ab77e8-01ea-436d-85f2-34e253944f11
# ╠═84192945-d4b6-4949-8f06-d94e04a7a56d
# ╠═63fe3259-7d79-40ec-98be-e0592e40ee6b
# ╠═4de6b5c9-4d2d-4bcb-bc88-c6f50a23f9b6
# ╟─2061e7fd-c740-4d4b-af5b-7a3a9444aafd
# ╠═f84beb4f-4136-4e5a-ba43-279b703fc75f
# ╠═c1278fb2-3e75-445f-893a-b8b8a7e931d3
# ╠═661531f7-f740-4dd4-9a59-89ddff06ba5c
# ╟─29fa4467-65ee-4dad-a660-5197864ddbdc
# ╠═9bb243cc-c69a-405b-bb35-6cddfde8fd30
# ╠═be369a01-a0c2-4a6b-831b-f716cc807240
# ╟─f4c78c61-19b1-4889-aa20-c6a3b157d435
# ╠═ed9b80e5-9678-4ba6-bb36-c2e0674ed9ba
# ╠═9ce4f63d-cd96-48d7-a637-07cb84fa88ab
# ╠═77eeefc7-e416-426b-8f87-1bc8439dae6d
# ╠═d503ee1e-1e1f-4235-b286-dc3137a2c96a
# ╟─ae1fe1ab-4a0e-4c80-bd6f-912201fb4bb4
# ╠═d5258595-60e4-406f-a71e-69111cdad8b9
# ╠═c98a90bf-1a3e-4681-a3b0-663c6844df6b
# ╠═a0c3067b-3aa5-493e-b132-89746483b5ce
# ╠═0998a9a7-d57a-476e-aacd-bee9396e9b8f
# ╟─cefb38c1-159e-42db-8088-294573fcece2
# ╠═a9f4f98f-ec2f-42d6-88da-4a8a6f727e93
# ╠═412c905f-050c-4b78-a66f-0d03978e7edf
# ╠═461481ef-f88b-4e4e-b57d-ce003abbfdf1
# ╟─6f1707ed-79ab-42dc-8ad8-d66a9e1a65b3