text:
#!/bin/bash
#=
exec julia --project=. -e 'include(popfirst!(ARGS))' \
    "${BASH_SOURCE[0]}" "$@"
=#

using MagneticResonanceSignals
using Statistics
using Unitful
using AxisArrays

# TODO: Move the signal processing parts (at least!) into src!
function twix_to_felix(twixname, felixname; repair=false)
    twix = load_twix(twixname)
    spectro = mr_load(twix; repair=repair)
    @info "Opened data file" twixname spectro

    signal = simple_averaging(spectro, downsample=2)

    # Extract some metadata required by the felix format
    frequency = standard_metadata(spectro).frequency
    f1_bandwidth = 1.0/step(AxisArrays.axes(signal, Axis{:time1}).val)
    f2_bandwidth = 1.0/step(AxisArrays.axes(signal, Axis{:time2}).val)

    save_felix(felixname, signal; bandwidth=(f2_bandwidth, f1_bandwidth), frequency=frequency)

    @info "Wrote to Felix file $felixname"
    signal
end

using ArgParse

argdef = ArgParseSettings(exc_handler=isinteractive() ? ArgParse.debug_handler : ArgParse.default_handler)
@add_arg_table argdef begin
    "--repair"
        help = "Attempt to repair a partially complete LCOSY sequence"
        action = :store_true
    "in_file"
        help = "path to the input twix file (Siemens meas_*.dat RAID format)"
        required = true
    "out_file"
        help = "path to the output file (Felix .dat format)."
        required = true
end

function process_args(argvec)
    parsed_args = parse_args(argvec, argdef)

    in_file = parsed_args["in_file"]
    out_file = parsed_args["out_file"]

    (in_file=in_file, out_file=out_file, repair=parsed_args["repair"])
end

args = process_args(ARGS)
twix_to_felix(args.in_file, args.out_file, repair=args.repair)

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
using Test

using Dice
using Dice: CuddMgr, DistBool

@testset "Cudd Backend Tests" begin

    mgr = CuddMgr()

    x = DistBool(mgr, 0.1)

    @test !!x == x
    @test Dice.isliteral(x)
    @test Dice.isliteral(!x)

    @test Dice.isposliteral(x)
    @test Dice.isnegliteral(!x)

    @test !Dice.isposliteral(!x)
    @test !Dice.isnegliteral(x)

    @test !Dice.Cudd_IsComplement(x.bit)
    @test Dice.Cudd_IsComplement((!x).bit)

    @test Dice.Cudd_Regular(x.bit) === x.bit
    @test Dice.Cudd_Regular((!x).bit) === x.bit

    # @test infer(x) ≈ 0.1
    # @test infer(!x) ≈ 0.9

end

ext: jl


text:
part1(filename::AbstractString) = length(visithouses(readchomp(filename)))

function visithouses(directions::AbstractString)
    x = 0
    y = 0
    visited = Set(((0, 0),))
    for dir in directions
        if dir == '<'
            x -= 1
        elseif dir == '>'
            x += 1
        elseif dir == '^'
            y -= 1
        elseif dir == 'v'
            y += 1
        end
        push!(visited, (x, y))
    end
    visited
end

function part2(filename::AbstractString)
    directions = readchomp(filename)
    human = directions[1:2:end]
    robot = directions[2:2:end]
    length(union(visithouses(human), visithouses(robot)))
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
#####
##### Global settings
#####
plot_opt = Dict(:titlefontsize=>28,
        :xticks=>false, :yticks=>false, :zticks=>false,
        :legend=>false, :fontfamily=>:serif)

################################################################################
#####
##### Define a function for library examples
#####

function plot_library(this_dat_name, this_fig_name, system_name="";
                        plot_opt=plot_opt,
                        plot_ind=501:1000,
                        plot_coordinate=1,
                        var_names=["X","Y","Z"])
        ### LOAD
        fname = this_dat_name*"uncontrolled.bson"
        @load fname dat_raw #numerical_grad_raw
        fname = this_dat_name*"ode_vars.bson";
        @load fname dat U_true
        fname = this_dat_name*"naive_vars_bayes.bson";
        @load fname noise_guess residual accepted_ind #ctr_guess
        fname = this_dat_name*"ctr_vars_sindy.bson";
        @load fname dat_ctr ctr_guess2 #sindy_grad_ctr

        ### PLOTS
        ## 1: Attractor view (uncontrolled)
        if size(dat_raw,1) == 3
                plot_raw = plot3d(dat_raw[1, :], dat_raw[2, :], dat_raw[3, :],
                        color=COLOR_DICT["data_uncontrolled"], lw=4,
                        #zlabel=system_name, guidefontsize=24, mirror=true
                        ;plot_opt...);
        elseif size(dat_raw,1) == 2
                plot_raw = plot(dat_raw[1, :], dat_raw[2, :],
                        color=COLOR_DICT["data_uncontrolled"], lw=4,
                        ;plot_opt...);
        end
        title!(system_name)

        ## 2: 1d time series (controlled)
        ind = plot_ind;
        plot_1d = plot(dat[plot_coordinate, ind],
                color=COLOR_DICT["data"], lw=4; plot_opt...)
                ylabel!("$(var_names[plot_coordinate])", guidefontsize=24)
        plot_ctr_true = plot(U_true[plot_coordinate,ind],
                color=COLOR_DICT["true"], lw=3;
                plot_opt...);

        ## 3: Residual with noise envelope (x coordinate)
        plot_residual1 = plot(residual[plot_coordinate,ind],
                            ribbon=mean(noise_guess),
                            fillalpha=0.5,
                            color=COLOR_DICT["residual"], lw=4;
                            plot_opt...)
            ylabel!("d$(var_names[plot_coordinate])/dt", guidefontsize=24);

        ## 4: Control signal guess (x coordinate)
        # Same ind as above
        plot_ctr_guess = plot(ctr_guess2[plot_coordinate,ind],
                        color=COLOR_DICT["control_time"], lw=6;
                        plot_opt...)

        ## 5: "Revealed underlying system"
        controlled_ind = 1:5000;
            ind = controlled_ind;
        let d = dat_ctr
            if size(dat_ctr,1) == 3
                    global plot_reconstruction = plot3d(
                        d[1, ind], d[2, ind], d[3, ind],
                        color=COLOR_DICT["model_controlled"], alpha=0.8, lw=3;
                        plot_opt...);
            elseif size(dat_ctr,1) == 2
                    global plot_reconstruction = plot(
                        d[1, ind], d[2, ind],
                        color=COLOR_DICT["model_controlled"], alpha=0.8, lw=3;
                        plot_opt...);
            end
        end;

        ##
        ## Put it all together
        ##
        lay = @layout [a b c [d; e] f]

        plot_final = plot(
                plot_raw, plot_1d, #plot_ctr_true,
                plot_residual1, plot_ctr_guess, plot_ctr_true, plot_reconstruction,
                layout=lay)
        plot!(size=(2000, 250))

        return plot_final
end

################################################################################
#####
##### Define a function for noise studies
#####

function plot_library_noise(this_dat_name, system_name="";
                        plot_opt=plot_opt)

        plot_opt[:xticks] = true
        plot_opt[:yticks] = true
        plot_opt[:tickfontsize] = 16
        # plot_opt[:colorbar] = true
        plot_opt[:color] = :Spectral#:RdYlBu

        ## Load
        this_dat_name = DAT_FOLDERNAME*this_dat_name;

        fname = this_dat_name*"metadata.bson"
        @load fname noise_vals control_signal_vals coef_norm

        fname = this_dat_name*"coefficients.bson"
        @load fname all_naive_err all_err

        fname = this_dat_name*"derivatives.bson";
        @load fname all_err_deriv_subsample all_naive_err_deriv
        # fname = this_dat_name*"coefficients.bson"
        # @load fname noise_vals vec_naive std_naive vec_err std_err
        #
        # fname = this_dat_name*"derivatives.bson";
        # @load fname vec_deriv std_deriv vec_naive_deriv std_naive_deriv
        # vec_naive, std_naive = mean_and_std(all_naive_err./coef_norm, 2)
        # vec_err, std_err = mean_and_std(all_err./coef_norm, 2)
        # vec_deriv, std_deriv = mean_and_std(all_err_deriv_subsample, 2)
        # vec_naive_deriv, std_naive_deriv = mean_and_std(all_naive_err_deriv, 2)

        ## Process
        # heatmap_coef_naive = mean(all_naive_err./coef_norm, dims=3)[:,:,1]
        # heatmap_coef_final = mean(all_err./coef_norm, dims=3)[:,:,1]

        heatmap_deriv_naive = mean(all_naive_err_deriv, dims=3)[:,:,1]
        heatmap_deriv_final = mean(all_err_deriv_subsample, dims=3)[:,:,1]
        clims = (0, 2*median(heatmap_deriv_naive))
        plot_opt[:clims] = clims

        plot_deriv = heatmap(control_signal_vals, noise_vals,
                heatmap_deriv_final; plot_opt...)
            xlabel!("Percent Data Perturbed", guidefontsize=20)
            xticks!(control_signal_vals)
            ylabel!("Noise", guidefontsize=20)
            title!(system_name*" Final Error")

        plot_opt[:yticks] = false
        plot_opt[:colorbar] = false
        plot_improve = heatmap(control_signal_vals, noise_vals,
                heatmap_deriv_naive .- heatmap_deriv_final; plot_opt...)
            xlabel!("Percent Data Perturbed", guidefontsize=20)
            xticks!(control_signal_vals)
            # ylabel!("Noise")
            title!("Improvement over SINDy")

        ## First panel: Error in coefficients
        # coef_plot = plot(noise_vals, vec_naive, ribbon=std_naive,
        #                 label="Intial SINDy", lw=3; plot_opt...)
        #     plot!(noise_vals, vec_err, ribbon=std_err, lw=3,
        #                 label="Final iteration")
        #     xlabel!("Noise", guidefontsize=24)
        #     ylabel!("Error", guidefontsize=20)
        #     # title!("Fractional Error in Coefficients")
        #     title!(system_name*" (coefficients)")

        ## Second panel: Error in derivatives
        # deriv_plot = plot(noise_vals, vec_naive_deriv, ribbon=std_naive_deriv,
        #                 label="Intial SINDy", lw=3; plot_opt...)
        #     plot!(noise_vals, vec_deriv, ribbon=std_deriv, lw=3,
        #                 label="Final iteration")
        #     xlabel!("Noise", guidefontsize=24)
        #     ylabel!("Error", guidefontsize=20)
        #     # title!("Error in Model Derivatives")
        #     title!(system_name*" (derivatives)")

        ##
        ## Put it all together
        ##
        lay = @layout [a b]

        # plot_final = plot(coef_plot, deriv_plot, layout=lay)
        plot_final = plot(plot_deriv, plot_improve, layout=lay)
        plot!(size=(2000, 250))

        ## Save
        # savefig(plot_final, this_dat_name * ".png")

        return plot_final
end

##
export plot_library, plot_library_noise


text:
"""
   vbmf_parameters

Compound type for vbmf computation. Contains the following fields:\n

    L::Int - number of rows of the original matrix
    M::Int - number of columns of the original matrix
    H::Int - internal product dimension
    H1::Int - number of columns of B that should belong to a contaminated class
    labels::Array{Int64, 1} - which columns of Y are labeled as non-contaminated
    AHat::Array{Float64, 2} - mean value of A, size (M, H)
    BHat::Array{Float64, 2} - mean value of B, size (L, H)
    SigmaA::Array{Float64, 2} - covariance of A, size (H, H)
    SigmaB::Array{Float64, 2} - covariance of B, size (H, H)
    CA::Array{Float64, 2} - prior covariance of A, size (H, H)
    CB::Array{Float64, 2} - prior covariance of B, size (H, H)
    invCA::Array{Float64, 2} - inverse of cA
    invCB::Array{Float64, 2} - inverse of cB
    sigma::Float64 - variance of data
    YHat::Array{Float64, 2} - estimate of Y, size (L, M)
"""
type vbmf_parameters
    L::Int
    M::Int
    H::Int
    H1::Int
    labels::Array{Int64,1}
    AHat::Array{Float64, 2}
    BHat::Array{Float64, 2}
    SigmaA::Array{Float64, 2}
    SigmaB::Array{Float64, 2}
    CA::Array{Float64, 2}
    CB::Array{Float64, 2}
    invCA::Array{Float64, 2}
    invCB::Array{Float64, 2}
    sigma2::Float64
    YHat::Array{Float64, 2}

    vbmf_parameters() = new()
end

"""
    vbmf_init(L::Int, M::Int, H::Int; ca::Float64, cb::Float64, sigma::Float64, H1::Int = 0,
    labels::Array{Int64,1} = Array{Int64,1}())

Returns an initialized structure of type vbmf_parameters.
"""
function vbmf_init(Y::Array{Float64,2}, H::Int; ca::Float64 = 1.0, cb::Float64 = 1.0, sigma2::Float64 = 1.0,
    H1::Int = 0, labels::Array{Int64,1} = Array{Int64,1}())
    params = vbmf_parameters()
    L, M = size(Y)

    params.L, params.M = L, M
    params.H = H
    params.H1 = H1
    params.labels = labels

    params.AHat = randn(M, H)
    # now set zeroes to where the columns of Y labeled as non-infected are
    # in A, these are rows
    params.AHat[labels, end-H1+1:end] = 0.0
    params.BHat = randn(L, H)
    params.SigmaA = zeros(H, H)
    params.SigmaB = zeros(H, H)
    params.CA = ca*eye(H)
    params.CB = cb*eye(H)
    params.invCA = inv(params.CA)
    params.invCB = inv(params.CB)
    params.sigma2 = sigma2
    params.YHat = params.BHat*params.AHat'

    return params
end

