text:
using Documenter
import NodeJS_16_jll
import LoggingExtras, LokiLogger, ProgressLogging, TerminalLoggers,
       SyslogLogging, LogRoller, LogCompose, Logging, Logging2, LoggingFormats,
       MiniLoggers

# "Backport" of JuliaLang/julia#40979
@eval Logging begin
    @doc """
    Alias for `LogLevel(-1000)`.
    """
    Debug
    @doc """
    Alias for `LogLevel(0)`.
    """
    Info
    @doc """
    Alias for `LogLevel(1000)`.
    """
    Warn
    @doc """
    Alias for `LogLevel(2000)`.
    """
    Error
end

# TODO: Upstream this
@eval LoggingFormats begin
    @doc """
        JSON

    Serialize log messages as JSON.
    """
    JSON

    @doc """
        LogFmt

    Format log messages as [logfmt](https://brandur.org/logfmt).
    """
    LogFmt

    @doc """
        Truncated

    Format log messages similar to `ConsoleLogger`, but truncate long messages.
    """
    Truncated
end

# Build the docs
makedocs(
    sitename = "JuliaLogging",
    pages = [
        "Home" => "index.md",
        "package-overview.md",
        "Tutorials" => [
            "tutorials/logging-basics.md",
            "tutorials/working-with-loggers.md",
            "tutorials/implement-a-new-logger.md",
        ],
        "How-to guides" => [
            "how-to/enable-debug.md",
            "how-to/log-to-file.md",
            "how-to/tee.md",
            "how-to/filter-messages.md",
            "how-to/rotate-log-files.md",
        ],
        "Reference" => [
            "reference/logcompose.md",
            "reference/logging.md",
            "reference/logging2.md",
            "reference/loggingextras.md",
            "reference/loggingformats.md",
            "reference/logroller.md",
            "reference/lokilogger.md",
            "reference/progresslogging.md",
            "reference/sysloglogging.md",
            "reference/terminalloggers.md",
            "reference/miniloggers.md",
        ],
        # "Background and discussion" => [
        # ],
    ],
    format = Documenter.HTML(
        assets = [
            "assets/julialogging.css",
            "assets/gruvbox-light-hard.css",
            "assets/gruvbox-dark-hard.css",
            "assets/favicon.ico",
        ],
        collapselevel = 1,
        ansicolor = true,
        prerender = get(ENV, "GITHUB_ACTIONS", nothing) == "true",
        node = NodeJS_16_jll.node(),
        highlightjs = joinpath(@__DIR__, "src", "assets", "julia.highlight.min.js"),
    ),
)

deploydocs(
    versions = nothing,
    push_preview = true,
    repo = "github.com/JuliaLogging/julialogging.github.io.git",
)


text:
@testset "TypeDefs" begin

using Base.Test
using DataValueArrays

x = DataValueArray(
    [1, 2, 3],
    [false, false, true]
)

y = DataValueArray(
    [
        1 2;
        3 4;
    ],
    [
        false false;
        true false;
    ],
)

@test isa(x, DataValueArray{Int, 1})
@test isa(x, DataValueVector{Int})

@test isa(y, DataValueArray{Int, 2})
@test isa(y, DataValueMatrix{Int})

end

ext: jl


text:
@testset "Examples" begin
    import Markdown: @md_str
    include("examples/quicktutorial.jl")
    include("examples/xorpruning.jl")
end


text:
# Functions to run pedigree models
function run_pedigree_model_multirecord(tbl, markerInfo, pedigree, Rvar, Gvar, Pvar; chainlength = 10000, outputfilename = "MCMC_samples", covariates = [])
    namesG11 = markerInfo[:,3] .* "_" .* markerInfo[:,1] .* "_" .* markerInfo[:,2] .* "_G11"
    namesG12 = markerInfo[:,3] .* "_" .* markerInfo[:,1] .* "_" .* markerInfo[:,2] .* "_G12"
    namesG22 = markerInfo[:,3] .* "_" .* markerInfo[:,1] .* "_" .* markerInfo[:,2] .* "_G22"
    covars = [namesG11[1]; namesG12[1]; namesG22[1]; covariates]
    printcovars = [covars; "anml_keyx"]
    eqcovars = [covars; "anml_key"; "anml_keyx"]
    model_equation = "yd_trait ="
    model_terms = join(eqcovars, " + ")
    model_equation = join([model_equation model_terms], " ")
    print("MODEL EQUATION:  ")
    println(model_equation)
    model = build_model(model_equation, Rvar)
    for term in covars
        set_covariate(model, term)
    end
    set_random(model,"anml_key",pedigree,Gvar)
    set_random(model, "anml_keyx", Pvar)
    for term in printcovars
        outputMCMCsamples(model, term)
    end
    out = runMCMC(model, tbl, output_samples_frequency = 10, chain_length = chainlength, output_folder = outputfilename)
end

function run_pedigree_model(tbl, markerInfo, pedigree, Rvar, Gvar; chainlength = 10000, outputfilename = "MCMC_samples", covariates = [])
    namesG11 = markerInfo[:,3] .* "_" .* markerInfo[:,1] .* "_" .* markerInfo[:,2] .* "_G11"
    namesG12 = markerInfo[:,3] .* "_" .* markerInfo[:,1] .* "_" .* markerInfo[:,2] .* "_G12"
    namesG22 = markerInfo[:,3] .* "_" .* markerInfo[:,1] .* "_" .* markerInfo[:,2] .* "_G22"
    covars = [namesG11[1]; namesG12[1]; namesG22[1]; covariates]
    printcovars = covars
    eqcovars = [covars; "anml_key"]
    model_equation = "yd_trait ="
    model_terms = join(eqcovars, " + ")
    model_equation = join([model_equation model_terms], " ")
    print("MODEL EQUATION:  ")
    println(model_equation)
    model = build_model(model_equation, Rvar)
    for term in covars
        set_covariate(model, term)
    end
    set_random(model,"anml_key",pedigree,Gvar)
    for term in printcovars
        outputMCMCsamples(model, term)
    end
    out = runMCMC(model, tbl, output_samples_frequency = 10, chain_length = chainlength, output_folder = outputfilename)
end

function run_pedigree_model_eQTL(tbl, pedigree, Rvar, Gvar; chainlength = 10000, freq = 10, outputfilename = "eQTL.MCMC_samples")
    pheno = String(names(tbl)[2])
    snp = String(names(tbl)[3])
    model_equation = pheno .* " = intercept + " .* snp .* " + anml_key"
    model = build_model(model_equation, Rvar)
    set_covariate(model, snp)
    set_random(model,"anml_key",pedigree, Gvar)
    outputMCMCsamples(model, snp)
    out = runMCMC(model, tbl, output_samples_frequency = freq, chain_length = chainlength, output_folder = outputfilename)
end

function run_spliceQTL(snp_data, pheno_data, pedigree, Rvar, Gvar; chainlength = 10000, freq = 10, outdir = "MCMC_samples")
    k = size(pheno_data.header, 1)
    m = size(snp_data.header, 1)
    for i in 1:k
        pheno_symbol = Symbol(pheno_data.header[i])
        phenoi = table(pheno_data.iid, pheno_data.data[:,i], names = [:anml_key, pheno_symbol], pkey=:anml_key)
        outputdir = outdir .* "/" .* pheno_data.header[i]
        run(`mkdir -p $outputdir`)
        for j in 1:m
            outputfilename = outputdir .* "/" .* snp_data.header[j]
            snp_symbol = Symbol(snp_data.header[j])
            if any(snp_data.data[:,j] .== "NA") # change NAs to mean(genotype)
                na_index = findall(snp_data.data[:,j] .== "NA")
                snp_data.data[na_index,j] = repeat([missing], size(na_index,1))
                snp_data.data[na_index,j] = repeat([mean(skipmissing(snp_data.data[:,j]))], size(na_index,1))
            end
            snpj = table(snp_data.iid, snp_data.data[:,j], names = [:anml_key, snp_symbol], pkey = :anml_key)
            tbl = join_tables(phenoi, snpj)
            meanA = mean(IndexedTables.select(tbl, snp_symbol))
            tbl = IndexedTables.transform(tbl, snp_symbol => snp_symbol => x -> x - meanA)
            df = DataFrame(tbl)
            run_pedigree_model_eQTL(df, pedigree, Rvar, Gvar; chainlength = chainlength, freq = freq, outputfilename = outputfilename)
        end
    end
end

max_stars_repo_name: egmreynolds/MANA

max_issues_repo_name: egmreynolds/MANA

max_forks_repo_name: egmreynolds/MANA


text:
"""
Library of mechanical models.
"""
module Mechanical

using ModelingToolkit

include("Rotational/Rotational.jl")

end


text:
function read_coordinates(exo_id::ExoID, num_dim::Int64, num_nodes::Int64)
    if num_dim == 1
        x_coords = Array{Float64}(undef, num_nodes)
        y_coords = Ref{Float64}(0.0)
        z_coords = Ref{Float64}(0.0)
    elseif num_dim == 2
        x_coords = Array{Float64}(undef, num_nodes)
        y_coords = Array{Float64}(undef, num_nodes)
        z_coords = Ref{Float64}(0.0)
    elseif num_dim == 3
        x_coords = Array{Float64}(undef, num_nodes)
        y_coords = Array{Float64}(undef, num_nodes)
        z_coords = Array{Float64}(undef, num_nodes)
    end
    error = ccall((:ex_get_coord, libexodus), ExodusError,
                  (ExoID, Ref{Float64}, Ref{Float64}, Ref{Float64}),
                  exo_id, x_coords, y_coords, z_coords)
    exodus_error_check(error, "read_coordinates")
    coords = Array{Float64}(undef, size(x_coords, 1), num_dim)
    if num_dim == 1
        coords = x_coords
    elseif num_dim == 2
        coords[:, 1] .= x_coords
        coords[:, 2] .= y_coords
    elseif num_dim == 3
        coords[:, 1] .= x_coords
        coords[:, 2] .= y_coords
        coords[:, 3] .= z_coords
    end
    return coords
end


text:
using BinaryBuilder

name = "OpenMPI"
version = v"4.1.1"
sources = [
    ArchiveSource("https://download.open-mpi.org/release/open-mpi/v$(version.major).$(version.minor)/openmpi-$(version).tar.gz",
                  "d80b9219e80ea1f8bcfe5ad921bd9014285c4948c5965f4156a3831e60776444"),
    DirectorySource("./bundled"),
]

