text:
using SnpArrays
using Knockoffs
using DelimitedFiles
using Random
using LinearAlgebra
using Distributions
using ProgressMeter
using MendelIHT
using VCFTools
using StatsBase
using CSV
using DataFrames
using Printf
using GLMNet
BLAS.set_num_threads(1)

global qctools_exe = "/scratch/users/bbchu/qctool/build/release/qctool_v2.0.7"
global snpknock2_exe = "/scratch/users/bbchu/knockoffgwas/snpknock2/bin/snpknock2"
global rapid_exe = "/scratch/users/bbchu/RaPID/RaPID_v.1.7"
global partition_exe = "/scratch/users/bbchu/knockoffgwas/knockoffgwas/utils/partition.R";

"""
    simulate_pop_structure(n, p)

Simulate genotypes with K = 2 populations. 1% of SNPs will have different allele
frequencies between the populations.

# Inputs
- `plinkfile`: Output plink file name.
- `n`: Number of samples
- `p`: Number of SNPs

# Output
- `x1`: n×p matrix of the 1st haplotype for each sample. Each row is a haplotype
- `x2`: n×p matrix of the 2nd haplotype for each sample. `x = x1 + x2`
- `populations`: Vector of length `n` indicating population membership for eachsample.
- `diff_markers`: Indices of the differentially expressed alleles.

# Reference
https://www.nature.com/articles/nrg2813
"""
function simulate_pop_structure(n::Int, p::Int)
    # first simulate genotypes treating all samples equally
    x1 = BitMatrix(undef, n, p)
    x2 = BitMatrix(undef, n, p)
    pmeter = Progress(p, 0.1, "Simulating genotypes...")
    @inbounds for j in 1:p
        d = Bernoulli(rand())
        for i in 1:n
            x1[i, j] = rand(d)
            x2[i, j] = rand(d)
        end
        next!(pmeter)
    end
    # assign populations and simulate 0.01p unually differentiated markers
    populations = rand(1:2, n)
    diff_markers = sample(1:p, Int(0.01p), replace=false)
    @inbounds for j in diff_markers
        pop1_allele_freq = 0.4rand()
        pop2_allele_freq = pop1_allele_freq + 0.6
        pop1_dist = Bernoulli(pop1_allele_freq)
        pop2_dist = Bernoulli(pop2_allele_freq)
        for i in 1:n
            d = isone(populations[i]) ? pop1_dist : pop2_dist
            x1[i, j] = rand(d)
            x2[i, j] = rand(d)
        end
    end
    return x1, x2, populations, diff_markers
end

"""
    simulate_IBD(h1::BitMatrix, h2::BitMatrix, populations::Vector{Int}, k::Int)

Simulate recombination events. Parent haplotypes `h1` and `h2` will be used to generate
`k` children, then both parent and children haplotypes will be returned.

In offspring simulation, we first randomly sample 2 parents from the same population.
Then generate offspring individuals by copying segments of the parents haplotype
directly to the offspring to represent IBD segments. The number of segments (i.e. places of
recombination) is 1 to 5 per sample chosen uniformly across all SNPs.

# Inputs
- `h1`: `n × p` matrix of the 1st haplotype for each parent. Each row is a haplotype
- `h2`: `n × p` matrix of the 2nd haplotype for each parent. `H = h1 + h2`
- `populations`: `populations[i]` is the population (represented as integer) of sample `i`.
- `k`: Total number of offsprings

# Output
- `H1`: `n+k × p` matrix of the 1st haplotype. The first `n` haplotypes are from parents
    and the next `k` haplotypes are the offsprings. Each row is a haplotype
- `H2`: `n+k × p` matrix of the 2nd haplotype. `x = x1 + x2`

# References
https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1003520
"""
function simulate_IBD(h1::BitMatrix, h2::BitMatrix, populations::Vector{Int}, k::Int)
    n, p = size(h1)
    unique_populations = unique(populations)
    # randomly designate gender for parents
    sex = bitrand(n)
    male_idx = findall(x -> x == true, sex)
    female_idx = findall(x -> x == false, sex)
    # simulate new samples
    x1 = falses(k, p)
    x2 = falses(k, p)
    fathers = Int[]
    mothers = Int[]
    pmeter = Progress(k, 0.1, "Simulating IBD segments...")
    for i in 1:k
        # assign parents (mom has to be from same population as dad)
        dad = rand(male_idx)
        mom = 0
        while true
            mom = rand(female_idx)
            populations[mom] == populations[dad] && break
        end
        push!(fathers, dad)
        push!(mothers, mom)
        # recombination
        recombine!(@view(x1[i, :]), @view(x2[i, :]), @view(h1[dad, :]),
                   @view(h2[dad, :]), @view(h1[mom, :]), @view(h2[mom, :]))
        # update progress
        next!(pmeter)
    end
    # combine offsprings and parents
    H1 = [h1; x1]
    H2 = [h2; x2]
    return H1, H2, fathers, mothers
end

function recombination_segments(breakpoints::Vector{Int}, snps::Int)
    start = 1
    result = UnitRange{Int}[]
    for bkpt in breakpoints
        push!(result, start:bkpt)
        start = bkpt + 1
    end
    push!(result, breakpoints[end]+1:snps)
    return result
end

function recombine!(child_h1, child_h2, dad_h1, dad_h2, mom_h1, mom_h2)
    p = length(child_h1)
    recombinations = rand(1:5)
    breakpoints = sort!(sample(1:p, recombinations, replace=false))
    segments = recombination_segments(breakpoints, p)
    for segment in segments
        dad_hap = rand() < 0.5 ? dad_h1 : dad_h2
        mom_hap = rand() < 0.5 ? mom_h1 : mom_h2
        copyto!(@view(child_h1[segment]), @view(dad_hap[segment]))
        copyto!(@view(child_h2[segment]), @view(mom_hap[segment]))
    end
end

function write_plink(outfile::AbstractString, x1::AbstractMatrix, x2::AbstractMatrix)
    n, p = size(x1)
    x = SnpArray(outfile * ".bed", n, p)
    for j in 1:p, i in 1:n
        c = x1[i, j] + x2[i, j]
        if c == 0
            x[i, j] = 0x00
        elseif c == 1
            x[i, j] = 0x02
        elseif c == 2
            x[i, j] = 0x03
        else
            error("matrix entries should be 0, 1, or 2 but was $c!")
        end
    end
    # create .bim file structure: https://www.cog-genomics.org/plink2/formats#bim
    open(outfile * ".bim", "w") do f
        for i in 1:p
            println(f, "1\tsnp$i\t0\t$(100i)\t1\t2")
        end
    end
    # create .fam file structure: https://www.cog-genomics.org/plink2/formats#fam
    open(outfile * ".fam", "w") do f
        for i in 1:n
            println(f, "$i\t1\t0\t0\t1\t-9")
        end
    end
    return nothing
end

function make_partition_mapfile(filename, p::Int)
    map_cM = LinRange(0.0, Int(p / 10000), p)
    open(filename, "w") do io
        println(io, "Chromosome\tPosition(bp)\tRate(cM/Mb)\tMap(cM)")
        for i in 1:p
            println(io, "chr1\t", 100i, '\t', 0.01rand(), '\t', map_cM[i])
        end
    end
end

function make_rapid_mapfile(filename, p::Int)
    map_cM = LinRange(0.0, Int(p / 10000), p)
    open(filename, "w") do io
        for i in 1:p
            println(io, i, '\t', map_cM[i])
        end
    end
end

function process_rapid_output(inputfile, outputfile)
    writer = open(outputfile, "w")
    df = readdlm(inputfile)
    println(writer, "CHR ID1 HID1 ID2 HID2 BP.start BP.end site.start site.end cM FAM1 FAM2")
    for r in eachrow(df)
        chr, id1, id2, hap1, hap2, start_pos, end_pos, genetic_len, start_site, end_site =
            Int(r[1]), Int(r[2]), Int(r[3]), Int(r[4]), Int(r[5]), Int(r[6]), Int(r[7]),
            r[8], Int(r[9]), Int(r[10])
        println(writer, chr, ' ', id1, ' ', hap1, ' ', id2, ' ', hap2, ' ',
            start_pos, ' ', end_pos, ' ', start_site, ' ', end_site, ' ',
            genetic_len, ' ', 1, ' ', 1)
    end
    close(writer)
end

function make_bgen_samplefile(filename, n)
    open(filename, "w") do io
        println(io, "ID_1 ID_2 missing sex")
        println(io, "0 0 0 D")
        for i in 1:n
            println(io, "$i 1 0 1")
        end
    end
end

function simulate_genotypes(n, p, offsprings, seed)
    Random.seed!(seed)
    h1, h2, populations, diff_markers = simulate_pop_structure(n, p)

    # simulate random mating to get IBD segments
    x1, x2 = simulate_IBD(h1, h2, populations, offsprings)

    # write phased genotypes to VCF format
    write_vcf("sim.phased.vcf.gz", x1, x2)

    # write unphased genotypes to PLINK binary format
    outfile = "sim"
    write_plink(outfile, x1, x2)

    # save pop1/pop2 index and unually differentiated marker indices
    writedlm("populations.txt", populations)
    writedlm("diff_markers.txt", diff_markers)

    # generate fake map file
    make_partition_mapfile("sim.partition.map", p)

    # also generate QC file that contains all SNPs and all samples
    snpdata = SnpData("sim")
    snpIDs = snpdata.snp_info[!, :snpid]
    sampleIDs = Matrix(snpdata.person_info[!, 1:2])
    writedlm("variants_qc.txt", snpIDs)
    writedlm("samples_qc.txt", sampleIDs)
end

function make_knockoffs(n, p, offsprings, seed)
    # fake map file for rapid
    make_rapid_mapfile("sim.rapid.map", p)

    # run rapid
    Random.seed!(seed)
    vcffile = "sim.phased.vcf.gz"
    mapfile = "sim.rapid.map"
    outfolder = "rapid"
    d = 3    # minimum IBD length in cM
    w = 3    # number of SNPs per window
    r = 10   # number of runs
    s = 2    # Minimum number of successes to consider a hit
    @time rapid(rapid_exe, vcffile, mapfile, d, outfolder, w, r, s)

    # unzip and postprocess rapid output to suit snpknock2
    run(pipeline(`gunzip -c ./rapid/results.max.gz`, stdout="./rapid/results.max"))
    process_rapid_output("./rapid/results.max", "sim.snpknock.ibdmap")

    # convert VCF to BGEN format
    outfile = "sim.bgen"
    run(`$qctools_exe -g $vcffile -og $outfile`)
    make_bgen_samplefile("sim.sample", n + offsprings)

    # snpknock2 arguments
    bgenfile = "sim"
    sample_qc = "samples_qc.txt"
    variant_qc = "variants_qc.txt"
    mapfile = "sim.partition.map"
    partfile = "sim.partition.txt"
    ibdfile = "sim.snpknock.ibdmap"
    K = 10
    cluster_size_min = 1000
    cluster_size_max = 10000
    hmm_rho = 1
    hmm_lambda = 1e-3
    windows = 0
    n_threads = Threads.nthreads()
    compute_references = true
    generate_knockoffs = true
    outfile = "sim.knockoffs"
    @time snpknock2(snpknock2_exe, bgenfile, sample_qc, variant_qc, mapfile, partfile, ibdfile,
        K, cluster_size_min, cluster_size_max, hmm_rho, hmm_lambda, windows, n_threads,
        seed, compute_references, generate_knockoffs, outfile)
end

# there is 100 causal snps, 50 on normally differentiated SNPs and 50 on abnormally diff snps
function simulate_beta_and_y(x::AbstractMatrix, diff_markers, seed)
    n, p = size(x)
    Random.seed!(seed)
    h2 = 0.5 # heritability
    k = 100 # number of causal snps
    d = Normal(0, sqrt(h2 / (2k))) # from paper: Efficient Implementation of Penalized Regression for Genetic Risk Prediction
    # simulate β
    β = zeros(p)
    norm_markers = setdiff(1:p, diff_markers)
    norm_causal_snps = sample(norm_markers, k >> 1, replace=false)
    β[norm_causal_snps] .= rand(d, k >> 1)
    shuffle!(diff_markers)
    diff_causal_markers = diff_markers[1:k >> 1]
    β[diff_causal_markers] .= rand(d, k >> 1)
    # simulate y
    ϵ = Normal(0, 1 - h2)
    y = x * β + rand(ϵ, n)
    # save results
    writedlm("y_true.txt", y)
    writedlm("beta_true.txt", β)
    writedlm("normal_markers.txt", norm_markers)
    writedlm("normal_causal_markers.txt", norm_causal_snps)
    writedlm("diff_markers.txt", diff_markers)
    writedlm("diff_causal_markers.txt", diff_causal_markers)
    return y, β
end

function R2(X::AbstractMatrix, y::AbstractVector, β̂::AbstractVector)
    μ = y - X * β̂
    tss = y .- mean(y)
    return 1 - dot(μ, μ) / dot(tss, tss)
end

function TP(correct_groups, signif_groups)
    return length(signif_groups ∩ correct_groups) / length(correct_groups)
end

function FDR(correct_groups, signif_groups)
    FP = length(signif_groups) - length(signif_groups ∩ correct_groups) # number of false positives
    # FPR = FP / (FP + TN) # https://en.wikipedia.org/wiki/False_positive_rate#Definition
    FDR = FP / length(signif_groups)
    return FDR
end

function tune_k(y::AbstractVector, xko_la::AbstractMatrix, original::Vector{Int},
    knockoff::Vector{Int}, groups::Vector{Int}, fdr::Float64, best_k::Int
    )
    # do a grid search for best sparsity level
    best_β = Float64[]
    best_err = Inf
    for cur_k in best_k:5:round(Int, 1.5best_k)
        result = fit_iht(y, xko_la, k=cur_k, init_beta=true, max_iter=500)
        W = coefficient_diff(result.beta, groups, original, knockoff)
        τ = threshold(W, fdr, :knockoff)
        detected = count(x -> x ≥ τ, W)
        if abs(detected - best_k) < best_err
            best_β = copy(result.beta)
            best_err = abs(detected - best_k)
        end
        println("wrapped CV says best_k = $best_k; using k = $cur_k detected $detected")
        GC.gc()
    end
    return best_β
end