"""
   copy(params_in::vbmf_parameters)

Copy function for vbmfa_parameters.
"""
function copy(params_in::vbmf_parameters)
    params = vbmf_parameters()

    for field in fieldnames(params_in)
        setfield!(params, field, getfield(params_in, field))
    end

    return params
end

"""
    updateA!(Y::Array{Float64,2}, params::vbmf_parameters)

Updates mean and covariance of the A matrix.
"""
function updateA!(Y::Array{Float64,2}, params::vbmf_parameters)
    params.SigmaA = params.sigma2*inv(params.BHat'*params.BHat +
        params.L*params.SigmaB + params.sigma2*params.invCA)
    params.AHat =  Y'*params.BHat*params.SigmaA/params.sigma2
    # now set zeroes to where the columns of Y labeled as non-infected are
    # in A, these are rows
    params.AHat[params.labels, end-params.H1+1:end] = 0.0
end

"""
    updateB!(Y::Array{Float64,2}, params::vbmf_parameters)

Updates mean and covariance of the B matrix.
"""
function updateB!(Y::Array{Float64,2}, params::vbmf_parameters)
    params.SigmaB = params.sigma2*inv(params.AHat'*params.AHat +
        params.M*params.SigmaA + params.sigma2*params.invCB)
    params.BHat =  Y*params.AHat*params.SigmaB/params.sigma2
end

"""
    updateYHat!(params::vbmf_parameters)

Updates estimate of Y.
"""
function updateYHat!(params::vbmf_parameters)
    params.YHat = params.BHat*params.AHat'
end

"""
    updateCA!(params::vbmf_parameters)

Updates the estimate of CA.
"""
function updateCA!(params::vbmf_parameters)
    for h in 1:params.H
        params.CA[h,h] = norm2(params.AHat[:,h])/params.M + params.SigmaA[h,h]
    end
    params.invCA = inv(params.CA)
end

"""
    updateCB!(params::vbmf_parameters)

Updates the estimate of CB.
"""
function updateCB!(params::vbmf_parameters)
    for h in 1:params.H
        params.CB[h,h] = norm2(params.BHat[:,h])/params.L + params.SigmaB[h,h]
    end
    params.invCB = inv(params.CB)
end