script = raw"""
# Enter the funzone
cd ${WORKSPACE}/srcdir/openmpi-*

if [[ "${target}" == *-freebsd* ]]; then
    # Help compiler find `complib/cl_types.h`.
    export CPPFLAGS="-I/opt/${target}/${target}/sys-root/include/infiniband"
fi

./configure --prefix=${prefix} \
    --build=${MACHTYPE} \
    --host=${target} \
    --enable-shared=yes \
    --enable-static=no \
    --without-cs-fs \
    --enable-mpi-fortran=usempif08 \
    --with-cross=${WORKSPACE}/srcdir/${target}

# Build the library
make -j${nproc}

# Install the library
make install

install_license LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.
#platforms = supported_platforms()
platforms = filter(p -> !Sys.iswindows(p) && !(arch(p) == "armv6l" && libc(p) == "glibc"), supported_platforms(; experimental=true))
platforms = expand_gfortran_versions(platforms)

products = [
    LibraryProduct("libmpi", :libmpi)
    ExecutableProduct("mpiexec", :mpiexec)
]

dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat="1.6", preferred_gcc_version=v"5")

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
using Makie: PointTrans, xyz_boundingbox, apply_transform

@testset "Basic transforms" begin
    function fpoint2(x::Point2)
        return Point2f(x[1] + 10, x[2] - 77)
    end

    function fpoint3(x::Point3)
        return Point3f(x[1] + 10, x[2] - 77, x[3] /  4)
    end
    trans2 = PointTrans{2}(fpoint2)
    trans3 = PointTrans{3}(fpoint3)
    points2 = [Point2f(0, 0), Point2f(0, 1)]
    bb = xyz_boundingbox(trans2, points2)
    @test bb == Rect(Vec3f(10, -77, 0), Vec3f(0, 1, 0))
    bb = xyz_boundingbox(trans3, points2)
    @test bb == Rect(Vec3f(10, -77, 0), Vec3f(0, 1, 0))

    points3 = [Point3f(0, 0, 4), Point3f(0, 1, -8)]
    bb = xyz_boundingbox(trans2, points3)
    @test bb == Rect(Vec3f(10, -77, -8), Vec3f(0, 1, 12))
    bb = xyz_boundingbox(trans3, points3)
    @test bb == Rect(Vec3f(10, -77, -2.0), Vec3f(0, 1, 3.0))

    @test apply_transform(trans2, points2) == fpoint2.(points2)
    @test apply_transform(trans3, points3) == fpoint3.(points3)

    @test_throws ErrorException PointTrans{2}(x::Int -> x)
    @test_throws ErrorException PointTrans{3}(x::Int -> x)
end

@testset "Tuple and identity transforms" begin
    t1 = sqrt
    t2 = (sqrt, log)
    t3 = (sqrt, log, log10)

    p2 = Point(2.0, 5.0)
    p3 = Point(2.0, 5.0, 4.0)

    @test apply_transform(identity, p2) == p2
    @test apply_transform(identity, p3) == p3

    @test apply_transform(t1, p2) == Point(sqrt(2.0), sqrt(5.0))
    @test apply_transform(t1, p3) == Point(sqrt(2.0), sqrt(5.0), sqrt(4.0))

    @test apply_transform(t2, p2) == Point2f(sqrt(2.0), log(5.0))
    @test apply_transform(t2, p3) == Point3f(sqrt(2.0), log(5.0), 4.0)

    @test apply_transform(t3, p3) == Point3f(sqrt(2.0), log(5.0), log10(4.0))

    i2 = (identity, identity)
    i3 = (identity, identity, identity)
    @test apply_transform(i2, p2) == p2
    @test apply_transform(i3, p3) == p3

    # test that identity gives back exact same arrays without copying
    p2s = Point2f[(1, 2), (3, 4)]
    @test apply_transform(identity, p2s) === p2s
    @test apply_transform(i2, p2s) === p2s
    @test apply_transform(i3, p2s) === p2s

    p3s = Point3f[(1, 2, 3), (3, 4, 5)]
    @test apply_transform(identity, p3s) === p3s
    @test apply_transform(i2, p3s) === p3s
    @test apply_transform(i3, p3s) === p3s

    @test apply_transform(identity, 1) == 1
    @test apply_transform(i2, 1) == 1
    @test apply_transform(i3, 1) == 1

    @test apply_transform(identity, 1..2) == 1..2
    @test apply_transform(i2, 1..2) == 1..2
    @test apply_transform(i3, 1..2) == 1..2

    pa = Point2f(1, 2)
    pb = Point2f(3, 4)
    r2 = Rect2f(pa, pb .- pa)
    @test apply_transform(t1, r2) == Rect2f(apply_transform(t1, pa), apply_transform(t1, pb) .- apply_transform(t1, pa) )
end


text:
"""
    findMinimum(p::Int, q::Int)
    findMinimum(p::Int, q::Int; optimizer::Symbol) where T <: Real

Finds a minimum of modelFreeEnergy in the box ```[-1,1]^modelDimension(q)``` at parameter (p,q).

Returns a Tuple with (p, q, temperature, optimizer, minimum configuration, starting configuration, timestamp , timespan of optimization, status).

Optional arguments:
-   ```optimizer```: ```:gradient``` for gradient descent, ```:ipopt``` for interior point method, ```:cplex``` for cplex optimizer. Default value: ```:gradient```.

```:gradient``` finds a local minimum.
Optional arguments for ```:gradient```.
-   ```startingCondition```: initialization for the Ipopt optimizer. Default value: random point in the box ```[-1,1]^modelDimension(q)```;
-   ```temperature```: temperature of the model. Default value: ```0```;
-   ```step```: step-size of gradient descent update. Default value: ```1e-3```;
-   ```maxIter```: maximum number of gradient descent steps allowed. Default value: ```10^6```;
-   ```threshold```: threshold for the norm of the gradient. If ```norm(grad) < threshold```, the gradient is considered null. Default value: 1e-3.

**CURRENTLY BROKEN!!**
```:ipopt``` finds a local minimum.
Optional arguments for ```:ipopt```:
-   ```verbose```: if ```true``` activates JuMP verbose mode. Default value: ```false```;
-   ```startingCondition```: initialization for the Ipopt optimizer. Default value: random point in the box ```[-1,1]^modelDimension(q)```;
-   ```temperature```: temperature of the model. Default value: ```0```.

```:cplex``` finds the global minimum at zero temperature.
Optional arguments for ```:cplex```:
-   ```verbose```: if ```true``` activates JuMP verbose mode. Default value: ```false```;
-   ```size```: to solve the problem, magnetizations are quantized in 1/size steps. Default value: ```100```.

"""
function modelFindMinimum(p::Int64,q::Int64;
                     optimizer::Symbol =:gradient,
                     kwargs...
                     ) where T <: Real

    id = Dates.format(now(), "yymmddHHMMSS")
    result = Tuple{Vector{Float64},Vector{Float64}}

    if optimizer == :gradient
        threshold = 1e-3
        if haskey(kwargs, :threshold)
            threshold = kwargs[:threshold]
        end
        maxIter = 1000000
        if haskey(kwargs, :maxIter)
            maxIter = kwargs[:maxIter]
        end
        step = 1e-4
        if haskey(kwargs, :step)
            step = kwargs[:step]
        end
        temperature = 0
        if haskey(kwargs, :temperature)
            temperature = kwargs[:temperature]
        end
        startingCondition = rand(Float64,modelDimension(q)) .|> u -> 2*u-1
        if haskey(kwargs, :startingCondition)
           startingCondition = kwargs[:startingCondition]
        end

        res = @timed localMinimum_gradientDescent(p,q,
                                                  startingCondition=startingCondition,
                                                  temperature=temperature,
                                                  threshold=threshold,
                                                  maxIter=maxIter,
                                                  step=step
                                                 )
        minimum,starting,status = res[:value]
        time = res[:time]

    elseif optimizer == :ipopt
        verbose = false
        if haskey(kwargs, :verbose)
            verbose = kwargs[:verbose]
        end
        temperature = 0
        if haskey(kwargs, :temperature)
            temperature = kwargs[:temperature]
        end
        startingCondition = rand(Float64,modelDimension(q)) .|> u -> 2*u-1
        if haskey(kwargs, :startingCondition)
           startingCondition = kwargs[:startingCondition]
        end

        minimum, starting, time, status = localMinimum_interiorPoint(p,q,
                                                                     startingCondition=startingCondition,
                                                                     temperature=temperature,
                                                                     verbose=verbose
                                                                    )
    elseif optimizer == :cplex
        size = 100
        if haskey(kwargs, :size)
            size = kwargs[:size]
        end
        verbose = false
        if haskey(kwargs, :verbose)
            verbose = kwargs[:verbose]
        end
        temperature = 0

        minimum, time, status = globalMinimum(p,q,verbose=verbose, size=size)
        starting = []

    else
        println("Wrong optimizer!")
        minimum = []
        starting = []
        time = -1
    end

    return ( p, q, temperature, string(optimizer), minimum, starting, id, time, status )
end

"""
    globalMinimum(p::Int,q::Int)
    globalMinimum(p::Int,q::Int; verbose::Bool, size::Int)

Finds the global minimum of modelFreeEnergy in the box ```[-1,1]^modelDimension(q)``` at zero temperature.

Returns a Tuple with (optimal magnetizations, time elapsed, JuMP status).

Optional arguments:
-   ```verbose```: if ```true``` activates JuMP verbose mode. Default value: ```false```;
-   ```size```: to solve the problem, magnetizations are quantized in 1/size steps. Default value: ```100```.
"""
function globalMinimum(p::Int,q::Int;
                       verbose::Bool=false,
                       size::Int=100)

    q2 = modelDimension(q)

    # initialize quadratic interaction matrix
    J = modelInteraction(p,q)

    # set optimizer to find global minima
    model = Model(CPLEX.Optimizer)
    set_optimizer_attribute(model, "CPXPARAM_OptimalityTarget", 3)

    # set verbosity of optimizer
    if !verbose
        set_silent(model)
    end

    # create quantized variables
    @variable(model, -size <= M[i=1:q2] <= size, Int)
    vars= JuMP.all_variables(model)

    # set quadratic objective function
    @objective(model, Min, dot(vars,J,vars))

    # optimize
    optimize!(model)

    # a priori there may be more optimal solutions
    # here we ignore this possibility

    GSconfiguration = (value.(vars))./size
    solveTime = solve_time(model)
    status = primal_status(model)                 |> string

    return (GSconfiguration,solveTime,status)
end

"""
    localMinimum_interiorPoint(p::Int,q::Int)
    localMinimum_interiorPoint(p::Int,q::Int; verbose::Bool, startingCondition::Vector{T}, temperature::Real) where T <: Real

CURRENTLY BROKEN AS IT EVALUATES THE ENTROPY OUTSIDE OF THE DOMAIN BY VERY SMALL AMOUNTS

Finds a local minimum of modelFreeEnergy in the box ```[-1,1]^modelDimension(q)``` at temperature t using the interior point method provided by Ipopt.

Returns a Tuple with (optimal magnetizations, starting condition, time elapsed, JuMP status).

Optional arguments:
-   ```verbose```: if ```true``` activates JuMP verbose mode. Default value: ```false```;
-   ```startingCondition```: initialization for the Ipopt optimizer. Default value: random point in the box ```[-1,1]^modelDimension(q)```;
-   ```temperature```: temperature of the model. Default value: ```0```.
"""
function localMinimum_interiorPoint(p::Int64, q::Int64;
                                    startingCondition::Vector{T} = rand(Float64,modelDimension(q)) .|> u -> 2*u-1,
                                    temperature::Real=0,
                                    verbose::Bool=false) where T <: Real

    q2 = modelDimension(q)
    if length(startingCondition) != q2
        throw(ArgumentError("Length of startingCondition incompatible with q!"))
    end

    # initialize quadratic interaction matrix
    J = modelInteraction(p,q)

    # set optimizer
    model = Model(Ipopt.Optimizer)

    register(model, :magnetizationEntropy, 1, magnetizationEntropy, autodiff=true)

    # set verbosity of optimizer
    if !verbose
        set_silent(model)
    end

    # create continuous variables and set starting condition
    @variable(model, -1 <= M[i=1:q2] <= 1)
    for i in 1:q2
        set_start_value(M[i],startingCondition[i])
    end
    vars= JuMP.all_variables(model)

    # initialize quadratic interaction matrix
    # and compute energy and entropy
    energy = @NLexpression(model,
                           sum(vars[i]*J[i,j]*vars[j] for i in 1:q2, j in 1:q2)
                          )
    entropy = @NLexpression(model,
                            2*magnetizationEntropy(vars[1])
                            + 4*sum(
                                    magnetizationEntropy(vars[i])
                                    for i in 2:q2-1)
                            + ( q%2==0 ? 2 : 4 )*magnetizationEntropy(vars[q2])
                           )

    # set objective function
    @NLobjective(model, Min, energy-temperature*entropy)

    # optimize
    optimize!(model)

    GSconfiguration = (value.(vars))
    solveTime = solve_time(model)
    status = primal_status(model) |> string

    return (GSconfiguration,startingCondition,solveTime,status)
end

"""
    localMinimum_gradientDescent(p::Int,q::Int)
    localMinimum_gradientDescent(p::Int,q::Int; startingCondition::Vector{T}, temperature::Real, step::Real, maxIter::Int, threshold::Real) where T <: Real

Finds a local minimum of modelFreeEnergy in the box ```[-1,1]^modelDimension(q)``` at temperature t using a custom gradient descent method that respects the box boundary conditions.