function get_signif_groups(β, groups)
    correct_groups = Int[]
    for i in findall(!iszero, β)
        g = groups[i]
        g ∈ correct_groups || push!(correct_groups, g)
    end
    return correct_groups
end

function summarize_result(groups, β_true, β_iht, β_iht_knockoff, β_iht_knockoff_cv, β_lasso, β_lasso_knockoff)
    length(groups) == length(β_true) || error("summarize_result: groups and β_true have different length!")
    df = DataFrame(Metric = String[], IHT = Float64[], IHT_ko = Float64[],
        IHT_ko_cv = Float64[], LASSO = Float64[], LASSO_ko = Float64[])

    # count non-zero entries of β returned from cross validation
    push!(df, hcat("beta_non_zero_count", count(!iszero, β_iht),
        count(!iszero, vec(readdlm("iht.knockoff.beta"))),
        count(!iszero, vec(readdlm("iht.knockoff.cv.beta"))),
        count(!iszero, β_lasso),
        count(!iszero, vec(readdlm("lasso.knockoff.beta")))
        ))
    # count non-zero entries after knockoff filter
    push!(df, hcat("beta_selected", count(!iszero, β_iht),
        count(!iszero, β_iht_knockoff),
        count(!iszero, β_iht_knockoff_cv), count(!iszero, β_lasso),
        count(!iszero, β_lasso_knockoff)
        ))
    # count TP proportion
    correct_snps = findall(!iszero, vec(readdlm("beta_true.txt")))
    correct_groups = get_signif_groups(β_true, groups)
    push!(df, hcat("TPP", TP(correct_snps, findall(!iszero, β_iht)),
        TP(correct_groups, get_signif_groups(β_iht_knockoff, groups)),
        TP(correct_groups, get_signif_groups(β_iht_knockoff_cv, groups)),
        TP(correct_groups, get_signif_groups(β_lasso, groups)),
        TP(correct_groups, get_signif_groups(β_lasso_knockoff, groups))
        ))
    # count FDR
    push!(df, hcat("FDR", FDR(correct_snps, findall(!iszero, β_iht)),
        FDR(correct_groups, get_signif_groups(β_iht_knockoff, groups)),
        FDR(correct_groups, get_signif_groups(β_iht_knockoff_cv, groups)),
        FDR(correct_groups, get_signif_groups(β_lasso, groups)),
        FDR(correct_groups, get_signif_groups(β_lasso_knockoff, groups))
        ))
    @show df
    CSV.write("summary.txt", df)
end

function one_simulation(n, p, offsprings, seed)
    cur_dir = pwd() * "/"

    # save all simulated data under data/sim$seed
    isdir("data") || mkdir("data")
    isdir("data/sim$seed") || mkdir("data/sim$seed")
    cd("data/sim$seed")
    data_dir = cur_dir * "data/sim$seed/"

    # simulate phased genotypes with pop structure and cryptic relatedness
    simulate_genotypes(n, p, offsprings, seed)

    # partition genotypes to groups (for making group knockoffs)
    partition(partition_exe, "sim", "sim.partition.map", "variants_qc.txt", "sim.partition.txt")

    # make knockoffs
    make_knockoffs(n, p, offsprings, seed)

    # loop over group resolution
    for group in 0:3
        x = SnpArray(data_dir * "knockoffs/sim.knockoffs_res$(group).bed")
        snpid = SnpData(data_dir * "knockoffs/sim.knockoffs_res$(group)").snp_info.snpid
        knockoff_idx = endswith.(snpid, ".k")
        original = findall(knockoff_idx .== false)
        knockoff = findall(knockoff_idx)
        groups = repeat(Vector{Int}(readdlm(data_dir * "knockoffs/sim.knockoffs_res$(group)_grp.txt",
            header=true)[1][:, 2]), inner=2) .+ 1
        xla = convert(Matrix{Float64}, @view(x[:, original]), center=true, scale=true, impute=true)
        xko_la = convert(Matrix{Float64}, x, center=true, scale=true, impute=true)
        diff_markers = Vector{Int}(vec(readdlm(data_dir * "diff_markers.txt")))
        group_dir = cur_dir * "res$group/"
        isdir(group_dir) || mkdir(group_dir)
        cd(group_dir)

        for fdr in [0.05, 0.1, 0.25, 0.5]
            top_dir = group_dir * "fdr$fdr/"
            new_dir = group_dir * "fdr$fdr/sim$seed"
            isdir(top_dir) || mkdir(top_dir)
            isdir(new_dir) || mkdir(new_dir)
            cd(new_dir)

            #
            # simulate beta and y
            #
            Random.seed!(seed)
            y, β = simulate_beta_and_y(xla, diff_markers, seed)

            #
            # Run standard IHT
            #
            Random.seed!(seed)
            path = 10:10:200
            mses = cv_iht(y, xla, path=path, init_beta=true)
            GC.gc()
            Random.seed!(seed)
            k_rough_guess = path[argmin(mses)]
            dense_path = (k_rough_guess - 9):(k_rough_guess + 9)
            mses_new = cv_iht(y, xla, path=dense_path, init_beta=true)
            GC.gc()
            Random.seed!(seed)
            result = fit_iht(y, xla, k=dense_path[argmin(mses_new)], init_beta=true, max_iter=500)
            @show result
            β_iht = result.beta
            writedlm("iht.beta", β_iht)
            GC.gc()

            #
            # Run standard lasso
            #
            Random.seed!(seed)
            cv = glmnetcv(xla, y, nfolds=5, parallel=true)
            β_lasso = GLMNet.coef(cv)
            writedlm("lasso.beta", β_lasso)

            #
            # run knockoff IHT
            #
            Random.seed!(seed)
            path = 10:10:200
            mses = cv_iht(y, xko_la, path=path, init_beta=true)
            GC.gc()
            Random.seed!(seed)
            k_rough_guess = path[argmin(mses)]
            dense_path = (k_rough_guess - 9):(k_rough_guess + 9)
            mses_new = cv_iht(y, xko_la, path=dense_path, init_beta=true)
            GC.gc()
            Random.seed!(seed)
            result = fit_iht(y, xko_la, k=dense_path[argmin(mses_new)],
                init_beta=true, max_iter=500)
            @show result
            β_iht_knockoff = extract_beta(result.beta, fdr, groups, original, knockoff)
            writedlm("iht.knockoff.beta", result.beta)
            writedlm("iht.knockoff.postfilter.beta", β_iht_knockoff)

            #
            # run knockoff IHT with wrapped cross validation
            #
            Random.seed!(seed)
            path = 10:10:200
            z = ones(Float64, size(xla, 1))
            mses = cv_iht_knockoff(y, xko_la, z, original, knockoff, fdr, path=path,
                init_beta=true, group_ko=groups)
            GC.gc()
            Random.seed!(seed)
            k_rough_guess = path[argmin(mses)]
            dense_path = (k_rough_guess - 9):(k_rough_guess + 9)
            mses_new = cv_iht_knockoff(y, xko_la, z, original, knockoff, fdr,
                path=dense_path, init_beta=true, group_ko=groups)
            GC.gc()
            # adjust sparsity level so it best matches sparsity chosen by ko filter
            Random.seed!(seed)
            best_k = dense_path[argmin(mses_new)]
            best_β = tune_k(y, xko_la, original, knockoff, groups, fdr, best_k)
            β_iht_knockoff_cv = extract_beta(best_β, fdr, groups, original, knockoff)
            writedlm("iht.knockoff.cv.beta", best_β)
            writedlm("iht.knockoff.cv.postfilter.beta", β_iht_knockoff_cv)

            #
            # Run knockoff lasso
            #
            Random.seed!(seed)
            cv = glmnetcv(xko_la, y, nfolds=5, parallel=true)
            β_lasso_knockoff = extract_beta(GLMNet.coef(cv), fdr, groups, original, knockoff)
            writedlm("lasso.knockoff.beta", GLMNet.coef(cv))
            writedlm("lasso.knockoff.postfilter.beta", β_lasso_knockoff)

            #
            # process results
            #
            summarize_result(groups[original], β, β_iht, β_iht_knockoff, β_iht_knockoff_cv, β_lasso, β_lasso_knockoff)
            GC.gc()
        end
        GC.gc()
    end
end

n = 5000
p = 50000
offsprings = 500
seed = parse(Int, ARGS[1])
one_simulation(n, p, offsprings, seed)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule yq_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("yq")
JLLWrappers.@generate_main_file("yq", UUID("174e927f-d384-57ac-82db-e88370a29c21"))
end  # module yq_jll

ext: jl


text:
# mutable struct KinemDef3D
#     alpha :: MotionDef
#     h :: MotionDef
#     u :: MotionDef
#     vartype :: String
#     vary ::Int8
#     add :: Array{Float64,2}

#     function KinemDef3D(alpha :: MotionDef, h::MotionDef, u::MotionDef, vartype = "Constant", vary = 0, add = zeros(1,1))
#         new(alpha, h, u, vartype, vary, add)
#     end
# end

struct ThreeDFieldSimple
    f2d :: Vector{TwoDFlowField}
    function ThreeDFieldSimple()
        f2d = TwoDFlowField[]
    new(f2d)
    end
end

struct ThreeDSurfSimple
    cref :: Float64
    AR :: Float64
    uref :: Float64
    pvt :: Float64
    lespcrit :: Vector{Float64}
    coord_file :: String
    ndiv :: Int8
    nspan :: Int8
    naterm :: Int8
    kindef :: KinemDef
    psi :: Vector{Float64}
    yle :: Vector{Float64}
    s2d :: Vector{TwoDSurf}
    a03d :: Vector{Float64}
    bc :: Vector{Float64}
    nshed :: Vector{Float64}
    bcoeff :: Vector{Float64}
    aterm3d :: Array{Float64}
    levstr :: Vector{Float64}

    function ThreeDSurfSimple(AR, kindef, coord_file, pvt, lespcrit = [10.;]; nspan = 12, cref = 1., uref=1., ndiv=70, naterm=35)

        bref = AR*cref

        psi = zeros(nspan)
        yle = zeros(nspan)

        s2d = TwoDSurf[]

	for i = 1:nspan
            psi[i] = real(i)*(pi/2)/nspan
            yle[i] = -bref*cos(psi[i])/2.
        end

        for i = 1:nspan
            kinem2d = KinemDef(kindef.alpha, kindef.h, kindef.u)
            lespc = lespcrit[1]
            push!(s2d, TwoDSurf(coord_file, pvt,  kinem2d, [lespc;], c=cref, uref=uref, ndiv=ndiv, naterm=naterm))
        end

        a03d = zeros(nspan)
        aterm3d = zeros(naterm, nspan)

        bc = zeros(nspan)
        nshed = [0.;]
        bcoeff = zeros(nspan)
        levstr = zeros(nspan)

        new(cref, AR, uref, pvt, lespcrit, coord_file,  ndiv, nspan, naterm, kindef, psi, yle, s2d, a03d, bc, nshed, bcoeff, aterm3d, levstr)
    end
end

struct KelvinConditionQSLLT
    surf :: ThreeDSurfSimple
    field :: ThreeDFieldSimple
end

function (kelv::KelvinConditionQSLLT)(tev_iter::Array{Float64})

    val = zeros(kelv.surf.nspan)

    for i = 1:kelv.surf.nspan
        uprev, wprev = ind_vel([kelv.field.f2d[i].tev[end]], kelv.surf.s2d[i].bnd_x, kelv.surf.s2d[i].bnd_z)

        kelv.field.f2d[i].tev[end].s = tev_iter[i]

        unow, wnow = ind_vel([kelv.field.f2d[i].tev[end]], kelv.surf.s2d[i].bnd_x, kelv.surf.s2d[i].bnd_z)

        kelv.surf.s2d[i].uind[:] = kelv.surf.s2d[i].uind[:] .- uprev .+ unow
        kelv.surf.s2d[i].wind[:] = kelv.surf.s2d[i].wind[:] .- wprev .+ wnow

        #Calculate downwash
        update_downwash(kelv.surf.s2d[i], [kelv.field.f2d[i].u[1], kelv.field.f2d[i].w[1]])

        #Calculate first two fourier coefficients
        update_a0anda1(kelv.surf.s2d[i])

        kelv.surf.bc[i] = kelv.surf.s2d[i].a0[1] + 0.5*kelv.surf.s2d[i].aterm[1]
    end

    calc_a0a13d(kelv.surf)

    for i = 1:kelv.surf.nspan
        val[i] = kelv.surf.s2d[i].uref*kelv.surf.s2d[i].c*pi*(kelv.surf.bc[i] + kelv.surf.a03d[i] + 0.5*kelv.surf.aterm3d[1,i] - kelv.surf.s2d[i].a0prev[1] - 0.5*kelv.surf.s2d[i].aprev[1]) + kelv.field.f2d[i].tev[end].s
    end

    return val
end

struct KelvinKuttaQSLLT
    surf :: ThreeDSurfSimple
    field :: ThreeDFieldSimple
    shed_ind :: Vector{Int}
end