"""
    updateSigma2!(Y::Array{Float64,2}, params::vbmf_parameters)

Updates estimate of sigma^2, the measurement noise.
"""
function updateSigma2!(Y::Array{Float64,2}, params::vbmf_parameters)
    params.sigma2 = (norm2(Y) - trace(2*Y'*params.BHat*params.AHat') +
        trace((params.AHat'*params.AHat + params.M*params.SigmaA)*
            (params.BHat'*params.BHat + params.L*params.SigmaB)))/(params.L*params.M)
end

"""
    vbmf!(Y::Array{Float64, 2}, params::vbmf_parameters, niter::Int, eps::Float64 = 1e-6, est_covs = false,
    est_var = false, logdir = "", desc = "", verb = false)

Computes variational bayes matrix factorization of Y = AB' + E. Independence of A and B is assumed.
Estimation of prior covariance cA and cB and of variance sigma can be turned on and off. Estimates
of cA and cB are empirical.
The prior model is following:

    p(Y|A,B) = N(Y|BA^T, sigma^2*I)
    p(A) = N(A|0, C_A), C_A = diag(c_a)
    p(B) = N(B|0, C_B), C_B = diag(c_b)

The params argument with initialized data is modified and contains the resulting estimates after the
algorithm stops.
"""
function vbmf!(Y::Array{Float64, 2}, params::vbmf_parameters, niter::Int; eps::Float64 = 1e-6, est_covs::Bool = false,
    est_var::Bool = false, logdir = "", desc = "", verb = false)
    priors = Dict()

    # create the log dictionary
    log = false
    if logdir !=""
        log = true
        logVar = create_log(params)
    end

    # choice of convergence control variable
    convergence_var = :BHat
    old = getfield(params, convergence_var)
    d = eps + 1.0 # delta
    i = 1

    # run the loop for a given number of iterations
    while ((i <= niter) && (d > eps))
        updateA!(Y, params)
        updateB!(Y, params)

        if est_covs
            updateCA!(params)
            updateCB!(params)
        end

        if est_var
            updateSigma2!(Y, params)
        end

        if log
            update_log!(logVar, params)
        end

        # check convergence
        d = delta(getfield(params, convergence_var), old)
        old = getfield(params, convergence_var)
        i += 1
    end

    # finally, compute the estimate of Y
    updateYHat!(params)

    # convergence info
    if verb
        print("Factorization finished after ", i-1, " iterations, eps = ", d, "\n")
    end

    # save inputs and outputs
    if log
        println("Saving outputs and inputs under ", logdir, "/")
        save_log(logVar, Y, priors, logdir, desc = desc)
    end

    return params
end

"""
    vbmf(Y::Array{Float64, 2}, params_in::vbmf_parameters, niter::Int, eps::Float64 = 1e-6, est_covs = false, est_var = false)

Calls vbmf_sparse!() but copies the params_in argument so that it is not modified and can be reused.
"""
function vbmf(Y::Array{Float64, 2}, params_in::vbmf_parameters, niter::Int; eps::Float64 = 1e-6, est_covs::Bool = false,
    est_var::Bool = false, logdir = "", desc = "", verb = false)
    # copy the input params
    params = copy(params_in)

    # run the algorithm
    vbmf!(Y, params, niter, eps = eps, est_covs = est_covs, est_var = est_var,
        logdir = logdir, desc = desc, verb = verb)

    return params
end


text:
TableTraits.isiterabletable(x::Dataset) = true

function IteratorInterfaceExtensions.getiterator(source::NDSparse)
    return rows(source)
end

function IteratorInterfaceExtensions.getiterator(source::IndexedTable)
    return TableTraitsUtils.create_tableiterator(values(columns(source)), collect(keys(columns(source))))
end

function ndsparse(x; idxcols=nothing, datacols=nothing, copy=false, kwargs...)
    if TableTraits.isiterable(x)
        source_data = collect_columns(IteratorInterfaceExtensions.getiterator(x))
        source_data isa Columns{<:Pair} && return ndsparse(source_data; copy=false, kwargs...)

        # For backward compatibility
        idxcols isa AbstractArray && (idxcols = Tuple(idxcols))
        datacols isa AbstractArray && (datacols = Tuple(datacols))

        if idxcols==nothing
            n = ncols(source_data)
            idxcols = (datacols==nothing) ? Between(1, n-1) : Not(datacols)
        end
        if datacols==nothing
            datacols = Not(idxcols)
        end

        hascolumns(source_data, idxcols) || error("Unknown idxcol")
        hascolumns(source_data, datacols) || error("Unknown datacol")

        idx_storage = rows(source_data, idxcols)
        data_storage = rows(source_data, datacols)

        return convert(NDSparse, idx_storage, data_storage; copy=false, kwargs...)
    elseif idxcols==nothing && datacols==nothing
        return convert(NDSparse, x, copy = copy, kwargs...)
    else
        throw(ArgumentError("x cannot be turned into an NDSparse."))
    end
end

NDSparse(x; kwargs...) = ndsparse(x; kwargs...)

function table(rows::AbstractArray{T}; copy=false, kwargs...) where {T<:Union{Tup, Pair}}
    table(collect_columns(rows); copy=false, kwargs...)
end

function table(iter; copy=false, kw...)
    if TableTraits.isiterable(iter)
        table(collect_columns(IteratorInterfaceExtensions.getiterator(iter)); copy=copy, kw...)
    elseif Tables.istable(typeof(iter))
        table(Tables.columntable(iter); copy=copy, kw...)
    else
        throw(ArgumentError("input satisfies neither IterableTables.jl nor Tables.jl"))
    end
end


text:
using Luxor
import Dates
using ProgressMeter

macro luxoremacs(body, width=600, height=600)
    quote
        path, io = Base.Filesystem.mktemp()
        close(io)

        path * ".png"
        lfname = path * ".png"
        # local lfname = Luxor._add_ext($(esc(fname)), :png)
        Drawing($width, $height, lfname)
        origin()
        background("white")
        sethue("black")
        $(esc(body))
        finish()
        # instead of preview, I'm printing out something
        println("$lfname")
        println("#<Image: $lfname>")
        # preview()
    end
end

function test_Luxor()
    @luxoremacs begin
        text("Hello world ...")
        circle(Luxor.Point(0, 0), 200, :stroke)
    end
    @luxoremacs begin
        radius=80
        # this change the current pen
        setdash("dot")
        sethue("gray30")
        # create points for side effect
        A, B = [Luxor.Point(x, 0) for x in [-radius, radius]]
        # draw a line between two points
        line(A, B, :stroke)
        # O is origin
        circle(O, radius, :stroke)
        # labels around points
        label("A", :NW, A)
        label("O", :N,  O)
        label("B", :NE, B)
        # circle.!!
        circle.([A, O, B], 2, :fill)
        circle.([A, B], 2radius, :stroke)
    end
end

function visualize(x, y, w, h, R)
    # shift all
    x = x .- R.xmax / 2
    y = y .- R.ymax / 2

    # apply scale
    scale = 500 / max(R.xmax, R.ymax)
    x = x .* scale
    y = y .* scale
    w = w .* scale
    h = h .* scale

    @svg begin
        setdash("dash")
        # the point is the center of box
        box(Luxor.Point(0,0), R.xmax * scale, R.ymax * scale, :stroke)
        setdash("solid")
        @showprogress 0.1 "drawing .." for i in 1:length(x)
            # this x y is center
            box(Luxor.Point(x[i],y[i]), w[i], h[i], :stroke)
        end
    end
end

function test()
    @luxoremacs begin
        box(Luxor.Point(0,0), 100, 100, :stroke)
        box(Luxor.Point(50,50), 20, 20, :stroke)
    end
end

max_stars_repo_name: lihebi/bhdl

max_issues_repo_name: lihebi/bhdl

max_forks_repo_name: lihebi/bhdl


text:
# randBall
@testset "randBall" begin
    res = FPRS.randBall(2,10)
    @test typeof(res) <: Array{Float64,2}
end

# randRotationMatrix(deflection=1.0)
@testset "randRotationMatrix" begin
           res = FPRS.randRotationMatrix()
           @test typeof(res) <: Array{Float64,2}
end

max_stars_repo_name: tsipenyuk/FPRS

max_issues_repo_name: tsipenyuk/FPRS

max_forks_repo_name: tsipenyuk/FPRS


text:
using SafeTestsets

@safetestset "test1" begin include("test1.jl") end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

# The following abstract types mimics the Concept type hierarchy:
#
# Concept
#   |__ Thing
#   |     |__ Relation
#   |     |__ Entity
#   |     |__ Attribute
#   |__ Type
#         |__ ThingType
#         |     |__ RelationType
#         |     |__ EntityType
#         |     |__ AttributeType{V}
#         |__ RoleType

abstract type AbstractConcept end
abstract type AbstractConceptManager end
abstract type AbstractExplainable end

abstract type AbstractThing <: AbstractConcept end
abstract type AbstractType  <: AbstractConcept end

# Things
abstract type AbstractRelation     <: AbstractThing end
abstract type AbstractEntity       <: AbstractThing end
abstract type AbstractAttribute    <: AbstractThing end

# Types
abstract type AbstractThingType     <: AbstractType end
abstract type AbstractRoleType      <: AbstractType end

# Thing types
abstract type AbstractRelationType  <: AbstractThingType end
abstract type AbstractEntityType    <: AbstractThingType end
abstract type AbstractAttributeType <: AbstractThingType end

function instantiate(concept::Proto.Concept)
    if hasproperty(concept, :thing)
        return instantiate(concept.thing)
    else
        return instantiate(concept._type)
    end
end

"""
    Remote

Wrapper type that encapsulates a concept and transaction.
"""
struct RemoteConcept{D <: AbstractConcept}
    concept::D
    transaction::AbstractCoreTransaction
end

"""
    as_remote(x, t)

Create a `Remote`(@ref) object for a concept `x` with a
transaction `t`.
"""
function as_remote(x::D, t::AbstractCoreTransaction) where {
    D<: AbstractConcept
}
    return RemoteConcept{D}(x, t)
end


text:
# This file is part of GenericSchur.jl, released under the MIT "Expat" license

"""
    triangularize(S::Schur{T}) => Schur{complex{T})

convert a (standard-form) quasi-triangular real Schur factorization into a
triangular complex Schur factorization.
"""
function triangularize(S::Schur{Ty}) where {Ty <: Real}
    CT = complex(Ty)
    Tr = S.T
    T = CT.(Tr)
    Z = CT.(S.Z)
    n = size(T,1)
    for j=n:-1:2
        if Tr[j,j-1] != zero(CT)
            # We want a unitary similarity transform from
            # ┌   ┐      ┌     ┐
            # │a b│      │w₁  x│
            # │c a│ into │0  w₂│ where bc < 0 (a,b,c real)
            # └   ┘      └     ┘
            # If we write it as
            # ┌     ┐
            # │u  v'│
            # │-v u'│
            # └     ┘
            # and make the Ansatz that u is real (so v is imaginary),
            # we arrive at:
            # θ = atan(sqrt(-Tr[j,j-1]/Tr[j-1,j]))
            # s = sin(θ)
            # c = cos(θ)
            s = sqrt(abs(Tr[j,j-1]))
            c = sqrt(abs(Tr[j-1,j]))
            r = hypot(s,c)
            G = Givens(j-1,j,complex(c/r),-im*(s/r))
            lmul!(G,T)
            rmul!(T,G')
            rmul!(Z,G')
        end
    end
    triu!(T)
    Schur(T,Z,diag(T))
end


text:
# Overview of how the ContainerManagers work:
#
# 1. Start a TCP server on the manager using a random port within the ephemeral range
# 2. Spawn additional jobs/tasks using definition overrides to spawn identical versions of
#    Julia which will connect to the manager over TCP. The Julia function `start_worker` is
#    run which sends the worker address and port to the manager via the TCP socket.
# 3. The manager receives all of the the workers addresses and stops the TCP server.
# 4. Using the reported addresses the manager connects to each of the worker like a typical
#    cluster manager.

# Determine the start of the ephemeral port range on this system. Used in `listenany` calls.
const PORT_HINT = if Sys.islinux()
    parse(Int, first(split(readchomp("/proc/sys/net/ipv4/ip_local_port_range"), '\t')))
elseif Sys.isapple()
    parse(Int, readchomp(`sysctl -n net.inet.ip.portrange.first`))
else
    49152  # IANA dynamic and/or private port range start (https://en.wikipedia.org/wiki/Ephemeral_port)
end

abstract type ContainerManager <: ClusterManager end

"""
    launch_timeout(mgr::ContainerManager) -> Int

The maximum duration (in seconds) a manager will wait for a worker to connect since the
manager initiated the spawning of the worker.
"""
launch_timeout(::ContainerManager)

"""
    desired_workers(mgr::ContainerManager) -> Tuple{Int, Int}

The minimum and maximum number of workers wanted by the manager.
"""
desired_workers(::ContainerManager)

function Distributed.launch(
    manager::ContainerManager, params::Dict, launched::Array, c::Condition
)
    min_workers, max_workers = desired_workers(manager)
    num_workers = 0

    # Determine the IP address of the current host within the specified range
    ips = filter!(getipaddrs()) do ip
        typeof(ip) === typeof(manager.min_ip) && manager.min_ip <= ip <= manager.max_ip
    end
    valid_ip = first(ips)

    # Only listen to the single IP address which the workers attempt to connect to.
    # TODO: Ideally should be using TLS connections.
    port, server = listenany(valid_ip, PORT_HINT)
    debug(LOGGER, "Manager accepting worker connections via: $valid_ip:$port")

    listen_task = @async begin
        while isopen(server) && num_workers < max_workers
            # TODO: Potential issue with a random connection consuming a worker slot?
            sock = accept(server)
            num_workers += 1

            # The worker will report it's own address through the socket. Eventually the
            # built in Julia cluster manager code will parse the stream and record the
            # address and port.
            config = WorkerConfig()
            config.io = sock

            # Note: `launched` is treated as a queue and will have elements removed from it
            # periodically.
            push!(launched, config)
            notify(c)
        end
    end

    # Generate command which starts a Julia worker and reports its information back to the
    # manager
    #
    # Typically Julia workers are started using the hidden `julia` flags `--bind-to` and
    # `--worker`. We won't use the `--bind-to` flag as we do not know where the container
    # will be started and what ports will be available. We don't want to use
    # `--worker COOKIE` as this essentially runs `start_worker(STDOUT, COOKIE)` which
    # reports the worker address and port to STDOUT. Instead we'll run the code ourselves
    # and report the connection information back to the manager over a socket.
    exec = """
        using Distributed
        using Sockets
        sock = connect(ip\"$valid_ip\", $port)
        start_worker(sock, \"$(cluster_cookie())\")
        """
    override_cmd = `julia -e $exec`

    # Non-blocking spawn of N-containers where N is equal to `max_workers`. Workers will
    # report back to the manager via the open port we just opened.
    spawn_containers(manager, override_cmd)

    # Wait up to the timeout for workers to inform the manager of their address.
    wait(listen_task, launch_timeout(manager))

    # TODO: Does stopping listening terminate the sockets from `accept`? If so, we could
    # potentially close the socket before we know the name of the connected worker. During
    # prototyping this has not been an issue.
    close(server)  # Causes `listen_task` to complete
    notify(c)

    if num_workers < max_workers
        if num_workers >= min_workers
            warn(LOGGER, "Only managed to launch $num_workers/$max_workers workers")
        else
            error("Unable to launch the minimum number of workers")
        end
    end
end

function Distributed.manage(
    manager::ContainerManager, id::Integer, config::WorkerConfig, op::Symbol
)
    # Note: Terminating the TCP connection from the master to the worker will cause the
    # worker to shutdown automatically.
end

# Waits for the `task` to complete. Stops waiting if duration waited exceeds the `timeout`.
function Base.wait(task::Task, timeout::Period)
    timeout_secs = Dates.value(Second(timeout))
    start = time()
    while !istaskdone(task) && (time() - start) < timeout_secs
        sleep(1)
    end
end

# The Docker container ID can be determined from inside a container by inspecting
# /proc/self/cgroup. The container ID is always a 64-character hexadecimal string.
#
# https://github.com/aws/amazon-ecs-agent/issues/1119
# Note: For AWS Batch array jobs the "ecs/<job_id>" does not include the array index
# Note: GitHub Actions uses the prefix "actions_job"
const CGROUP_REGEX = r"\/(?:docker|ecs\/[0-9a-f]{32}|actions_job)\/(?<container_id>[0-9a-f]{64})\b"

# Determine the container ID of the currently running container
function container_id()
    id = ""
    isfile("/proc/self/cgroup") || return id
    open("/proc/self/cgroup") do fp
        while !eof(fp)
            line = chomp(readline(fp))
            value = split(line, ':')[3]
            m = match(CGROUP_REGEX, value)
            if m !== nothing
                id = m[:container_id]
                break
            end
        end
    end
    return String(id)
end


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie, RDatasets
let
    dset = dataset("datasets", "iris")
    byCat = dset.Species
    categ = unique(byCat)
    markers = [:circle, :diamond, :utriangle]

    fig = Figure(resolution = (700, 450))
    ax = Axis(fig, xlabel = "Sepal Length", ylabel = "Sepal Width")
    for (idx,c) in enumerate(categ)
        indc = findall(x->x == c, byCat)
        scatter!(dset.SepalLength[indc], dset.SepalWidth[indc],
            marker = markers[idx], markersize = 15, label = "$(c)")
    end
    axislegend("Species")
    fig[1,1] = ax
    fig
    save(joinpath(@__DIR__, "output", "irisDataSet.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "RDatasets"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using RCall
R"""
memory.limit(40000)
require(data.table)
N=1e7; K=100
set.seed(1)
DT <- data.table(
  id1 = sample(sprintf("id%03d",1:K), N, TRUE),      # large groups (char)
  id2 = sample(sprintf("id%03d",1:K), N, TRUE),      # large groups (char)
  id3 = sample(sprintf("id%010d",1:(N/K)), N, TRUE), # small groups (char)
  id4 = sample(K, N, TRUE),                          # large groups (int)
  id5 = sample(K, N, TRUE),                          # large groups (int)
  id6 = sample(N/K, N, TRUE),                        # small groups (int)
  v1 =  sample(5, N, TRUE),                          # int in range [1,5]
  v2 =  sample(5, N, TRUE),                          # int in range [1,5]
  v3 =  sample(round(runif(100,max=100),4), N, TRUE) # numeric e.g. 23.5749
)
# fwrite(DT,"DT.csv")
cat("GB =", round(sum(gc()[,2])/1024, 3), "\n")
system.time( DT[, sum(v1), keyby="id1,id2"] )
"""

using DataBench, DataFrames, CSV, TextParse, IterableTables
# df = DataBench.createSynDataFrame(10_000_000, 100)
# @time df = CSV.read("DT.csv");
@time df = TextParse.csvread("DT.csv");
df = DataFrame(df)

fn = sum
byvec = [:id1, :id2]
@time ids = fastby(sum, df, byvec);

aggrb(df) = aggregate(df[[:id1,:id2,:v1]], [:id1,:id2], sum);

using BenchmarkTools
@btime aggrb($df)
@time by(df[[:id1,:id2,:v1]], [:id1,:id2], subdf -> sum(subdf[:v1]));
@time by(df, [:id1,:id2], subdf -> sum(subdf[:v1]));

using FastGroupBy
@time fastby(sum, df, :id1, :v1);

showcols(df)

using SortingAlgorithms
using DataStructures
idstrvec = rand("id".*dec.(1:1_000_000,10), 100_000_000)

using FastGroupBy
@time fastby(sum, idstrvec, fcollect(length(idstrvec)))

sort(idstrvec[1:10])            # compile
@time sort(idstrvec)
#  80.119106 seconds (88 allocations: 1.118 GiB, 5.72% gc time)

"Return `key=>count` pairs for elements in the argument sorted by the `key`."
function sorted_counts(v)
    c = counter(v)
    sort(collect(c.map), by = first)
end
using SortingAlgorithms
sorted_counts(idstrvec[1:10])   # compile
@time sorted_counts(idstrvec)
# 49.105053 seconds (1.00 M allocations: 114.762 MiB, 0.35% gc time)


text:
mountaincar = QuickMDP(
    function (s, a, rng)
        x, v = s
        vp = clamp(v + a*0.001 + cos(3*x)*-0.0025, -0.07, 0.07)
        xp = x + vp
        if xp > 0.5
            r = 100.0
        else
            r = -1.0
        end
        return (sp=(xp, vp), r=r)
    end,
    actions = [-1., 0., 1.],
    initialstate = Deterministic((-0.5, 0.0)),
    discount = 0.95,
    isterminal = s -> s[1] > 0.5
)


text:
# Module to communicate with GW instek PSP series power supplies
#=
note:
1)  The serial port on the PSP-603 (and possibly others) is non standard.
    Pin 4 needs +12V to power the line driver.  The transmitter will swing from
    Vpin4 - 4v to 0v.  This worked with my RS232 adapter.
2)  The equipment and probably your computer are both DTE.  Null modem
    is required.
3)  Configure COM port for 2400,8,1,N flow control none.
=#

# not including margining commands

module GWinstekPSP

# functions that read values
include("psp_status.jl") # L
include("psp_readvoltage.jl") # V
include("psp_readcurrent.jl") # A
include("psp_readpower.jl") # W
include("psp_readvlimit.jl") # U
include("psp_readilimit.jl") # I
include("psp_readplimit.jl") # P
include("psp_readflags.jl") # F

# functions that set values
include("psp_incrementvoltage.jl")  # SV+
include("psp_decrementvoltage.jl")  # SV-
include("psp_incrementvlimit.jl")   # SU+
include("psp_decrementvlimit.jl")   # SU-

include("psp_incrementilimit.jl")   # SI+
include("psp_decrementilimit.jl")   # SI-

include("psp_incrementplimit.jl")   # SP+
include("psp_decrementplimit.jl")   # SP-

include("psp_setvlimitmax.jl") # SUM
include("psp_setilimitmax.jl") # SIM
include("psp_setplimitmax.jl") # SPM

include("psp_setknobfine.jl") # KF
include("psp_setknobnormal.jl") # KN

include("psp_togglerelay.jl") # KO
include("psp_setrelayon.jl") # KOE
include("psp_setrelayoff.jl") # KOD

include("psp_savepreset.jl") # EEP

include("psp_setvoltage.jl") # SV
include("psp_setvlimit.jl") # SU
include("psp_setilimit.jl") # SI
include("psp_setplimit.jl") # SP

end # module


text:
# Julia Module Manipulate the Coefficient Matrix in Standard LP or MILP Problem
# Version: 2.0
# Auther: Edward J. Xu
# Date: Aprial 7th, 2019
module CoeffMatrix
# Function to initiate the coefficient matrix, update the last row, or add a row in the last.
# 1. default double-script order:
#   |     n(j)|
#   |m(i)  *  |
#
# 2. default single-script order: k = (j - 1) * m + i (To bottom, then right)
#    to convert back: i = k % m
#                     j = Int8((k - i) / m)
#
# 3. if inverse == TURE, use inversed single-script order: k = (i - 1) * n + j (To right, then bottom)
#    to convert back: j = k % n
#                     i = Int8((k - j) / n)
#
function updateRow(; mat_coeff = 0, add = false, m, n, mat_para,
    vec_i_x, vec_j_x, vec_i_a = vec_i_x, vec_j_a = vec_j_x, inverse = false)
    num_x = m * n
    # 1. Check if m, n and mat_coeff are coherent.
    if mat_coeff != 0
        (num_row_coeff, num_col_coeff) = size(mat_coeff)
        if num_x != num_col_coeff
            print("Wrong m or n, of which the product should be same with col number of mat_coeff!")
        end
    end
    # 2. Generate mat_para_cal with same row and col numbers with mat_x, if mat_para is not a matrix
    if isa(mat_para, Number)  # If mat_para is a number or a vector, generate a m*n matrix for later use
        mat_para_cal = mat_para * ones(m, n)
    end
    # If row vector, check if col number equals that of mat_x.
    # If col vector, check if row number equals that of mat_x.
    (num_row_a, num_col_a) = size(mat_para)
    if ((num_row_a != 1) && (num_row_a != m)) || ((num_col_a != 1) && (num_col_a != n))
        print("Wrong mat_para, which should be vector or a matix in same size with mat_x!")
        return (mat_coeff, mat_para)
    elseif num_row_a < m  # if mat_para is a row vector, duplicate it to m-row matrix
        mat_para_cal = ones(m, n)
        for i = 1: m
            mat_para_cal[i, :] = mat_para[1, :]
        end
        mat_para = mat_para_new
    elseif num_col_a < n  # if mat_para is a strict column vector, duplicate it to n-column matrix
        mat_para_cal = ones(m, n)
        for j = 1: n
            mat_para_cal[:, j] = mat_para[:, 1]
        end
    else
        mat_para_cal = mat_para
    end
    # 3. Calculate the new row vector for coefficient matrix
    vecRow_coeff_new = zeros(1, num_x)
    for p = 1: length(vec_i_x)
        for q = 1: length(vec_j_x)
            i = vec_i_x[p]
            j = vec_j_x[q]
            if inverse  # inverse single-script order
                k = (i - 1) * n + j  # To right, then bottom
                vecRow_coeff_new[(i - 1) * n + vec_j_x[j]] = mat_para_cal[vec_i_a[p], vec_j_a[q]]
            else  # default single-script order
                k = (j - 1) * m + i  # To bottom, then right
                vecRow_coeff_new[k] = mat_para_cal[vec_i_a[p], vec_j_a[q]]
            end
        end
    end
    # 4. Update the last row of coefficient matrix or add a row in the last, if "add = 1".
    if mat_coeff == 0  # Initiate the coefficient matrix, if no coefficient matrix input.
        mat_coeff = vecRow_coeff_new
    elseif add  # Add a row in the last
        mat_coeff = vcat(mat_coeff, vecRow_coeff_new)
    else  # Update the last row
        mat_coeff[num_row_coeff, :] = mat_coeff[num_row_coeff, :] + transpose(vecRow_coeff_new)
    end
    return (mat_coeff, mat_para_cal)
end
# Example:
# mat_para = [99.74  99.21  100.21  99.76  100.48  100.7  99.42  99.11  97.69  98.94  97.22  98.99;
#          1.89   2.02   1.95    2.16   2.16    2.08   2.08   2.09   1.97   1.99   1.89   1.95;
#          15.58  15.11  14.93   14.86  15.06   15.51  14.41  14.96  15.62  14.4   15.64  14.52]
# (mat_coeff, mat_para_cal) = updateRow(m = 3, n = 12, mat_para = - hcat(mat_para[:,1]),
#     vec_i_x = [1 2 3], vec_j_x = [1])
# (mat_coeff, mat_para_cal) = updateRow(add = 1, mat_coeff = mat_coeff, m = 3, n = 12,
#     mat_para = - mat_para, vec_i_x = [1 2 3], vec_j_x = collect(2: 1: 12))
# (mat_coeff, mat_para_cal) = updateRow(add = 1, mat_coeff = mat_coeff, m = 3, n = 12,
#     mat_para = mat_para, vec_i_x = [1 2 3], vec_j_x = collect(1: 1: 11), vec_j_a = collect(1: 1: 11) .+ 1)
end  # module CoeffMatrix


text:
using NBodyIPs, JuLIP, LinearAlgebra
using NBodyIPs: combinebasis, BASIS

e1 = rand()
e2 = rand()
b1 = OneBody(e1)
b2 = OneBody(e2)
at = bulk(:Cu, cubic=true) * 3
println(@test(energy(b1, at) ≈ length(at) * b1.E0))
println(@test(maximum(norm.(forces(b1, at))) == 0.0))
println(@test(OneBody(2*e1+3*e2) == combinebasis([b1, b2], [2,3])))
println(@test(hash(BASIS(), b1) == hash(BASIS(), b2)))

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
"""
    module Vision

Data blocks, encodings and more for computer vision.

The most important [`Block`] is [`Image`](#)`{N}`.

Blocks:

- [`Image`](#)`{N}`: an `N`-dimensional color image
- [`Mask`](#)`{N}`: an `N`-dimensional categorical mask
- [`Keypoints`](#)`{N}`: a fixed number of `N`-dimensional keypoints

Encodings:

- [`OneHot`](#) is implemented for `Mask`s
- [`ImagePreprocessing`](#) prepares an `Image` for training by
    normalizing and expanding the color channels
- [`KeypointPreprocessing`](#) prepares `Keypoints` for training by
    normalizing them.

"""
module Vision

using ..FastAI
using ..FastAI:
    # blocks
    Block, WrapperBlock, AbstractBlock, OneHotTensor, OneHotTensorMulti, Label,
    LabelMulti, wrapped,
    # encodings
    Encoding, StatefulEncoding, OneHot,
    # visualization
    ShowText,
    # other
    FASTAI_METHOD_REGISTRY, registerlearningmethod!, Datasets
import FastAI.Datasets

# for tests
using ..FastAI: testencoding

# extending
import ..FastAI:
    blockmodel, blockbackbone, blocklossfn, encode, decode, checkblock,
    encodedblock, decodedblock, showblock!, mockblock, setup, encodestate,
    decodestate

import Colors: colormaps_sequential, Colorant, Color, Gray, Normed, RGB,
    alphacolor, deuteranopic, distinguishable_colors
import FixedPointNumbers: N0f8
import DataAugmentation
import DataAugmentation: apply, Identity, ToEltype, ImageToTensor, Normalize,
    BufferedThreadsafe, ScaleKeepAspect, PinOrigin, RandomCrop, CenterResizeCrop,
    AdjustBrightness, AdjustContrast, Maybe,
    ResizePadDivisible, itemdata
import ImageInTerminal
import IndirectArrays: IndirectArray
import Requires: @require
import StaticArrays: SVector
import Statistics: mean, std
import UnicodePlots

using InlineTest
using ShowCases

# Blocks
include("blocks/bounded.jl")
include("blocks/image.jl")
include("blocks/mask.jl")
include("blocks/keypoints.jl")

include("encodings/onehot.jl")
include("encodings/imagepreprocessing.jl")
include("encodings/keypointpreprocessing.jl")
include("encodings/projective.jl")

include("models/Models.jl")
include("models.jl")
include("learningmethods/utils.jl")
include("learningmethods/classification.jl")
include("learningmethods/segmentation.jl")
include("learningmethods/keypointregression.jl")
include("recipes.jl")

include("tests.jl")

function __init__()
    _registerrecipes()
    @require Makie="ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" begin
        import .Makie
        import .Makie: @recipe, @lift
        import .FastAI: ShowMakie
        include("makie.jl")
    end
end

export Image, Mask, Keypoints, Bounded,
    # encodings
    ImagePreprocessing, KeypointPreprocessing, ProjectiveTransforms,
    # learning methods
    ImageClassificationSingle, ImageClassificationMulti,
    ImageKeypointRegression, ImageSegmentation

end


text:
include("main.jl")

function tester()
    n     = nworkers()
    out   = open("data.log", "w")
    cd("../instances")
    files = readdir()
    for file in files
        name = split(splitext(file)[1], '-')[1]
        if name == "uf20"
            print("$file \t")
            info = main(file)
            println(info)
            write(out, "$file $info\n")
            flush(out)
        end
    end
    cd("../src")
    close(out)
end

function tester_article()
    names = ["../instances/uf20-01.cnf",
             "../instances/uf50-01.cnf",
             "../instances/uf75-01.cnf",
             "../instances/uf100-01.cnf",
             "../instances/uf250-01.cnf"]

    out   = open("data.dat", "w")
    mmm = 10

    func = "hcos"

    for name in names
        println("$name")
        times   = 0
        iters   = 0
        solveds = 0
        xs      = []
        ys      = []
        for i in 1:mmm
            println("$i")
            time, iter, solved, x, y = main(name)
            times   += time
            iters   += iter
            solveds += solved
            if i == 1
                xs = copy(x)
                ys = copy(y)
            else
                for j in 1:length(xs)
                    xs[j] += x[j]
                    ys[j] += y[j]
                end
            end
        end

        write(out, "$name $(times/mmm) $(iters/mmm) $(solveds/mmm)\n")
        flush(out)

        name = @sprintf("out_%s_%s.dat", splitext(basename(name))[1], func)
        f    = open(name, "w")

        for j in 1:length(xs)
            write(f, "$(xs[j]/mmm) $(ys[j]/mmm)\n")
        end
        close(f)
    end

    close(out)
end


text:
@testset "Setting the puzzle" begin
    @testset "Puzzle is ORDPUSSLE; Setting ORDPUSSLE again; Rejected" begin
        # Arrange
        dictionary = SwedishDictionary(["ORDPUSSEL"])
        game = NiancatGame(dictionary, TestingGameService())
        team = Team(1, "defaultteam", "")
        user = User(1, "name", "displayname", team)
        gamecommand(game, user, SetPuzzle("ORDPUSSLE"))

        # Act
        response = gamecommand(game, user, SetPuzzle("ORDPUSSLE"))

        # Assert
        @test response isa Rejected
    end

    @testset "Puzzle is ORDPUSSLE; Setting the puzzle to DATORSPLE; Allowed" begin
        # Arrange
        dictionary = SwedishDictionary(["ORDPUSSEL", "DATORSPEL"])
        game = NiancatGame(dictionary, TestingGameService())
        team = Team(1, "defaultteam", "")
        user = User(1, "name", "displayname", team)
        gamecommand(game, user, SetPuzzle("ORDPUSSLE"))

        # Act
        response = gamecommand(game, user, SetPuzzle("DATORSPLE"))

        # Assert
        @test response == NewPuzzle("DATORSPLE")
    end

    @testset "Puzzle is ORDPUSSLE; Setting the anagram ORDPUSLES; Rejected" begin
        # Arrange
        dictionary = SwedishDictionary(["ORDPUSSEL"])
        game = NiancatGame(dictionary, TestingGameService())
        team = Team(1, "defaultteam", "")
        user = User(1, "name", "displayname", team)
        gamecommand(game, user, SetPuzzle("ORDPUSSLE"))

        # Act
        response = gamecommand(game, user, SetPuzzle("ORDPUSLES"))

        # Assert
        @test response isa Rejected
    end

    @testset "Puzzle is ORDPUSSLE; Setting the non-normal ordpusles; Rejected" begin
        # Arrange
        dictionary = SwedishDictionary(["ORDPUSSEL"])
        game = NiancatGame(dictionary, TestingGameService())
        team = Team(1, "defaultteam", "")
        user = User(1, "name", "displayname", team)
        gamecommand(game, user, SetPuzzle("ORDPUSSLE"))

        # Act
        response = gamecommand(game, user, SetPuzzle("ordpusles"))

        # Assert
        @test response isa Rejected
    end

    @testset "Puzzle has no solutions; Rejected" begin
        # Arrange
        dictionary = SwedishDictionary(["ORDPUSSEL"])
        game = NiancatGame(dictionary, TestingGameService())
        team = Team(1, "defaultteam", "")
        user = User(1, "name", "displayname", team)

        # Act
        response = gamecommand(game, user, SetPuzzle("DATORSPLE"))

        # Assert
        @test response isa Rejected
    end

    @testset "Puzzle is rejected because of no solutions; The round is still the same" begin
        # Arrange
        dictionary = SwedishDictionary(["ORDPUSSEL"])
        game = NiancatGame(dictionary, TestingGameService())
        team = Team(1, "defaultteam", "")
        user = User(1, "name", "displayname", team)
        gamecommand(game, user, SetPuzzle("ORDPUSSLE"))
        roundbefore = gameround(game)

        # Act
        gamecommand(game, user, SetPuzzle("NOSUCHWORD"))
        roundafter = gameround(game)

        # Assert
        @test roundbefore == roundafter
    end

    @testset "Puzzle is ORDPUSSLE; Setting the anagram ORDPUSLES; The round is unchanged" begin
        # Arrange
        dictionary = SwedishDictionary(["ORDPUSSEL"])
        game = NiancatGame(dictionary, TestingGameService())
        team = Team(1, "defaultteam", "")
        user = User(1, "name", "displayname", team)
        gamecommand(game, user, SetPuzzle("ORDPUSSLE"))
        roundbefore = gameround(game)

        # Act
        response = gamecommand(game, user, SetPuzzle("ORDPUSLES"))
        roundafter = gameround(game)

        # Assert
        @test roundbefore == roundafter
    end

    @testset "Puzzle is ORDPUSSLE; Setting the puzzle to DATORSPLE; The round is changed" begin
        # Arrange
        dictionary = SwedishDictionary(["ORDPUSSEL", "DATORSPEL"])
        game = NiancatGame(dictionary, TestingGameService())
        team = Team(1, "defaultteam", "")
        user = User(1, "name", "displayname", team)
        gamecommand(game, user, SetPuzzle("ORDPUSSLE"))
        roundbefore = gameround(game)

        # Act
        response = gamecommand(game, user, SetPuzzle("DATORSPLE"))
        roundafter = gameround(game)

        # Assert
        @test roundbefore != roundafter
    end

    @testset "Setting the puzzle to PUSSGUKRA; One solution; One solution event sent" begin
        # Arrange
        dictionary = SwedishDictionary(["PUSSGURKA"])
        gameservice = TestingGameService()
        game = NiancatGame(dictionary, gameservice)
        team = Team(1, "defaultteam", "")
        user = User(1, "name", "displayname", team)

        # Act
        response = gamecommand(game, user, SetPuzzle("PUSSGUKRA"))

        # Assert
        @test length(gameservice.gameevents) != 0
        event = lastgameevent(gameservice)
        @test event.eventtype == Gameface.GameEvent_Solution
        @test event.round == gameround(game)
        @test event.data == normalizedword(dictionary, "PUSSGURKA")
    end

    @testset "Setting the puzzle to DATORSPLE; Two solutions; Two solution events sent" begin
        # Arrange
        dictionary = SwedishDictionary(["DATORSPEL", "LEDARPOST"])
        gameservice = TestingGameService()
        game = NiancatGame(dictionary, gameservice)
        team = Team(1, "defaultteam", "")
        user = User(1, "name", "displayname", team)

        # Act
        response = gamecommand(game, user, SetPuzzle("DATORSPLE"))

        # Assert
        @test length(gameservice.gameevents) >= 2
        solutions = [event.data
                     for event in gameservice.gameevents
                     if event.eventtype == Gameface.GameEvent_Solution]
        sortedsolutions = sort(solutions)
        @test sortedsolutions == ["datorspel", "ledarpost"]
    end

    @testset "Puzzle is not set; Setting ORDPUSSLE; NewPuzzle response sent" begin
        # Arrange
        dictionary = SwedishDictionary(["ORDPUSSEL"])
        game = NiancatGame(dictionary, TestingGameService())
        team = Team(1, "defaultteam", "")
        user = User(1, "name", "displayname", team)

        # Act
        response = gamecommand(game, user, SetPuzzle("ORDPUSSLE"))

        # Assert
        @test response isa NewPuzzle
    end

    @testset "Puzzle is PUSSGUKRA; Setting a new puzzle DATORSPLE; A solution board notification is sent" begin
        # Arrange
        gameservice = TestingGameService()
        dictionary = SwedishDictionary(["PUSSGURKA", "DATORSPEL"])
        game = NiancatGame(dictionary, gameservice)
        team = Team(1, "defaultteam", "")
        user = User(1, "name", "displayname", team)
        gamecommand(game, user, SetPuzzle("PUSSGRUKA"))
        instanceid = gameinstanceid(game)
        round = gameround(game)

        # Act
        response = gamecommand(game, user, SetPuzzle("DATORSPLE"))
        # Precondition: Response is a NewPuzzle
        @test response isa NewPuzzle

        # Assert
        @test length(gameservice.notifications) != 0
        @test SolutionboardNotification(instanceid, round) in gameservice.notifications
    end

    @testset "No puzzle is set; Setting a new puzzle DATORSPLE; No solution board notification is sent" begin
        # Arrange
        gameservice = TestingGameService()
        dictionary = SwedishDictionary(["PUSSGURKA", "DATORSPEL"])
        game = NiancatGame(dictionary, gameservice)
        team = Team(1, "defaultteam", "")
        user = User(1, "name", "displayname", team)

        # Act
        gamecommand(game, user, SetPuzzle("DATORSPLE"))

        # Assert
        @test !any(x -> x isa SolutionboardNotification, gameservice.notifications)
    end

    @testset "Puzzle is PUSSGUKRA; Setting a new puzzle DATORSPLE; A score board notification is sent" begin
        # Arrange
        gameservice = TestingGameService()
        dictionary = SwedishDictionary(["PUSSGURKA", "DATORSPEL"])
        game = NiancatGame(dictionary, gameservice)
        team = Team(1, "defaultteam", "")
        user = User(1, "name", "displayname", team)
        gamecommand(game, user, SetPuzzle("PUSSGRUKA"))
        instanceid = gameinstanceid(game)
        round = gameround(game)

        # Act
        response = gamecommand(game, user, SetPuzzle("DATORSPLE"))
        # Precondition: Response is a NewPuzzle
        @test response isa NewPuzzle

        # Assert
        @test length(gameservice.notifications) != 0
        @test ScoreboardNotification(instanceid, round) in gameservice.notifications
    end

    @testset "No puzzle is set; Setting a new puzzle DATORSPLE; No score board notification is sent" begin
        # Arrange
        gameservice = TestingGameService()
        dictionary = SwedishDictionary(["PUSSGURKA", "DATORSPEL"])
        game = NiancatGame(dictionary, gameservice)
        team = Team(1, "defaultteam", "")
        user = User(1, "name", "displayname", team)

        # Act
        gamecommand(game, user, SetPuzzle("DATORSPLE"))

        # Assert
        @test !any(x -> x isa ScoreboardNotification, gameservice.notifications)
    end
end


text:
## Partitions of Sample Space Ω
## ============================

"""A partition of a set containing both an under approximation, and the rest
  The rest is an overapproximation - rest"""
type ApproxPartition{D<:Domain}
  under::Vector{D}
  rest::Vector{D}
end

"Lower and upper bounds for measure of partition"
function measure(p::ApproxPartition)
  lb = isempty(p.under) ? 0 : sum([measure(event) for event in p.under])
  rest = isempty(p.rest) ? 0 : sum([measure(event) for event in p.rest])
  ub = lb + rest
  Interval(lb,ub)
end

"Get an over approximation"
collect{D}(p::ApproxPartition{D}) = vcat(p.under, p.rest)

## Sampling
## ========

"A partition which is efficient for drawing many samples."
type SampleablePartition{D<:Domain}
  over::Vector{D}
  last_under::Int
  cat::Categorical

  function SampleablePartition(under::Vector, rest::Vector)
    over = vcat(under, rest)
    vols = Float64[measure(box) for box in over]
    pnormalize!(vols)
    cat = Categorical(vols, Distributions.NoArgCheck())
    new(over,length(over),cat)
  end
end

SampleablePartition{D}(p::ApproxPartition{D}) = SampleablePartition{D}(p.under, p.rest)

"Point sample from preimage - may be invalid point due to approximations"
abstract_sample(p::SampleablePartition) = p.over[rand(p.cat)]

"`n` Point samples from preimage - may be invalid point due to approximations"
function point_sample(p::SampleablePartition, n::Integer)
  rand_indices = rand(p.cat, n)
  events = [p.over[i] for i in rand_indices]
  [rand(event) for event in events]
end

point_sample(p::SampleablePartition) = point_sample(p,1)[1]

"Do refined rejection sampling from preimage"
function point_sample_exact(p::SampleablePartition, Y::RandVar{Bool}; maxtries = 1E7)
  for i = 1:maxtries
    sample = rand(p)
    if call(Y,sample) return sample end
  end
  #TODO: Better error handling
  error("Could not get sample in $maxtries tries")
end

# Preimage of Y under F, unioned with X
function pre_partition{T}(
    Y::SymbolicRandVar{Bool},
    ::Type{T};
    RandVarType::Type = default_randvar(),
    precision = default_precision(),
    args...)

  init_box = unit_box(LazyBox{Float64}, dims(Y))
  Y_conv = convert(RandVarType{Bool}, Y)
  set_precision!(Y_conv, precision)
  pre_partition(Y_conv, init_box, T; precision=precision, args...)
end


text:
HEADER = "// Generated by [BasicBSplineExporter.jl](https://github.com/hyrodium/BasicBSplineExporter.jl)\n"

"""
2 -> "2"
"""
function _scriptpov(x::Integer)
    if isfinite(x)
        return repr(x)
        # return @sprintf "%1.24f" x
    else
        error("numerical value must be finite")
    end
end

"""
2.3 -> "2.3"
"""
function _scriptpov(x::Real)
    if isfinite(x)
        return repr(convert(Float64,x))
        # return @sprintf "%1.24f" x
    else
        error("numerical value must be finite")
    end
end

"""
[2.3, -5.2] -> "<2.3,-5.2>"
"""
function _scriptpov(x::AbstractVector{<:Real})
    return "<"*join(_scriptpov.(x), ", ")*">"
end

"""
[[2.3, -5.2], [-1.1, 8.2]] -> "<2.3,-5.2>, <-1.1,8.2>"
"""
function _scriptpov(x::AbstractVector{<:AbstractVector{<:Real}})
    return join(_scriptpov.(x), ", ")
end

"""
RGB(1,2,3) -> "rgb<1,2,3>"
"""
function _scriptpov(c::AbstractRGB)
    _c = RGB(c)
    v = [_c.r,_c.g,_c.b]
    return "rgb"*_scriptpov(v)
end

"""
RGB(1,2,3) -> "rgb<1,2,3>"
"""
function _scriptpov(c::ColorAlpha)
    _c = RGBA(c)
    v = [_c.r,_c.g,_c.b,1-_c.alpha]
    return "rgbt"*_scriptpov(v)
end

function _spheres(M::BSplineManifold{Dim, Deg, <:StaticVector{3}}; preindent=0, radius_name="radius_sphere") where {Dim,Deg}
    _spheres(controlpoints(M), preindent=preindent, radius_name=radius_name)
end

function _cylinders(M::BSplineManifold{Dim, Deg, <:StaticVector{3}}; preindent=0, radius_name="radius_cylinder") where {Dim,Deg}
    _cylinders(controlpoints(M), preindent=preindent, radius_name=radius_name)
end

function _spheres(a::AbstractArray{<:AbstractVector{<:Real}}; preindent=0, radius_name="radius_sphere")
    script = "  "^(preindent)
    script *= "union{\n" * "  "^(preindent)
    for ai in a
        script *= "  sphere{$(_scriptpov(ai)), $(radius_name)}\n" * "  "^(preindent)
    end
    script *= "}\n"
end

function _cylinders(a::AbstractArray{<:AbstractVector{<:Real},1}; preindent=0, radius_name="radius_cylinder")
    n = length(a)
    script = "  "^(preindent)
    script *= "union{\n" * "  "^(preindent)
    for i in 1:n-1
        script *= "  cylinder{$(_scriptpov(a[i])), $(_scriptpov(a[i+1])), $(radius_name)}\n" * "  "^(preindent)
    end
    script *= "}\n"
end

function _cylinders(a::AbstractArray{<:AbstractVector{<:Real},2}; preindent=0, radius_name="radius_cylinder")
    n1, n2 = size(a)
    script = "  "^(preindent)
    script *= "union{\n"
    for i1 in 1:n1
        script *= _cylinders(a[i1,:], preindent=preindent+1, radius_name=radius_name)
    end
    for i2 in 1:n2
        script *= _cylinders(a[:,i2], preindent=preindent+1, radius_name=radius_name)
    end
    script *= "  "^(preindent)
    script *= "}\n"
end

function _cylinders(a::AbstractArray{<:AbstractVector{<:Real},3}; preindent=0, radius_name="radius_cylinder")
    n1, n2, n3 = size(a)
    script = "  "^(preindent)
    script *= "union{\n"
    for i2 in 1:n2, i3 in 1:n3
        script *= _cylinders(a[:,i2,i3], preindent=preindent+1, radius_name=radius_name)
    end
    for i1 in 1:n1, i3 in 1:n3
        script *= _cylinders(a[i1,:,i3], preindent=preindent+1, radius_name=radius_name)
    end
    for i1 in 1:n1, i2 in 1:n2
        script *= _cylinders(a[i1,i2,:], preindent=preindent+1, radius_name=radius_name)
    end
    script *= "  "^(preindent)
    script *= "}\n"
end

function _curve(M::BSplineManifold{1, Deg, <:StaticVector{3}}; mesh=10, preindent=0, radius_name="radius_curve") where Deg
    P = bsplinespaces(M)[1]
    p = degree(P)
    k = knotvector(P)
    l = length(k)
    m = mesh

    ts = unique!(vcat([range(k[i],k[i+1],length = m+1) for i in p+1:l-p-1]...))

    𝒑s = M.(ts)

    script = "  "^(preindent)
    script *= "union{\n" * "  "^(preindent)
    script *= "  object{\n"
    script *= _spheres(𝒑s, preindent=preindent+2, radius_name=radius_name) * "  "^(preindent)
    script *= "  }\n" * "  "^(preindent)
    script *= "  object{\n"
    script *= _cylinders(𝒑s, preindent=preindent+2, radius_name=radius_name) * "  "^(preindent)
    script *= "  }\n" * "  "^(preindent)
    script *= "}\n"

    return script
end

"""
Generate POV-Ray script of `mesh2`
http://povray.org/documentation/3.7.0/r3_4.html#r3_4_5_2_4
"""
function _surface(M::BSplineManifold{2, Deg, <:StaticVector{3}}; mesh::Int=10, smooth=true, preindent=0) where Deg
    P = bsplinespaces(M)
    p1, p2 = p = degree.(P)
    k1, k2 = k = knotvector.(P)
    l1, l2 = l = length.(k)
    m1, m2 = m = (mesh, mesh)

    ts1 = unique!(vcat([range(k1[i],k1[i+1],length = m1+1) for i in p1+1:l1-p1-1]...))
    ts2 = unique!(vcat([range(k2[i],k2[i+1],length = m2+1) for i in p2+1:l2-p2-1]...))
    N1 = length(ts1)-1
    N2 = length(ts2)-1
    cs1 = [(ts1[i]+ts1[i+1])/2 for i in 1:N1]
    cs2 = [(ts2[i]+ts2[i+1])/2 for i in 1:N1]

    # 𝒆(t) = normalize(cross(BasicBSpline.tangentvectors(M,t...)...))

    𝒑s = [M(t1,t2) for t1 in ts1, t2 in ts2]
    𝒑c = [M(t1,t2) for t1 in cs1, t2 in cs2]

    # 𝒆s = 𝒆.(ts)
    # 𝒆c = 𝒆.(tc)

    Ns(i1, i2) = i1 + (N1+1) * (i2-1)
    Nc(i1, i2) = i1 + N1 * (i2-1)

    F1 = [[Ns(i1,i2)-1, Ns(i1+1,i2)-1, (N1+1)*(N2+1)+Nc(i1,i2)-1] for i1 in 1:N1, i2 in 1:N2]
    F2 = [[Ns(i1,i2)-1, Ns(i1,i2+1)-1, (N1+1)*(N2+1)+Nc(i1,i2)-1] for i1 in 1:N1, i2 in 1:N2]
    F3 = [[Ns(i1+1,i2+1)-1, Ns(i1+1,i2)-1, (N1+1)*(N2+1)+Nc(i1,i2)-1] for i1 in 1:N1, i2 in 1:N2]
    F4 = [[Ns(i1+1,i2+1)-1, Ns(i1,i2+1)-1, (N1+1)*(N2+1)+Nc(i1,i2)-1] for i1 in 1:N1, i2 in 1:N2]

    np = (N1+1)*(N2+1) + N1*N2
    nf = 4*N1*N2

    script = "  "^(preindent)
    script *= "mesh2{\n" * "  "^(preindent)
    script *= "  vertex_vectors{\n" * "  "^(preindent)
    script *= "    " * _scriptpov(np) * ", \n" * "  "^(preindent)
    script *= "    " * _scriptpov([𝒑s...]) * ", \n" * "  "^(preindent)
    script *= "    " * _scriptpov([𝒑c...]) * "\n" * "  "^(preindent)
    script *= "  }\n" * "  "^(preindent)

    # if smooth
    #     script *= "  normal_vectors{\n" * "  "^(preindent)
    #     script *= "    " * _scriptpov(np) * ", \n" * "  "^(preindent)
    #     script *= "    " * _scriptpov([𝒆s...]) * ", \n" * "  "^(preindent)
    #     script *= "    " * _scriptpov([𝒆c...]) * "\n" * "  "^(preindent)
    #     script *= "  }\n" * "  "^(preindent)
    # end

    script *= "  face_indices{\n" * "  "^(preindent)
    script *= "    " * _scriptpov(nf) * ", \n" * "  "^(preindent)
    script *= "    " * _scriptpov([F1...]) * "\n" * "  "^(preindent)
    script *= "    " * _scriptpov([F2...]) * "\n" * "  "^(preindent)
    script *= "    " * _scriptpov([F3...]) * "\n" * "  "^(preindent)
    script *= "    " * _scriptpov([F4...]) * "\n" * "  "^(preindent)
    script *= "  }\n" * "  "^(preindent)
    script *= "}\n"

    return script
end

function _solid(M::BSplineManifold{3, Deg, <:StaticVector{3}}; mesh=10, smooth=true, preindent=0) where Deg
    surfaces = _bsplinesurfaces(M)

    script = "  "^(preindent)
    script *= "union{\n"
    for M in surfaces
        script *= _surface(M, mesh=mesh, preindent=preindent+1)
    end
    script *= "  "^(preindent)
    script *= "}\n"

    return script
end

function _save_povray_1d3d(name::String, M::BSplineManifold{1, Deg, <:StaticVector{3}}; mesh::Int=10, points=true, thickness=0.1, maincolor=RGB(1,0,0), subcolor=RGB(.5,.5,.5)) where Deg
    radius_curve = thickness
    radius_cylinder = thickness
    radius_sphere = 2*radius_cylinder
    color_curve = maincolor
    color_cylinder = subcolor
    color_sphere = weighted_color_mean(0.5, subcolor, colorant"black")
    script = HEADER
    script *= "#local radius_curve = $(radius_curve);\n"
    script *= "#local radius_sphere = $(radius_sphere);\n"
    script *= "#local radius_cylinder = $(radius_cylinder);\n"
    script *= "#local color_curve = $(_scriptpov(color_curve));\n"
    script *= "#local color_sphere = $(_scriptpov(color_sphere));\n"
    script *= "#local color_cylinder = $(_scriptpov(color_cylinder));\n"
    script *= "union{\n"
    script *= "  object{\n"
    script *= _curve(M, mesh=mesh, preindent=2, radius_name=radius_curve)
    script *= "    pigment{color_curve}\n"
    script *= "  }\n"
    if points
        script *= "  object{\n"
        script *= _spheres(M, preindent=2)
        script *= "    pigment{color_sphere}\n"
        script *= "  }\n"
        script *= "  object{\n"
        script *= _cylinders(M, preindent=2)
        script *= "    pigment{color_cylinder}\n"
        script *= "  }\n"
    end
    script *= "}"

    open(name, "w") do f
        write(f,script)
    end

    return nothing
end

function _save_povray_2d3d(name::String, M::BSplineManifold{2, Deg, <:StaticVector{3}}; mesh::Int=10, points=true, thickness=0.1, maincolor=RGB(1,0,0), subcolor=RGB(.5,.5,.5)) where Deg
    radius_cylinder = thickness
    radius_sphere = 2*radius_cylinder
    color_cylinder = subcolor
    color_sphere = weighted_color_mean(0.5, subcolor, colorant"black")
    color_surface = maincolor
    script = HEADER
    script *= "#local radius_sphere = $(radius_sphere);\n"
    script *= "#local radius_cylinder = $(radius_cylinder);\n"
    script *= "#local color_sphere = $(_scriptpov(color_sphere));\n"
    script *= "#local color_cylinder = $(_scriptpov(color_cylinder));\n"
    script *= "#local color_surface = $(_scriptpov(color_surface));\n"
    script *= "union{\n"
    script *= "  object{\n"
    script *= _surface(M, mesh=mesh, preindent=2)
    script *= "    pigment{color_surface}\n"
    script *= "  }\n"
    if points
        script *= "  object{\n"
        script *= _spheres(M, preindent=2)
        script *= "    pigment{color_sphere}\n"
        script *= "  }\n"
        script *= "  object{\n"
        script *= _cylinders(M, preindent=2)
        script *= "    pigment{color_cylinder}\n"
        script *= "  }\n"
    end
    script *= "}"

    open(name, "w") do f
        write(f,script)
    end

    return nothing
end

function _save_povray_3d3d(name::String, M::BSplineManifold{3, Deg, <:StaticVector{3}}; mesh::Int=10, points=true, thickness=0.1, maincolor=RGB(1,0,0), subcolor=RGB(.5,.5,.5)) where Deg
    radius_cylinder = thickness
    radius_sphere = 2*radius_cylinder
    color_cylinder = subcolor
    color_sphere = weighted_color_mean(0.5, subcolor, colorant"black")
    color_solid = maincolor
    script = HEADER
    script *= "#local radius_sphere = $(radius_sphere);\n"
    script *= "#local radius_cylinder = $(radius_cylinder);\n"
    script *= "#local color_sphere = $(_scriptpov(color_sphere));\n"
    script *= "#local color_cylinder = $(_scriptpov(color_cylinder));\n"
    script *= "#local color_solid = $(_scriptpov(color_solid));\n"
    script *= "union{\n"
    script *= "  object{\n"
    script *= _solid(M, mesh=mesh, preindent=2)
    script *= "    pigment{color_solid}\n"
    script *= "  }\n"
    if points
        script *= "  object{\n"
        script *= _spheres(M, preindent=2)
        script *= "    pigment{color_sphere}\n"
        script *= "  }\n"
        script *= "  object{\n"
        script *= _cylinders(M, preindent=2)
        script *= "    pigment{color_cylinder}\n"
        script *= "  }\n"
    end
    script *= "}"

    open(name, "w") do f
        write(f,script)
    end

    return nothing
end

function save_pov(name::String, M::AbstractBSplineManifold{1}; mesh::Int=10, points=true, thickness=0.1, maincolor=RGB(1,0,0), subcolor=RGB(.5,.5,.5))
    _save_povray_1d3d(name,M,mesh=mesh,points=points,thickness=thickness,maincolor=maincolor,subcolor=subcolor)
end

function save_pov(name::String, M::AbstractBSplineManifold{2}; mesh::Int=10, points=true, thickness=0.1, maincolor=RGB(1,0,0), subcolor=RGB(.5,.5,.5))
    _save_povray_2d3d(name,M,mesh=mesh,points=points,thickness=thickness,maincolor=maincolor,subcolor=subcolor)
end

function save_pov(name::String, M::AbstractBSplineManifold{3}; mesh::Int=10, points=true, thickness=0.1, maincolor=RGB(1,0,0), subcolor=RGB(.5,.5,.5))
    _save_povray_3d3d(name,M,mesh=mesh,points=points,thickness=thickness,maincolor=maincolor,subcolor=subcolor)
end


text:
using AeroMDAO
using Test

@testset "2D Panel Method - Doublet-Source Kulfan CST" begin
    # Define Kulfan CST coefficients
    alpha_u = [0.2, 0.3, 0.2, 0.15, 0.2]
    alpha_l = [-0.2, -0.1, -0.1, -0.001]
    dzs     = (0., 0.)

    # Define airfoil
    airfoil = (Foil ∘ kulfan_CST)(alpha_u, alpha_l, dzs, 0.0, 60);

    # Define uniform flow
    uniform = Uniform2D(1., 5.)

    # Evaluate case
    cl, cls, cms, cps, panels = solve_case(airfoil, uniform; num_panels = 80)

    @test cl       ≈  0.84188988 atol = 1e-6
    @test sum(cls) ≈  0.84073703 atol = 1e-6
    @test sum(cms) ≈ -0.26104277 atol = 1e-6
end

@testset "Geometry - Airfoil Processing" begin
    # Import and read airfoil coordinates
    foilpath = joinpath((dirname ∘ dirname ∘ pathof)(AeroMDAO), "test/CRM.dat")
    coords   = read_foil(foilpath)

    # Cosine spacing
    cos_foil = cosine_foil(coords, 51)

    # Split airfoil
    up, low  = split_foil(cos_foil)

    # Convert coordinates to Kulfan CST variables
    num_dv   = 4
    alpha_u, alpha_l = coords_to_CST(up, num_dv), coords_to_CST(low, num_dv)

    # Generate same airfoil using Kulfan CST parametrisation
    cst_foil = (Foil ∘ kulfan_CST)(alpha_u, alpha_l, (0., 0.), 0.0)

    uniform  = Uniform2D(1., 5.)
    cl, cls, cms, cps, panels = solve_case(cst_foil, uniform; num_panels = 80)

    @test cl       ≈  0.85736965 atol = 1e-6
    @test sum(cls) ≈  0.85976886 atol = 1e-6
    @test sum(cms) ≈ -0.29766116 atol = 1e-6
end

@testset "Geometry - Two-Section Trapezoidal Wing" begin
    # Define wing
    wing_right = HalfWing(chords    = [1.0, 0.6, 0.2],
                          twists    = [2.0, 0.0, -0.2],
                          spans     = [5.0, 0.5],
                          dihedrals = [5., 5.],
                          sweep_LEs = [5., 5.]);

    # Get wing info
    b, S, c, AR = info(wing_right)
    λ           = taper_ratio(wing_right)
    wing_mac    = mean_aerodynamic_center(wing_right)

    @test b        ≈ 5.50000000                    atol = 1e-6
    @test S        ≈ 4.19939047                    atol = 1e-6
    @test c        ≈ 0.79841269                    atol = 1e-6
    @test AR       ≈ 7.20342634                    atol = 1e-6
    @test λ        ≈ 0.20000000                    atol = 1e-6
    @test wing_mac ≈ [0.42092866, 1.33432539, 0.0] atol = 1e-6
end

@testset "Vortex Lattice Method - NACA 0012 Rectangular Wing" begin
    # Define wing
    wing = Wing(foils     = Foil.(naca4((0,0,1,2)) for i ∈ 1:2),
                chords    = [0.18, 0.16],
                twists    = [0., 0.],
                spans     = [0.5,],
                dihedrals = [5.],
                sweep_LEs = [1.14])

    # Define reference values
    ρ   = 1.225
    ref = [0.25 * mean_aerodynamic_chord(wing), 0., 0.]
    Ω   = [0.0, 0.0, 0.0]

    # Define freestream condition
    uniform = Freestream(10.0, 2.0, 2.0, Ω)

    # Evaluate stability case
    nf_coeffs, ff_coeffs, dv_coeffs = solve_stability_case(wing, uniform; rho_ref = ρ, r_ref = ref, span_num = 20, chord_num = 5)

    # Test values
    nf_tests = [0.001189, -0.000228, 0.152203, -0.000242, -0.003486, -8.1e-5, 0.0, 0.0, 0.0]
    ff_tests = [0.00123,  -0.000271, 0.152198]
    dv_tests = [ 0.068444 -0.000046 -0.000711  0.023607  0.000337;
                 0.010867 -0.007536  0.129968  0.021929 -0.012086;
                 4.402229 -0.012973 -0.070654  6.833903  0.001999;
                 0.031877 -0.013083  0.460035  0.091216 -0.039146;
                -0.112285 -0.004631  0.105695 -0.852395 -0.007696;
                -0.002218 -0.002115  0.008263 -0.003817  0.001079]

    # Nearfield coefficients test
    [ @test nf_c ≈ nf_t atol = 1e-6 for (nf_c, nf_t) in zip(nf_coeffs, nf_tests) ]
    # Farfield coefficients test
    [ @test ff_c ≈ ff_t atol = 1e-6 for (ff_c, ff_t) in zip(ff_coeffs, ff_tests) ]
    # Stability derivatives' coefficients test
    [ @test dv_c ≈ dv_t atol = 1e-6 for (dv_c, dv_t) in zip(dv_coeffs, dv_tests) ]
end

@testset "Vortex Lattice Method - Vanilla Aircraft" begin
    ## Wing
    wing = Wing(foils = Foil.(fill(naca4((0,0,1,2)), 2)),
                chords    = [1.0, 0.6],
                twists    = [0.0, 0.0],
                spans     = [5.0],
                dihedrals = [11.39],
                sweep_LEs = [0.]);

    # Horizontal tail
    htail_foils = Foil.(fill(naca4((0,0,1,2)), 2))
    htail = Wing(foils     = htail_foils,
                 chords    = [0.7, 0.42],
                 twists    = [0.0, 0.0],
                 spans     = [1.25],
                 dihedrals = [0.],
                 sweep_LEs = [6.39])

    # Vertical tail
    vtail_foils = Foil.(fill(naca4((0,0,0,9)), 2))
    vtail = HalfWing(foils     = vtail_foils,
                     chords    = [0.7, 0.42],
                     twists    = [0.0, 0.0],
                     spans     = [1.0],
                     dihedrals = [0.],
                     sweep_LEs = [7.97])

    ## Assembly
    wing_panels  = panel_wing(wing, 16, 10;
                              spacing = "cosine")
    htail_panels = panel_wing(htail, 6, 6;
                              position = [4., 0, 0],
                              angle    = deg2rad(-2.),
                              axis     = [0., 1., 0.],
                              spacing  = "cosine")
    vtail_panels = panel_wing(vtail, 5, 6;
                              position  = [4., 0, 0],
                              angle     = π/2,
                              axis      = [1., 0., 0.],
                              spacing   = "cosine")

    aircraft = Dict("Wing"            => wing_panels,
                    "Horizontal Tail" => htail_panels,
                    "Vertical Tail"   => vtail_panels)

    ## Reference quantities
    ac_name = "My Aircraft"
    S, b, c = projected_area(wing), span(wing), mean_aerodynamic_chord(wing)
    ρ       = 1.225
    ref     = [0.25c, 0., 0.]
    V, α, β = 1.0, 1.0, 1.0
    Ω       = [0.0, 0.0, 0.0]
    fs      = Freestream(V, α, β, Ω)

    ## Stability case
    dv_data = solve_stability_case(aircraft, fs;
                                   rho_ref   = ρ,
                                   r_ref     = ref,
                                   area_ref  = S,
                                   span_ref  = b,
                                   chord_ref = c,
                                   name      = ac_name);

    nfs, ffs, dvs = dv_data[ac_name]

    nf_tests = [0.000258, -0.006642, 0.074301, -0.003435, 0.075511, 0.001563, 0.0, 0.0, 0.0]
    ff_tests = [0.000375, -0.006685, 0.074281]
    dv_tests = [ 0.016795  0.003460  0.003761   0.093303 -0.000674;
                -0.000863 -0.374410  0.403476   0.000630 -0.253848;
                 5.749765  0.046649 -0.01346   15.571205  0.020396;
                 0.022674 -0.196605  0.660392   0.099065 -0.039688;
                -2.70367  -0.132928  0.070111 -37.372278 -0.064439;
                 0.002034  0.087382  0.014991   0.005840  0.091088]

    # Nearfield coefficients test
    [ @test nf_c ≈ nf_t atol = 1e-6 for (nf_c, nf_t) in zip(nfs, nf_tests) ]
    # Farfield coefficients test
    [ @test ff_c ≈ ff_t atol = 1e-6 for (ff_c, ff_t) in zip(ffs, ff_tests) ]
    # Stability derivatives' coefficients test
    [ @test dv_c ≈ dv_t atol = 1e-6 for (dv_c, dv_t) in zip(dvs, dv_tests) ]
end;

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
using Test
using LogicCircuits

@testset "BitCircuit test" begin

    function test_integrity(bc)
        for el in 1:size(bc.elements,2)
            d = bc.elements[1,el]
            @test bc.nodes[1,d] <= el
            @test bc.nodes[2,d] >= el
            p = bc.elements[2,el]
            @test el ∈ bc.parents[bc.nodes[3,p]:bc.nodes[4,p]]
            s = bc.elements[3,el]
            @test el ∈ bc.parents[bc.nodes[3,s]:bc.nodes[4,s]]
        end
        for node in 1:size(bc.nodes,2)
            first_el = bc.nodes[1,node]
            last_el = bc.nodes[2,node]
            if first_el != 0
                for i = first_el:last_el
                    @test bc.elements[1,i] == node
                end
            end
            first_par = bc.nodes[3,node]
            last_par = bc.nodes[3,node]
            if first_par != 0
                for i = first_par:last_par
                    par = bc.parents[i]
                    @test bc.elements[2,par] == node || bc.elements[3,par] == node
                end
            else
                @test node == num_nodes(bc) || node <= num_leafs(bc)
            end
        end
        @test sum(length, bc.layers) == size(bc.nodes,2)
    end

    r = fully_factorized_circuit(PlainLogicCircuit, 10)
    test_integrity(BitCircuit(r, 10))

    vtree = PlainVtree(10, :balanced)
    r = fully_factorized_circuit(StructLogicCircuit, vtree)
    test_integrity(BitCircuit(r, 10))

    num_vars = 7
    mgr = SddMgr(num_vars, :balanced)
    v = Dict([(i => compile(mgr, Lit(i))) for i=1:num_vars])
    c = (v[1] | !v[2] | v[3]) &
        (v[2] | !v[7] | v[6]) &
        (v[3] | !v[4] | v[5]) &
        (v[1] | !v[4] | v[6])

    r = smooth(PlainLogicCircuit(c)) # flows don't make sense unless the circuit is smooth; cannot smooth trimmed SDDs

    test_integrity(BitCircuit(r, 7))

    l1 = LogicCircuit(Lit(1))
    l2 = LogicCircuit(Lit(2))
    l3 = LogicCircuit(Lit(-1))
    l4 = LogicCircuit(Lit(-2))
    r = (l1 & l2) | (l3 & l4)
    test_integrity(BitCircuit(r,2))

end


text:
function create_param_dual_cum_sum(model::Optimizer)
    return zeros(model.number_of_parameters_in_model)
end

function calculate_dual_of_parameters(model::Optimizer)
    param_dual_cum_sum = create_param_dual_cum_sum(model)
    update_duals_with_affine_constraint_cache!(param_dual_cum_sum, model)
    update_duals_with_quadratic_constraint_cache!(param_dual_cum_sum, model)
    update_duals_in_affine_objective!(
        param_dual_cum_sum,
        model.affine_objective_cache,
    )
    update_duals_in_quadratic_objective!(
        param_dual_cum_sum,
        model.quadratic_objective_cache_pc,
    )
    empty_and_feed_duals_to_model(model, param_dual_cum_sum)
    return model
end

function update_duals_with_affine_constraint_cache!(
    param_dual_cum_sum::Vector{Float64},
    model::Optimizer,
)
    for (F, S) in keys(model.affine_constraint_cache.dict)
        affine_constraint_cache_inner = model.affine_constraint_cache[F, S]
        affine_added_cache_inner = model.affine_added_cache[F, S]
        if !isempty(affine_constraint_cache_inner)
            update_duals_with_affine_constraint_cache!(
                param_dual_cum_sum,
                model.optimizer,
                affine_constraint_cache_inner,
                affine_added_cache_inner,
            )
        end
    end
    return
end

function update_duals_with_affine_constraint_cache!(
    param_dual_cum_sum::Vector{Float64},
    optimizer::OT,
    affine_constraint_cache_inner::MOI.Utilities.DoubleDicts.DoubleDictInner{
        F,
        S,
        V1,
    },
    affine_added_cache_inner::MOI.Utilities.DoubleDicts.DoubleDictInner{F,S,V2},
) where {OT,F,S,V1,V2}
    for (ci, param_array) in affine_constraint_cache_inner
        calculate_parameters_in_ci!(
            param_dual_cum_sum,
            optimizer,
            param_array,
            affine_added_cache_inner[ci],
        )
    end
    return
end

function update_duals_with_quadratic_constraint_cache!(
    param_dual_cum_sum::Vector{Float64},
    model::Optimizer,
)
    for (F, S) in keys(model.quadratic_constraint_cache_pc.dict)
        quadratic_constraint_cache_pc_inner =
            model.quadratic_constraint_cache_pc[F, S]
        if !isempty(quadratic_constraint_cache_pc_inner)
            update_duals_with_quadratic_constraint_cache!(
                param_dual_cum_sum,
                model,
                quadratic_constraint_cache_pc_inner,
            )
        end
    end
    return
end

function update_duals_with_quadratic_constraint_cache!(
    param_dual_cum_sum::Vector{Float64},
    model::Optimizer,
    quadratic_constraint_cache_pc_inner::MOI.Utilities.DoubleDicts.DoubleDictInner{
        F,
        S,
        V,
    },
) where {F,S,V}
    for (poi_ci, param_array) in quadratic_constraint_cache_pc_inner
        moi_ci = model.quadratic_added_cache[poi_ci]
        calculate_parameters_in_ci!(
            param_dual_cum_sum,
            model.optimizer,
            param_array,
            moi_ci,
        )
    end
    return
end

function calculate_parameters_in_ci!(
    param_dual_cum_sum::Vector{Float64},
    optimizer::OT,
    param_array::Vector{MOI.ScalarAffineTerm{T}},
    ci::CI,
) where {OT,CI,T}
    cons_dual = MOI.get(optimizer, MOI.ConstraintDual(), ci)

    for param in param_array
        param_dual_cum_sum[param.variable.value-PARAMETER_INDEX_THRESHOLD] -=
            cons_dual * param.coefficient
    end
    return
end

function update_duals_in_affine_objective!(
    param_dual_cum_sum::Vector{Float64},
    affine_objective_cache::Vector{MOI.ScalarAffineTerm{T}},
) where {T}
    for param in affine_objective_cache
        param_dual_cum_sum[param.variable.value-PARAMETER_INDEX_THRESHOLD] -=
            param.coefficient
    end
    return
end

function update_duals_in_quadratic_objective!(
    param_dual_cum_sum::Vector{Float64},
    quadratic_objective_cache_pc::Vector{MOI.ScalarAffineTerm{T}},
) where {T}
    for param in quadratic_objective_cache_pc
        param_dual_cum_sum[param.variable.value-PARAMETER_INDEX_THRESHOLD] -=
            param.coefficient
    end
    return
end

function empty_and_feed_duals_to_model(
    model::Optimizer,
    param_dual_cum_sum::Vector{Float64},
)
    empty!(model.dual_value_of_parameters)
    model.dual_value_of_parameters = zeros(model.number_of_parameters_in_model)
    for (vi_val_minus_threshold, param_dual) in enumerate(param_dual_cum_sum)
        model.dual_value_of_parameters[vi_val_minus_threshold] = param_dual
    end
    return
end

"""
    ParameterDual <: MOI.AbstractVariableAttribute

Attribute defined to get the dual values associated to parameters

# Example

```julia
MOI.get(model, POI.ParameterValue(), p)
```
"""
struct ParameterDual <: MOI.AbstractVariableAttribute end

MOI.is_set_by_optimize(::ParametricOptInterface.ParameterDual) = true

function MOI.get(model::Optimizer, ::ParameterDual, v::MOI.VariableIndex)
    if !is_additive(
        model,
        MOI.ConstraintIndex{MOI.VariableIndex,Parameter}(v.value),
    )
        error("Cannot calculate the dual of a multiplicative parameter")
    end
    return model.dual_value_of_parameters[v.value-PARAMETER_INDEX_THRESHOLD]
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintDual,
    cp::MOI.ConstraintIndex{MOI.VariableIndex,Parameter},
)
    if !is_additive(model, cp)
        error("Cannot calculate the dual of a multiplicative parameter")
    end
    return model.dual_value_of_parameters[cp.value-PARAMETER_INDEX_THRESHOLD]
end

function is_additive(model::Optimizer, cp::MOI.ConstraintIndex)
    if cp.value in model.multiplicative_parameters
        return false
    end
    return true
end


text:
"""
The base type used by all koza trees. A new, concrete type is dynamically
created for each problem, tailored to its specifics, giving it a higher level
of performance. This has the effect of essentially creating a highly optimised
interpreter for each problem.
"""
abstract KozaTree <: KozaParent

"""
Constructs the koza tree type for a given problem, based on the inputs to that
tree.
"""
function compose_tree_type(inputs::Vector{KozaInput})
  # Generate an anonymous type for the tree.
  t = anonymous_type(koza, "type <: KozaTree
    root::KozaNode

    constructor() = new()
    constructor(n::KozaNode) = new(n)
  end")

  # Implement the execution method for this tree type.
  src = join(vcat(["_t::$(t)"], ["$(a.label)::$(a.ty)" for a in inputs]), ",")
  src = "execute($(src)) = execute($(join(vcat(["_t.root"], [
  a.label for a in inputs]), ",")))"
  eval(koza, Base.parse(src))

  # Implement the cloning method for this tree type.
  return t