Returns a Tuple with (optimal magnetizations, starting condition, status).
Status can be:
-   ```"good"```: optimization ended due to null gradient;
-   ```"max iter reached": no convergence was reached before the maximum number of allower iterations;
-   ```"stuck in position"```: the descent is stuck in a position while the gradient is non null.
Only ```"good"``` should be considered as a successful optimization.

Optional arguments:
-   ```startingCondition```: initialization for the Ipopt optimizer. Default value: random point in the box ```[-1,1]^modelDimension(q)```;
-   ```temperature```: temperature of the model. Default value: ```0```;
-   ```step```: step-size of gradient descent update. Default value: ```1e-3```;
-   ```maxIter```: maximum number of gradient descent steps allowed. Default value: ```10^6```;
-   ```threshold```: threshold for the norm of the gradient. If ```norm(grad) < threshold```, the gradient is considered null. Default value: 1e-3.
"""
function localMinimum_gradientDescent(p::Int,q::Int;
                                      startingCondition::Vector{T} = rand(Float64,modelDimension(q)) .|> u -> 2*u-1,
                                      temperature::Real =0.,
                                      threshold::Real = 1e-3,
                                      maxIter::Int = 1000000,
                                      step::Real = 1e-3
                                     ) where T <: Real

    q2 = modelDimension(q)
    if length(startingCondition) != q2
        throw(ArgumentError("Length of startingCondition incompatible with q!"))
    end
    # set starting point and starting gradient
    old_position = startingCondition.+0.1 # add 0.1 to make it different from actual_postition before starting the while cycle
    actual_position = startingCondition
    grad = gradModelFreeEnergy(actual_position,p,q,temperature)
    iter = 0

    # while gradient not too small, keep following the steepest descent
    while norm(grad) > threshold && iter < maxIter && actual_position != old_position
        iter += 1
        # update position with gradient
        old_position = actual_position
        actual_position -= step * grad
        # make sure that the new position is in the box
        # otherwise push it back
        for i in 1:q2
            if actual_position[i] > 1
                actual_position[i] = 1
            elseif actual_position[i] < -1
                actual_position[i] = -1
            end
        end
        # compute new gradient
        grad = gradModelFreeEnergy(actual_position,p,q,temperature)
        # if gradient descent direction (i.e. -gradient) points out of the box, put that component to zero
        for i in 1:q2
            if ( actual_position[i] == 1 && grad[i] < 0 ) || ( actual_position[i] == -1 && grad[i] > 0 )
                grad[i] = 0
            end
        end
    end

    if iter == maxIter
        status = "max iter reached"
    elseif actual_position == old_position
        status = "stuck in position"
    else
        status = "good"
    end

    return (actual_position, startingCondition, status)
end


text:
### Abstract parent types for Metropolis-Hastings samplers
abstract AbstractMetropolisHastingsSampler <: AbstractSampler
abstract AbstractMetropolisHastingsSamplerState <: AbstractSamplerState{BaseSample}

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

### Type holding the parameters for a scalable Metropolis-Hastings sampler (i.e., any proposal density with a scale parameter)
immutable ScalableMetropolisHastingsSampler{T<:AbstractFloat,M<:AbstractArray} <: AbstractMetropolisHastingsSampler
    density::Symbol
    initialscalefactor::T
    scaleparameters::M
    extraargs::Tuple
    ScalableMetropolisHastingsSampler(d::Symbol,i::T,p::AbstractArray{T},args...) = new(d,i,p,args)
end

###Factory functions
@inline _sampler(::Type{Val{:mh}},den::Symbol,i::AbstractFloat,m::AbstractArray,args...) = ScalableMetropolisHastingsSampler{typeof(i),typeof(m)}(den,i,m,args...)
@inline _sampler(::Type{Val{:mh}},den::Symbol,i::AbstractFloat,d::Integer,args...) = _sampler(Val{:mh},den,i,ones(typeof(i),d),args...)

###Size functions
@inline numparas(sampler::ScalableMetropolisHastingsSampler) = size(sampler.scaleparameters,1)

###Name functions (used in println for generic AbstractSampler)
@inline samplername(sampler::ScalableMetropolisHastingsSampler) = "Scalable Metropolis-Hastings"

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

### Type holding the state of the Markov Chain for a Generalized Metropolis-Hastings sampler
type MetropolisHastingsSamplerState{D<:AbstractProposalDensity} <: AbstractMetropolisHastingsSamplerState
    density::D
    from::BaseSample
    proposals::BaseSample
    acceptance::Vector
    MetropolisHastingsSamplerState(d::D,f::BaseSample,p::BaseSample,a::Vector) = new(d,f,p,a)
end

###Factory functions
@inline function _samplerstate{N<:Number,T<:AbstractFloat}(s::ScalableMetropolisHastingsSampler,nsamples::Integer,::Type{N},::Type{T},::Bool)
    nparas = numparas(s)
    f = _samples(Val{:base},nparas,1,N,T)
    p = _samples(Val{:base},nparas,nsamples,N,T)
    d = _density(Val{s.density},zeros(N,nparas),deepcopy(s.scaleparameters),deepcopy(s.extraargs)...)
    s.initialscalefactor!=one(T)?scale!(d,s.initialscalefactor):nothing #rescale the distribution for a scale factor different from 1.0
    MetropolisHastingsSamplerState{typeof(d)}(d,f,p,zeros(T,nsamples))
end

###Name functions (used in println for AbstractSamplerState)
@inline samplerstatename{D<:SymmetricDensity}(state::MetropolisHastingsSamplerState{D}) = "Metropolis"
@inline samplerstatename{D<:ASymmetricDensity}(state::MetropolisHastingsSamplerState{D}) = "Metropolis-Hastings"

###update the sampler state, called from SMHRunner
@inline function prepare!(state::MetropolisHastingsSamplerState,updatefrom::Bool=false)
    if updatefrom
        #set the current point around which we are sampling
        setfrom!(state,state.proposals)
        #condition the density on the current point
        condition!(state.density,state.from.values)
    end
    #return the state
    state
end

###update the auxiliary samplerstate with information from the indicator samplerstate
@inline function prepareauxiliary!(indicatorstate::MetropolisHastingsSamplerState,auxiliarystate::MetropolisHastingsSamplerState)
    #set the auxiliary point as the from point in the auxiliarystate
    setfrom!(auxiliarystate,indicatorstate.proposals)
    #condition the density on the current point
    condition!(auxiliarystate.density,auxiliarystate.from.values)
    #return the auxiliary state
    auxiliarystate
end

###when the "from" field does not need updating
@inline prepareindicator!(indicatorstate::MetropolisHastingsSamplerState) = indicatorstate

###when the "from" field has been updated
@inline function prepareindicator!(indicatorstate::MetropolisHastingsSamplerState,auxiliarystate::MetropolisHastingsSamplerState,i::Int)
    #set the current point around which we are sampling
    setfrom!(indicatorstate,auxiliarystate.proposals,i)
    #condition the density on the current point
    condition!(indicatorstate.density,indicatorstate.from.values)
    #return the indicator state
    indicatorstate
end

###propose new samples from the current point
@inline function propose!(state::MetropolisHastingsSamplerState)
    #propose n values from the current point
    propose!(state.density,state.proposals.values)
    #return the proposals
    state
end

###for a symmetric proposal density, no adjustment of the acceptance rate is needed
function acceptance!{D<:SymmetricDensity}(state::MetropolisHastingsSamplerState{D})
    fl = state.from.loglikelihood[1] + state.from.logprior[1]
    for i=1:numsamples(state.proposals)
        @inbounds state.acceptance[i] = state.proposals.loglikelihood[i] + state.proposals.logprior[i] - fl
    end
    state.acceptance
end

###for an asymmetric proposal density, an adjustment needs to be performed
###based on the difference between K(from,val)/K(val,from)
function acceptance!{D<:ASymmetricDensity}(state::MetropolisHastingsSamplerState{D})
    #calculate the log-probability for each of the proposed points K(from,val)
    logprobability!(state.acceptance,state.density,state.proposals.values)
    #now calculate the inverse probability K(val,from)
    temp1 = similar(state.acceptance,1)
    fl = state.from.loglikelihood[1] + state.from.logprior[1]
    for i = 1:numsamples(state.proposals)
        #condition the proposal density on each proposed point
        @inbounds condition!(state.density,state.proposals.values[:,i])
        #calculate the logprobability K(val,from)
        logprobability!(temp1,state.density,state.from.values)
        #adjust the result accordingly
        @inbounds state.acceptance[i] += (state.proposals.loglikelihood[i] + state.proposals.logprior[i] - temp1[1] - fl)
    end
    #restore the previous state of the proposal density
    condition!(state.density,state.from.values)
    #return the result
    state.acceptance
end

@inline function tune!(state::MetropolisHastingsSamplerState,scalefactor::AbstractFloat)
    scale!(state.density,scalefactor)
    state
end


text:
using RawFile

sizes = [(100,),(75,100),(10,20,50)]
fname = "testfile.raw"
for t in [UInt8,Int32,Int64,Float32,Float64,ComplexF16,ComplexF32]
    for s in sizes
        d = rand(t,s)
        saveraw(d,fname)
        dd = readraw(fname)
        dd_batch = t[]
        readraw(fname,10) do d
            append!(dd_batch,d[:])
        end
        dd_batch = reshape(dd_batch,s)

        @assert d == dd
        @assert d == dd_batch
        @assert s == rawsize(fname)

        saveraw(fname) do f
           for i=0:Int(s[end]/10)-1
                write(f,view(d,[Colon() for j=1:ndims(d)-1]...,i*10+1:i*10+10))
            end
        end

        dd = readraw(fname)
        @info size(d)
        @info size(dd)
        @info d[70:80]
        @info dd[70:80]
        @assert d == dd

        saveraw(view(d,[Colon() for j=1:ndims(d)-1]...,1:10),fname)
        for i=1:Int(s[end]/10)-1
            appendraw(view(d,[Colon() for j=1:ndims(d)-1]...,i*10+1:i*10+10),fname)
        end
        dd = readraw(fname)
        @assert d == dd

        if length(s)==1
            saveraw(view(d,1:10),fname)
            for i=11:s[1]
                appendraw(d[i],fname)
            end
            dd = readraw(fname)
            @assert d == dd
        end

        isfile(fname) && rm(fname)
    end
end


text:
@reexport module _Series

include("binarysplitting.jl")

include("reciprocal.jl")
include("sqrt.jl")
include("arctan.jl")
include("pi.jl")
include("euler.jl")
include("ln2.jl")

end


text:
import Base.reinterpret

type DAF
    filename::ASCIIString
    array::Vector{UInt8}
    little::Bool
    id::ASCIIString
    nd::Int32
    ni::Int32
    name::ASCIIString
    first::Int32
    last::Int32
    ss::Int32
    nc::Int32
end

const RECORD_LENGTH = 1024

function getrecord(array, n)
    array[1+RECORD_LENGTH*(n-1):RECORD_LENGTH*n]
end

getrecord(daf::DAF, n) = getrecord(daf.array, n)

function DAF(filename)
    if !isfile(filename)
        error("'$filename' does not exist.")
    end
    array = Mmap.mmap(filename, Vector{UInt8})
    fr = readfilerecord(getrecord(array, 1))
    DAF(filename, array, fr...)
end

function reinterpret(t::Type, array::Vector{UInt8}, littleendian::Bool)
    out = reinterpret(t, array)
    if littleendian
        return out
    else
        return map!(ntoh, out)
    end
end

function readint(record, address, littleendian=true)
    reinterpret(Int32, record[address+1:address+4], littleendian)[1]
end

function readascii(record, address, len)
    rstrip(ASCIIString(record[address+1:address+len]))
end

function islittleendian(record, legacy)
    if legacy
        nd = readint(record, 8)
        if nd == 2
            return true
        end
        nd = readint(record, 8, false)
        if nd == 2
            return false
        else
            error("Endianess could not be detected.")
        end
    else
        endianess = readascii(record, 88, 8)
        return endianess == "LTL-IEEE"
    end
end

function readfilerecord(record)
    ftpstr = b"FTPSTR:\r:\n:\r\n:\r\x00:\x81:\x10\xce:ENDFTP"
    id = readascii(record,0,8)
    legacy = id == "NAIF/DAF"
    little = islittleendian(record, legacy)
    nd = readint(record, 8, little)
    ni = readint(record, 12, little)
    name = readascii(record, 16, 60)
    first = readint(record, 76, little)
    last = readint(record, 80, little)
    ss = nd + div(ni+1, 2)
    nc = 8*ss
    if ~legacy
        if ftpstr != record[699+1:699+28]
            error("This DAF file is damaged.")
        end
    end
    return little, id, nd, ni, name, first, last, ss, nc
end

function summaryheader(record, little)
    next, _, nsum = round(Int32, reinterpret(Float64, record[1:24], little))
    return next, nsum
end

function addsummaries!(summaries, record, names, nsum, ss, nc)
    for i = 1:nsum
        push!(summaries, (
            rstrip(ASCIIString(names[1+(i-1)*nc:i*nc])),
            record[25+(i-1)*ss*8:25+i*ss*8]
        ))
    end
end

function getsummaries(daf::DAF)
    summaries = Tuple{ASCIIString,Vector{UInt8}}[]
    record = getrecord(daf, daf.first)
    names = getrecord(daf, daf.first+1)
    next, nsum = summaryheader(record, daf.little)
    addsummaries!(summaries, record, names, nsum, daf.ss, daf.nc)
    while next != 0
        record = getrecord(daf, next)
        names = getrecord(daf, next+1)
        next, nsum = summaryheader(record, daf.little)
        addsummaries!(summaries, record, names, nsum, daf.ss, daf.nc)
    end
    return summaries
end


text:
# ## Data

n = 9
k = 6;

# ## Model

using Turing

@model function globe_toss(n, k)
    θ ~ Beta(1, 1)
    k ~ Binomial(n, θ)
    return k, θ
end;

# ## Output

using Random

Random.seed!(1)
chains = sample(globe_toss(n, k), NUTS(0.65), 1000)

# \defaultoutput{}


text:
using Dates

struct Route
    lon1::Float64
    lon2::Float64
    lat1::Float64
    lat2::Float64
    x_nodes::Int64
    y_nodes::Int64
end

"""Identify the shortest path given arrays of locations and an array of the earliest time at each point."""
function shortest_path(indx, pindx, sp)
    ix = findfirst(isequal(sp[end]), indx)
    if ix == nothing
        return sp
    end
    pix = pindx[ix]
    append!(sp, pix)
    if pix == 0.0
        return Array(sp)
    else
        return shortest_path(indx, pindx, sp)
    end
end

"""Get locations of the shortest path"""
function get_locs(indx, sp, x_locs, y_locs)
    X = []
    Y = []
    for k in sp[1:end-1]
        idx = findfirst(isequal(k), indx)
        append!(X, x_locs[idx])
        append!(Y, y_locs[idx])
    end
    return hcat(X, Y)
end

"Find index of initial start time in array"
@inline function time_to_index(time, time_values)
    idx = findfirst(isequal(time), time_values)
    if isequal(idx, nothing) == true
        return size(time_values)[1]
    else
        return idx
    end
end

"""Convert hours in float to an index value to be added to the start index. This needs to be changed based on the weather data used."""
@inline function convert_time(time::Float64)
    mm = round(time)
    return Dates.Hour(mm)
end

"Time dependent shortest path."
function route_solve(route::Route, performance::Performance, start_time::DateTime, times, x, y, wisp, widi, wadi, wahi, cusp, cudi)
    start_time_idx = time_to_index(start_time, times)
    earliest_times = fill(Inf, size(x))
    prev_node = zero(x)
    node_indices = reshape(1:length(x), size(x))
    arrival_time = Inf
    final_node = 0
    earliest_times = fill(Inf, size(x))
    idx_range = size(x)[2]
    idy_range = size(x)[1]
    @simd for idx in 1:idx_range
        d, b = haversine(route.lon1, route.lat1, x[1, idx], y[1, idx])
        wd_int = widi[start_time_idx, idx, 1]
        ws_int = wisp[start_time_idx, idx, 1]
        cs_int = cusp[start_time_idx, idx, 1]
        cd_int = cudi[start_time_idx, idx, 1]
        wadi_int = wadi[start_time_idx, idx, 1]
        if isnan(wadi_int) == true
            wadi_int = cd_int
        end
        wahi_int = wahi[start_time_idx, idx, 1]
        if isnan(wahi_int) == true
            wahi_int = 0.0
        end
        speed = cost_function(performance, cd_int, cs_int, wd_int, ws_int,
                              wadi_int, wahi_int, b)
        if speed >= 0.0
            earliest_times[1, idx] = d/speed
        end
    end
    for idy in 1:idy_range-1
        for idx1 in 1:idx_range
            if isinf(earliest_times[idy, idx1]) == false
                t = start_time + convert_time(earliest_times[idy, idx1])
                t_idx = time_to_index(t, times)
                wd_int = widi[t_idx, idx1, idy]
                ws_int = wisp[t_idx, idx1, idy]
                wadi_int = wadi[t_idx, idx1, idy]
                wahi_int = wahi[t_idx, idx1, idy]
                cs_int = cusp[t_idx, idx1, idy]
                cd_int = cudi[t_idx, idx1, idy]
                @simd for idx2 in 1:idx_range
                    d, b = haversine(x[idy, idx1], y[idy, idx1],
                                        x[idy+1, idx2], y[idy+1, idx2])
                    speed = cost_function(performance, cd_int, cs_int,
                                                    wd_int, ws_int, wadi_int, wahi_int, b)
                    if speed >= 0.0
                        tt = earliest_times[idy, idx1] + d/speed
                        if earliest_times[idy+1, idx2] > tt
                            earliest_times[idy+1, idx2] = tt
                            prev_node[idy+1, idx2] = node_indices[idy, idx1]
                        end
                    end
                end
            end
        end
    end

    @simd for idx in 1:idx_range
        if isinf(earliest_times[end, idx]) == false
            d, b = haversine(x[end, idx], y[end, idx], route.lon2, route.lat2)
            t = start_time + convert_time(earliest_times[end, idx])
            t_idx = time_to_index(t, times)
            wd_int = widi[t_idx, idx, end]
            ws_int = wisp[t_idx, idx, end]
            wadi_int = wadi[t_idx, idx, end]
            wahi_int = wahi[t_idx, idx, end]
            cs_int = cusp[t_idx, idx, end]
            cd_int = cudi[t_idx, idx, end]
            speed = cost_function(performance, cd_int, cs_int,
                                wd_int, ws_int, wadi_int, wahi_int, b)
            if speed >= 0.0
                tt = earliest_times[end, idx] + d/speed
                if arrival_time > tt
                    arrival_time = tt
                    final_node = node_indices[end, idx]
                end
            end
        end
    end
    sp = shortest_path(node_indices, prev_node, [final_node])
    locs = get_locs(node_indices, sp, x, y)
    x_path = vcat([route.lon1], locs[:, 1], [route.lon2])
    y_path = vcat([route.lat1], locs[:, 2], [route.lat2])
    locs = hcat(x_path, y_path)
    return arrival_time, locs, earliest_times
end


text:
#!/usr/bin/env julia

# CUDAnative.jl version, using the fancy profiler instead of manual events

using CUDAdrv, CUDAnative

function kernel_dummy(ptr)
    Base.pointerset(ptr, Float32(0), Int(blockIdx().x), 8)
    return nothing
end

const len = 1000

const ITERATIONS = 5000

function main()
    dev = CuDevice(0)
    ctx = CuContext(dev)

    cpu_time = Vector{Float64}(ITERATIONS)

    gpu_arr = CuArray{Float32}(len)
    for i in 1:ITERATIONS
        i == ITERATIONS-4 && CUDAdrv.start_profiler()

        cpu_tic = time_ns()
        CUDAnative.@profile begin
            @cuda (len,1) kernel_dummy(pointer(gpu_arr))
        end
        cpu_toc = time_ns()

        cpu_time[i] = (cpu_toc-cpu_tic)/1000
    end
    CUDAdrv.stop_profiler()

    @printf("CPU time: %.2fus\n", median(cpu_time))
    CUDAnative.Profile.print()

    destroy(ctx)
end

main()


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Bus <: Topology
        number::Int
        name::String
        bustype::Union{Nothing, BusTypes.BusType}
        angle::Union{Nothing, Float64}
        magnitude::Union{Nothing, Float64}
        voltage_limits::Union{Nothing, Min_Max}
        base_voltage::Union{Nothing, Float64}
        area::Union{Nothing, Area}
        load_zone::Union{Nothing, LoadZone}
        ext::Dict{String, Any}
        internal::InfrastructureSystemsInternal
    end

A power-system bus.

# Arguments
- `number::Int`: number associated with the bus
- `name::String`: the name of the bus
- `bustype::Union{Nothing, BusTypes.BusType}`: bus type
- `angle::Union{Nothing, Float64}`: angle of the bus in radians, validation range: `(-1.571, 1.571)`, action if invalid: `error`
- `magnitude::Union{Nothing, Float64}`: voltage as a multiple of basevoltage, validation range: `voltage_limits`, action if invalid: `warn`
- `voltage_limits::Union{Nothing, Min_Max}`: limits on the voltage variation as multiples of basevoltage
- `base_voltage::Union{Nothing, Float64}`: the base voltage in kV, validation range: `(0, nothing)`, action if invalid: `error`
- `area::Union{Nothing, Area}`: the area containing the bus
- `load_zone::Union{Nothing, LoadZone}`: the load zone containing the bus
- `ext::Dict{String, Any}`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct Bus <: Topology
    "number associated with the bus"
    number::Int
    "the name of the bus"
    name::String
    "bus type"
    bustype::Union{Nothing, BusTypes.BusType}
    "angle of the bus in radians"
    angle::Union{Nothing, Float64}
    "voltage as a multiple of basevoltage"
    magnitude::Union{Nothing, Float64}
    "limits on the voltage variation as multiples of basevoltage"
    voltage_limits::Union{Nothing, Min_Max}
    "the base voltage in kV"
    base_voltage::Union{Nothing, Float64}
    "the area containing the bus"
    area::Union{Nothing, Area}
    "the load zone containing the bus"
    load_zone::Union{Nothing, LoadZone}
    ext::Dict{String, Any}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal

    function Bus(number, name, bustype, angle, magnitude, voltage_limits, base_voltage, area, load_zone, ext, internal, )
        (number, name, bustype, angle, magnitude, voltage_limits, base_voltage, area, load_zone, ext, internal, ) = check_bus_params(
            number,
            name,
            bustype,
            angle,
            magnitude,
            voltage_limits,
            base_voltage,
            area,
            load_zone,
            ext,
            internal,
        )
        new(number, name, bustype, angle, magnitude, voltage_limits, base_voltage, area, load_zone, ext, internal, )
    end
end

function Bus(number, name, bustype, angle, magnitude, voltage_limits, base_voltage, area=nothing, load_zone=nothing, ext=Dict{String, Any}(), )
    Bus(number, name, bustype, angle, magnitude, voltage_limits, base_voltage, area, load_zone, ext, InfrastructureSystemsInternal(), )
end

function Bus(; number, name, bustype, angle, magnitude, voltage_limits, base_voltage, area=nothing, load_zone=nothing, ext=Dict{String, Any}(), internal=InfrastructureSystemsInternal(), )
    Bus(number, name, bustype, angle, magnitude, voltage_limits, base_voltage, area, load_zone, ext, internal, )
end

# Constructor for demo purposes; non-functional.
function Bus(::Nothing)
    Bus(;
        number=0,
        name="init",
        bustype=nothing,
        angle=0.0,
        magnitude=0.0,
        voltage_limits=(min=0.0, max=0.0),
        base_voltage=nothing,
        area=nothing,
        load_zone=nothing,
        ext=Dict{String, Any}(),
    )
end

"""Get [`Bus`](@ref) `number`."""
get_number(value::Bus) = value.number

InfrastructureSystems.get_name(value::Bus) = value.name
"""Get [`Bus`](@ref) `bustype`."""
get_bustype(value::Bus) = value.bustype
"""Get [`Bus`](@ref) `angle`."""
get_angle(value::Bus) = value.angle
"""Get [`Bus`](@ref) `magnitude`."""
get_magnitude(value::Bus) = value.magnitude
"""Get [`Bus`](@ref) `voltage_limits`."""
get_voltage_limits(value::Bus) = value.voltage_limits
"""Get [`Bus`](@ref) `base_voltage`."""
get_base_voltage(value::Bus) = value.base_voltage
"""Get [`Bus`](@ref) `area`."""
get_area(value::Bus) = value.area
"""Get [`Bus`](@ref) `load_zone`."""
get_load_zone(value::Bus) = value.load_zone
"""Get [`Bus`](@ref) `ext`."""
get_ext(value::Bus) = value.ext
"""Get [`Bus`](@ref) `internal`."""
get_internal(value::Bus) = value.internal

"""Set [`Bus`](@ref) `number`."""
set_number!(value::Bus, val) = value.number = val

InfrastructureSystems.set_name!(value::Bus, val) = value.name = val
"""Set [`Bus`](@ref) `bustype`."""
set_bustype!(value::Bus, val) = value.bustype = val
"""Set [`Bus`](@ref) `angle`."""
set_angle!(value::Bus, val) = value.angle = val
"""Set [`Bus`](@ref) `magnitude`."""
set_magnitude!(value::Bus, val) = value.magnitude = val
"""Set [`Bus`](@ref) `voltage_limits`."""
set_voltage_limits!(value::Bus, val) = value.voltage_limits = val
"""Set [`Bus`](@ref) `base_voltage`."""
set_base_voltage!(value::Bus, val) = value.base_voltage = val
"""Set [`Bus`](@ref) `area`."""
set_area!(value::Bus, val) = value.area = val
"""Set [`Bus`](@ref) `load_zone`."""
set_load_zone!(value::Bus, val) = value.load_zone = val
"""Set [`Bus`](@ref) `ext`."""
set_ext!(value::Bus, val) = value.ext = val


text:
using Printf
using LinearAlgebra

include("../common/PWGrid_v01.jl")
include("../common/wrappers_fft.jl")
include("../common/gen_lattice.jl")

include("gen_dr.jl")
include("gen_rho.jl")
include("Poisson_solve.jl")

function test_main()
    #
    Ns = [64, 64, 64]
    LatVecs = gen_lattice_sc(16.0)
    #
    pw = PWGrid( Ns, LatVecs )
    #
    Npoints = pw.Npoints
    Ω = pw.Ω
    r = pw.r
    Ns = pw.Ns
    #
    # Generate array of distances
    #
    center = sum(LatVecs,dims=2)/2
    println("center = ", center)
    dr = gen_dr( r, center )
    #
    # Generate charge density
    #
    σ1 = 0.75
    σ2 = 0.50
    rho = gen_rho( dr, σ1, σ2 )
    #
    # Solve Poisson equation and calculate Hartree energy
    #
    phi = Poisson_solve( pw, rho )
    Ehartree = 0.5*dot( phi, rho ) * Ω/Npoints
    #
    Uanal = ( (1/σ1 + 1/σ2)/2 - sqrt(2) / sqrt( σ1^2 + σ2^2 ) ) / sqrt(pi)
    @printf("Num, ana, diff = %18.10f %18.10f %18.10e\n", Ehartree, Uanal, abs(Ehartree-Uanal))
end

@time test_main()


text:
# functions for plotting results

"""
    plot_temperature(result::MolecularDynamicsResult, stride::Integer)::Plot

Plot the temperature of a `MolecularDynamicsResult` against time, sampling every `stride` points.
"""
function plot_temperature(result::MolecularDynamicsResult, stride::Integer)
    N = length(get_system(result))
    p = plot(
        title = "Temperature during Simulation [n = $(N)]",
        xlab = "Time",
        ylab = "Temperature",
    )
    plot_temperature!(p, result, stride, true)
end
"""
    plot_temperature!(p::Plot, result::MolecularDynamicsResult, stride::Integer, first_plot::Bool = false)::Plot

Plot the temperature of a `MolecularDynamicsResult` against time, sampling every `stride` points.
Add the new line to an existing plot and update the legend only if it is the first plot.
If it is not the first plot, add a vertical line to differentiate the segments of the simulation in the plot.
"""
function plot_temperature!(p::Plot, result::MolecularDynamicsResult, stride::Integer, first_plot::Bool = false)
    time_range = [(t, auconvert(u"ps", time)) for (t, time) ∈ enumerate(get_time_range(result)) if (t - 1) % stride == 0]
    if (!first_plot)
        vline!(
            p,
            [time_range[1][2]],
            label = false,
            color = :black,
            linestyle = :dot,
            lw = 2
        )
    end
    plot!(
        p,
        [time for (t, time) ∈ time_range],
        [auconvert(u"K", temperature(result, austrip(t))) for (t, time) ∈ time_range],
        label = first_plot ? "Simulation Temperature" : nothing,
        color = 1,
    )
    reference_temp = reference_temperature(result)
    if (!ismissing(reference_temp))
        plot!(
            p,
            [time for (t, time) ∈ time_range],
            [auconvert(u"K", reference_temp) for (t, time) ∈ time_range],
            label = first_plot ? "Reference Temperature" : nothing,
            color = 2,
            linestyle = :dash,
            lw = 2
        )
    end
    p
end

"""
    plot_energy(result::MolecularDynamicsResult, stride::Integer)::Plot

Plot the kinetic, potential, and total energy of a `MolecularDynamicsResult` against time, sampling every `stride` points.
"""
function plot_energy(result::MolecularDynamicsResult, stride::Integer)
    N = length(get_system(result))
    p = plot(
        title = "Energy during Simulation [n = $(N)]",
        xlab = "Time",
        ylab = "Energy",
        legend = :right
    )
    plot_energy!(p, result, stride, true)
end
"""
    plot_energy!(p::Plot, result::MolecularDynamicsResult, stride::Integer, first_plot::Bool = false)::Plot

Plot the kinetic, potential, and total energy of a `MolecularDynamicsResult` against time, sampling every `stride` points.
Add the new lines to an existing plot and update the legend only if it is the first plot.
If it is not the first plot, add a vertical line to differentiate the segments of the simulation in the plot.
"""
function plot_energy!(p::Plot, result::MolecularDynamicsResult, stride::Integer, first_plot::Bool = false)
    time_range = [(t, auconvert(u"ps", time)) for (t, time) ∈ enumerate(get_time_range(result)) if (t - 1) % stride == 0]
    if (!first_plot)
        vline!(
            p,
            [time_range[1][2]],
            label = false,
            color = :black,
            linestyle = :dot,
            lw = 2
        )
    end
    plot!(
        p,
        [time for (t, time) ∈ time_range],
        [kinetic_energy(result, austrip(t))u"hartree" for (t, time) ∈ time_range],
        label = first_plot ? "Kinetic Energy" : nothing,
        color = 2
    )
    plot!(
        p,
        [time for (t, time) ∈ time_range],
        [potential_energy(result, austrip(t))u"hartree" for (t, time) ∈ time_range],
        label = first_plot ? "Potential Energy" : nothing,
        color = 1
    )
    plot!(
        p,
        [time for (t, time) ∈ time_range],
        [total_energy(result, austrip(t))u"hartree" for (t, time) ∈ time_range],
        label = first_plot ? "Total Energy" : nothing,
        color = 3
    )
end

"""
    plot_rdf(result::MolecularDynamicsResult, σ::Real, sample_fraction::Float64 = 1.0)::Plot

Plot the radial distribution function of a `MolecularDynamicsResult` sampling only the trailing `sample_fraction` of the timesteps.
Use σ (from Lennard Jones) as a normalization factor for the radius (assumed to be in atomic units).
"""
function plot_rdf(result::MolecularDynamicsResult, σ::Real, sample_fraction::Float64 = 1.0)
    plot_rdf(result, σ * u"bohr", sample_fraction)
end
"""
    plot_rdf(result::MolecularDynamicsResult, σ::Unitful.Length, sample_fraction::Float64 = 1.0)::Plot

Plot the radial distribution function of a `MolecularDynamicsResult` sampling only the trailing `sample_fraction` of the timesteps.
Use σ (from Lennard Jones) as a normalization factor for the radius.
"""
function plot_rdf(result::MolecularDynamicsResult, σ::Unitful.Length, sample_fraction::Float64 = 1.0)
    @assert 0 < sample_fraction ≤ 1
    N = length(get_system(result))
    T = length(get_time_range(result)) - 1
    rs, grf = rdf(result, sample_fraction)
    @assert length(rs) == length(grf)
    plot(
        title = "Radial Distribution Function [n = $(N)] [T = $(T)]",
        xlab = "Distance r/σ",
        ylab = "Radial Distribution g(r)",
        legend = false
    )
    plot!(
        rs * u"bohr" / auconvert(σ),
        grf
    )
end

max_stars_repo_name: jrdegreeff/cesmix-julia

max_issues_repo_name: jrdegreeff/cesmix-julia

max_forks_repo_name: jrdegreeff/cesmix-julia


text:
using Revise
using Pkg
Pkg.activate("C:/Users/cmmenzen/.julia/dev/BigMat")
using BigMat
using Test
using Random
using LinearAlgebra
using Optim

@testset "BigMat.jl" begin
    include("MPT.jl")
    include("ComputeTN.jl")
    include("MatrixAlgebra.jl") # inversion
end

max_stars_repo_name: clarazen/BigMat

max_issues_repo_name: clarazen/BigMat

max_forks_repo_name: clarazen/BigMat


text:
# routines to determine correct matching of landmarks (cyclic shift detection)

"""
    update_cyclicmatching(X, ll,obs_info, Xᵒ, W, Q)

## writes into / modifies
- `X`
- `ll`

## Returns
`Q`, `obsinfo`, `accept`
"""
function update_cyclicmatching(X, ll,obsinfo, Xᵒ, W, Q, At)
    # shift each element in xobsT, this means
        # - creating xobsTᵒ (shift each element of xobsT)
        # - creating obsinfoᵒ
        # - creating  Qᵒ from Q by using xobsTᵒ and obsinfoᵒ
        # - recomputing backward recursion

    direction = 2 * rand(Bool) -1
    #    xobsTᵒ = Base.circshift.(Q.xobsT,direction) # this would be for shifting all shapes
    xobsTᵒ = copy(Q.xobsT)
    k = sample(1:Q.nshapes) # randomly pick a shape to which we cyclically shift indices
    xobsTᵒ[k] = Base.circshift(Q.xobsT[k],direction)
    oi = obsinfo
    obsinfoᵒ = ObsInfo(oi.L0, oi.LT, oi.Σ0, oi.ΣT, oi.xobs0, xobsTᵒ, oi.obs_atzero, oi.fixinitmomentato0, oi.n, oi.nshapes)
    Qᵒ = construct_gp_xobsT(Q, xobsTᵒ)
    Qᵒ = update_guidrec!(Qᵒ, obsinfoᵒ)

    x0 = X[1].yy[1]
    llᵒ, Xᵒ = gp!(LeftRule(), Xᵒ, x0, W, Qᵒ; skip=sk)
    if log(rand()) <= (sum(llᵒ) - sum(ll))
        ll .= llᵒ
        Q = Qᵒ
        obsinfo = obsinfoᵒ
                # for k ∈ 1:Q.nshapes
                #     for i ∈ eachindex(X[1].yy)
                #         X[k].yy[i] .= Xᵒ[k].yy[i]
                #     end
                # end
        X = copypaths!(X,Xᵒ)
        accept = 1
        println("Cyclic shift for shape $k in direction $direction.")
    else
        accept = 0
    end
    obsinfo, accept, Q, X, ll
end


text:
using NumericalIntegration
using Test
using InteractiveUtils # for subtypes
using StaticArrays
using HCubature # for testing n-dimensional integration

@testset "compare with analytic result" begin
    x = collect(-π : 2*π/2048 : π)
    y = sin.(x)
    p = collect(0 : π/1024 : π)
    q = sin.(p)
    for M in subtypes(IntegrationMethod)
        for T in [Float32, Float64, BigFloat]
            for (xs,ys,val,atol) in [
                                     (x,y,0,1e-4),
                                     (p,q,2,1e-4),
                                  ]
                result = @inferred integrate(T.(xs), T.(ys),M())
                @test isapprox(result, val, atol=atol)
                @test typeof(result) == T
            end
        end
    end
end

@testset "cumulative integration" begin
    x = collect(-π : 2*π/2048 : π)
    y = sin.(x)
    exact = @. -cos(x) + cos(-π)
    for M in subtypes(IntegrationMethod)
        hasmethod(cumul_integrate, Tuple{AbstractVector, AbstractVector, M}) || continue
        for T in [Float32, Float64, BigFloat]
            result = @inferred cumul_integrate(T.(x), T.(y), M())
            @test typeof(result) == Vector{T}
            @test all(isapprox.(result, exact, atol=1e-4))
        end
    end
end

@testset "n-dimensional integration testing" begin
    X = range(0,stop=2π,length=10)
    Y = range(-π,stop=π,length=10)
    Z = range(0,stop=2,length=10)

    A = Array{Float64}(undef,length(X),length(Y),length(Z))
    f(x) = sin(x[1]) + cos(x[2]) + 2x[3]

    for (k,z) in enumerate(Z)
        for (j,y) in enumerate(Y)
            for (i,x) in enumerate(X)
                A[i,j,k] = f([x,y,z])
            end
        end
    end

    @test isapprox(integrate((X,Y,Z),A), hcubature(f,[0,-π,0],[2π,π,2])[1], atol=1e-4)
end

@testset "SVector" begin
    xs = range(0, stop=1, length=9)
    ys1 = randn(9)
    ys2 = randn(9)
    ys = map(SVector, ys1, ys2)
    for M in subtypes(IntegrationMethod)
        m = M()
        res1 = integrate(xs,ys1,m)
        res2 = integrate(xs,ys2,m)
        res = @inferred integrate(xs,ys,m)
        @test res ≈ @SVector [res1, res2]
    end
end

@testset "Raising Warnings" begin
    xs = collect(-1.0 : 0.5 : 1.0)
    ys = xs.^2
    m = RombergEven()
    expwarn = "RombergEven :: final step reached, but accuracy not: 1.3333333333333335 > 1.0e-12"
    @test_logs (:warn, expwarn) integrate(xs, ys, m)
end

using Unitful

@testset "Unitful compatibility" begin
    x = LinRange(1u"s", 10u"s", 9)
    y = rand(9)*u"m/s"
    for M in subtypes(IntegrationMethod)
        @test typeof(integrate(x,y,M())) == typeof(1.0u"m")
        if hasmethod(cumul_integrate, Tuple{AbstractVector, AbstractVector, M})
            @test typeof(cumul_integrate(x,y,M())) == Vector{typeof(1.0u"m")}
        end
    end
end


text:
using ManifoldsBase

using Test
import Base: *
struct NonManifold <: AbstractManifold{ManifoldsBase.ℝ} end
struct NonMPoint <: AbstractManifoldPoint end
struct NonTVector <: TVector end
struct NonCoTVector <: CoTVector end
struct NotImplementedRetraction <: AbstractRetractionMethod end
struct NotImplementedInverseRetraction <: AbstractInverseRetractionMethod end
*(t::Float64, X::NonTVector) = X
@testset "AbstractManifold with empty implementation" begin
    M = NonManifold()
    p = NonMPoint()
    v = NonTVector()
    @test base_manifold(M) === M
    @test number_system(M) === ℝ
    @test representation_size(M) === nothing

    @test_throws MethodError manifold_dimension(M)

    # by default isapprox compares given points or vectors
    @test isapprox(M, [0], [0])
    @test isapprox(M, [0], [0]; atol = 1e-6)
    @test !isapprox(M, [0], [1])
    @test !isapprox(M, [0], [1]; atol = 1e-6)
    @test isapprox(M, [0], [0], [0])
    @test isapprox(M, [0], [0], [0]; atol = 1e-6)
    @test !isapprox(M, [0], [0], [1])
    @test !isapprox(M, [0], [0], [1]; atol = 1e-6)

    exp_retr = ManifoldsBase.ExponentialRetraction()

    @test_throws MethodError retract!(M, p, p, v)
    @test_throws MethodError retract!(M, p, p, v, exp_retr)
    @test_throws MethodError retract!(M, p, p, [0.0], 0.0)
    @test_throws MethodError retract!(M, p, p, [0.0], 0.0, exp_retr)
    @test_throws MethodError retract!(M, [0], [0], [0])
    @test_throws MethodError retract!(M, [0], [0], [0], exp_retr)
    @test_throws MethodError retract!(M, [0], [0], [0], 0.0)
    @test_throws MethodError retract!(M, [0], [0], [0], 0.0, exp_retr)
    @test_throws MethodError retract(M, [0], [0])
    @test_throws MethodError retract(M, [0], [0], exp_retr)
    @test_throws MethodError retract(M, [0], [0], 0.0)
    @test_throws MethodError retract(M, [0], [0], 0.0, exp_retr)
    @test_throws MethodError retract(M, [0.0], [0.0])
    @test_throws MethodError retract(M, [0.0], [0.0], exp_retr)
    @test_throws MethodError retract(M, [0.0], [0.0], 0.0)
    @test_throws MethodError retract(M, [0.0], [0.0], 0.0, exp_retr)
    @test_throws MethodError retract(M, [0.0], [0.0], NotImplementedRetraction())

    log_invretr = ManifoldsBase.LogarithmicInverseRetraction()

    @test_throws MethodError inverse_retract!(M, p, p, p)
    @test_throws MethodError inverse_retract!(M, p, p, p, log_invretr)
    @test_throws MethodError inverse_retract!(M, [0], [0], [0])
    @test_throws MethodError inverse_retract!(M, [0], [0], [0], log_invretr)
    @test_throws MethodError inverse_retract(M, [0], [0])
    @test_throws MethodError inverse_retract(M, [0], [0], log_invretr)
    @test_throws MethodError inverse_retract(M, [0.0], [0.0])
    @test_throws MethodError inverse_retract(M, [0.0], [0.0], log_invretr)
    @test_throws MethodError inverse_retract(
        M,
        [0.0],
        [0.0],
        NotImplementedInverseRetraction(),
    )

    @test_throws MethodError project!(M, p, [0])
    @test_throws MethodError project!(M, [0], [0])
    @test_throws MethodError project(M, [0])

    @test_throws MethodError project!(M, v, p, [0.0])
    @test_throws MethodError project!(M, [0], [0], [0])
    @test_throws MethodError project(M, [0], [0])
    @test_throws MethodError project(M, [0.0], [0.0])

    @test_throws MethodError inner(M, p, v, v)
    @test_throws MethodError inner(M, [0], [0], [0])
    @test_throws MethodError norm(M, p, v)
    @test_throws MethodError norm(M, [0], [0])
    @test_throws MethodError angle(M, p, v, v)
    @test_throws MethodError angle(M, [0], [0], [0])

    @test_throws MethodError distance(M, [0.0], [0.0])

    @test_throws MethodError exp!(M, p, p, v)
    @test_throws MethodError exp!(M, p, p, v, 0.0)
    @test_throws MethodError exp!(M, [0], [0], [0])
    @test_throws MethodError exp!(M, [0], [0], [0], 0.0)
    @test_throws MethodError exp(M, [0], [0])
    @test_throws MethodError exp(M, [0], [0], 0.0)
    @test_throws MethodError exp(M, [0.0], [0.0])
    @test_throws MethodError exp(M, [0.0], [0.0], 0.0)

    @test_throws MethodError embed!(M, p, [0]) # no copy for NoPoint p
    @test embed!(M, [0], [0]) == [0]
    @test embed(M, [0]) == [0]

    # Identity
    @test_throws MethodError embed!(M, v, p, [0.0]) # no copyto
    @test embed!(M, [0], [0], [0]) == [0]
    @test_throws MethodError embed(M, [0], v) # no copyto
    @test embed(M, [0.0], [0.0]) == [0.0]

    @test_throws MethodError log!(M, v, p, p)
    @test_throws MethodError log!(M, [0], [0], [0])
    @test_throws MethodError log(M, [0.0], [0.0])

    @test_throws MethodError vector_transport_to!(M, [0], [0], [0], [0])
    @test_throws MethodError vector_transport_to(M, [0], [0], [0])
    @test_throws MethodError vector_transport_to!(M, [0], [0], [0], ProjectionTransport())

    @test_throws MethodError vector_transport_direction!(M, [0], [0], [0], [0])
    @test_throws MethodError vector_transport_direction(M, [0], [0], [0])

    @test_throws MethodError ManifoldsBase.vector_transport_along!(M, [0], [0], [0], x -> x)
    @test_throws MethodError vector_transport_along(M, [0], [0], x -> x)

    @test_throws MethodError injectivity_radius(M)
    @test_throws MethodError injectivity_radius(M, [0])
    @test_throws MethodError injectivity_radius(M, [0], exp_retr)
    @test_throws MethodError injectivity_radius(M, exp_retr)

    @test_throws MethodError zero_vector!(M, [0], [0])
    @test_throws MethodError zero_vector(M, [0])

    @test ManifoldsBase.check_point(M, [0]) === nothing
    @test ManifoldsBase.check_point(M, p) === nothing
    @test is_point(M, [0])
    @test ManifoldsBase.check_point(M, [0]) === nothing

    @test ManifoldsBase.check_vector(M, [0], [0]) === nothing
    @test ManifoldsBase.check_vector(M, p, v) === nothing
    @test is_vector(M, [0], [0])
    @test ManifoldsBase.check_vector(M, [0], [0]) === nothing

    @test_throws MethodError hat!(M, [0], [0], [0])
    @test_throws MethodError vee!(M, [0], [0], [0])
end


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

JLLWrappers.@generate_main_file_header("gb")
JLLWrappers.@generate_main_file("gb", UUID("006bdf2e-00d8-555c-bee2-993406487cbf"))
end  # module gb_jll


text:
# support for device-side exceptions (from CUDAnative/src/exceptions.jl)

## exception type

struct KernelException <: Exception
    dev::RuntimeDevice
    exstr::Union{String,Nothing}
end
KernelException(dev) = KernelException(dev, nothing)

function Base.showerror(io::IO, err::KernelException)
    print(io, "KernelException: exception(s) thrown during kernel execution on device $(err.dev.device)")
    if err.exstr !== nothing
        println(io, ":")
        print(io, err.exstr)
    end
end

## exception ring buffer

struct ExceptionEntry
    kern::UInt64
    ptr::DevicePtr{UInt8,AS.Global}
end
ExceptionEntry(kern) = ExceptionEntry(kern, DevicePtr{UInt8,AS.Global}(0))
ExceptionEntry() = ExceptionEntry(0)

## exception codegen

# emit a global variable for storing the current exception status
function emit_exception_user!(mod::LLVM.Module)
    # add a fake user for __ockl_hsa_signal_store and __ockl_hsa_signal_load
    if !haskey(LLVM.functions(mod), "__fake_global_exception_flag_user")
        JuliaContext() do ctx
            ft = LLVM.FunctionType(LLVM.VoidType(ctx))
            fn = LLVM.Function(mod, "__fake_global_exception_flag_user", ft)
            Builder(ctx) do builder
                entry = BasicBlock(fn, "entry", ctx)
                position!(builder, entry)
                T_nothing = LLVM.VoidType(ctx)
                T_i32 = LLVM.Int32Type(ctx)
                T_i64 = LLVM.Int64Type(ctx)
                T_signal_store = LLVM.FunctionType(T_nothing, [T_i64, T_i64, T_i32])
                signal_store = LLVM.Function(mod, "__ockl_hsa_signal_store", T_signal_store)
                call!(builder, signal_store, [ConstantInt(0,ctx),
                                              ConstantInt(0,ctx),
                                              # __ATOMIC_RELEASE == 3
                                              ConstantInt(Int32(3), ctx)])
                T_signal_load = LLVM.FunctionType(T_i64, [T_i64, T_i32])
                signal_load = LLVM.Function(mod, "__ockl_hsa_signal_load", T_signal_load)
                loaded_value = call!(builder, signal_load, [ConstantInt(0,ctx),
                                                            # __ATOMIC_ACQUIRE == 2
                                                            ConstantInt(Int32(2), ctx)])
                ret!(builder)
            end
        end
    end
    @assert haskey(LLVM.functions(mod), "__fake_global_exception_flag_user")
end
function delete_exception_user!(mod::LLVM.Module)
    if haskey(LLVM.functions(mod), "__fake_global_exception_flag_user")
        delete!(mod, LLVM.functions(mod)["__fake_global_exception_flag_user"])
    end
    @assert !haskey(LLVM.functions(mod), "__fake_global_exception_flag_user")
end


text:
# Robust Statistics

#############################
#
#   Trimming outliers
#
#############################

# Trimmed set
"""
    trim(x; prop=0.0, count=0)

Return a copy of `x` with either `count` or proportion `prop` of the highest
and lowest elements removed.  To compute the trimmed mean of `x` use
`mean(trim(x))`; to compute the variance use `trimvar(x)` (see [`trimvar`](@ref)).

# Example
```julia
julia> trim([1,2,3,4,5], prop=0.2)
3-element Array{Int64,1}:
 2
 3
 4