function (kelv::KelvinKuttaQSLLT)(v_iter::Array{Float64})

    nshed = length(kelv.shed_ind)
    val = zeros(kelv.surf.nspan+nshed)

    levcount = 0
    for i = 1:kelv.surf.nspan
        if i in kelv.shed_ind
            uprev, wprev = ind_vel([kelv.field.f2d[i].tev[end]; kelv.field.f2d[i].lev[end]], kelv.surf.s2d[i].bnd_x, kelv.surf.s2d[i].bnd_z)

            kelv.surf.levstr[i] -= kelv.field.f2d[i].lev[end].s

            kelv.field.f2d[i].tev[end].s = v_iter[i]
            levcount += 1
            kelv.field.f2d[i].lev[end].s = v_iter[kelv.surf.nspan+levcount]

            kelv.surf.levstr[i] += kelv.field.f2d[i].lev[end].s

            unow, wnow = ind_vel([kelv.field.f2d[i].tev[end]; kelv.field.f2d[i].lev[end]], kelv.surf.s2d[i].bnd_x, kelv.surf.s2d[i].bnd_z)

        else
            uprev, wprev = ind_vel([kelv.field.f2d[i].tev[end]], kelv.surf.s2d[i].bnd_x, kelv.surf.s2d[i].bnd_z)

            kelv.field.f2d[i].tev[end].s = v_iter[i]

            unow, wnow = ind_vel([kelv.field.f2d[i].tev[end]], kelv.surf.s2d[i].bnd_x, kelv.surf.s2d[i].bnd_z)
        end

        kelv.surf.s2d[i].uind[:] = kelv.surf.s2d[i].uind[:] .- uprev .+ unow
        kelv.surf.s2d[i].wind[:] = kelv.surf.s2d[i].wind[:] .- wprev .+ wnow

        #Calculate downwash
        update_downwash(kelv.surf.s2d[i], [kelv.field.f2d[i].u[1], kelv.field.f2d[i].w[1]])

        #Calculate first two fourier coefficients
        update_a0anda1(kelv.surf.s2d[i])

        kelv.surf.bc[i] = kelv.surf.s2d[i].a0[1] + 0.5*kelv.surf.s2d[i].aterm[1]
    end

    #calc_a0a13d_wlev(kelv.surf, kelv.shed_ind, v_iter[kelv.surf.nspan+1:end])
    calc_a0a13d_wlev(kelv.surf)

    levcount = 0
    for i = 1:kelv.surf.nspan
        if i in kelv.shed_ind
            val[i] = kelv.surf.s2d[i].uref*kelv.surf.s2d[i].c*pi*(kelv.surf.bc[i] + kelv.surf.a03d[i] + 0.5*kelv.surf.aterm3d[1,i] - kelv.surf.s2d[i].a0prev[1] - 0.5*kelv.surf.s2d[i].aprev[1]) + kelv.field.f2d[i].tev[end].s + kelv.field.f2d[i].lev[end].s
            levcount += 1
            if kelv.surf.s2d[i].a0[1] > 0
                lesp_cond = kelv.surf.s2d[i].lespcrit[1]
            else
                lesp_cond = -kelv.surf.s2d[i].lespcrit[1]
            end
            val[kelv.surf.nspan+levcount] = kelv.surf.s2d[i].a0[1] + kelv.surf.a03d[i] - lesp_cond

        else
            val[i] = kelv.surf.s2d[i].uref*kelv.surf.s2d[i].c*pi*(kelv.surf.bc[i] + kelv.surf.a03d[i] + 0.5*kelv.surf.aterm3d[1,i] - kelv.surf.s2d[i].a0prev[1] - 0.5*kelv.surf.s2d[i].aprev[1]) + kelv.field.f2d[i].tev[end].s
        end
    end

    return val
end

struct KelvinConditionQSLLT_sep
    surf :: ThreeDSurfSimple
    field :: ThreeDFieldSimple
    ispan :: Int
end

function (kelv::KelvinConditionQSLLT_sep)(tev_iter::Array{Float64})

    i = kelv.ispan

    uprev, wprev = ind_vel([kelv.field.f2d[i].tev[end]], kelv.surf.s2d[i].bnd_x, kelv.surf.s2d[i].bnd_z)

    kelv.field.f2d[i].tev[end].s = tev_iter[1]

    unow, wnow = ind_vel([kelv.field.f2d[i].tev[end]], kelv.surf.s2d[i].bnd_x, kelv.surf.s2d[i].bnd_z)

    kelv.surf.s2d[i].uind[:] = kelv.surf.s2d[i].uind[:] .- uprev .+ unow
    kelv.surf.s2d[i].wind[:] = kelv.surf.s2d[i].wind[:] .- wprev .+ wnow

    #Calculate downwash
    update_downwash(kelv.surf.s2d[i], [kelv.field.f2d[i].u[1], kelv.field.f2d[i].w[1]])

    #Calculate first two fourier coefficients
    update_a0anda1(kelv.surf.s2d[i])

    kelv.surf.bc[i] = kelv.surf.s2d[i].a0[1] + 0.5*kelv.surf.s2d[i].aterm[1]

    val = kelv.surf.s2d[i].uref*kelv.surf.s2d[i].c*pi*(kelv.surf.bc[i] + kelv.surf.a03d[i] + 0.5*kelv.surf.aterm3d[1,i] - kelv.surf.s2d[i].a0prev[1] - 0.5*kelv.surf.s2d[i].aprev[1]) + kelv.field.f2d[i].tev[end].s

    return val
end

struct KelvinKuttaQSLLT_sep
    surf :: ThreeDSurfSimple
    field :: ThreeDFieldSimple
    ispan :: Int
end

function (kelv::KelvinKuttaQSLLT_sep)(v_iter::Array{Float64})

    i = kelv.ispan
    val = zeros(2)

    uprev, wprev = ind_vel([kelv.field.f2d[i].tev[end]; kelv.field.f2d[i].lev[end]], kelv.surf.s2d[i].bnd_x, kelv.surf.s2d[i].bnd_z)

    kelv.field.f2d[i].tev[end].s = v_iter[1]
    kelv.field.f2d[i].lev[end].s = v_iter[2]

    unow, wnow = ind_vel([kelv.field.f2d[i].tev[end]; kelv.field.f2d[i].lev[end]], kelv.surf.s2d[i].bnd_x, kelv.surf.s2d[i].bnd_z)

    kelv.surf.s2d[i].uind[:] = kelv.surf.s2d[i].uind[:] .- uprev .+ unow
    kelv.surf.s2d[i].wind[:] = kelv.surf.s2d[i].wind[:] .- wprev .+ wnow

    #Calculate downwash
    update_downwash(kelv.surf.s2d[i], [kelv.field.f2d[i].u[1], kelv.field.f2d[i].w[1]])

    #Calculate first two fourier coefficients
    update_a0anda1(kelv.surf.s2d[i])

    kelv.surf.bc[i] = kelv.surf.s2d[i].a0[1] + 0.5*kelv.surf.s2d[i].aterm[1]

    val[1] = kelv.surf.s2d[i].uref*kelv.surf.s2d[i].c*pi*(kelv.surf.bc[i] + kelv.surf.a03d[i] + 0.5*kelv.surf.aterm3d[1,i] - kelv.surf.s2d[i].a0prev[1] - 0.5*kelv.surf.s2d[i].aprev[1]) + kelv.field.f2d[i].tev[end].s + kelv.field.f2d[i].lev[end].s

    if kelv.surf.s2d[i].a0[1] > 0
        lesp_cond = kelv.surf.s2d[i].lespcrit[1]
    else
        lesp_cond = -kelv.surf.s2d[i].lespcrit[1]
    end

    val[2] = kelv.surf.s2d[i].a0[1] + kelv.surf.a03d[i] - lesp_cond

    return val
end

mutable struct KinemCantilever
    eta :: Array{Float64}
    etad :: Array{Float64}
    etadd :: Array{Float64}
    eta_pr :: Array{Float64}
    etad_pr :: Array{Float64}
    etadd_pr :: Array{Float64}
    eta_pr2 :: Array{Float64}
    etad_pr2 :: Array{Float64}
    etadd_pr2 :: Array{Float64}
    eta_pr3 :: Array{Float64}
    etad_pr3 :: Array{Float64}
    etadd_pr3 :: Array{Float64}

    function KinemCantilever(eta, etad)
        z = zeros(length(eta))
        new(eta, etad, z, eta, etad, z, eta, etad, z, eta, etad, z)
    end
end


text:
using Test, bad

import Distributions

p = 0:0.1:1
for n in 0:15, x = 0:n, pp in p
    @test Distributions.pdf.(Distributions.Binomial(n, pp), x) ≈ pmf(x, n, pp)
    @test Distributions.cdf.(Distributions.Binomial(n, pp), x) ≈ cdf(x, n, pp)
end

p = Beta(5, 7)

# check that marginal distribution (of x1) sums to 1
for xpartial in 0:10, npartial in xpartial:10
    pmff = pmf.(0:10, 10, p; partial = (xpartial, npartial) )
    @test sum(pmff) ≈ 1.0
    @test cumsum(pmff) ≈ cdf.(0:10, 10, p; partial = (xpartial, npartial) )
end

function simon(r1::Int, n1::Int, r::Int, n::Int)
    nn2 = zeros(Int, n1 + 1)
    nn2[(r1 + 2):end] .= n - n1
    cc = repeat([Inf], n1 + 1)
    cc[(r1 + 2):end] .= r
    cc2 = cc .- collect(0:n1)
    cc2[cc2 .< 0]     .= -Inf
    cc2[cc2 .>= nn2]  .= Inf
    nn2[cc2 .== -Inf] .= 0
    Design(nn2, cc2)
end

design = simon(4, 18, 10, 33)
XX     = sample_space(design)
x1, x2 = XX[:,1], XX[:,2]

# check that joint distribution sums to one for all possible partial x1 observations
for x1partial in 0:n1(design)
    pmf_x1x2 = pmf.(x2, n2.(design, x1), p) .*
        pmf.(x1, n1(design), p; partial = (x1partial, n1(design)) )
    @test 1.0 ≈ sum( pmf_x1x2 )
end

# check that conditional distribution (of x2 given x2) sums to 1
for x1 in 0:n1(design)
    pmff = pmf.(0:n2(design, x1), n2(design, x1), p)
    @test sum(pmff) ≈ 1.0
    @test all( cumsum(pmff) .≈ cdf.(0:n2(design, x1), n2(design, x1), p))
end

@test true


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dcc_dropdown