end

"""
Produces a deep-clone of a provided Koza Tree.
"""
clone{T <: KozaTree}(t::T) = T(clone(t.root), t)

"""
Selects a node from a provided Koza Tree at random uniformly.
"""
sample(t::KozaTree) = sample(t.root)

"""
Determines whether a given node is the root of the tree to which it belongs.
"""
isroot(n::KozaNode) = isa(n.parent, KozaTree)

"""
Returns the number of nodes in a given Koza Tree.
"""
num_nodes(t::KozaTree) = num_nodes(t.root)
height(t::KozaTree) = height(t.root)

"""
Produces a string-based description of a given Koza Tree.
"""
describe(t::KozaTree) = describe(t.root)

"""
Destructively swaps the position of two given nodes within the same Koza Tree.
"""
function swap!(n1::KozaNode, n2::KozaNode)
  t = n1.parent
  swap!(n2.parent, n2, n1)
  swap!(t, n1, n2)
end
function swap!(parent::KozaNode, child::KozaNode, with::KozaNode)
  parent.children[findfirst(parent.children, child)] = with
  with.parent = parent
end
function swap!(parent::KozaTree, child::KozaNode, with::KozaNode)
  parent.root = with
  with.parent = parent
end


text:
using Base: Unicode, SecretBuffer, SecretBuffer!, getpass, shred!