```
"""
function trim(x::AbstractVector; prop::Real=0.0, count::Integer=0)
    trim!(copy(x); prop=prop, count=count)
end

"""
    trim!(x; prop=0.0, count=0)

A variant of [`trim`](@ref) that modifies `x` in place.
"""
function trim!(x::AbstractVector; prop::Real=0.0, count::Integer=0)
    n = length(x)
    n > 0 || throw(ArgumentError("x can not be empty."))

    if count == 0
        0 <= prop < 0.5 || throw(ArgumentError("prop must satisfy 0 ≤ prop < 0.5."))
        count = floor(Int, n * prop)
    else
        prop == 0 || throw(ArgumentError("prop and count can not both be > 0."))
        0 <= count < n/2 || throw(ArgumentError("count must satisfy 0 ≤ count < length(x)/2."))
    end

    partialsort!(x, (n-count+1):n)
    partialsort!(x, 1:count)
    deleteat!(x, (n-count+1):n)
    deleteat!(x, 1:count)

    return x
end

"""
    winsor(x; prop=0.0, count=0)

Return a copy of `x` with either `count` or proportion `prop` of the lowest
elements of `x` replaced with the next-lowest, and an equal number of the
highest elements replaced with the previous-highest.  To compute the Winsorized
mean of `x` use `mean(winsor(x))`.

# Example
```julia
julia> winsor([1,2,3,4,5], prop=0.2)
5-element Array{Int64,1}:
 2
 2
 3
 4
 4