"""
    dcc_dropdown(;kwargs...)

A Dropdown component.
Dropdown is an interactive dropdown element for selecting one or more
items.
The values and labels of the dropdown items are specified in the `options`
property and the selected item(s) are specified with the `value` property.

Use a dropdown when you have many options (more than 5) or when you are
constrained for space. Otherwise, you can use RadioItems or a Checklist,
which have the benefit of showing the users all of the items at once.
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique across all of the
components in an app.
- `options` (optional): An array of options {label: [string|number], value: [string|number]},
an optional disabled field can be used for each option. options has the following type: Array of lists containing elements 'label', 'value', 'disabled', 'title'.
Those elements have the following types:
  - `label` (String | Real; required): The dropdown's label
  - `value` (String | Real; required): The value of the dropdown. This value
corresponds to the items specified in the
`value` property.
  - `disabled` (Bool; optional): If true, this option is disabled and cannot be selected.
  - `title` (String; optional): The HTML 'title' attribute for the option. Allows for
information on hover. For more information on this attribute,
see https://developer.mozilla.org/en-US/docs/Web/HTML/Global_attributes/titles
- `value` (String | Real | Array of String | Reals; optional): The value of the input. If `multi` is false (the default)
then value is just a string that corresponds to the values
provided in the `options` property. If `multi` is true, then
multiple values can be selected at once, and `value` is an
array of items with values corresponding to those in the
`options` prop.
- `optionHeight` (Real; optional): height of each option. Can be increased when label lengths would wrap around
- `className` (String; optional): className of the dropdown element
- `clearable` (Bool; optional): Whether or not the dropdown is "clearable", that is, whether or
not a small "x" appears on the right of the dropdown that removes
the selected value.
- `disabled` (Bool; optional): If true, this dropdown is disabled and the selection cannot be changed.
- `multi` (Bool; optional): If true, the user can select multiple values
- `placeholder` (String; optional): The grey, default text shown when no option is selected
- `searchable` (Bool; optional): Whether to enable the searching feature or not
- `search_value` (String; optional): The value typed in the DropDown for searching.
- `style` (Dict; optional): Defines CSS styles which will override styles previously set.
- `loading_state` (optional): Object that holds the loading state object coming from dash-renderer. loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `persistence` (Bool | String | Real; optional): Used to allow user interactions in this component to be persisted when
the component - or the page - is refreshed. If `persisted` is truthy and
hasn't changed from its previous value, a `value` that the user has
changed while using the app will keep that change, as long as
the new `value` also matches what was given originally.
Used in conjunction with `persistence_type`.
- `persisted_props` (Array of a value equal to: 'value's; optional): Properties whose user interactions will persist after refreshing the
component or the page. Since only `value` is allowed this prop can
normally be ignored.
- `persistence_type` (a value equal to: 'local', 'session', 'memory'; optional): Where persisted user changes will be stored:
memory: only kept in memory, reset on page refresh.
local: window.localStorage, data is kept after the browser quit.
session: window.sessionStorage, data is cleared once the browser quit.
"""
function dcc_dropdown(; kwargs...)
        available_props = Symbol[:id, :options, :value, :optionHeight, :className, :clearable, :disabled, :multi, :placeholder, :searchable, :search_value, :style, :loading_state, :persistence, :persisted_props, :persistence_type]
        wild_props = Symbol[]
        return Component("dcc_dropdown", "Dropdown", "dash_core_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: asc-csa/dash-core-components

max_issues_repo_name: asc-csa/dash-core-components

max_forks_repo_name: asc-csa/dash-core-components


text:
@testset "@cstruct" begin
	@cstruct OpaqueStruct
	@cstruct OpaqueStruct
	@test_throws ErrorException sizeof(OpaqueStruct)
	@test sizeof(Ptr{OpaqueStruct}) == sizeof(Ptr{Cvoid})

	@cstruct EmptyStruct {
	} __packed__
	@test sizeof(EmptyStruct) == 0

	@cstruct CcharStruct {
		c::Cchar
	} __packed__
	@test sizeof(CcharStruct) == sizeof(Cchar)
	@test :c in propertynames(CcharStruct)

	@cstruct CintStruct {
		i::Cint
	} __packed__
	@test sizeof(CintStruct) == sizeof(Cint)
	@test :i in propertynames(CintStruct)

	@cstruct CstructStruct {
		i::CintStruct
	} __packed__
	@test sizeof(CstructStruct) == sizeof(CintStruct)
	@test :i in propertynames(CstructStruct)

	@cstruct CintCcharStruct {
		i::CintStruct
		c::CcharStruct
	} __packed__
	@test sizeof(CintCcharStruct) == sizeof(CintStruct)+sizeof(CcharStruct)
	@test :i in propertynames(CintCcharStruct)
	@test :c in propertynames(CintCcharStruct)

	@cstruct CcharCintPackedStruct {
		c::CcharStruct
		i::CintStruct
	} __packed__
	@test sizeof(CcharCintPackedStruct) == sizeof(CintStruct)+sizeof(CcharStruct)
	@test :c in propertynames(CcharCintPackedStruct)
	@test :i in propertynames(CcharCintPackedStruct)

	@cstruct CcharCintAlignedStruct {
		c::CcharStruct
		@calign sizeof(CintStruct)
		i::CintStruct
	} __packed__
	@test sizeof(CcharCintAlignedStruct) == max(sizeof(CintStruct), sizeof(Cint))+max(sizeof(CcharStruct), sizeof(Cint))
	@test :c in propertynames(CcharCintAlignedStruct)
	@test :i in propertynames(CcharCintAlignedStruct)

	@cstruct CunionStruct {
		@cunion {
			c::Cuchar[4]
			s::Cushort[2]
			i::Cuint
		} __packed__
	} __packed__
	@test :c in propertynames(CunionStruct)
	@test :s in propertynames(CunionStruct)
	@test :i in propertynames(CunionStruct)

	@cstruct PtrToStructStruct
	@cstruct PtrToPtrToStructStruct {
		p::Ptr{PtrToStructStruct}
	} __packed__
	@cstruct PtrToStructStruct {
		p::Ptr{PtrToStructStruct}
	} __packed__
	@test sizeof(PtrToPtrToStructStruct) == sizeof(Ptr)
	@test :p in propertynames(PtrToPtrToStructStruct)
	@test sizeof(PtrToStructStruct) == sizeof(Ptr)
	@test :p in propertynames(PtrToStructStruct)

	cus = CunionStruct()
	@test cus.i == 0
	@test cus.c[1] == 0
	@test cus.c[4] == 0
	@test cus.s[1] == 0
	@test cus.s[2] == 0
	cus.i = 0xff0000ff
	@test cus.i == 0xff0000ff
	@test cus.c[1] == 0xff
	@test cus.c[4] == 0xff
	@test cus.s[1] != 0
	@test cus.s[2] != 0

	@cstruct Cint32Bitfield {
		(i:32)::Cint
	} __packed__
	@test sizeof(Cint32Bitfield) == sizeof(Cint)
	@test :i in propertynames(Cint32Bitfield)

	@cstruct Cuint32Bitfield {
		(u:32)::Cuint
	} __packed__
	@test sizeof(Cuint32Bitfield) == sizeof(Cuint)
	@test :u in propertynames(Cuint32Bitfield)

	@cstruct Cint2Bitfield {
		(i:2)::Cint
	} __packed__
	@test sizeof(Cint2Bitfield) == sizeof(Cchar)
	@test :i in propertynames(Cint2Bitfield)

	@cstruct Cuint2Bitfield {
		(u:2)::Cuint
	} __packed__
	@test sizeof(Cuint2Bitfield) == sizeof(Cchar)
	@test :u in propertynames(Cuint2Bitfield)

	@cstruct Cuint32Cint32Bitfields {
		(u:32)::Cuint
		(i:32)::Cint
	} __packed__
	@test sizeof(Cuint32Cint32Bitfields) == sizeof(Cuint)+sizeof(Cint)
	@test :u in propertynames(Cuint32Cint32Bitfields)
	@test :i in propertynames(Cuint32Cint32Bitfields)

	@cstruct Cuint16Cint16Bitfields {
		(u:16)::Cuint
		(i:16)::Cint
	} __packed__
	@test sizeof(Cuint16Cint16Bitfields) == sizeof(Cuint)
	@test :u in propertynames(Cuint16Cint16Bitfields)
	@test :i in propertynames(Cuint16Cint16Bitfields)

	@cstruct Cint16Cint16Bitfields {
		(i1:16)::Cint
		(i2:16)::Cint
	} __packed__
	@test sizeof(Cint16Cint16Bitfields) == sizeof(Cuint)

	@cstruct Cint16CalignCint16Bitfields {
		(i1:16)::Cint
		@calign sizeof(Cint)
		(i2:16)::Cint
	} __packed__
	@test sizeof(Cint16CalignCint16Bitfields) == sizeof(Cint)*2

	@cstruct Cint16CintBitfields {
		(u:16)::Cuint
		i::Cint
	} __packed__
	@test sizeof(Cint16CintBitfields) == 2+sizeof(Cint)
	@test :u in propertynames(Cint16CintBitfields)
	@test :i in propertynames(Cint16CintBitfields)

	bf = Cuint16Cint16Bitfields()
	@test bf.u == 0
	@test bf.i == 0
	bf.i = Cint(-1)
	@test bf.u == 0
	@test bf.i == -1
	bf.u = Cuint(0x0001)
	@test bf.u == 0x0001
	@test bf.i == -1

	bf = Cint16CalignCint16Bitfields()
	@test bf.i1 == 0
	@test bf.i2 == 0
	bf.i1 = -1
	@test bf.i1 == -1
	@test bf.i2 == 0
	bf.i2 = -1
	@test bf.i1 == -1
	@test bf.i2 == -1
	bf.i2 = 32768
	@test bf.i1 == -1
	@test bf.i2 < 0

	# https://github.com/analytech-solutions/CBinding.jl/issues/6
	@cstruct AnonymousAggregateField {
		u::@cunion {
			x::Cint
			y::Cchar
		}
	}
	aaf = AnonymousAggregateField(u = (x = -1,))
	@test aaf.u.x == -1
	aaf.u = (x = 0,)
	@test aaf.u.x == 0

	@cstruct NestedAnonymousAggregateField {
		u::@cunion {
			x::@cunion {
				z::Cint
			}
			y::Cchar
		}
	}
	aaf = NestedAnonymousAggregateField(u = (x = (z = -1,),))
	@test aaf.u.x.z == -1
	aaf.u = (x = (z = 0,),)
	@test aaf.u.x.z == 0

	# https://github.com/analytech-solutions/CBinding.jl/issues/8
	@cstruct BrokenLayout {
		x::@cstruct {
			y::Cint
			z::Cint
		}
	}
	bl = BrokenLayout()
	@test bl.x.y == 0 && bl.x.z == 0
	bl.x.z = 123
	@test bl.x.y == 0 && bl.x.z == 123

	@cstruct ArrayOfStruct {
		x::BrokenLayout[2]
	}
	@test sizeof(ArrayOfStruct) == 2*2*sizeof(Cint)
	aos = ArrayOfStruct()
	@test typeof(aos.x) <: CBinding.Caccessor
	@test typeof(aos.x[]) <: Carray
	@test typeof(aos.x[2]) <: CBinding.Caccessor
	@test typeof(aos.x[2][]) <: BrokenLayout
	@test typeof(aos.x[2].x) <: CBinding.Caccessor
	@test typeof(aos.x[2].x.y) <: Cint
end


text:
using Literate
Literate.markdown("nextjournal.jl", "."; documenter=false)
Literate.notebook("nextjournal.jl", ".")


text:
using Documenter
using PolynomialGTM

makedocs(
    sitename = "PolynomialGTM",
    format = Documenter.HTML(),
    modules = [PolynomialGTM],
    pages = [
        "Quick Start" => "index.md",
        "Docstrings" => "docstrings.md"
    ]
)

deploydocs(
    target    = "build",
    repo      = "github.com/cadojo/PolynomialGTM.jl.git",
    branch    = "gh-pages",
    devbranch = "main",
    versions  = ["stable" => "v^", "manual", "v#.#", "v#.#.#"],
)

max_stars_repo_name: cadojo/PolynomialGTMModel

max_issues_repo_name: cadojo/PolynomialGTMModel

max_forks_repo_name: cadojo/PolynomialGTMModel


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.1
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# [quote="BambOoxX, post:12, topic:64568"]
# you define the `SCI` package such that it patches the `COM` package is this right ?
# [/quote]
#
# Yes.  It is a standard pattern.
#
# [quote="BambOoxX, post:12, topic:64568"]
# If so : What would happen if `COM` is not loaded prior to `SCI` .
# [/quote]
#
# Assume that, unlike the MWE I have shown above, both COM and SCI are packaged and can be loaded with `using COM` and `using SCI`.  Then the result of `using SCI; using COM` and the result of `using COM; using SCI` will be the same.
#
# In the following example, `COM` package has a function `prettify(x)` whose methods can be defined for user-defined types, and a function `prettyprint(x)` that uses `prettify(x)`.   `SCI` package makes use of them.
#
# __src/COM.jl of COM package:__
# ```julia
# """
# `COM.prettyprint(x)` prints nicely the object `x` for which `COM.prettify(x)` method is defined.
# """
# module COM
# prettify(x) = sprint(io -> show(io, "text/plain", x))
# prettyprint(io::IO, x) = print(io, prettify(x))
# prettyprint(x) = prettyprint(stdout, x)
# end
# ```
#
# __Example of COM__
# ```
# using COM
# COM.prettyprint([π, 2π])
# ```
# __Output:__
# ```
# 2-element Vector{Float64}:
#  3.141592653589793
#  6.283185307179586
# ```
#
# __src/SCI.jl of SCI package (with deps COM):__
# ```julia
# """A scirntific module (calculate exp(x))"""
# module SCI
#
# abstract type AbstractProblem end
# struct Problem{T} <: AbstractProblem x::T end
#
# abstract type AbstractAlgorithm end
# struct Builtin <: AbstractAlgorithm end
# Base.@kwdef struct Taylor <: AbstractAlgorithm n::Int = 10 end
# default_algorithm(prob::Problem) = Builtin()
#
# struct Solution{R, P<:AbstractProblem, A<:AbstractAlgorithm} result::R; prob::P; alg::A end
# solve(prob::AbstractProblem) = solve(prob, default_algorithm(prob))
# solve(prob::AbstractProblem, alg::Builtin) = Solution(exp(prob.x), prob, alg)
# solve(prob::AbstractProblem, alg::Taylor) = Solution(sum(prob.x^k/factorial(k) for k in 0:alg.n), prob, alg)
#
# using COM
#
# COM.prettify(sol::Solution{R, P, A}) where {R, P<:AbstractProblem, A<:Builtin} = """
# Problem:   x = $(sol.prob.x)
# Algorithm: builtin exp(x)
# Result:    $(sol.result)
# """
#
# COM.prettify(sol::Solution{R, P, A}) where {R, P<:AbstractProblem, A<:Taylor} = """
# Problem:   x = $(sol.prob.x)
# Algorithm: Taylor series of exp(x) upto degree $(sol.alg.n)
# Result:    $(sol.result)
# """
#
# end
# ```
#
# __Code to run:__
# ```julia
# # The order of the following two lines may be reversed.
# using SCI
# using COM
#
# prob = SCI.Problem(1)
# COM.prettyprint(SCI.solve(prob))
# println()
# COM.prettyprint(SCI.solve(prob, SCI.Taylor()))
# ```
#
# __Output:__
# ```
# Problem:   x = 1
# Algorithm: builtin exp(x)
# Result:    2.718281828459045
#
# Problem:   x = 1
# Algorithm: Taylor series of exp(x) upto degree 10
# Result:    2.7182818011463845
# ```
#
# In this way, you can use the pretty printing functions provided by `COM` package in any other package without changing any code in `COM` package.

# %%
# The order of the following two lines may be reversed.
using SCI
using COM

prob = SCI.Problem(1)
COM.prettyprint(SCI.solve(prob))
println()
COM.prettyprint(SCI.solve(prob, SCI.Taylor()))

max_stars_repo_path: 0011/Both SCI module and SCI2 module use COM module Part 2.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
# Julia program to be called from the command line to run continous variable simulation
# Remaining code is in the module ContVarEvolution.
# Example run:  julia run.jl examples/example1
# Example run:  julia run.jl examples/example1  1     # runs with random seed = 1
# Example run:  julia -p 4 run.jl examples/example1
using DataStructures
using Distributed
using Random
using Dates
using Statistics
@everywhere include("ContVarEvolution.jl")
#include("ContVarEvolution.jl")

@doc """function my_isdefined()
  S  should be a symbol that corresponds to a global variable.
  For example:  global A
  The corresponding symbol is :A
  my_isdefined(:A)  returns true if variable A has a value assigned to it, and false otherwise.
  A replacement for the  julia v6 function "isdefined()" which is deprecated in v7
  This definition works in both v6 and v7
"""
function my_isdefined( S::Symbol )
  try
    eval(S)
    return true
  catch
    return false
  end
end

@doc """ function save_params()
  Save the parameters specified in the configuration file into ContVarEvolution.cont_var_result() data structure.
  The parameters are global variables.
"""
function save_params()
  global mutation_stddev_list
  global N_mut_list
  # sim_record  is a record containing both the parameters and the results for a trial
  if my_isdefined(:mutation_stddev_list)
    sim_record = cont_var_result( N_list, num_attributes_list, mutation_stddev_list, Float64[], num_trials,N_list[1],
       num_subpops,num_attributes_list[1], ngens, burn_in,
       mutation_stddev_list[1], ideal, fit_slope, neutral )
  elseif my_isdefined( :N_mut_list )
    sim_record = cont_var_result( N_list, num_attributes_list, Float64[], N_mut_list, num_trials,N_list[1],
       num_subpops,num_attributes_list[1], ngens, burn_in,
       N_mut_list[1]/N_list[1], ideal, fit_slope, neutral )
  else
    error("Either mutation_stddev_list or N_mut_list must be defined in the configuration file.")
  end
  return sim_record
end

@doc """ function check_parameters()
  Check that all necessary parameters have been read from the parameter file.
  With julia v. 6.2 on CentOS, missing parameters can cause a segmentation fault.
"""
function check_parameters()
  global num_trials, N_list, num_subpops, num_attributes_list, ngens, burn_in, mutation_stddev_list, N_mut_list, ideal, fit_slope, neutral
  param_list = [:num_trials, :N_list, :num_subpops, :num_attributes_list, :ngens, :burn_in, :ideal, :fit_slope, :neutral]
  for p in param_list
    if !my_isdefined(p)
      error("The parameter $(String(p)) is not defined in the parameter file: $(simname).jl.")
    end
  end
end

if length(ARGS) == 0
  simname = "examples/example0"
else
  simname = ARGS[1]
  if length(ARGS) >= 2   # second command-line argument is random number seed
    seed = parse(Int,ARGS[2])
    println("seed: ",seed)
    srand(seed)
  end
end
include("$(simname).jl")
#println("simname: ",simname)
println("simtype: ",simtype)
check_parameters()
sim_record = save_params()
run_trials( simname, sim_record )


text:
# Please run from the `basic_gppp` directory.
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using Stheno, Random, Plots, Statistics

###########################  Define and inspect our model  ###########################

rng = MersenneTwister(123456);

#=
In this example, `f` is an unknown real-valued function that we wish to infer. To achieve
this, we have access to two sensors. The first returns noisy estimates of `f`, where we have
been reliably informed by whoever designed the sensor that the mean of the noise is given by
`sin(x) - 5 + sqrt(abs(x))`, and that it's variance is low (1e-2). How the designer
estimated this function, and why a sensor might possibly have such a strange mean error, is
beyond the scope of this example. The second returns biased measurements of `f`, where the
bias is known to be 3.5. The model below specifies a model for this scenario.
=#
model = @gppp let

    # Define a smooth latent process that we wish to infer.
    f = GP(SEKernel())

    # Define the two noise processes described.
    noise1 = sqrt(1e-2) * GP(WhiteKernel()) + (x->sin.(x) .- 5.0 .+ sqrt.(abs.(x)))
    noise2 = sqrt(1e-1) * GP(3.5, WhiteKernel())

    # Define the processes that we get to observe.
    y1 = f + noise1
    y2 = f + noise2
end;

# Generate some toy observations of `y₁` and `y₂`.
x1 = GPPPInput(:y1, sort(rand(rng, 3) * 10));
x2 = GPPPInput(:y2, sort(rand(rng, 10) * 10));
x = BlockData(x1, x2);
ŷ = rand(rng, model(x));
ŷ1, ŷ2 = split(x, ŷ);

# Compute the posterior processes.
model′ = posterior(model(x), ŷ);

# Sample jointly from the posterior processes and compute posterior marginals.
xp_ = range(-2.5, stop=12.5, length=500);
xp_f = GPPPInput(:f, xp_);
xp_y1 = GPPPInput(:y1, xp_);
xp_y2 = GPPPInput(:y2, xp_);
xp = BlockData(xp_f, xp_y1, xp_y2);
model′_xp = rand(rng, model′(xp, 1e-9));
f′xp, y1′xp, y2′xp = split(xp, model′_xp);

###########################  Plot results  ###########################

gr();

posterior_plot = plot();

# Plot posterior over y1.
plot!(posterior_plot, xp_, model′(xp_y1); color=:red, label="y1");
plot!(posterior_plot, xp_, y1′xp; color=:red, label="", linewidth=1, linealpha=0.2);

# Plot posterior over y2.
plot!(posterior_plot, xp_, model′(xp_y2); color=:green, label="y2");
plot!(posterior_plot, xp_, y2′xp; color=:green, label="", linewidth=1, linealpha=0.2);

# Plot posterior over f.
plot!(posterior_plot, xp_, model′(xp_f); color=:blue, label="Latent Function");
plot!(posterior_plot, xp_, f′xp; color=:blue, label="", linewidth=1, linealpha=0.2);

# Plot samples on which we conditioned.
scatter!(posterior_plot, x1.x, ŷ1;
    markercolor=:red,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.8,
    label="Sensor 1",
);
scatter!(posterior_plot, x2.x, ŷ2;
    markercolor=:green,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.8,
    label="Sensor 2",
);

display(posterior_plot);


text:
# example with constraint in struct
# struct Vector2D{T<:Number}
#    x::T
#    y::T
# end

# function dot(u::Vector2D{T}, v::Vector2D{T}) where T
#     u.x*v.x + u.y*v.y
# end

# function cross(u::Vector2D{T}, v::Vector2D{T}) where T
#     u.x*v.y - u.y*v.x
# end

# julia> a = Vector2D(1, 2)
# Vector2D{Int64}(1, 2)

# julia> b = Vector2D(3, 4)
# Vector2D{Int64}(3, 4)

# julia> c =Vector2D{Char}('A', 'B')
# ERROR: TypeError: in Vector2D, in T, expected T<:Number, got Type{Char}
# Stacktrace:
#  [1] top-level scope
#    @ REPL[11]:1

# julia> dot(a,b)
# 11

# julia> cross(a,b)
# -2

# dot(a,c)

# same example with simpler method signature
# in this version  u and v can be of different Number subtype and still works
# struct Vector2D{T<:Number}
#    x::T
#    y::T
# end

# function dot(u::Vector2D, v::Vector2D)
#     u.x*v.x + u.y*v.y
# end

# function cross(u::Vector2D, v::Vector2D)
#     u.x*v.y - u.y*v.x
# end

# julia> a = Vector2D(1, 2)
# Vector2D{Int64}(1, 2)

# julia> b = Vector2D(3.0, 4.0)
# Vector2D{Float64}(3.0, 4.0)

# julia> c =Vector2D{Char}('A', 'B')
# ERROR: TypeError: in Vector2D, in T, expected T<:Number, got Type{Char}
# Stacktrace:
#  [1] top-level scope
#    @ REPL[7]:1

# julia> dot(a,b)
# 11.0

# julia> cross(a,b)
# -2.0

# example with constraint in method - allows non Number Vecttor2
# but methods are type protected
# struct Vector2D{T}
#    x::T
#    y::T
# end

# function dot(u::Vector2D{T}, v::Vector2D{T}) where T <: Number
#     u.x*v.x + u.y*v.y
# end

# function cross(u::Vector2D{T}, v::Vector2D{T}) where T <: Number
#     u.x*v.y - u.y*v.x
# end

# julia> a = Vector2D(1, 2)
# Vector2D{Int64}(1, 2)

# julia> b = Vector2D(3, 4)
# Vector2D{Int64}(3, 4)

# julia> c =Vector2D{Char}('A', 'B')
# Vector2D{Char}('A', 'B')

# julia> a = Vector2D(1, 2)
# Vector2D{Int64}(1, 2)

# julia> b = Vector2D(3, 4)
# Vector2D{Int64}(3, 4)

# julia> c =Vector2D{Char}('A', 'B')
# Vector2D{Char}('A', 'B')

# julia> dot(a,b)
# 11

# julia> cross(a,b)
# -2

# julia> dot(a,c)
# ERROR: MethodError: no method matching dot(::Vector2D{Int64}, ::Vector2D{Char})
# Closest candidates are:
#   dot(::Vector2D{T}, ::Vector2D{T}) where T<:Number at REPL[4]:1
# Stacktrace:
#  [1] top-level scope
#    @ REPL[14]:1

# julia> cross(a,c)
# ERROR: MethodError: no method matching cross(::Vector2D{Int64}, ::Vector2D{Char})
# Closest candidates are:
#   cross(::Vector2D{T}, ::Vector2D{T}) where T<:Number at REPL[5]:1
# Stacktrace:
#  [1] top-level scope
#    @ REPL[15]:1

max_stars_repo_name: FourMInfo/JuliaforBeginners

max_issues_repo_name: FourMInfo/JuliaforBeginners

max_forks_repo_name: FourMInfo/JuliaforBeginners


text:
using JuMP, Cbc
using ExampleModeling

m, n, k = 3, 1, 1
a = [-8]
b = [-1]
c = [14, -33, 20]
A = reshape([-1, -4, 2], m, n)
B = reshape([-2, -1, 1], m, k)

specs = Specs(relax_integer=false)
indices = Indices(m, n, k)
params = Params(a, b, c, A, B)
model = ExampleModel(specs, indices, params)

set_optimizer(model, Cbc.Optimizer)
optimize!(model)

if termination_status(model) == MOI.INFEASIBLE
    println("Model is infeasible.")
    exit()
end

variables = Variables(model)
objectives = Objectives(model)

println(indices)
println(params)
println(variables)
println(objectives)

max_stars_repo_name: nnhjy/JuMP-project-template

max_issues_repo_name: nnhjy/JuMP-project-template

max_forks_repo_name: nnhjy/JuMP-project-template


text:
# CXXFFI - The Julia C++ FFI
#
# This file contains the julia parts of the C++ FFI.
# For bootstrapping purposes, there is a small C++ shim
# that we will call out to. In general I try to keep the amount of
# code duplication on the julia side to a minimum, even if this
# means having more code on the C++ side (the original version
# had a messy 4 step bootstrap process).
#
#
# The main interface provided by this file is the @cpp macro.
# The following usages are currently supportd.
#
#   - Stactic function call
#       @cpp mynamespace::func(args...)
#   - Membercall (where m is a CppPtr, CppRef or CppValue)
#       @cp m->foo(args...)
#
# Note that unary * inside a call, e.g. @cpp foo(*a) is treated as
# a (C++ side) derefence of a.

using Base.Meta
import Base: ==
import Base.Intrinsics.llvmcall

# # # Bootstap initialization

const libcxxffi = joinpath(dirname(Base.source_path()),string("../deps/usr/lib/libcxxffi",
    ccall(:jl_is_debugbuild, Cint, ()) != 0 ? "-debug" : ""))

function init()
    ccall((:init_julia_clang_env,libcxxffi),Void,())
    ccall((:init_header,libcxxffi),Void,(Ptr{Uint8},),joinpath(dirname(Base.source_path()),"test.h"))
end
init()

# # # Types we will use to represent C++ values

# The equiavlent of a C++ T<targs...> *
immutable CppPtr{T,targs}
    ptr::Ptr{Void}
end

==(p1::CppPtr,p2::Ptr) = p1.ptr == p2
==(p1::Ptr,p2::CppPtr) = p1 == p2.ptr

# The equiavlent of a C++ T<targs...> &
immutable CppRef{T,targs}
    ptr::Ptr{Void}
end

# The equiavlent of a C++ T<targs...>
immutable CppValue{T,targs}
    data::Vector{Uint8}
end

# Represents a forced cast form the value T
# (which may be any C++ compatible value)
# to the the C++ type To
immutable CppCast{T,To}
    from::T
end
CppCast{T,To}(from::T,::Type{To}) = CppCast{T,To}(from)
cast{T,To}(from::T,::Type{To}) = CppCast{T,To}(from)

# Represents a C++ Deference
immutable CppDeref{T}
    val::T
end
CppDeref{T}(val::T) = CppDeref{T}(val)

macro pcpp_str(s)
    CppPtr{symbol(s),()}
end

macro rcpp_str(s)
    CppRef{symbol(s),()}
end

macro vcpp_str(s)
    CppValue{symbol(s),()}
end

import Base: cconvert

cconvert(::Type{Ptr{Void}},p::CppPtr) = p.ptr
cconvert(::Type{Ptr{Void}},p::CppRef) = p.ptr

# Bootstrap definitions
pointerTo(t::pcpp"clang::Type") = pcpp"clang::Type"(ccall((:getPointerTo,libcxxffi),Ptr{Void},(Ptr{Void},),t.ptr))
referenceTo(t::pcpp"clang::Type") = pcpp"clang::Type"(ccall((:getReferenceTo,libcxxffi),Ptr{Void},(Ptr{Void},),t.ptr))

tovdecl(p::pcpp"clang::Decl") = pcpp"clang::ValueDecl"(ccall((:tovdecl,libcxxffi),Ptr{Void},(Ptr{Void},),p.ptr))
tovdecl(p::pcpp"clang::ParmVarDecl") = pcpp"clang::ValueDecl"(p.ptr)

CreateDeclRefExpr(p::pcpp"clang::ValueDecl",ty::pcpp"clang::Type") = pcpp"clang::Expr"(ccall((:CreateDeclRefExpr,libcxxffi),Ptr{Void},(Ptr{Void},Ptr{Void}),p.ptr,ty.ptr))
CreateDeclRefExpr(p,ty::pcpp"clang::Type") = CreateDeclRefExpr(tovdecl(p),ty)

CreateParmVarDecl(p::pcpp"clang::Type") = pcpp"clang::ParmVarDecl"(ccall((:CreateParmVarDecl,libcxxffi),Ptr{Void},(Ptr{Void},),p.ptr))

CreateMemberExpr(base::pcpp"clang::Expr",isarrow::Bool,member::pcpp"clang::ValueDecl") = pcpp"clang::Expr"(ccall((:CreateMemberExpr,libcxxffi),Ptr{Void},(Ptr{Void},Cint,Ptr{Void}),base.ptr,isarrow,member.ptr))
BuildCallToMemberFunction(me::pcpp"clang::Expr", args::Vector{pcpp"clang::Expr"}) = pcpp"clang::Expr"(ccall((:build_call_to_member,libcxxffi),Ptr{Void},(Ptr{Void},Ptr{Ptr{Void}},Csize_t),
        me.ptr,[arg.ptr for arg in args],length(args)))

BuildMemberReference(base, t, IsArrow, name) =
    pcpp"clang::Expr"(ccall((:BuildMemberReference,libcxxffi),Ptr{Void},(Ptr{Void},Ptr{Void},Cint,Ptr{Uint8}),base, t, IsArrow, name))

DeduceReturnType(expr) = pcpp"clang::Type"(ccall((:DeduceReturnType,libcxxffi),Ptr{Void},(Ptr{Void},),expr))

CreateFunction(rt,argt) = pcpp"llvm::Function"(ccall((:CreateFunction,libcxxffi),Ptr{Void},(Ptr{Void},Ptr{Void},Csize_t),rt,argt,length(argt)))

ExtractValue(v::pcpp"llvm::Value",idx) = pcpp"llvm::Value"(ccall((:create_extract_value,libcxxffi),Ptr{Void},(Ptr{Void},Csize_t),v.ptr,idx))

tollvmty(t::pcpp"clang::Type") = pcpp"llvm::Type"(ccall((:tollvmty,libcxxffi),Ptr{Void},(Ptr{Void},),t.ptr))

getTemplateArgs(tmplt::pcpp"clang::ClassTemplateSpecializationDecl") =
    rcpp"clang::TemplateArgumentList"(ccall((:getTemplateArgs,libcxxffi),Ptr{Void},(Ptr{Void},),tmplt))

getTargsSize(targs) =
 ccall((:getTargsSize,libcxxffi),Csize_t,(Ptr{Void},),targs)

getTargTypeAtIdx(targs, i) =
    pcpp"clang::Type"(ccall((:getTargTypeAtIdx,libcxxffi),Ptr{Void},(Ptr{Void},Csize_t),targs,i))

getUndefValue(t::pcpp"llvm::Type") =
    pcpp"llvm::Value"(ccall((:getUndefValue,libcxxffi),Ptr{Void},(Ptr{Void},),t))

getStructElementType(t::pcpp"llvm::Type",i) =
    pcpp"llvm::Type"(ccall((:getStructElementType,libcxxffi),Ptr{Void},(Ptr{Void},Uint32),t,i))

CreateRet(builder,val::pcpp"llvm::Value") =
    pcpp"llvm::Value"(ccall((:CreateRet,libcxxffi),Ptr{Void},(Ptr{Void},Ptr{Void}),builder,val))

CreateRetVoid(builder) =
    pcpp"llvm::Value"(ccall((:CreateRetVoid,libcxxffi),Ptr{Void},(Ptr{Void},),builder))

CreateBitCast(builder,val,ty) =
    pcpp"llvm::Value"(ccall((:CreateBitCast,libcxxffi),Ptr{Void},(Ptr{Void},Ptr{Void},Ptr{Void}),builder,val,ty))

BuildCXXTypeConstructExpr(t::pcpp"clang::Type", exprs::Vector{pcpp"clang::Expr"}) =
    pcpp"clang::Expr"(ccall((:typeconstruct,libcxxffi),Ptr{Void},(Ptr{Void},Ptr{Ptr{Void}},Csize_t),
        t,[expr.ptr for expr in exprs],length(exprs)))

cxxsizeof(t::pcpp"clang::CXXRecordDecl") = ccall((:cxxsizeof,libcxxffi),Csize_t,(Ptr{Void},),t)
cxxsizeof(t::pcpp"clang::Type") = ccall((:cxxsizeofType,libcxxffi),Csize_t,(Ptr{Void},),t)

createDerefExpr(e::pcpp"clang::Expr") = pcpp"clang::Expr"(ccall((:createDerefExpr,libcxxffi),Ptr{Void},(Ptr{Void},),e.ptr))

getType(v::pcpp"llvm::Value") = pcpp"llvm::Type"(ccall((:getValueType,libcxxffi),Ptr{Void},(Ptr{Void},),v))

isPointerType(t::pcpp"llvm::Type") = ccall((:isLLVMPointerType,libcxxffi),Cint,(Ptr{Void},),t) > 0

CreateConstGEP1_32(builder,x,idx) = pcpp"llvm::Value"(ccall((:CreateConstGEP1_32,libcxxffi),Ptr{Void},(Ptr{Void},Ptr{Void},Uint32),builder,x,uint32(idx)))

getPointerTo(t::pcpp"llvm::Type") = pcpp"llvm::Type"(ccall((:getLLVMPointerTo,libcxxffi),Ptr{Void},(Ptr{Void},),t))

CreateLoad(builder,val::pcpp"llvm::Value") = pcpp"llvm::Value"(ccall((:createLoad,libcxxffi),Ptr{Void},(Ptr{Void},Ptr{Void}),builder.ptr,val.ptr))

BuildDeclarationNameExpr(name, ctx::pcpp"clang::DeclContext") =
    pcpp"clang::Expr"(ccall((:BuildDeclarationNameExpr,libcxxffi),Ptr{Void},(Ptr{Uint8},Ptr{Void}),name,ctx))

getContext(decl::pcpp"clang::Decl") = pcpp"clang::DeclContext"(ccall((:getContext,libcxxffi),Ptr{Void},(Ptr{Void},),decl))

CreateCallExpr(Fn::pcpp"clang::Expr",args::Vector{pcpp"clang::Expr"}) = pcpp"clang::Expr"(
    ccall((:CreateCallExpr,libcxxffi),Ptr{Void},(Ptr{Void},Ptr{Ptr{Void}},Csize_t),
        Fn,[arg.ptr for arg in args],length(args)))

toLLVM(t::pcpp"clang::Type") = pcpp"llvm::Type"(ccall((:ConvertTypeForMem,libcxxffi),Ptr{Void},(Ptr{Void},),t))

getDirectCallee(t::pcpp"clang::CallExpr") = pcpp"clang::FunctionDecl"(ccall((:getDirectCallee,libcxxffi),Ptr{Void},(Ptr{Void},),t))
getCalleeReturnType(t::pcpp"clang::CallExpr") = pcpp"clang::Type"(ccall((:getCalleeReturnType,libcxxffi),Ptr{Void},(Ptr{Void},),t))

const CK_Dependent      = 0
const CK_BitCast        = 1
const CK_LValueBitCast  = 2
const CK_LValueToRValue = 3
const CK_NoOp           = 4
const CK_BaseToDerived  = 5
const CK_DerivedToBase  = 6
const CK_UncheckedDerivedToBase = 7
const CK_Dynamic = 8
const CK_ToUnion = 9
const CK_ArrayToPointerDecay = 10
const CK_FunctionToPointerDecay = 11
const CK_NullToPointer = 12
const CK_NullToMemberPointer = 13
const CK_BaseToDerivedMemberPointer = 14
const CK_DerivedToBaseMemberPointer = 15
const CK_MemberPointerToBoolean = 16
const CK_ReinterpretMemberPointer = 17
const CK_UserDefinedConversion = 18
const CK_ConstructorConversion = 19
const CK_IntegralToPointer = 20
const CK_PointerToIntegral = 21
const CK_PointerToBoolean = 22
const CK_ToVoid = 23
const CK_VectorSplat = 24
const CK_IntegralCast = 25

createCast(arg,t,kind) = pcpp"clang::Expr"(ccall((:createCast,libcxxffi),Ptr{Void},(Ptr{Void},Ptr{Void},Cint),arg,t,kind))
# Built-in clang types
chartype() = pcpp"clang::Type"(unsafe_load(cglobal((:cT_cchar,libcxxffi),Ptr{Void})))
cpptype(::Type{Uint8}) = chartype()#pcpp"clang::Type"(unsafe_load(cglobal(:cT_uint8,Ptr{Void})))
cpptype(::Type{Int8}) = pcpp"clang::Type"(unsafe_load(cglobal((:cT_int8,libcxxffi),Ptr{Void})))
cpptype(::Type{Uint32}) = pcpp"clang::Type"(unsafe_load(cglobal((:cT_uint32,libcxxffi),Ptr{Void})))
cpptype(::Type{Int32}) = pcpp"clang::Type"(unsafe_load(cglobal((:cT_int32,libcxxffi),Ptr{Void})))
cpptype(::Type{Uint64}) = pcpp"clang::Type"(unsafe_load(cglobal((:cT_uint64,libcxxffi),Ptr{Void})))
cpptype(::Type{Int64}) = pcpp"clang::Type"(unsafe_load(cglobal((:cT_int64,libcxxffi),Ptr{Void})))
cpptype(::Type{Bool}) = pcpp"clang::Type"(unsafe_load(cglobal((:cT_int1,libcxxffi),Ptr{Void})))
cpptype(::Type{Void}) = pcpp"clang::Type"(unsafe_load(cglobal((:cT_void,libcxxffi),Ptr{Void})))

# CXX Level Casting

for (rt,argt) in ((pcpp"clang::ClassTemplateSpecializationDecl",pcpp"clang::Decl"),
                  (pcpp"clang::CXXRecordDecl",pcpp"clang::Decl"))
    s = split(string(rt.parameters[1]),"::")[end]
    isas = symbol(string("isa",s))
    ds = symbol(string("dcast",s))
    # @cxx llvm::isa{$rt}(t)
    @eval $(isas)(t::$(argt)) = ccall($(quot(isas)),Int,(Ptr{Void},),t) != 0
    @eval $(ds)(t::$(argt)) = ($rt)(ccall($(quot(ds)),Ptr{Void},(Ptr{Void},),t))
end

# Clang Type* Bootstrap

for s in (:isVoidType,:isBooleanType,:isPointerType,:isReferenceType,
    :isCharType,:isIntegerType)
    @eval ($s)(t::pcpp"clang::Type") = ccall($(quot(s)),Int,(Ptr{Void},),t) != 0
end

for (r,s) in ((pcpp"clang::CXXRecordDecl",:getPointeeCXXRecordDecl),
              (pcpp"clang::CXXRecordDecl",:getAsCXXRecordDecl))
    @eval ($s)(t::pcpp"clang::Type") = ($(quot(r)))(ccall($(quot(s)),Ptr{Void},(Ptr{Void},),t))
end

# Main interface

immutable CppExpr{T,targs}; end

# Checks the arguments to make sure we have concrete C++ compatible
# types in the signature. This is used both in cases where type inference
# cannot infer the appropriate argument types during staging (in which case
# the error here will cause it to fall back to runtime) as well as when the
# user has bad types (in which case you'll get a compile-time (but not stage-time)
# error.
function check_args(argt,f)
    for (i,t) in enumerate(argt)
        if isa(t,UnionType) || (isa(t,DataType) && t.abstract) ||
            (!(t <: CppPtr) && !(t <: CppRef) && !(t <: CppValue) && !(t <: CppCast) && !in(t,
                [Bool, Uint8, Uint32, Uint64, Int64, Ptr{Void}, Ptr{Uint8}]))
            error("Got bad type information while compiling $f (got $t for argument $i)")
        end
    end
end

# # # Clang Level Decl Access

# Do C++ name lookup
translation_unit() = pcpp"clang::Decl"(ccall((:tu_decl,libcxxffi),Ptr{Void},()))

primary_ctx(p::pcpp"clang::DeclContext") =
    pcpp"clang::DeclContext"(p == C_NULL ? C_NULL : ccall((:get_primary_dc,libcxxffi),Ptr{Void},(Ptr{Void},),p))

toctx(p::pcpp"clang::Decl") =
    pcpp"clang::DeclContext"(p == C_NULL ? C_NULL : ccall((:decl_context,libcxxffi),Ptr{Void},(Ptr{Void},),p))

to_decl(p::pcpp"clang::DeclContext") =
    pcpp"clang::Decl"(p == C_NULL ? C_NULL : ccall((:to_decl,libcxxffi),Ptr{Void},(Ptr{Void},),p))

function lookup_name(fname::String, ctx::pcpp"clang::DeclContext")
    @assert ctx != C_NULL
    pcpp"clang::Decl"(
        ccall((:lookup_name,libcxxffi),Ptr{Void},(Ptr{Uint8},Ptr{Void}),bytestring(fname),ctx))
end

function lookup_name(parts)
    cur = translation_unit()
    for fpart in parts
        cur = lookup_name(fpart,primary_ctx(toctx(cur)))
        cur == C_NULL && error("Could not find $fpart as part of $(join(parts,"::")) lookup")
    end
    cur
end

lookup_ctx(fname::String) = lookup_name(split(fname,"::"))
lookup_ctx(fname::Symbol) = lookup_ctx(string(fname))

function specialize_template(cxxt::pcpp"clang::ClassTemplateDecl",targs,cpptype)
    @assert cxxt != C_NULL
    ts = pcpp"clang::Type"[cpptype(t) for t in targs]
    d = pcpp"clang::ClassTemplateSpecializationDecl"(ccall((:SpecializeClass,libcxxffi),Ptr{Void},
            (Ptr{Void},Ptr{Void},Uint32),cxxt.ptr,[p.ptr for p in ts],length(ts)))
    d
end

function cppdecl(fname,targs)
    # Let's get a clang level representation of this type
    ctx = lookup_ctx(fname)

    # If this is a templated class, we need to do template
    # resolution
    cxxt = cxxtmplt(ctx)
    if cxxt != C_NULL
        ts = map(cpptype,targs)
        deduced_class = specialize_template(cxxt,targs,cpptype)
        ctx = deduced_class
    end

    ctx
end
cppdecl{fname,targs}(::Union(Type{CppPtr{fname,targs}}, Type{CppValue{fname,targs}},
    Type{CppRef{fname,targs}})) = cppdecl(fname,targs)

# @cpp (@cpp dyn_cast{vcpp"clang::TypeDecl"}(d))->getTypeForDecl()
typeForDecl(d::pcpp"clang::Decl") = pcpp"clang::Type"(ccall((:typeForDecl,libcxxffi),Ptr{Void},(Ptr{Void},),d))
typeForDecl(d::pcpp"clang::CXXRecordDecl") = typeForDecl(pcpp"clang::Decl"(d.ptr))
typeForDecl(d::pcpp"clang::ClassTemplateSpecializationDecl") = typeForDecl(pcpp"clang::Decl"(d.ptr))

cpptype{s,t}(p::Type{CppPtr{s,t}}) = pointerTo(typeForDecl(cppdecl(p)))
cpptype{T<:CppPtr}(p::Type{CppPtr{T,()}}) = pointerTo(cpptype(T))
cpptype{s,t}(p::Type{CppValue{s,t}}) = typeForDecl(cppdecl(p))
cpptype{s,t}(p::Type{CppRef{s,t}}) = referenceTo(typeForDecl(cppdecl(p)))
cpptype{T}(p::Type{Ptr{T}}) = pointerTo(cpptype(T))

function _decl_name(d)
    @assert d != C_NULL
    s = ccall((:decl_name,libcxxffi),Ptr{Uint8},(Ptr{Void},),d)
    ret = bytestring(s)
    c_free(s)
    ret
end

function simple_decl_name(d)
    s = ccall((:simple_decl_name,libcxxffi),Ptr{Uint8},(Ptr{Void},),d)
    ret = bytestring(s)
    c_free(s)
    ret
end

get_pointee_name(t) = _decl_name(getPointeeCXXRecordDecl(t))
get_name(t) = _decl_name(getAsCXXRecordDecl(t))

function juliatype(t::pcpp"clang::Type")
    if isVoidType(t)
        return Void
    elseif isBooleanType(t)
        return Bool
    elseif isPointerType(t)
        cxxd = getPointeeCXXRecordDecl(t)
        if cxxd != C_NULL
            return CppPtr{symbol(get_pointee_name(t)),()}
        else
            pt = pcpp"clang::Type"(ccall((:referenced_type,libcxxffi),Ptr{Void},(Ptr{Void},),t.ptr))
            tt = juliatype(pt)
            if tt <: CppPtr
                return CppPtr{tt,()}
            else
                return Ptr{tt}
            end
        end
    elseif isReferenceType(t)
        t = pcpp"clang::Type"(ccall((:referenced_type,libcxxffi),Ptr{Void},(Ptr{Void},),t.ptr))
        return CppRef{symbol(get_name(t)),()}
    elseif isCharType(t)
        return Uint8
    elseif isIntegerType(t)
        if t == cpptype(Int64)
            return Int64
        elseif t == cpptype(Uint64)
            return Uint64
        elseif t == cpptype(Uint32)
            return Uint32
        elseif t == cpptype(Int32)
            return Int32
        elseif t == cpptype(Uint8) || t == chartype()
            return Uint8
        elseif t == cpptype(Int8)
            return Int8
        end
        # This is wrong. Might be any int. Need to access the BuiltinType::Kind Enum
        return Int
    else
        rd = getAsCXXRecordDecl(t)
        if rd.ptr != C_NULL
            targt = ()
            if isaClassTemplateSpecializationDecl(pcpp"clang::Decl"(rd.ptr))
                tmplt = dcastClassTemplateSpecializationDecl(pcpp"clang::Decl"(rd.ptr))
                targs = getTemplateArgs(tmplt)
                args = pcpp"clang::Type"[]
                for i = 0:(getTargsSize(targs)-1)
                    push!(args,getTargTypeAtIdx(targs,i))
                end
                targt = tuple(map(juliatype,args)...)
            end
            return CppValue{symbol(get_name(t)),targt}
        end
    end
    return Ptr{Void}
end

julia_to_llvm(x::ANY) = pcpp"llvm::Type"(ccall(:julia_type_to_llvm,Ptr{Void},(Any,),x))
# Various clang conversions (bootstrap definitions)

# @cxx llvm::dyn_cast{vcpp"clang::ClassTemplateDecl"}
cxxtmplt(p::pcpp"clang::Decl") = pcpp"clang::ClassTemplateDecl"(ccall((:cxxtmplt,libcxxffi),Ptr{Void},(Ptr{Void},),p))

stripmodifier{s,targs}(p::Union(Type{CppPtr{s,targs}},
    Type{CppRef{s,targs}}, Type{CppValue{s,targs}})) = p
stripmodifier{T,To}(p::Type{CppCast{T,To}}) = T
stripmodifier{T}(p::Type{CppDeref{T}}) = T
stripmodifier{T}(p::Type{Ptr{T}}) = Ptr{T}
stripmodifier(p::Union(Type{Bool},Type{Int64},Type{Uint32})) = p

resolvemodifier{s,targs}(p::Union(Type{CppPtr{s,targs}}, Type{CppRef{s,targs}},
    Type{CppValue{s,targs}}), e::pcpp"clang::Expr") = e
resolvemodifier(p::Union(Type{Bool},Type{Int64}, Type{Uint32}), e::pcpp"clang::Expr") = e
resolvemodifier{T}(p::Type{Ptr{T}}, e::pcpp"clang::Expr") = e
resolvemodifier{T,To}(p::Type{CppCast{T,To}}, e::pcpp"clang::Expr") =
    createCast(e,cpptype(To),CK_BitCast)
resolvemodifier{T}(p::Type{CppDeref{T}}, e::pcpp"clang::Expr") =
    createDerefExpr(e)

resolvemodifier_llvm{s,targs}(builder, t::Type{CppCast{s,targs}}, v::pcpp"llvm::Value") =
    resolvemodifier_llvm(builder, t.parameters[1], ExtractValue(v,0))

resolvemodifier_llvm{s,targs}(builder, t::Union(Type{CppPtr{s,targs}}, Type{CppRef{s,targs}}),
        v::pcpp"llvm::Value") = ExtractValue(v,0)

resolvemodifier_llvm{ptr}(builder, t::Type{Ptr{ptr}}, v::pcpp"llvm::Value") = v
resolvemodifier_llvm(builder, t::Union(Type{Int64},Type{Uint32},Type{Bool}), v::pcpp"llvm::Value") = v
#resolvemodifier_llvm(builder, t::Type{Uint8}, v::pcpp"llvm::Value") = v

function resolvemodifier_llvm{s,targs}(builder, t::Type{CppValue{s,targs}}, v::pcpp"llvm::Value")
    ty = cpptype(t)
    @assert isPointerType(getType(v))
    # Get the array
    array = CreateConstGEP1_32(builder,v,1)
    arrayp = CreateLoad(builder,CreateBitCast(builder,array,getPointerTo(getType(array))))
    # Get the data pointer
    data = CreateConstGEP1_32(builder,arrayp,1)
    dp = CreateBitCast(builder,data,getPointerTo(getPointerTo(tollvmty(ty))))
    # A pointer to the actual data
    CreateLoad(builder,dp)
end

# LLVM-level manipulation
function llvmargs(builder, f, argt)
    args = Array(pcpp"llvm::Value", length(argt))
    for i in 1:length(argt)
        t = argt[i]
        args[i] = pcpp"llvm::Value"(ccall((:get_nth_argument,libcxxffi),Ptr{Void},(Ptr{Void},Csize_t),f,i-1))
        args[i] = resolvemodifier_llvm(builder, t, args[i])
        if args[i] == C_NULL
            error("Failed to process argument")
        end
    end
    args
end

# C++ expression manipulation
function buildargexprs(argt)
    callargs = pcpp"clang::Expr"[]
    pvds = pcpp"clang::ParmVarDecl"[]
    for i in 1:length(argt)
        #@show argt[i]
        t = argt[i]
        st = stripmodifier(t)
        argit = cpptype(st)
        st <: CppValue && (argit = pointerTo(argit))
        argpvd = CreateParmVarDecl(argit)
        push!(pvds, argpvd)
        expr = CreateDeclRefExpr(argpvd,argit)
        st <: CppValue && (expr = createDerefExpr(expr))
        expr = resolvemodifier(t, expr)
        push!(callargs,expr)
    end
    callargs, pvds
end

function associateargs(builder,argt,args,pvds)
    for i = 1:length(args)
        t = stripmodifier(argt[i])
        argit = cpptype(t)
        if t <: CppValue
            argit = pointerTo(argit)
        end
        AssociateValue(pvds[i],argit,args[i])
    end
end

AssociateValue(d::pcpp"clang::ParmVarDecl", ty::pcpp"clang::Type", V::pcpp"llvm::Value") = ccall((:AssociateValue,libcxxffi),Void,(Ptr{Void},Ptr{Void},Ptr{Void}),d,ty,V)

irbuilder() = pcpp"clang::CodeGen::CGBuilderTy"(ccall((:clang_get_builder,libcxxffi),Ptr{Void},()))

# # # Staging

# Main implementation.
# The two staged functions cppcall and cppcall_member below just change
# the thiscall flag depending on which ones is called.

setup_cpp_env(f::pcpp"llvm::Function") = ccall((:setup_cpp_env,libcxxffi),Ptr{Void},(Ptr{Void},),f)
cleanup_cpp_env(state) = ccall((:cleanup_cpp_env,libcxxffi),Void,(Ptr{Void},),state)

dump(d::pcpp"clang::Decl") = ccall((:cdump,libcxxffi),Void,(Ptr{Void},),d)
dump(d::pcpp"clang::FunctionDecl") = ccall((:cdump,libcxxffi),Void,(Ptr{Void},),d)
dump(expr::pcpp"clang::Expr") = ccall((:exprdump,libcxxffi),Void,(Ptr{Void},),expr)
dump(t::pcpp"clang::Type") = ccall((:typedump,libcxxffi),Void,(Ptr{Void},),t)

function _cppcall(expr, thiscall, argt)
    check_args(argt, expr)

    pvds = pcpp"clang::ParmVarDecl"[]

    rslot = C_NULL
    fname = string(expr.parameters[1])

    if thiscall
        # Create an expression to refer to the `this` parameter.
        ct = cpptype(argt[1])
        if argt[1] <: CppValue
            # We want the this argument to be modified in place,
            # so we make a deref of the pointer to the
            # data.
            pct = pointerTo(ct)
            pvd = CreateParmVarDecl(pct)
            push!(pvds, pvd)
            dre = createDerefExpr(CreateDeclRefExpr(pvd,pct))
        else
            pvd = CreateParmVarDecl(ct)
            push!(pvds, pvd)
            dre = CreateDeclRefExpr(pvd,ct)
        end
        me = BuildMemberReference(dre, ct, argt[1] <: CppPtr, fname)

        if me == C_NULL
            error("Could not find member function $fname")
        end

        # And now the expressions for all the arguments
        callargs, callpvds = buildargexprs(argt[2:end])
        append!(pvds, callpvds)
        # The actual call expression
        mce = BuildCallToMemberFunction(me,callargs)

        # At this point we're done with creating the clang AST,
        # and can move on to llvm code generation

        # First we need to get the return type of the C++ expression
        rt = getCalleeReturnType(pcpp"clang::CallExpr"(mce.ptr))
        rett = juliatype(rt)
        llvmrt = julia_to_llvm(rett)
        llvmargt = argt

        @show rett
        @show llvmrt

        # Let's create an LLVM fuction
        f = CreateFunction(llvmrt,
            map(julia_to_llvm,[argt...]))

        # Clang's code emitter needs some extra information about the function, so let's
        # initialize that as well
        state = setup_cpp_env(f)

        builder = irbuilder()

        # First compute the llvm arguments (unpacking them from their julia wrappers),
        # then associate them with the clang level variables
        args = llvmargs(builder, f, argt)
        associateargs(builder,argt,args,pvds)

        ret = pcpp"llvm::Value"(ccall((:emitcppmembercallexpr,libcxxffi),Ptr{Void},(Ptr{Void},Ptr{Void}),mce.ptr,rslot))

        #if jlrslot != C_NULL && rslot != C_NULL && (length(args) == 0 || rslot != args[1])
        #    @cpp1 builder->CreateRet(pcpp"llvm::Value"(jlrslot.ptr))
        #else

        #return createReturn(builder,f,argt,llvmrt,rt,ret,state)
    else
        d = lookup_ctx(fname)
        @assert d.ptr != C_NULL
        # If this is a typedef or something we'll try to get the primary one
        primary_decl = to_decl(primary_ctx(toctx(d)))
        if primary_decl != C_NULL
            d = primary_decl
        end
        # Let's see if we're constructing something. And, if so, let's
        # setup an array to accept the result
        cxxd = dcastCXXRecordDecl(d)
        cxxt = cxxtmplt(d)
        if cxxd != C_NULL || cxxt != C_NULL
            targs = expr.parameters[2]
            T = CppValue{symbol(fname),tuple(targs...)}
            if cxxd == C_NULL
                cxxd = specialize_template(cxxt,targs,cpptype)
            end
            juliart = T

            # The arguments to llvmcall will have an extra
            # argument for the return slot
            llvmargt = [Ptr{Uint8},argt...]

            # And now the expressions for all the arguments
            callargs, callpvds = buildargexprs(argt)
            append!(pvds, callpvds)

            # Let's create an LLVM fuction
            f = CreateFunction(julia_to_llvm(Void),
                map(julia_to_llvm,map(stripmodifier,llvmargt)))

            state = setup_cpp_env(f)
            builder = irbuilder()

            # First compute the llvm arguments (unpacking them from their julia wrappers),
            # then associate them with the clang level variables
            args = llvmargs(builder, f, llvmargt)
            associateargs(builder,argt,args[2:end],pvds)

            ctce = BuildCXXTypeConstructExpr(typeForDecl(cxxd),callargs)

            ccall((:emitexprtomem,libcxxffi),Void,(Ptr{Void},Ptr{Void},Cint),ctce,args[1],1)
            ret = Void

            CreateRetVoid(builder)

            cleanup_cpp_env(state)

            arguments = [:(pointer(r.data)), [:(args[$i]) for i = 1:length(argt)]]

            size = cxxsizeof(cxxd)
            rr = Expr(:block,
                :( r = ($(T))(Array(Uint8,$size)) ),
                Expr(:call,:llvmcall,f.ptr,Void,tuple(llvmargt...),arguments...),
                :r)
            return rr
        else
            myctx = getContext(d)
            dne = BuildDeclarationNameExpr(split(fname,"::")[end],myctx)

            # And now the expressions for all the arguments
            callargs, callpvds = buildargexprs(argt)
            append!(pvds, callpvds)

            ce = CreateCallExpr(dne,callargs)

            # First we need to get the return type of the C++ expression
            rt = DeduceReturnType(ce)
            rett = juliatype(rt)

            llvmargt = [argt...]

            issret = (rett != None) && rett <: CppValue

            if issret
                llvmargt = [Ptr{Uint8},llvmargt]
            end

            llvmrt = julia_to_llvm(rett)

            # Let's create an LLVM fuction
            f = CreateFunction(issret ? julia_to_llvm(Void) : llvmrt,
                map(julia_to_llvm,llvmargt))

            # Clang's code emitter needs some extra information about the function, so let's
            # initialize that as well
            state = setup_cpp_env(f)

            builder = irbuilder()

            args = llvmargs(builder, f, llvmargt)
            associateargs(builder,argt,issret ? args[2:end] : args,pvds)

            if issret
                rslot = CreateBitCast(builder,args[1],getPointerTo(toLLVM(rt)))
            end

            ret = pcpp"llvm::Value"(ccall((:emitcallexpr,libcxxffi),Ptr{Void},(Ptr{Void},Ptr{Void}),ce,rslot))

            if rett <: CppValue
                ret = C_NULL
            end
            #return createReturn(builder,f,argt,llvmrt,rt,ret,state)
        end
    end

    # Common return path for everything that's calling a normal function
    # (i.e. everything but constructors)
    jlrt = rett
    if ret == C_NULL
        jlrt = Void
        CreateRetVoid(builder)
    else
        #@show rett
        if rett == Void
            CreateRetVoid(builder)
        else
            if rett <: CppPtr || rett <: CppRef
                undef = getUndefValue(llvmrt)
                elty = getStructElementType(llvmrt,0)
                ret = CreateBitCast(builder,ret,elty)
                ret = pcpp"llvm::Value"(ccall((:create_insert_value,libcxxffi),Ptr{Void},(Ptr{Void},Ptr{Void},Csize_t),undef.ptr,ret.ptr,0))
            end
            CreateRet(builder,ret)
        end
    end

    cleanup_cpp_env(state)

    if (rett != None) && rett <: CppValue
        arguments = [:(pointer(r.data)), [:(args[$i]) for i = 1:length(argt)]]
        size = cxxsizeof(rt)
        return Expr(:block,
            :( r = ($(rett))(Array(Uint8,$size)) ),
            Expr(:call,:llvmcall,f.ptr,Void,tuple(llvmargt...),arguments...),
            :r)
    else
        return Expr(:call,:llvmcall,f.ptr,rett,argt,
            [:(args[$i]) for i = 1:length(argt)]...)
    end

end

stagedfunction cppcall(expr, args...)
    _cppcall(expr, false, args)
end

stagedfunction cppcall_member(expr, args...)
    _cppcall(expr, true, args)
end

# Builds a call to the cppcall staged functions that represents a
# call to a C++ function.
# Arguments:
#   - cexpr:    The (julia-side) call expression
#   - this:     For a member call the expression corresponding to
#               object of which the function is being called
#   - prefix:   For a call namespaced call, all namespace qualifiers
#
# E.g.:
#   - @cpp foo::bar::baz(a,b,c)
#       - cexpr == :( baz(a,bc) )
#       - this === nothing
#       - prefix == "foo::bar"
#   - @cpp m->DoSomething(a,b,c)
#       - cexpr == :( DoSomething(a,b,c) )
#       - this == :( m )
#       - prefix = ""
#
function build_cpp_call(cexpr, this, prefix = "")
    @assert isexpr(cexpr,:call)
    targs = []

    # Turn prefix and call expression, back into a fully qualified name
    # (and optionally type arguments)
    if isexpr(cexpr.args[1],:curly)
        fname = string(prefix, cexpr.args[1].args[1])
        targs = map(macroexpand, copy(cexpr.args[1].args[2:end]))
    else
        fname = string(prefix, cexpr.args[1])
        targs = ()
    end

    arguments = cexpr.args[2:end]

    # Unary * is treated as a deref
    for (i, arg) in enumerate(arguments)
        if isexpr(arg,:call)
            # is unary *
            if length(arg.args) == 2 && arg.args[1] == :*
                arguments[i] = :( CppDeref($arg) )
            end
        end
    end

    # Add `this` as the first argument
    this !== nothing && unshift!(arguments, this)

    # The actual call to the staged function
    ret = Expr(:call, this === nothing ? :cppcall : :cppcall_member,
        :(CppExpr{$(quot(symbol(fname))),$(quot(targs))}()), arguments...)

    ret
end

function to_prefix(expr)
    if isa(expr,Symbol)
        return string(expr)
    elseif isexpr(expr,:(::))
        return string(to_prefix(expr.args[1]),"::",to_prefix(expr.args[2]))
    end
    error("Invalid NSS $expr")
end

function cpps_impl(expr,prefix="")
    # Expands a->b to
    # llvmcall((cpp_member,typeof(a),"b"))
    if isa(expr,Symbol)
        error("Unfinished")
        return cpp_ref(expr,prefix,uuid)
    elseif expr.head == :(->)
        a = expr.args[1]
        b = expr.args[2]
        i = 1
        while !(isexpr(b,:call) || isa(b,Symbol))
            b = expr.args[2].args[i]
            if !(isexpr(b,:call) || isexpr(b,:line) || isa(b,Symbol))
                error("Malformed C++ call. Expected member not $b")
            end
            i += 1
        end
        if isexpr(b,:call)
            return build_cpp_call(b,a)
        else
            error("Unsupported!")
        end
    elseif isexpr(expr,:(=))
        error("Unimplemented")
    elseif isexpr(expr,:(::))
        return cpps_impl(expr.args[2],string(to_prefix(expr.args[1]),"::"))
    elseif isexpr(expr,:call)
        return build_cpp_call(expr,nothing,prefix)
    end
    error("Unrecognized CPP Expression ",expr," (",expr.head,")")
end

macro cxx(expr)
    cpps_impl(expr)
end


text:
struct SparseValueIterationSolver <: Solver
    max_iterations::Int64
    belres::Float64 # the Bellman Residual
    include_Q::Bool
    verbose::Bool
    init_util::Vector{Float64}
end

function SparseValueIterationSolver(;max_iterations=500,
                                    belres::Float64=1e-3,
                                    include_Q::Bool=true,
                                    verbose::Bool=false,
                                    init_util::Vector{Float64}=Vector{Float64}(undef, 0))
    return SparseValueIterationSolver(max_iterations, belres, include_Q, verbose, init_util)
end

@POMDP_require solve(solver::SparseValueIterationSolver, mdp::Union{MDP,POMDP}) begin
    P = typeof(mdp)
    S = statetype(P)
    A = actiontype(P)
    @req discount(::P)
    @req n_states(::P)
    @req n_actions(::P)
    @subreq ordered_states(mdp)
    @subreq ordered_actions(mdp)
    @req transition(::P,::S,::A)
    @req reward(::P,::S,::A,::S)
    @req stateindex(::P,::S)
    @req actionindex(::P, ::A)
    @req actions(::P, ::S)
    as = actions(mdp)
    ss = states(mdp)
    a = first(as)
    s = first(ss)
    dist = transition(mdp, s, a)
    D = typeof(dist)
    @req support(::D)
    @req pdf(::D,::S)
end

function qvalue!(m::Union{MDP,POMDP}, transition_A_S_S2, reward_S_A::AbstractMatrix{F}, value_S::AbstractVector{F}, out_qvals_S_A) where {F}
    @assert size(out_qvals_S_A) == (n_states(m), n_actions(m))
    for a in 1:n_actions(m)
        out_qvals_S_A[:, a] = view(reward_S_A, :, a) + discount(m) * transition_A_S_S2[a] * value_S
    end
end

function transition_matrix_a_s_sp(mdp::MDP)
    # Thanks to zach
    na = n_actions(mdp)
    ns = n_states(mdp)
    transmat_row_A = [Float64[] for _ in 1:n_actions(mdp)]
    transmat_col_A = [Float64[] for _ in 1:n_actions(mdp)]
    transmat_data_A = [Float64[] for _ in 1:n_actions(mdp)]

    for s in states(mdp)
        si = stateindex(mdp, s)
        for a in actions(mdp, s)
            ai = actionindex(mdp, a)
            if !isterminal(mdp, s, a) # if terminal, the transition probabilities are all just zero
                td = transition(mdp, s, a)
                for (sp, p) in weighted_iterator(td)
                    if p > 0.0
                        spi = stateindex(mdp, sp)
                        push!(transmat_row_A[ai], si)
                        push!(transmat_col_A[ai], spi)
                        push!(transmat_data_A[ai], p)
                    end
                end
            end
        end
    end
    transmats_A_S_S2 = [sparse(transmat_row_A[a], transmat_col_A[a], transmat_data_A[a], n_states(mdp), n_states(mdp)) for a in 1:n_actions(mdp)]
    # Note: not valid for terminal states
    # @assert all(all(sum(transmats_A_S_S2[a], dims=2) .≈ ones(n_states(mdp))) for a in 1:n_actions(mdp)) "Transition probabilities must sum to 1"
    return transmats_A_S_S2
end

function reward_s_a(mdp::MDP)
    reward_S_A = fill(-Inf, (n_states(mdp), n_actions(mdp))) # set reward for all actions to -Inf unless they are in actions(mdp, s)
    for s in states(mdp)
        for a in actions(mdp, s)
            td = transition(mdp, s, a)
            r = 0.0
            for (sp, p) in weighted_iterator(td)
                if p > 0.0
                    r += p*reward(mdp, s, a, sp)
                end
            end
            reward_S_A[stateindex(mdp, s), actionindex(mdp, a)] = r
        end
    end
    return reward_S_A
end

function solve(solver::SparseValueIterationSolver, mdp::Union{MDP, POMDP})
    nS = n_states(mdp)
    nA = n_actions(mdp)
    if isempty(solver.init_util)
        v_S = zeros(nS)
    else
        @assert length(solver.init_util) == nS "Input utility dimension mismatch"
        v_S = solver.init_util
    end

    transition_A_S_S2 = transition_matrix_a_s_sp(mdp)
    reward_S_A = reward_s_a(mdp)
    qvals_S_A = zeros(nS, nA)
    maxchanges_T = zeros(solver.max_iterations)

    total_time = 0.0
    for i in 1:solver.max_iterations
        iter_time = @elapsed begin
            qvalue!(mdp, transition_A_S_S2, reward_S_A, v_S, qvals_S_A)
            new_v_S = dropdims(maximum(qvals_S_A, dims=2), dims=2)
            @assert size(v_S) == size(new_v_S)
            maxchanges_T[i] = maximum(abs.(new_v_S .- v_S))
            v_S = new_v_S
        end
        total_time += iter_time
        if solver.verbose
            @info "residual: $(maxchanges_T[i]), time: $(iter_time), total time: $(total_time) " i
        end
        maxchanges_T[i] < solver.belres ? break : nothing
    end
    qvalue!(mdp, transition_A_S_S2, reward_S_A, v_S, qvals_S_A)
    # Rounding to avoid floating point error noise
    policy_S = dropdims(getindex.(argmax(round.(qvals_S_A, digits=20), dims=2), 2), dims=2)
    if solver.include_Q
        policy = ValueIterationPolicy(mdp, qvals_S_A, v_S, policy_S)
    else
        policy = ValueIterationPolicy(mdp, utility=v_S, policy=policy_S, include_Q=false)
    end
    return policy
end


text:
using Jieba
using Test

@testset "Jieba.jl" begin
    # 分词
    @time cut("我是真的大帅哥") |> println
    @time cut_for_search("卡梅隆安东尼戴维斯") |> println

    # 特征提取
    @time analyse.extract_tags("提取TF-IDF特征向量") |> println
    @time analyse.textrank("提取TextRank特征向量") |> println

    # 词性标注
    @time pseg.lcut("我爱北京天安门") |> println
end

max_stars_repo_name: lotcher/Jieba

max_issues_repo_name: lotcher/Jieba

max_forks_repo_name: lotcher/Jieba


text:
function eplotsetup(eroot::AbstractDict)

    pdir  = eroot["plot"]
    ppys  = joinpath(eroot["plot"],"eplot-lsm.py")

    fID = open(ppys,"w");

    write(fID,"#!/usr/bin/env python\n");
    write(fID,"import cdsapi\n");
    write(fID,"c = cdsapi.Client()\n\n");

    write(fID,"c.retrieve(\"reanalysis-era5-single-levels\",\n");
    write(fID,"    {\n");
    write(fID,"        \"product_type\": \"reanalysis\",\n");
    write(fID,"        \"variable\": \"land_sea_mask\",\n");
    write(fID,"        \"year\": \"$(year(Dates.now()))\",\n");
    write(fID,"        \"month\": \"01\",\n");
    write(fID,"        \"day\": \"01\",\n");
    write(fID,"        \"time\": \"00:00\",\n");
    write(fID,"        \"format\": \"netcdf\",\n");
    write(fID,"    },\n");
    write(fID,"    \"eplot-GLBx0.25-lsm.nc\")\n\n");
    write(fID,"})\n\n");

end

function eplotsubregion(ID::AbstractString,pdir::AbstractString)

    fnc = joinpath(pdir,"eplot-GLBx0.25-lsm.nc");
    if !isfile(fnc)
        error("$(Dates.now()) - The Global Land-Sea Mask file \"eplot-GLBx0.25-lsm.nc\" does not exist in the folder $(eroot["plot"]). Please run the script \"eplot-lsm.py\" in the folder to download the ERA5 Land-Sea Mask file.")
    end

    ds = NCDataset(fnc)
    lsm = ds["lsm"][:]*1; lon = ds["longitude"][:]*1; lat = ds["latitude"][:]*1
    close(ds)

    rlon,rlat,rinfo = gregiongridvec(ID,lon,lat)
    rlsm = regionextractgrid(lsm,rinfo);
    if ndims(rlsm) == 3; rlsm = dropdims(rlsm,dims=3); end
    rfnc = joinpath(pdir,"eplot-$(ID)x0.25-lsm.nc")
    if isfile(rfnc)
        @info "$(Dates.now()) - Stale NetCDF file $(rfnc) detected.  Overwriting ..."
        rm(rfnc);
    end

    ds = NCDataset(rfnc,"c",attrib = Dict(
        "Conventions" => "CF-1.6",
        "history"     => "Created on $(Dates.now())"
    ))

    ds.dim["longitude"] = length(rlon)
    ds.dim["latitude"]  = length(rlat)

    nclon = defVar(ds,"longitude",Float32,("longitude",),attrib = Dict(
        "units"     => "degrees_east",
        "long_name" => "longitude",
    ))

    nclat = defVar(ds,"latitude",Float32,("latitude",),attrib = Dict(
        "units"     => "degrees_north",
        "long_name" => "latitude",
    ))

    nclsm = defVar(ds,"lsm",Float32,("longitude","latitude"),attrib = Dict(
        "units"         => "(0 - 1)",
        "long_name"     => "Land-sea mask",
        "standard_name" => "land_binary_mask",
    ))

    nclon[:] = rlon; nclat[:] = rlat; nclsm[:] = rlsm

    close(ds)

end

function eplotlsm(ID::AbstractString,pdir::AbstractString;overwrite=false)

    fnc = joinpath(pdir,"eplot-$(ID)x0.25-lsm.nc")
    if !isfile(fnc) || overwrite; eplotsubregion(ID,pdir) end
    ds  = NCDataset(fnc);
    lsm = ds["lsm"][:]*1; lon = ds["longitude"][:]*1; lat = ds["latitude"][:]*1
    close(ds)

    return lsm,lon,lat

end

function eplotlsmgrid(
    gridbounds::Array, pdir::AbstractString;
    pID::AbstractString
)

    fnc = joinpath(pdir,"eplot-$(pID)x0.25-lsm.nc"); ds = NCDataset(fnc);
    lsm = ds["lsm"][:]*1; lon = ds["longitude"][:]*1; lat = ds["latitude"][:]*1
    close(ds)

    rlon,rlat,rinfo = regiongridvec(gridbounds,lon,lat)
    rlsm = regionextractgrid(lsm,rinfo);

    return rlsm,rlon,rlat

end


text:
# This file was generated, do not modify it. # hide
clf = machine(KNNClassifier(K=3), Xs, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)

accuracy(ŷ, y[test]) |> r3


text:
# Adapted from:
# SOSDEMO4 --- Matrix Copositivity
# Section 3.4 of SOSTOOLS User's Manual

@testset "SOSDEMO4 with $(factory.constructor)" for factory in sdp_factories
    @polyvar x[1:5]

    # The matrix under consideration
    J = [1 -1  1  1 -1;
        -1  1 -1  1  1;
         1 -1  1 -1  1;
         1  1 -1  1 -1;
        -1  1  1 -1  1]

    xs = vec(x).^2
    xsJxs = dot(xs, J*xs)
    r = sum(xs)

    m0 = SOSModel(factory)
    @constraint m0 xsJxs >= 0
    JuMP.optimize!(m0)
    @test JuMP.dual_status(m0) == MOI.InfeasibilityCertificate

    m1 = SOSModel(factory)
    @constraint m1 r*xsJxs >= 0
    JuMP.optimize!(m1)
    @test JuMP.primal_status(m1) == MOI.FeasiblePoint
    # Program is feasible. The matrix J is copositive
end