const PASSWD_PADDING =
    [0x28, 0xBF, 0x4E, 0x5E, 0x4E, 0x75, 0x8A, 0x41,
     0x64, 0x00, 0x4E, 0x56, 0xFF, 0xFA, 0x01, 0x08,
     0x2E, 0x2E, 0x00, 0xB6, 0xD0, 0x68, 0x3E, 0x80,
     0x2F, 0x0C, 0xA9, 0xFE, 0x64, 0x53, 0x69, 0x7A]

const AES_SUFFIX = [0x73, 0x41, 0x6C, 0x54] # sAlT

function UTF8ToPDFEncoding!(s::SecretBuffer)
    try
        r = SecretBuffer()
        while !eof(s)
            write(r, UnicodeToPDFEncoding(read(s, Char)))
        end
        return seekstart(r)
    finally
        shred!(s)
    end
end

struct StdSecHandler <: SecHandler
    v::Int
    version::Tuple{Int, Int}
    length::Int
    r::Int
    o::Vector{UInt8}
    u::Vector{UInt8}
    oe::Vector{UInt8}
    ue::Vector{UInt8}
    p::UInt32
    perms::Vector{UInt8}
    encMetadata::Bool
    id::Vector{UInt8}
    cf::Union{CosDict, CosNullType}
    stmf::CosName
    strf::CosName
    eff::CosName
    access::Function
    skey_path::String
    iv_path::String
    keys::Dict{CosName, Tuple{UInt32, Vector{UInt8}}}