```
"""
function winsor(x::AbstractVector; prop::Real=0.0, count::Integer=0)
    winsor!(copy(x); prop=prop, count=count)
end

"""
    winsor!(x; prop=0.0, count=0)

A variant of [`winsor`](@ref) that modifies vector `x` in place.
"""
function winsor!(x::AbstractVector; prop::Real=0.0, count::Integer=0)
    n = length(x)
    n > 0 || throw(ArgumentError("x can not be empty."))

    if count == 0
        0 <= prop < 0.5 || throw(ArgumentError("prop must satisfy 0 ≤ prop < 0.5."))
        count = floor(Int, n * prop)
    else
        prop == 0 || throw(ArgumentError("prop and count can not both be > 0."))
        0 <= count < n/2 || throw(ArgumentError("count must satisfy 0 ≤ count < length(x)/2."))
    end

    partialsort!(x, (n-count+1):n)
    partialsort!(x, 1:count)
    x[1:count] .= x[count+1]
    x[n-count+1:end] .= x[n-count]

    return x
end

#############################
#
#   Other
#
#############################

# Variance of a trimmed set.
"""
    trimvar(x; prop=0.0, count=0)

Compute the variance of the trimmed mean of `x`. This function uses
the Winsorized variance, as described in Wilcox (2010).
"""
function trimvar(x::AbstractVector; prop::Real=0.0, count::Integer=0)
    n = length(x)
    n > 0 || throw(ArgumentError("x can not be empty."))

    if count == 0
        0 <= prop < 0.5 || throw(ArgumentError("prop must satisfy 0 ≤ prop < 0.5."))
        count = floor(Int, n * prop)
    else
        0 <= count < n/2 || throw(ArgumentError("count must satisfy 0 ≤ count < length(x)/2."))
        prop = count/n
    end

    return var(winsor(x, count=count)) / (n * (1 - 2prop)^2)
end


text:
using Pkg.Registry
Registry.add(Registry.RegistrySpec(url = "https://github.com/JuliaRegistries/General.git"))
Registry.add(Registry.RegistrySpec(url = "https://github.com/BioJulia/BioJuliaRegistry.git"))


text:
using BCVUMPS
using Plots
using JLD2
using FileIO

# 1x1 Ising magnetisation
begin
    magplot = plot()
    model = Ising()
    for D = [2,4,8,16]
        mag = []
        for β = 0.41:0.0025:0.48
            printstyled("D = $(D), β = $(β) \n"; bold=true, color=:red)

            env = bcvumps_env(model, β,D)
            mag = [mag;magnetisation(env, model, β)]
            chkp_file = "./data/$(model)_β$(round(β+0.0025,digits=4))_D$(D).jld2"
            if isfile(chkp_file) == false
                save(chkp_file, "env", env)
            end
        end
        β = 0.41:0.0025:0.48
        magplot = plot!(β,mag,seriestype = :scatter,title = "magnetisation", label = "BCVUMPS D = $(D)", lw = 3)
    end

    tmag = []
    for β = 0.41:0.001:0.48
        tmag = [tmag;magofβ(Ising(), β)]
    end
    β = 0.41:0.001:0.48
    magplot = plot!(β,tmag,label = "exact", lw = 2)
    xlabel!("β")
    ylabel!("mag")
    savefig(magplot,"./plot/2Disingmag.svg")
end


text:
import HTTP.Parsers
import HTTP.URIs

@testset "utils.jl" begin
    @test HTTP.Strings.escapehtml("&\"'<>") == "&"'<>"

    @test HTTP.Cookies.isurlchar('\u81')
    @test !HTTP.Cookies.isurlchar('\0')

    @test HTTP.Strings.tocameldash("accept") == "Accept"
    @test HTTP.Strings.tocameldash("Accept") == "Accept"
    @test HTTP.Strings.tocameldash("eXcept-this") == "Except-This"
    @test HTTP.Strings.tocameldash("exCept-This") == "Except-This"
    @test HTTP.Strings.tocameldash("not-valid") == "Not-Valid"
    @test HTTP.Strings.tocameldash("♇") == "♇"
    @test HTTP.Strings.tocameldash("bλ-a") == "Bλ-A"
    @test HTTP.Strings.tocameldash("not fixable") == "Not fixable"
    @test HTTP.Strings.tocameldash("aaaaaaaaaaaaa") == "Aaaaaaaaaaaaa"
    @test HTTP.Strings.tocameldash("conTENT-Length") == "Content-Length"
    @test HTTP.Strings.tocameldash("Sec-WebSocket-Key2") == "Sec-Websocket-Key2"
    @test HTTP.Strings.tocameldash("User-agent") == "User-Agent"
    @test HTTP.Strings.tocameldash("Proxy-authorization") == "Proxy-Authorization"
    @test HTTP.Strings.tocameldash("HOST") == "Host"
    @test HTTP.Strings.tocameldash("ST") == "St"
    @test HTTP.Strings.tocameldash("X-\$PrototypeBI-Version") == "X-\$prototypebi-Version"
    @test HTTP.Strings.tocameldash("DCLK_imp") == "Dclk_imp"

    for (bytes, utf8) in (
            (UInt8[], ""),
            (UInt8[0x00], "\0"),
            (UInt8[0x61], "a"),
            (UInt8[0x63, 0x61, 0x66, 0xe9, 0x20, 0x63, 0x72, 0xe8, 0x6d, 0x65], "café crème"),
            (UInt8[0x6e, 0x6f, 0xeb, 0x6c], "noël"),
            (UInt8[0xc4, 0xc6, 0xe4], "ÄÆä"),
        )
        @test HTTP.Strings.iso8859_1_to_utf8(bytes) == utf8
    end

    withenv("HTTPS_PROXY"=>nothing, "https_proxy"=>nothing) do
        @test HTTP.ConnectionRequest.getproxy("https", "https://julialang.org/") == nothing
    end
    withenv("HTTPS_PROXY"=>"") do
        # to be compatible with Julia 1.0
        @test HTTP.ConnectionRequest.getproxy("https", "https://julialang.org/") == nothing
    end
    withenv("https_proxy"=>"") do
        @test HTTP.ConnectionRequest.getproxy("https", "https://julialang.org/") == nothing
    end
    withenv("HTTPS_PROXY"=>"https://user:pass@server:80") do
        @test HTTP.ConnectionRequest.getproxy("https", "https://julialang.org/") == "https://user:pass@server:80"
    end
    withenv("https_proxy"=>"https://user:pass@server:80") do
        @test HTTP.ConnectionRequest.getproxy("https", "https://julialang.org/") == "https://user:pass@server:80"
    end

    withenv("HTTP_PROXY"=>nothing, "http_proxy"=>nothing) do
        @test HTTP.ConnectionRequest.getproxy("http", "http://julialang.org/") == nothing
    end
    withenv("HTTP_PROXY"=>"") do
        @test HTTP.ConnectionRequest.getproxy("http", "http://julialang.org/") == nothing
    end
    withenv("http_proxy"=>"") do
        @test HTTP.ConnectionRequest.getproxy("http", "http://julialang.org/") == nothing
    end
    withenv("HTTP_PROXY"=>"http://user:pass@server:80") do
        @test HTTP.ConnectionRequest.getproxy("http", "http://julialang.org/") == "http://user:pass@server:80"
    end
    withenv("http_proxy"=>"http://user:pass@server:80") do
        @test HTTP.ConnectionRequest.getproxy("http", "http://julialang.org/") == "http://user:pass@server:80"
    end
end # testset


text:
## Exercise 6-6
## A palindrome is a word that is spelled the same backward and forward, like “noon” and “redivider”. Recursively, a word is a palindrome if the first and last letters are the same and the middle is a palindrome.

## The following are functions that take a string argument and return the first, last, and middle letters:

function first(word)
    first = firstindex(word)
    word[first]
end

function last(word)
    last = lastindex(word)
    word[last]
end

function middle(word)
    first = firstindex(word)
    last = lastindex(word)
    word[nextind(word, first) : prevind(word, last)]
end

## We’ll see how they work in Strings

## 1. Test these functions out. What happens if you call middle with a string with two letters? One letter? What about the empty string, which is written "" and contains no letters?
println("Ans 1: ")

println(first("first")) # returns 'f'
println(last("last")) # returns 't'
println(middle("middle")) # returns "iddl"

println(middle("as")) # returns ''
println(middle("A")) # returns ''
# println(middle("")) # error, empty string

## 2. Write a function called ispalindrome that takes a string argument and returns true if it is a palindrome and false otherwise. Remember that you can use the built-in function length to check the length of a string.
println("Ans 2: ")

function ispalindrome(word)
    if length(word) == 0
        return true
    end

    if first(word) == last(word)
        return ispalindrome(middle(word))
    else
        return false
    end
end

println(ispalindrome("noon"))
println(ispalindrome("redivider"))
println(ispalindrome("name"))

println("End.")


text:
mutable struct LKJPrecisionXMLElement <: MyXMLElement
    el::XMLOrNothing
    var::AbstractMatrix{Float64}
    diag_param::Parameter
    off_diag_parm::Parameter
    id::String

    function LKJPrecisionXMLElement(var_mat::AbstractArray{Float64},
                                    id::AbstractString)
        check_posdef(var_mat)
        vars, corrs = var_corrs(var_mat)
        diag_param = Parameter(vars, "$id.diagonalVar")
        diag_param.lower = 0.0
        off_diag_param = Parameter(corrs, "$id.correlation")
        off_diag_param.lower = -1.0
        off_diag_param.upper = 1.0
        return new(nothing, var_mat, diag_param, off_diag_param, id)
    end
end

function LKJPrecisionXMLElement(p::Int)
    return LKJPrecisionXMLElement(Diagonal(ones(p)), bn.DIFF_PREC_ID)
end

function make_xml(lkj::LKJPrecisionXMLElement)
    var, corrs = var_corrs(lkj.var)

    set_value(lkj.diag_param, var)
    set_value(lkj.off_diag_parm, corrs)

    el = new_element(bn.CACHED_MATRIX_INVERSE)
    set_id!(el, lkj.id)

    cs_el = new_child(el, bn.COMPOUND_SYMMETRIC_MATRIX)
    set_attributes(cs_el,
                   Dict(bn.AS_CORRELATION => bn.TRUE, bn.IS_CHOLESKY => bn.TRUE)
                  )

    diag_el = new_child(cs_el, bn.DIAGONAL)
    dp_el = make_xml(lkj.diag_param)
    add_child(diag_el, dp_el)
    od_el = new_child(cs_el, bn.OFF_DIAGONAL)
    op_el = make_xml(lkj.off_diag_parm)
    add_child(od_el, op_el)

    lkj.el = el
    return el
end

function var_corrs(var::AbstractArray{Float64})
    vars = diag(var)
    corr_mat = cov2cor(var, sqrt.(diag(var)))
    p = size(var, 1)
    n_corrs = div(p * (p - 1), 2)
    corrs = zeros(n_corrs)
    ind = 0
    for i = 1:p
        for j = (i + 1):p
            ind += 1
            corrs[ind] = corr_mat[i, j]
        end
    end
    return vars, corrs
end

################################################################################
## Priors
################################################################################

mutable struct LKJCorrelationPrior <: MyXMLElement
    el::XMLOrNothing
    corr_param::Parameter
    shape::Float64

    function LKJCorrelationPrior(corr_param::Parameter)
        return new(nothing, corr_param, 1.0)
    end
end

function make_xml(lkj::LKJCorrelationPrior)
    make_xml(lkj.corr_param)
    el = new_element(bn.LJK_CORRELATION_PRIOR)
    set_id!(el, get_id(lkj.corr_param.el) * ".prior")
    set_attribute(el, bn.SHAPE_PARAMETER, lkj.shape)
    dim = div(1 + isqrt(1 + 8 * length(lkj.corr_param)), 2)
    set_attribute(el, bn.DIMENSION, dim)
    set_attribute(el, bn.CHOLESKY, bn.TRUE)
    data_el = new_child(el, bn.DATA)
    add_ref_el(data_el, lkj.corr_param.el)
    lkj.el = el
    return el
end

mutable struct LogNormalPriorXMLElement <: MyXMLElement
    el::XMLOrNothing
    param::Parameter
    μ::Float64
    σ::Float64
    offset::Float64

    function LogNormalPriorXMLElement(param::Parameter)
        return new(nothing, param, 0.0, 1.0, 0.0)
    end
end

function make_xml(lnp::LogNormalPriorXMLElement)
    make_xml(lnp.param)
    el = new_element(bn.LOG_NORMAL_PRIOR)
    set_id!(el, get_id(lnp.param.el) * ".prior")
    set_attributes(el, Dict(
                            bn.MEAN => string(lnp.μ),
                            bn.STDEV => string(lnp.σ),
                            bn.OFFSET => string(lnp.offset),
                            bn.MEAN_IN_REAL_SPACE => bn.FALSE
                           )
                  )
    add_ref_el(el, lnp.param.el)
    lnp.el = el
    return el
end

function LKJPrecisionPriors(lkj::LKJPrecisionXMLElement)
    corr_prior = LKJCorrelationPrior(lkj.off_diag_parm)
    var_prior = LogNormalPriorXMLElement(lkj.diag_param)
    return CompoundXMLElement([corr_prior, var_prior])
end


text:
module Viruses
export meltrate

"""
    meltrate(temp, ddf)

Calculates the melt rate using temperature

Input:
- temp -- temperature [C]
- ddf  -- melt factor [mm/h/C]

Output:
- meltrate [mm/h]
"""
function meltrate(temp, ddf)
    if temp<=0
        # it's freezing, no melt
        0.0
    elseif temp>0
        temp*ddf
    end
end

end # module


text:
"""
    MatHeatDiffModule

Module for linear heat diffusion material models.
"""
module MatHeatDiffModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.MatModule: AbstractMat
using FinEtools.MatrixUtilityModule: mulCAB!

"""
    MatHeatDiff{MTAN<:Function, MUPD<:Function} <: AbstractMat

Type of material model for heat diffusion.
"""
struct MatHeatDiff{MTAN<:Function, MUPD<:Function} <: AbstractMat
	thermal_conductivity::Array{FFlt, 2};# Thermal conductivity
	specific_heat::FFlt;# Specific heat per unit volume
	mass_density::FFlt # mass density
	tangentmoduli!::MTAN
	update!::MUPD
end

"""
    MatHeatDiff(thermal_conductivity)

Construct material model for heat diffusion.

Supply the matrix of thermal conductivity constants.
"""
function MatHeatDiff(thermal_conductivity)
    return MatHeatDiff(thermal_conductivity, 0.0, 0.0, tangentmoduli!, update!)
end

"""
    MatHeatDiff(thermal_conductivity, specific_heat)

Construct material model for heat diffusion.

Supply the matrix of thermal conductivity constants.
"""
function MatHeatDiff(thermal_conductivity, specific_heat)
    return MatHeatDiff(thermal_conductivity, specific_heat, 0.0, tangentmoduli!, update!)
end

"""
    tangentmoduli!(self::MatHeatDiff, kappabar::FFltMat, t::FFlt, dt::FFlt, loc::FFltMat, label::FInt)

Calculate the thermal conductivity matrix.

- `kappabar` = matrix of thermal conductivity (tangent moduli) in material
  coordinate system, supplied as a buffer and overwritten.
"""
function tangentmoduli!(self::MatHeatDiff, kappabar::FFltMat, t::FFlt = 0.0, dt::FFlt = 0.0, loc::FFltMat = reshape(FFlt[],0,0), label::FInt = 0)
    copyto!(kappabar, self.thermal_conductivity);
    return kappabar
end

"""
    update!(self::MatHeatDiff, heatflux::FFltVec, output::FFltVec, gradT::FFltVec, t::FFlt= 0.0, dt::FFlt= 0.0,  loc::FFltMat=FFltMat[], label::FInt=0, quantity=:nothing)

Update material state.

# Arguments
- `gradT` = thermal gradient vector,
- `t` = current time,
- `dt` = current time step,
- `loc` = location of the quadrature point in global Cartesian coordinates,
- `label` = label of the finite element in which the quadrature point is located.
- `quantity` = quantity to be output (`:heatflux`)

# Output
- `heatflux` = heat flux vector, allocated by the caller with a size of
  the embedding space. The components of the heat flux vector are
  calculated and stored in the `heatflux` vector.
- `output` =  array which is (if necessary) allocated  in an appropriate size, filled with the output quantity, and returned.
"""
function update!(self::MatHeatDiff, heatflux::FFltVec, output::FFltVec, gradT::FFltVec, t::FFlt= 0.0, dt::FFlt= 0.0, loc::FFltMat=reshape(FFlt[],0,0), label::FInt=0, quantity=:nothing)
	sdim = size(self.thermal_conductivity, 2)
    @assert length(heatflux) == sdim
    mulCAB!(heatflux, self.thermal_conductivity, -gradT);
    if quantity == :nothing
        #Nothing to be copied to the output array
    elseif quantity == :heatflux
        (length(output) >= sdim) || (output = zeros(sdim)) # make sure we can store it
        copyto!(output, heatflux);
    end
    return output
end

end

max_stars_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff

max_issues_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff

max_forks_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff


text:
TEST_PARAMETERS = Dict(
    "data" => rand(2, 50),
    "labels" => fill(:U, 50),
    "params" => Dict(
        "C" => 0.1,
        "gamma" => 4,
        "classifier" => "SVDDneg",
        "query_strategy" => "RandomPQs"),
    "query_history" => [],
    "subspaces" => [[1, 2]],
    "subspace_grids" => [hcat([[x,y] for x in 1:10 for y in 1:10]...)])

RAW_PARAMETERS = JSON.parse(JSON.json(TEST_PARAMETERS))


text:
using ParallelPOMCPOW
using Base.Test

using POMCPOW
using POMDPs
using POMDPModels
using POMDPToolbox

m = BabyPOMDP()

ps = POMCPOWSolver(tree_queries=100_000, max_depth=10, enable_action_pw=false)
solver = ParallelPOMCPOWSolver(ps, 50)

planner = solve(solver, m)

action(planner, initial_state_distribution(m))
action(planner, initial_state_distribution(m))
@time a, info = action_info(planner, initial_state_distribution(m))
@show info[:search_time_us]/1e6
@show info[:tree_queries]

splanner = solve(ps, m)
action(splanner, initial_state_distribution(m))
@time a, info = action_info(splanner, initial_state_distribution(m))
@show info[:search_time_us]/1e6


text:
using Test
using TestSetExtensions

@testset "All the tests" begin
    @includetests ARGS
end


text:
abstract type RDataException <: Exception end

struct CodecMissingError <: RDataException
    formatName::Symbol
end

function Base.showerror(io::IO, e::CodecMissingError)
    print(io, string(
        "$(typeof(e)): Codec$(e.formatName) package is required to read ",
        "$(e.formatName)-compressed RData files. Run ",
        "Pkg.add(\"Codec$(e.formatName)\") to install it. Then in your code call ",
        "\"using Codec$(e.formatName)\" before \"using RData\"."
    ))
end


text:
struct SubProblem{T <: AbstractFloat, A <: AbstractVector, PT <: AbstractPenaltyterm}
    id::Int
    probability::T
    optimizer::MOI.AbstractOptimizer
    objective::MOI.AbstractScalarFunction

    decisions::Decisions
    projection_targets::Vector{MOI.VariableIndex}
    ξ::Vector{Decision{T}}

    x::A
    ρ::A

    penaltyterm::PT

    function SubProblem(model::JuMP.Model,
                        id::Integer,
                        π::AbstractFloat,
                        penaltyterm::AbstractPenaltyterm)
        T = typeof(π)
        # Get optimizer backend and initial objective
        optimizer = backend(model)
        F = MOI.get(optimizer, MOI.ObjectiveFunctionType())
        objective = MOI.get(optimizer, MOI.ObjectiveFunction{F}())
        # Ensure that we use AffineDecisionFunction
        if !(F <: AffineDecisionFunction)
            F = AffineDecisionFunction{Float64}
            MOI.set(optimizer, MOI.ObjectiveFunction{F}(), convert(F, objective))
        end
        # Get decisions
        decisions = get_decisions(model)::Decisions
        # Optimize wait-and-see model to generate
        # initial decision
        MOI.optimize!(optimizer)
        status = MOI.get(optimizer, MOI.TerminationStatus())
        if !(status ∈ StochasticPrograms.AcceptableTermination)
            error("Initial wait-and-see problem could not be solved in subproblem $id, returned status $status.")
        end
        x₀ = map(decisions.undecided) do vi
            T(MOI.get(optimizer, MOI.VariablePrimal(), vi))
        end
        A = typeof(x₀)
        ξ = map(x₀) do val
            Decision(val, T)
        end
        # Penalty term
        PT = typeof(penaltyterm)
        subproblem = new{T,A,PT}(id,
                                 π,
                                 optimizer,
                                 objective,
                                 decisions,
                                 Vector{MOI.VariableIndex}(undef, length(x₀)),
                                 ξ,
                                 x₀,
                                 zero(x₀),
                                 penaltyterm)
        return subproblem
    end
end

struct SubproblemSolution{T}
    status::MOI.TerminationStatusCode
    value::T
end
function Base.:+(lhs::SubproblemSolution{T}, rhs::SubproblemSolution{T}) where T
    val = lhs.value + rhs.value
    if lhs.status == rhs.status
        return SubproblemSolution(lhs.status, val)
    end
    # Ensure that non-optimal status is propagated
    if lhs.status in AcceptableTermination
        return SubproblemSolution(rhs.status, val)
    end
    if rhs.status in AcceptableTermination
        return SubproblemSolution(lhs.status, val)
    end
    # Let lhs dictate end status
    return SubproblemSolution(lhs.status, val)
end
Base.zero(::Type{SubproblemSolution{T}}) where T = SubproblemSolution(MOI.OPTIMAL, zero(T))

function initialize!(subproblem::SubProblem, penalty::AbstractFloat)
    # Add projection targets
    add_projection_targets!(subproblem)
    # Initialize penalty
    initialize_penaltyterm!(subproblem.penaltyterm,
                            subproblem.optimizer,
                            penalty / 2,
                            subproblem.decisions.undecided,
                            subproblem.projection_targets)
end

function add_projection_targets!(subproblem::SubProblem)
    ξ = subproblem.ξ
    model = subproblem.optimizer
    for i in eachindex(ξ)
        name = add_subscript(:ξ, i)
        var_index, _ = MOI.add_constrained_variable(model, StochasticPrograms.SingleKnownSet(ξ[i]))
        set_known_decision!(subproblem.decisions, var_index, ξ[i])
        MOI.set(model, MOI.VariableName(), var_index, name)
        subproblem.projection_targets[i] = var_index
    end
    return nothing
end

function update_subproblem!(subproblem::SubProblem, ξ::AbstractVector, r::AbstractFloat)
    x = subproblem.x
    ρ = subproblem.ρ
    ρ .= ρ + r * (x - ξ)
    return nothing
end
update_subproblems!(subproblems::Vector{<:SubProblem}, ξ::AbstractVector, r::AbstractFloat) =
    map(prob -> update_subproblem!(prob, ξ, r), subproblems)

function reformulate_subproblem!(subproblem::SubProblem, ξ::AbstractVector, r::AbstractFloat)
    model = subproblem.optimizer
    f = subproblem.objective
    F = AffineDecisionFunction{Float64}
    # Update dual penalty
    for (i,vi) in enumerate(subproblem.decisions.undecided)
        j = if typeof(f) <: AffineDecisionFunction
            j = something(findfirst(t -> t.variable_index == vi,
                                    f.decision_part.terms), 0)
        else
            j = 0
        end
        coefficient = iszero(j) ? 0.0 : f.decision_part.terms[j].coefficient
        MOI.modify(model, MOI.ObjectiveFunction{F}(),
                   StochasticPrograms.DecisionCoefficientChange(vi, coefficient + subproblem.ρ[i]))
    end
    # Update projection targets
    for i in eachindex(ξ)
        subproblem.ξ[i].value = ξ[i]
    end
    # Update penalty
    update_penaltyterm!(subproblem.penaltyterm,
                        model,
                        r / 2,
                        subproblem.decisions.undecided,
                        subproblem.projection_targets)
    return nothing
end

function restore_subproblem!(subproblem::SubProblem)
    model = subproblem.optimizer
    # Delete penalty-term
    remove_penalty!(subproblem.penaltyterm, model)
    # Delete projection targets
    for var in subproblem.projection_targets
        StochasticPrograms.remove_known_decision!(subproblem.decisions, var)
        MOI.delete(model, var)
    end
    empty!(subproblem.projection_targets)
    # Restore objective
    f = subproblem.objective
    F = typeof(f)
    MOI.set(model, MOI.ObjectiveFunction{F}(), f)
    return nothing
end

function (subproblem::SubProblem{T})(ξ::AbstractVector) where T <: AbstractFloat
    MOI.optimize!(subproblem.optimizer)
    status = MOI.get(subproblem.optimizer, MOI.TerminationStatus())
    if status ∈ StochasticPrograms.AcceptableTermination
        subproblem.x .= _get_iterate(subproblem)
        return SubproblemSolution(status, T(_objective_value(subproblem)))
    elseif status == MOI.INFEASIBLE
        val = MOI.get(subproblem.optimizer, MOI.ObjectiveSense()) == MOI.MAX_SENSE ? -Inf : Inf
        return SubproblemSolution(status, T(val))
    elseif status == MOI.DUAL_INFEASIBLE
        val = MOI.get(subproblem.optimizer, MOI.ObjectiveSense()) == MOI.MAX_SENSE ? Inf : -Inf
        return SubproblemSolution(status, T(val))
    else
        return SubproblemSolution(status, NaN)
    end
end

function _objective_value(subproblem::SubProblem)
    objective = subproblem.objective
    obj_val = MOIU.eval_variables(objective) do vi
        MOI.get(subproblem.optimizer, MOI.VariablePrimal(), vi)
    end
    return subproblem.probability * obj_val
end

function _get_iterate(subproblem::SubProblem)
    return map(subproblem.decisions.undecided) do vi
        MOI.get(subproblem.optimizer, MOI.VariablePrimal(), vi)
    end
end