end

# Output 32-bytes padded password
function get_padded_pw(password::SecretBuffer)
    up, l = SecretBuffer(), password.size
    if l < 32
        write(up, password)
        write(up, PASSWD_PADDING[1:32-l])
    else
        for i = 1:32
            write(up, password.data[i])
        end
    end
    return seekstart(up)
end

function algo02(h::StdSecHandler, password::SecretBuffer)
    # step a, b
    mdctx = DigestContext("md5")
    shred!(get_padded_pw(password)) do pw
        update!(mdctx, pw)
    end
    # c
    update!(mdctx, h.o)
    # d
    p = copy(reinterpret(UInt8, [h.p]))
    ENDIAN_BOM == 0x01020304 && reverse!(p)
    update!(mdctx, p)
    # e, ID
    update!(mdctx, h.id)
    # f
    h.r >= 4 && !h.encMetadata && update!(mdctx, fill(0xff, 4))
    # g
    md = close(mdctx)
    # h
    n = 5
    if h.r >= 3
        n = div(h.length, 8)
        for i = 1:50
            reset(mdctx)
            init(mdctx)
            md = md[1:n]
            update!(mdctx, md)
            md = close(mdctx)
        end
    end
    md = md[1:n]
    return md
end

function algo02a(h::StdSecHandler, pw::SecretBuffer, isowner::Bool)
    # b
    seekstart(pw)
    pw.size > 127 && (pw.size = 127)

    o, u, oe, ue = h.o, h.u, h.oe, h.ue
    ov, ovs, oks = (@view o[1:32]), (@view o[33:40]), (@view o[41:48])
    uv, uvs, uks = (@view u[1:32]), (@view u[33:40]), (@view u[41:48])
    # c
    if isowner
        or = shred!(SecretBuffer()) do pv
            write(pv, pw); write(pv, ovs); write(pv, u)
            return algo02b(h, pw, seekstart(pv), true)
        end
        or != ov && return nothing
    end
    # d, e
    ik = shred!(SecretBuffer()) do pv
        if isowner
            write(pv, pw); write(pv, oks); write(pv, u)
        else
            write(pv, pw); write(pv, uks)
        end
        return SecretBuffer!(algo02b(h, pw, seekstart(pv), isowner))
    end
    iv = SecretBuffer!(fill(0x00, 16))
    try
        cctx = CipherContext("aes_256_cbc", ik, iv, false)
        set_padding(cctx, 0)
        de = isowner ? oe : ue
        fek = SecretBuffer!(update!(cctx, de))
        write(seekend(fek), close(cctx))
        # f
        return algo13(h, seekstart(fek)) ? seekstart(fek) : nothing
    finally
        shred!(ik); shred!(iv)
    end
end

function algo02b(h::StdSecHandler, password::SecretBuffer,
                 input::SecretBuffer, isOwner::Bool)
    dgst_algo = "sha256"
    mdctx = DigestContext(dgst_algo)
    update!(mdctx, input)
    k = close(mdctx)
    if h.r > 5
        round = 0
        for i = 1:100 # The spec says by 80 most cases should converge
            # a
            e = shred!(SecretBuffer()) do k1
                write(k1, password); write(k1, k); (isOwner && write(k1, h.u))
                for i = 1:6 # 2^6 = 64
                    write(k1, k1)
                end
                # b
                key = SecretBuffer!(k[1:16]); iv = SecretBuffer!(k[17:32])
                cctx = CipherContext("aes_128_cbc", key, iv, true)
                shred!(key); shred!(iv)
                set_padding(cctx, 0)
                e = update!(cctx, k1)
                append!(e, close(cctx))
                return e
            end
            # c
            bn = BigNum(e[1:16])
            m = mod(bn, UInt(3))
            dgst_algo = m == 0 ? "sha256" : m == 1 ? "sha384" : "sha512"
            # d
            mdctx = DigestContext(dgst_algo)
            update!(mdctx, e)
            k = close(mdctx)
            i <= 64 && continue
            # e
            e[end] > (i - 33) && continue
            # f
            round = i - 1
            break
        end
        round == 0 && error("Unable to hash the password")
    end
    return k[1:32]
end

function algo03p_a_d(h::StdSecHandler, password::SecretBuffer)
    # a, b
    mdctx = DigestContext("md5")
    shred!(get_padded_pw(password)) do pw
        update!(mdctx, pw)
    end
    md = close(mdctx)
    # c
    if h.r >= 3
        for i = 1:50
            reset(mdctx); init(mdctx)
            update!(mdctx, md)
            md = close(mdctx)
        end
    end
    # d
    n = h.r == 2 ? 5 : div(h.length, 8)
    return md[1:n]
end

algo03(h::StdSecHandler, password::Vector{UInt8}) = error(E_NOT_IMPLEMENTED)

function algo04(h::StdSecHandler, password::SecretBuffer)
    key = SecretBuffer!(algo02(h, password))
    iv = SecretBuffer!(UInt8[])
    cctx = CipherContext("rc4", key, iv, true)
    shred!(iv)
    u = update!(cctx, PASSWD_PADDING)
    c = close(cctx)
    reset(cctx)
    return append!(u, c), key
end

function modkey(key::SecretBuffer, f::Int)
    mk = SecretBuffer()
    for i = 1:key.size
        write(mk, xor(key.data[i], UInt8(f)))
    end
    return mk
end

function algo05p_a_e(h::StdSecHandler, password::SecretBuffer)
    key = SecretBuffer!(algo02(h, password))
    iv = SecretBuffer!(UInt8[])
    try
        mdctx = DigestContext("md5")
        update!(mdctx, PASSWD_PADDING)
        update!(mdctx, h.id)
        md = close(mdctx)
        cctx = CipherContext("rc4", key, iv, true)
        c = update!(cctx, md)
        append!(c, close(cctx))

        for i = 1:19
            reset(cctx)
            shred!(modkey(key, i)) do tkey
                init(cctx, tkey, iv, true)
            end
            c = update!(cctx, c)
            append!(c, close(cctx))
        end
        return c, key
    finally
        shred!(iv)
    end
end

algo05(h::StdSecHandler, password::Vector{UInt8}) = error(E_NOT_IMPLEMENTED)

function algo06(h::StdSecHandler, password::SecretBuffer)
    (u, key) = h.r == 2 ? algo04(h, password) : algo05p_a_e(h, password)
    u == h.u[1:16] && return key
    shred!(key)
    return nothing
end

function algo07(h::StdSecHandler, password::SecretBuffer)
    # a
    key = SecretBuffer!(algo03p_a_d(h, password))
    # b
    iv = SecretBuffer!(UInt8[])
    try
        cctx = CipherContext("rc4", key, iv, false)
        od, count = h.o, (h.r == 2 ? 0 : 19)
        for i = count:-1:0
            reset(cctx)
            shred!(modkey(key, i)) do tkey
                init(cctx, tkey, iv, false)
            end
            od = update!(cctx, od)
            append!(od, close(cctx))
        end
        return shred!(upw -> algo06(h, upw), SecretBuffer!(od))
    finally
        shred!(iv); shred!(key)
    end
end

algo08(h::StdSecHandler, password::Vector{UInt8}) = error(E_NOT_IMPLEMENTED)
algo09(h::StdSecHandler, password::Vector{UInt8}) = error(E_NOT_IMPLEMENTED)
algo10(h::StdSecHandler, password::Vector{UInt8}) = error(E_NOT_IMPLEMENTED)

function algo11(h::StdSecHandler, password::SecretBuffer)
    u = h.u
    uv, uvs, uks = (@view u[1:32]), (@view u[33:40]), (@view u[41:48])
    uc = shred!(SecretBuffer()) do upw
        write(upw, password); write(upw, uvs)
        return algo02b(h, password, seekstart(upw), false)
    end
    return uc == u[1:32] ? algo02a(h, password, false) : nothing
end

function algo12(h::StdSecHandler, password::SecretBuffer)
    o = h.o
    ov, ovs, oks = (@view o[1:32]), (@view o[33:40]), (@view o[41:48])
    oc = shred!(SecretBuffer()) do opw
        write(opw, password); write(opw, ovs); write(opw, h.u)
        return algo02b(h, password, seekstart(opw), true)
    end
    return oc == o[1:32] ? algo02a(h, password, true) : nothing
end

function algo13(h::StdSecHandler, fek::SecretBuffer)
    perms = h.perms
    cctx = shred!(SecretBuffer!(fill(0x00, 16))) do iv
        CipherContext("aes_256_ecb", fek, iv, false)
    end
    set_padding(cctx, 0)
    dperm = update!(cctx, perms)
    append!(dperm, close(cctx))
    dperm[10:12] != b"adb" && return false
    dpermp = dperm[1:4]
    ENDIAN_BOM == 0x01020304 && reverse!(dpermp)
    return (h.p == reinterpret(UInt32, dpermp)[1])
end

function SecHandler(doc::CosDoc, access::Function)
    enc = doc.encrypt
    if get(enc, cn"Filter") === cn"Standard"
        v      = get(get(enc, cn"V", CosInt(0)))
        v >= 1 || return nothing
        return StdSecHandler(doc, v, access)
    end
    subfilter = get(enc, cn"SubFilter")
    subfilter in [cn"adbe.pkcs7.s3", cn"adbe.pkcs7.s4", cn"adbe.pkcs7.s5"] &&
        return PKISecHandler(doc, subfilter, access)
    error("Incompatible security handler used in the document")
end

function StdSecHandler(doc::CosDoc, v::Int, access::Function)
    enc = doc.encrypt
    length = get(get(enc, cn"Length", CosInt(40)))
    r      = get(get(enc, cn"R"))
    vrsn   = doc.version
    r == 1 && error("Unsupported release of security handler")
    o      = Vector{UInt8}(get(enc, cn"O"))
    u      = Vector{UInt8}(get(enc, cn"U"))
    oe     = Vector{UInt8}(get(enc, cn"OE"))
    ue     = Vector{UInt8}(get(enc, cn"UE"))
    p      = reinterpret(UInt32, Int32(get(get(enc, cn"P"))))
    perms  = Vector{UInt8}(get(enc, cn"Perms"))
    encMetadata = get(get(enc, cn"EncryptMetadata", CosTrue))
    ids = cosDocGetID(doc)
    id = Vector{UInt8}(ids[1])

    cf   = get(enc, cn"CF", CosDict())
    if r < 4
        d = CosDict()
        set!(cf, cn"StdCF", d)
        set!(d, cn"Type", cn"CryptFilter")
        set!(d, cn"CFM", cn"V2")
        set!(d, cn"AuthEvent", cn"DocOpen")
        stmf = get(enc, cn"StmF", cn"StdCF")
        strf = get(enc, cn"StrF", cn"StdCF")
    else
        stmf = get(enc, cn"StmF", cn"Identity")
        strf = get(enc, cn"StrF", cn"Identity")
    end
    eff  = get(enc, cn"EFF", stmf)
    eff !== stmf &&
        Base.@warn("Embedded file streams without Crypt filters may not decrypt")
    access === identity && (access = doUI)
    skey_path, io = get_tempfilepath()
    # Ensure the files have enough bytes in them.
    shred!(SecretBuffer!(crypto_random(rand(1000:10000)))) do s
        write(io, s)
        close(io)
    end
    iv_path, io = get_tempfilepath()
    # Ensure the files have enough bytes in them.
    shred!(SecretBuffer!(crypto_random(rand(1000:10000)))) do s
        write(io, s)
        close(io)
    end
    keys = Dict(cn"Identity" => (0xffffffff, Vector{UInt8}()))
    return StdSecHandler(v, vrsn, length, r, o, u, oe, ue,
                         p, perms, encMetadata, id,
                         cf, stmf, strf, eff, access,
                         skey_path, iv_path, keys)
end

function validateUserPW(h::StdSecHandler, cfn::CosName, password::SecretBuffer)
    cfm = get_cfm(h, cfn)
    cfm === cn"None"  && return SecretBuffer!(Vector{UInt8}(undef, 0))
    cfm === cn"V2"    && return algo06(h, password)
    cfm === cn"AESV2" && return algo06(h, password)
    cfm === cn"AESV3" && return algo11(h, password)
    error(E_INVALID_CRYPT)
end

function validateOwnerPW(h::StdSecHandler, cfn::CosName, password::SecretBuffer)
    cfm = get_cfm(h, cfn)
    cfm === cn"None"  && return SecretBuffer!(Vector{UInt8}(undef, 0))
    cfm === cn"V2"    && return algo07(h, password)
    cfm === cn"AESV2" && return algo07(h, password)
    cfm === cn"AESV3" && return algo12(h, password)
    error(E_INVALID_CRYPT)
end

doUI() = Base.getpass("Enter the password to open the document")

function validatePW(h::StdSecHandler, cfn::CosName, s::SecretBuffer)
    key = validateOwnerPW(h, cfn, seekstart(s))
    key !== nothing && return 0xffffffff, key
    key = validateUserPW(h, cfn, seekstart(s))
    key !== nothing && return h.p, key
    return nothing
end

funicode!(h::StdSecHandler, cfn::CosName, pw::SecretBuffer) =
    shred!(UTF8ToPDFEncoding!, pw)

function get_key(h::StdSecHandler, cfn::CosName, password::Function)
    vpw = shred!(pw -> validatePW(h, cfn, pw), SecretBuffer(""))
    vpw === nothing || return store_key(h, cfn, vpw)
    i = 1
    while vpw === nothing && i <= 3
        vpw = shred!(pw -> validatePW(h, cfn, pw), funicode!(h, cfn, password()))
        i += 1
    end
    vpw === nothing && error(E_INVALID_PASSWORD)
    return store_key(h, cfn, vpw)
end