text:
export x̂ , ŷ , ẑ
export HIGH , LOW
export EVERYWHERE

struct Direction{T} end
const x̂ = Direction{1}()
const ŷ = Direction{2}()
const ẑ = Direction{3}()

struct GridType{T} end
const p̂ = GridType{:Primal}()
const d̂ = GridType{:Dual}()

abstract type Side end
struct High <: Side end
struct Low  <: Side end

const HIGH = High()
const LOW = Low()

const Float = typeof(0.0)
const CFloat = Complex{Float}
const NT = NTuple

const EVERYWHERE = x -> true


text:
module UpdateJulia

export update_julia

using JSON, Suppressor

## Version fetching ##
const last_fetched = Ref(0.0)
const versions = Ref{Dict{VersionNumber, Dict{String, Any}}}()
const nightly_version = Ref{VersionNumber}()
function fetch()
    download_delete("https://julialang-s3.julialang.org/bin/versions.json") do file
        open(file) do io
            versions[] = Dict(VersionNumber(k)=>v for (k, v) in JSON.parse(read(io, String)))
        end
    end
    try # theoretically this is typically wasteful, but it is simple, and it is empirically fast
        download_delete("https://raw.githubusercontent.com/JuliaLang/julia/master/VERSION") do file
            open(file) do io
                nightly_version[] = VersionNumber(read(io, String))
            end
        end
    catch
        if !isassigned(nightly_version)
            v = maximum(keys(UpdateJulia.versions[]))
            nightly_version[] = VersionNumber("$(v.major).$(v.minor+1).0-DEV")
        end
    end
    versions[][nightly_version[]] = Dict()
    last_fetched[] = time()
    nothing
end

"""
    prefer(v1, v2)

Whether to prefer v1 over v2.

Not part of the public API.
"""
function prefer(v1::Union{VersionNumber, Missing}, v2::Union{VersionNumber, Missing})
    ismissing(v1) && return false
    ismissing(v2) && return true
    !is_stable(v1) && is_stable(v2) && return false
    is_stable(v1) && !is_stable(v2) && return true
    is_nightly(v1) && !is_nightly(v2) && return false
    !is_nightly(v1) && is_nightly(v2) && return true
    return v1 > v2
end
is_stable(v::VersionNumber) = isempty(v.prerelease)
is_nightly(v::VersionNumber) = v.prerelease == ("DEV",)

function latest(prefix="")
    kys = collect(filter(v->startswith(string(v), prefix), keys(versions[])))
    isempty(kys) && throw(ArgumentError("No released versions starting with \"$prefix\""))
    partialsort!(kys, 1, lt=prefer)
end

## Keyword argument helpers ##
macro os(windows, apple=windows, freebsd=apple, linux=freebsd, other=linux)
    @static if Sys.iswindows()
        windows
    elseif Sys.isapple()
        apple
    elseif @static VERSION >= v"1.1" && Sys.isfreebsd()
        freebsd
    elseif Sys.islinux()
        linux
    else
        other
    end
end

function default_install_location(systemwide, v)
    default = if systemwide
        @os "\\Program Files" "/Applications" "/opt"
    else
        (joinpath(homedir(), (@os "AppData\\Local\\Programs" "Applications" ".local")))
    end
    current = dirname(dirname(
        @static Sys.isapple() ? dirname(dirname(dirname(Base.Sys.BINDIR))) : Base.Sys.BINDIR))

    if systemwide == startswith(Base.Sys.BINDIR, homedir())
        default # installing at a different systemwide level
    elseif current == default || (@static Sys.isunix() && (@static !Sys.isapple() && (
                                ( systemwide && startswith(current, "/opt/julia")) ||
                                (!systemwide && startswith(current, homedir()   )) )))
        current # current is already a conventional location
    elseif v === nothing
        "$default\" or \"$current"
    else
        println("julia-$VERSION is currently installed in $current. Install julia-$v in $default instead? Y/n")
        response = readline(stdin)
        if isempty(response) || uppercase(first(response)) == 'Y'
            default
        elseif uppercase(first(response)) == 'N'
            current
        else
            error("Unknown input: $response. Expected 'y' or 'n'.")
        end
    end
end

function v_url(version_str, os_str, arch_str, prefer_gui)
    if version_str == "nightly"
        arch_dir = arch_str == "aarch64" ? "aarch64" : "x$(Sys.WORD_SIZE)"
        arch_append = arch_str == "aarch64" ? "aarch64" : "$(Sys.WORD_SIZE)"
        os_append = os_str == "winnt" ? "win" : os_str
        extension = @static Sys.iswindows() ?  (prefer_gui ? "exe" : "zip") : (@static Sys.isapple() ?  "dmg" : "tar.gz")

        nightly_version[], "https://julialangnightlies-s3.julialang.org/bin/$os_str/$arch_dir/julia-latest-$os_append$arch_append.$extension"
    else
        v = latest(version_str)
        v == nightly_version[] && return v_url("nightly", os_str, arch_str, prefer_gui)

        options = filter(x -> x["os"] == os_str && x["arch"] == arch_str, versions[][v]["files"])
        isempty(options) && error("No valid download for \"$version_str\" matching os=\"$os_str\" and arch=\"$arch_str\"")
        sort!(options, by = x->x["kind"], rev=prefer_gui)

        v, first(options)["url"]
    end
end

## Main ##
"""
    update_julia(version::AbstractString="")

Install the latest version of Julia from https://julialang.org

If `version` is provided, installs the latest version that starts with `version`.
If `version == "nightly"`, then installs the bleeding-edge nightly version.

# Keyword Arguments
Behavior flags
- `dry_run = false` skip the actual download and installation
- `verbose = dry_run` print the final value of all arguments
$(Sys.iswindows() ? "- `prefer_gui = false` if true, prefer using the \"installer\" version rather than downloading the \"archive\" version and letting UpdateJulia automatically install it. Incompatible with `migrate_packages`." : "")
- `migrate_packages = <upgrading to a later version of Julia without an existing global environment>` whether to migrate packages in the default global environment. May be `true`, `false`, or `:force`. Only `:force` will replace an existing Project.toml

Destination
- `aliases = ["julia", "julia-\$(v.major).\$(v.minor)", "julia-\$v"]` which aliases to attempt to create for the installed version of Julia. Regardless, will not replace stable versions with less stable versions or newer versions with older versions of the same stability.
- `systemwide = $(!startswith(Base.Sys.BINDIR, homedir()))` install for all users, `false` only installs for current user.
- `install_location = systemwide ? "$(default_install_location(true, nothing))" : "$(default_install_location(false, nothing))"` directory to put installed binaries
$(Sys.iswindows() ? "" : "- `bin = systemwide ? \"/usr/local/bin\" : \"$(joinpath(homedir(), ".local/bin"))\"` directory to store links to the binaries")

Source
- `os_str = "$(@os "winnt" "mac" "freebsd" "linux")"` string representation of the operating system: "linux", "mac", "winnt", or "freebsd".
- `arch = "$(string(Base.Sys.ARCH))"` string representation of the CPU architecture: "x86_64", "i686", "aarch64", "armv7l", or "powerpc64le".

- `v = ...` the `VersionNumber` to install
- `url = ...` URL to download that version from, if you explicitly set `url`, also explicitly set `v` lest they differ
"""
function update_julia(version::AbstractString="";
    os_str = (@os "winnt" "mac" "freebsd" "linux"),
    arch = string(Base.Sys.ARCH),
    prefer_gui = false,
    fetch = time() > last_fetched[] + 60 * 15, # 15 minutes
    _v_url = ((fetch && UpdateJulia.fetch()); v_url(version, os_str, arch, prefer_gui)),
    v = first(_v_url),
    migrate_packages = (VERSION.major, VERSION.minor) < (v.major, v.minor) && !ispath(joinpath(first(Base.DEPOT_PATH), "environments", "v$(v.major).$(v.minor)", "Project.toml")),
    url = last(_v_url),
    aliases = ["julia", "julia-$(v.major).$(v.minor)", "julia-$v"],
    systemwide = !startswith(Base.Sys.BINDIR, homedir()),
    install_location = default_install_location(systemwide, v),
    bin = (@static Sys.iswindows() ? nothing : (systemwide ? "/usr/local/bin" : joinpath(homedir(), ".local/bin"))),
    dry_run = false,
    verbose = dry_run)

    @static VERSION >= v"1.1" && verbose && display(Base.@locals)
    @static Sys.iswindows() && bin !== nothing && (println("bin not supported for windows"); bin=nothing)

    prereport(v) #TODO should this report more info?

    if @static Sys.iswindows() && endswith(url, ".exe")
        prefer_gui || printstyled("A GUI installer was available but not an archive.\n", color=Base.warn_color())
        dry_run && (println("aborting before download & install"); return v)
        download_delete(url) do file
            mv(file, file*".exe")
            try
                printstyled("Launching GUI installer now:\n", color=:green)
                run(`$file.exe`)
            finally
                mv(file*".exe", file)
            end
        end
        return v
    elseif prefer_gui
        printstyled("An archive was available but not a GUI. Installing the archive now.\n", color=Base.warn_color())
    end

    dry_run && (println("aborting before download & install"); return v)
    executable = try
        download_delete(url) do file
            extract(install_location, file, v)
        end
    catch x
        if x isa Base.IOError && systemwide && occursin("permission denied", x.msg)
            printstyled("Permission denied attempting to perform a systemwide installation."*
                "Try again with `systemwide=false` or run with elevated permissions.\n",
                color=Base.error_color())
        end
        rethrow()
    end

    @static if Sys.iswindows()
        # Windows doesn't use the bin system, instead adding each individual julia
        # installation to path. This approach does a worse job of handling multiple versions
        # overwriting each other, but Windows doesn't support symlinks for ordinary users,
        # and hardlinks to julia from different executables don't run, so while possible, it
        # would be much more work to use the more effective unix approach. For now, we
        # create version specific executables, add everything to path, and let the user deal
        # with ordering path entries if they want to.
        bin = dirname(executable)
    end

    isdir(bin) || (println("Making path to $bin"); mkpath(bin))
    ensure_on_path(bin, systemwide, v)

    for command in aliases
        link(executable, bin, command * (@os ".exe" ""), systemwide, v)
    end

    #Try these standard commands to see if they work, even if we didn't create them just now
    report(union(aliases, ["julia", "julia-$(v.major).$(v.minor)", "julia-$v"]), v)

    # Try to migrate after reporting successful julia installation because migrate may fail
    # even if the install succeeds.
    (migrate_packages == :force || migrate_packages) &&
        UpdateJulia.migrate_packages(v, migrate_packages == :force)

    v
end

## Prereport ##
function prereport(v)
    if v == latest()
        printstyled("installing the latest version of julia: $v\n", color = :green)
    elseif "DEV" ∈ v.prerelease
        printstyled("installing julia $v\n"*
        "This version is an experimental development build not recommended for most users. "*
        "The latest official release is $(latest())\n", color = :red)
    else
        printstyled("installing julia $v\n", color = :yellow)
        printstyled("This version is $(v > latest() ? "un-released" : "out of date"). "*
        "The latest official release is $(latest())\n", color = :yellow)
    end
end

## Download ##
function download_delete(f, url)
    # use download instead of Downloads.download for backwards compatibility
    file = @suppress_err download(url)
    try
        f(file)
    finally
        rm("$file")
    end
end

## Extract ##
function extract(install_location, download_file, v)
    isdir(install_location) || (println("Making path to $install_location"); mkpath(install_location))

    @static if Sys.isapple()
        run(`hdiutil attach $download_file`)
        volumes = filter(x->startswith(x, "Julia-$v"), readdir("/Volumes"))
        folder = last(volumes)
        try
            cp("/Volumes/$folder/Julia-$(v.major).$(v.minor).app", "$install_location/$folder.app", force=true)
        finally
            for volume in volumes
                run(`hdiutil detach /Volumes/$volume`)
            end
        end
        "$install_location/$folder.app/Contents/Resources/julia/bin/julia"
    else
        # We must extract to a temporary location instead of directly into install_location
        # because we don't know what the name of the extracted folder is. Specifically, on
        # nightlies, it is (today) julia-db1d2f5891. We need to return the executable
        # location which entails determining this extension.
        extract_location = mktempdir()
        @static if Sys.iswindows()
            run(`powershell.exe -nologo -noprofile -command "& { Add-Type -A 'System.IO.Compression.FileSystem'; [IO.Compression.ZipFile]::ExtractToDirectory('$download_file', '$extract_location'); }"`)
        else
            run(`tar zxf $download_file -C $extract_location`)
        end
        folders = readdir(extract_location)
        @assert length(folders) == 1
        folder = first(folders)
        folder == "julia-$v" || startswith(folder, "julia-") && v.prerelease == ("DEV",) ||
            @warn "Unexpected install folder: $folder"

        mv(joinpath(extract_location, folder), joinpath(install_location, folder), force=true)
        rm(extract_location)

        joinpath(install_location, folder, "bin", @os "julia.exe" "julia")
    end
end

## Link ##
function ensure_on_path(bin, systemwide, v)
    @static if Sys.iswindows()
        # Long term solution
        path = strip(open(io -> read(io, String), `powershell.exe -nologo -noprofile -command "[Environment]::GetEnvironmentVariable(\"PATH\"$(systemwide ? "" : ", \"User\""))"`))
        new_path = insert_path(path, bin, v)
        if path != new_path
            run(`powershell.exe -nologo -noprofile -command "[Environment]::SetEnvironmentVariable(\"PATH\", \"$new_path\"$(systemwide ? "" : ", \"User\""))"`)
            println("Adding $bin to $(systemwide ? "system" : "user") path. Shell/PowerShell restart may be required.")
        end

        # Short term solution
        ENV["PATH"] = insert_path(ENV["PATH"], bin, v)
    else
        # Long term solution
        if !occursin(bin, ENV["PATH"])
            printstyled("Please add $bin to path\n", color=Base.warn_color())
        end

        # Short term solution
        occursin(bin, ENV["PATH"]) || (ENV["PATH"] *= ":$bin")
    end
end

"""
    insert_path(path, entry, v)

Insert path entry after valid & preferred julia but before unknown & unpreferred julia

Insert entry into path following these guidelines
- after versions `prefer`red over `v`
- before versions `v` is `prefer`red over (including unknown versions)
- skip operation if `entry` already meets above guidelines
- before existing entries for `v`
- as late as possible

Not part of the public API
"""
function insert_path(path, entry, v)
    @assert Sys.iswindows()
    entries = split(path, ";")
    keys = map(entries) do entry
        l = skipmissing(try VersionNumber(m[2]) catch; missing end
            for m in eachmatch(r"julias?(-|\\)([0-9.a-zA-Z\-+]*)\\\\", entry))
        (isempty(l) ? missing : maximum(l), occursin("julia", lowercase(entry)))
    end

    # after versions `prefer`red over `v`
    last_better = findlast(k->prefer(k[1], v), keys)
    last_better === nothing && (last_better = 0)

    # before versions `v` is `prefer`red over (including unknown versions) & before existing entries for `v`
    first_worse_or_eq = findfirst(k -> k[2] && !prefer(k[1], v), keys[last_better+1:end])
    first_worse_or_eq === nothing && (first_worse_or_eq = lastindex(entries)+1-last_better)
    first_worse_or_eq += last_better

    # skip operation if `entry` already meets above guidelines
    entry ∈ entries[last_better+1:min(end, first_worse_or_eq)] && return path

    join(insert!(entries, first_worse_or_eq, entry), ";")
end

function link(executable, bin, command, systemwide, v)
    link = joinpath(bin, command)
    @static if Sys.iswindows()
        # Make a hard link from julia.exe to julia-1.6.4.exe within the same directory
        # because that hardlink won't work across directories and a symlink requires
        # privileges
        isfile(link) || run(`cmd.exe -nologo -noprofile /c mklink /H $link $executable`)
    else
        # Make a link from the executable in the install location to a bin shared with other julia versions
        old = version_of(command)
        if !prefer(old, v) # If v is as good or better than old, a symlink is warranted
            run(`ln -sf $executable $link`) # Because force is not available via Base.symlink

            old = version_of(command)
            if prefer(v, old) # A worse symlink has higher precedence

                link = strip(open(x -> read(x, String), `$(@os "which.exe" "which") $command`))
                if !systemwide && !startswith(link, homedir())
                    printstyled("`$command` points to $link, which points to version $old. Not editing $link because this is not a systemwide installation.\n", color=Base.warn_color())
                else
                    printstyled("Replacing $link with a symlink to this installation\n", color=Base.info_color())
                    run(`ln -sf $executable $link`) # Because force is not available via Base.symlink
                end

            end
        end
    end
end

## Copy Packages ##
function migrate_packages(v, force=false)
    envpath = joinpath(first(Base.DEPOT_PATH), "environments")
    source = joinpath(envpath, "v$(VERSION.major).$(VERSION.minor)", "Project.toml")
    target = joinpath(envpath, "v$(v.major).$(v.minor)", "Project.toml")

    if !isfile(source)
        printstyled("Package migration failed because there was no Project.toml to copy from at $source\n", color=Base.warn_color())
    elseif !force && ispath(target)
        printstyled("Package migration failed because $target already exists. To overwrite that file, try again with `migrate_packages = :force` or run `UpdateJulia.migrate_packages(v\"$v\", true)`.\n", color=Base.warn_color())
    else
        mkpath(dirname(target))
        source == target || cp(source, target, force=force)
        run(`julia-$(v.major).$(v.minor) -e "using Pkg; Pkg.activate(\"$(escape_string(dirname(target)))\"); Pkg.update()"`)
        printstyled("Migrated packages from $(VERSION.major).$(VERSION.minor) to $(v.major).$(v.minor)\n", color=:green)
    end
end

## Test ##
function report(commands, v)
    successes = filter(c->version_of(c)==v, commands)
    @assert !isempty(successes)
    printstyled("Success! \`$(join(successes, "\` & \`"))\` now to point to $v\n", color=:green)
end

function version_of(command)
    try
        str = open(f->read(f, String), `$command -v`)
        @assert startswith(str, "julia version ")
        @assert endswith(str, "\n")
        VersionNumber(str[15:end-1])
    catch x
        missing
    end
end

end

ext: jl


text:
## Author: Victor H. Aguiar

## questions: slack

## Import packages, make sure to install them.
## Julia version: v"1.5.2"
import JuMP
import Ipopt

## define a model
robinson=JuMP.Model(Ipopt.Optimizer)

# define variables
#coconuts
@JuMP.variable(robinson,c>=0)
#labor
@JuMP.variable(robinson,l>=0)
#wages
@JuMP.variable(robinson,w>=0)

# define parameters
α1=.4
α2=(1-α1)
#@JuMP.NLparameter(robinson, β == 1/2)
β = .1
A=50
Lbar=100.0

# function prodfunction(x)
#        try
#              A*x^(β)
#       catch
#             NaN
#       end
# end
#
# function prodfunction_prime(x)
#        try
#              A*β*(x)^(β-1)
#       catch
#             NaN
#       end
# end
#
# function prodfunction_prime_prime(x)
#        try
#               A*β^2*x^(β-2)
#       catch
#             NaN
#       end
# end
#
# function prodfunction_prime_prime_prime(x)
#        try
#               A*β^3*x^(β-3)
#       catch
#             NaN
#       end
# end
prodfunction_prime(x) = A*x^β
prodfunction_prime(x) = A*β*x^(β-1)
prodfunction_prime_prime(x) = A*β^2*x^(β-2)
prodfunction_prime_prime_prime(x) = A*β^3*x^(β-3)
# define equations, here in JuMP we must define them as (nonlinear) constraint
JuMP.register(robinson, :prodfunction, 1, prodfunction, prodfunction_prime,
         prodfunction_prime_prime)

JuMP.register(robinson, :prodfunction_prime, 1, prodfunction_prime, prodfunction_prime_prime,
                  prodfunction_prime_prime_prime)

# FOC of the consumer
@JuMP.NLconstraint(robinson,α1/(α2)*(Lbar-l)*w==c)
#@JuMP.NLconstraint(robinson,α1/(α2)*(Lbar-l)*prodfunction_prime(l+ 1e-7)==c)

# Marginal productivity of labor equal to wage
@JuMP.NLconstraint(robinson,prodfunction_prime(l)==w)

#@JuMP.constraint(robinson,c>=0)
#@JuMP.constraint(robinson,l>=0)
#@JuMP.constraint(robinson,w>=0)
# Budget constraint in equilibrium
@JuMP.NLconstraint(robinson,prodfunction(l)==c)

# Solve model
@JuMP.NLobjective(robinson,Min,1)
JuMP.optimize!(robinson)

JuMP.value(c)
JuMP.value(l)
JuMP.value(w)


text:
using ChaosTools
using Test
using StatsBase
using Statistics

println("\nTesting Takens' best estimate")
@testset "Takens best" begin
  @testset "Henon map" begin
    ds = Systems.henon()
    ts = trajectory(ds, 5000)
    x = ts[:, 1]
    X = embed(x, 2, 1)
    D_C, D_C_95u, D_C_95l = takens_best_estimate(X, std(x)/4)
    @test 1.15 < D_C < 1.25
    @test D_C_95u < 1.05*D_C
    @test D_C_95l > 0.95*D_C
  end
  @testset "Lorenz System" begin
    ds = Systems.lorenz()
    tr = trajectory(ds, 2000; dt = 0.1)
    x = tr[:, 1]
    τ = estimate_delay(x, "mi_min", 1:20)
    X = embed(x, 4, τ)
    D_C, D_C_95u, D_C_95l = takens_best_estimate(X, std(x)/4)
    @test 1.85 < D_C < 2.1
    @test D_C_95u < 1.05*D_C
    @test D_C_95l > 0.95*D_C
  end
end


text:
# using Alexya

@testset "Controls" begin

widgets = (
    inspect = Button("Inspect", @hexpand),
    close = Button("Close", @hexpand)
)

inputs = (
    text_field = "Lorem ipsum dolor",
    number = 10,
    δ = (10, -50:50)
)

numbers = (
    complex = 3 + 4im,
    point = Point(-5, 5)
)

others = (
    slider = -20:20,
    φ = (-100:100, 50),
    color = colorant"#000",
    switch = true
)

options = Controls(
    widgets = widgets,
    inputs = inputs,
    numbers = numbers,
    others = others
)

onevent(:clicked, options.widgets.inspect) do
    @show options
end

app = Window("Controsl test", 400, 800) do
    ScrolledWindow() do
        options
    end
end

showall(app)

@test options.widgets.inspect[]     == widgets.inspect
@test options.widgets.close[]       == widgets.close
@test options.inputs.text_field[]   == inputs.text_field
@test options.inputs.number[]       == inputs.number
@test options.inputs.δ[]            == 10.0
@test options.numbers.complex[]     == numbers.complex
@test options.numbers.point[]       == numbers.point
@test options.others.slider[]       == 0.0
@test options.others.φ[]            == 50.0
@test options.others.color[]        == convert(RGBA, others.color)
@test options.others.switch[]       == others.switch

options.inputs.text_field[] = "Okay"
options.inputs.number[] = 5
options.inputs.δ[] = 10
options.numbers.complex[] = 3.14 + 1im
options.numbers.point[] = Point(10, 10)
options.others.slider[] = -9/2
options.others.φ[] = -1/12
options.others.color[] = colorant"#fff"
options.others.switch[] = false

@test options.inputs.text_field[]   == "Okay"
@test options.inputs.number[]       == 5
@test options.inputs.δ[]            == 10
@test options.numbers.complex[]     == 3.14 + 1im
@test options.numbers.point[]       == Point(10, 10)
@test options.others.slider[]       == -9/2
@test options.others.φ[]            == -1/12
@test options.others.color[]        == RGBA(colorant"#fff")
@test options.others.switch[]       == false

destroy(app)

end # testset

@testset "Options Macro" begin

@init "Options test" 400 800

@layout aside(:v, 250)

options = @create @options begin
    widgets = (
        inspect = Button("Inspect", @hexpand),
        close = Button("Close", @hexpand)
    )

    inputs = (
        text_field = "Lorem ipsum dolor",
        number = 10,
        δ = (10, -50:50)
    )

    numbers = (
        complex = 3 + 4im,
        point = Point(-5, 5)
    )

    others = (
        slider = -20:20,
        φ = (-100:100, 50),
        color = colorant"#000",
        switch = true
    )
end

@use function setup()
    #@show options
    noLoop()
end

@test inputs.text_field[]   == "Lorem ipsum dolor"
@test inputs.number[]       == 10.0
@test inputs.δ[]            == 10.0
@test numbers.complex[]     == 3 + 4im
@test numbers.point[]       == Point(-5, 5)
@test others.slider[]       == 0.0
@test others.φ[]            == 50.0
@test others.color[]        == RGBA(0, 0, 0, 1)
@test others.switch[]       == true

start()

end # testset


text:
using CSV, JuMP, MathProgBase, Gurobi, Random, LinearAlgebra, StatsBase, Distributions, Plots, BenchmarkTools, Suppressor

# Formulate optimization problem for Ordinary Least Squares Regression

function OrdinaryLeastSquares(x,y)
#     OLS_model = Model(solver = GurobiSolver(LogFile = append_name,OutputFlag = 0))
    OLS_model = Model(solver = GurobiSolver(OutputFlag = 0))

    # Define variables
    n = size(x,1) # observations
    p = size(x,2) # features

    # Define decision variables
    @variable(OLS_model, β[1:p+1])

    # Define constraints - none

    # Define objective function
    @objective(OLS_model, Min, sum((y[i] - β[1] - β[2:p+1]' * x[i,:])^2 for i = 1:n))
#     return(OLS_model)

    # Run the optimization
    solve(OLS_model)
    betas = getvalue(β)
    objective = getobjectivevalue(OLS_model)
    return(betas)
end

# Formulate optimization problem for Least Quantile Squares Regression
# *********** q: 0, 1, 2, 3, 4

function LeastQuantileSquares(x,y,quantile)
    LQS_model = Model(solver = GurobiSolver(OutputFlag = 0,TimeLimit = 180))
#     LQS_model = Model(solver = GurobiSolver(OutputFlag = 0))

    # Define variables
    M = 10000
    n = size(x,1) # observations
    p = size(x,2) # features
    q = ceil(n*(quantile/4))

    # Define decision variables
    @variable(LQS_model, γ)
    @variable(LQS_model, z[1:n], Bin)
    @variable(LQS_model, μ[1:n])
    @variable(LQS_model, β[1:p+1])
    @variable(LQS_model, r[1:n])
    @variable(LQS_model, a[1:n])

    # Define constraints
    @constraint(LQS_model, [i=1:n], r[i] == y[i] - β[1] - β[2:p+1]'*x[i,:])
    @constraint(LQS_model, [i=1:n], γ >= a[i] - μ[i])
    @constraint(LQS_model, [i=1:n], a[i] >= r[i])
    @constraint(LQS_model, [i=1:n], a[i] >= -r[i])
    @constraint(LQS_model, [i=1:n], μ[i] <= M*(1 - z[i]))
    @constraint(LQS_model, sum(z[i] for i=1:n) == q)
    @constraint(LQS_model, [i=1:n], μ[i] >= 0)

     # Define objective function
    @objective(LQS_model, Min, γ)
#     return(LQS_model)

    # Run the optimization
    solve(LQS_model)
    betas = getvalue(β)
    qth_residual = getvalue(γ)
    return(betas,qth_residual)
end

# Formulate optimization problem for Least Trimmed Squares Regression
# Note: k is a user specified parameter, k <= (n - p - 1)/2

function LeastTrimmedSquares(k,x,y)
    LTS_model = Model(solver = GurobiSolver(OutputFlag = 0,TimeLimit = 500))

    # Define variables
    M = 10000
    n = size(x,1) # observations
    p = size(x,2) # features

    # Define decision variables
    @variable(LTS_model, z[1:n], Bin)
    @variable(LTS_model, β[1:p+1])
    @variable(LTS_model, r[1:n])

    # Define constraints
    @constraint(LTS_model, [i=1:n], β[1] + β[2:p+1]'*x[i,:] - y[i] <= r[i] + M*z[i])
    @constraint(LTS_model, [i=1:n], -β[1] - β[2:p+1]'*x[i,:] + y[i] <= r[i] + M*z[i])
    @constraint(LTS_model, sum(z[i] for i=1:n) <= k)
    @constraint(LTS_model, [i=1:n], r[i] >= 0)

    # Define objective function as QMIP ******** (minimizing over β, z)
    @objective(LTS_model, Min, sum(r[i]^2 for i=1:n))
#     return(LTS_model)

    # Run the optimization
    solve(LTS_model)
    betas = getvalue(β)
    outlier_indicators = getvalue(z)
    return(betas, outlier_indicators)
end

# Formulate optimization problem for Least Trimmed Squares Regression
# Note: k is a user specified parameter, k <= (n - p - 1)/2

function LeastTrimmedSumResiduals(k,x,y)
    LTSR_model = Model(solver = GurobiSolver(OutputFlag = 0,TimeLimit = 500))

    # Define variables
    M = 10000
    n = size(x,1) # observations
    p = size(x,2) # features

    # Define decision variables
    @variable(LTSR_model, z[1:n], Bin)
    @variable(LTSR_model, β[1:p+1])
    @variable(LTSR_model, r[1:n])

    # Define constraints
    @constraint(LTSR_model, [i=1:n], β[1] + β[2:p+1]'*x[i,:] - y[i] <= r[i] + M*z[i])
    @constraint(LTSR_model, [i=1:n], -β[1] - β[2:p+1]'*x[i,:] + y[i] <= r[i] + M*z[i])
    @constraint(LTSR_model, sum(z[i] for i=1:n) <= k)
    @constraint(LTSR_model, [i=1:n], r[i] >= 0)

    # Define objective function as QMIP ******** (minimizing over β, z)
    @objective(LTSR_model, Min, sum(r[i] for i=1:n))
#     return(LTS_model)

    # Run the optimization
    solve(LTSR_model)
    betas = getvalue(β)
    outlier_indicators = getvalue(z)
    return(betas, outlier_indicators)
end

# Formulate optimization problem for Least Trimmed Squares Regression
# Note: k is a user specified parameter, k <= (n - p - 1)/2

function BoundedInfluence(α,x,y)
    BI_model = Model(solver = GurobiSolver(OutputFlag = 0,TimeLimit = 500))

    # Define variables
    M = 10000
    n = size(x,1) # observations
    p = size(x,2) # features
    bkdn = ceil((n-p-1)/2)

    # Define parameters
    x_int = hcat(ones(n),x)
    β_ols = (x_int'*x_int)^-1*x_int'*y
    e_residuals = y-x_int*β_ols
    σ_squared = e_residuals'*e_residuals/(n-p+1)
    H = x*(x'*x)^-1*x'
    h = diag(H)
    δ = ((e_residuals).^2/(σ_squared*(1 .- h)))*(1/(p+1))*(h./(1 .- h))
    σ_h = (e_residuals./(1 .- h).^0.5)
    σ_h_squared = (e_residuals./(1 .- h).^0.5).^2
#     println(δ)

    # Define decision variables
    @variable(BI_model, z[1:n], Bin)
    @variable(BI_model, β[1:p+1])
    @variable(BI_model, r[1:n])

    # Define constraints
    @constraint(BI_model, [i=1:n], β[1] + β[2:p+1]'*x[i,:] - y[i] <= r[i] + M*z[i])
    @constraint(BI_model, [i=1:n], -β[1] - β[2:p+1]'*x[i,:] + y[i] <= r[i] + M*z[i])
    @constraint(BI_model, sum(z[i] for i=1:n) <= bkdn)
    @constraint(BI_model, [i=1:n], r[i] >= 0)

    # Define objective function
    @objective(BI_model, Min, α*sum(r[i]^2 + z[i]*δ[i]^2 for i=1:n) + (1-α)*sum(z[i] for i=1:n))

    # Run the optimization
    solve(BI_model)
    betas = getvalue(β)
    outlier_indicators = getvalue(z)
    return(h, betas, outlier_indicators)
end

star_data =  [  1  4.37  5.23
                2  4.56  5.74
                3  4.26  4.93
                4  4.56  5.74
                5  4.30  5.19
                6  4.46  5.46
                7  3.84  4.65
                8  4.57  5.27
                9  4.26  5.57
                10  4.37  5.12
                11  3.49  5.73
                12  4.43  5.45
                13  4.48  5.42
                14  4.01  4.05
                15  4.29  4.26
                16  4.42  4.58
                17  4.23  3.94
                18  4.42  4.18
                19  4.23  4.18
                20  3.49  5.89
                21  4.29  4.38
                22  4.29  4.22
                23  4.42  4.42
                24  4.49  4.85
                25  4.38  5.02
                26  4.42  4.66
                27  4.29  4.66
                28  4.38  4.90
                29  4.22  4.39
                30  3.48  6.05
                31  4.38  4.42
                32  4.56  5.10
                33  4.45  5.22
                34  3.49  6.29
                35  4.23  4.34
                36  4.62  5.62
                37  4.53  5.10
                38  4.45  5.22
                39  4.53  5.18
                40  4.43  5.57
                41  4.38  4.62
                42  4.45  5.06
                43  4.50  5.34
                44  4.45  5.34
                45  4.55  5.54
                46  4.45  4.98
                47  4.42  4.50 ];

x_star_data = star_data[:,2]
y_star_data = star_data[:,3];

star_betas_OLS = OrdinaryLeastSquares(x_star_data,y_star_data)

star_betas_LQS, star_qth_resid_LQS = LeastQuantileSquares(x_star_data,y_star_data,2)

star_betas_LTS, star_outlier_indicators_LTS = LeastTrimmedSquares(4,x_star_data,y_star_data)

findall(x->x==1, star_outlier_indicators_LTS)

star_betas_LTSR, star_outlier_indicators_LTSR = LeastTrimmedSumResiduals(4,x_star_data,y_star_data)

findall(x->x==1, star_outlier_indicators_LTSR)

h_star, star_betas_BI, star_outlier_indicators_BI = BoundedInfluence(0.4,x_star_data,y_star_data)

findall(x->x==1, star_outlier_indicators_BI)

# Plot all the fits
# gr()
star_data_plot = scatter(x_star_data,y_star_data,series_annotations = text.(1:47, :bottom),
#     title="Fitted Regressions on Hertzsprung-Russell Star Data",
    label="",
    legend=:outerbottom,
    markersize = 2,
    markercolor = :black,
    dpi=300,
    fontfamily="Times",
    size = (600, 500))

f_ols(x_ols) = star_betas_OLS[1] + x_ols.*star_betas_OLS[2]
plot!(f_ols, 3.5, 5,label="Ordinary Least Squares")

f_lqs(x_lqs) = star_betas_LQS[1] + x_lqs.*star_betas_LQS[2]
plot!(f_lqs, 3.5, 5,label="Least Median Squares")

f_lts(x_lts) = star_betas_LTS[1] + x_lts.*star_betas_LTS[2]
plot!(f_lts, 3.5, 5,label="Least Trimmed Squares")

f_ltsr(x_ltsr) = star_betas_LTSR[1] + x_ltsr.*star_betas_LTSR[2]
plot!(f_ltsr, 3.5, 5,label="Least Trimmed Summed Residuals")

f_bi(x_bi) = star_betas_BI[1] + x_bi.*star_betas_BI[2]
plot!(f_bi, 3.5, 5,label="Bounded Influence Least Squares")

savefig("star_data_plot2")

star_data_plot2

stack_loss = [ 1  80  27  89  42
               1  80  27  88  37
               1  75  25  90  37
               1  62  24  87  28
               1  62  22  87  18
               1  62  23  87  18
               1  62  24  93  19
               1  62  24  93  20
               1  58  23  87  15
               1  58  18  80  14
               1  58  18  89  14
               1  58  17  88  13
               1  58  18  82  11
               1  58  19  93  12
               1  50  18  89   8
               1  50  18  86   7
               1  50  19  72   8
               1  50  19  79   8
               1  50  20  80   9
               1  56  20  82  15
               1  70  20  91  15 ];

x_stack_loss = stack_loss[:,2:4];

y_stack_loss = stack_loss[:,5];

stack_loss_betas_OLS = OrdinaryLeastSquares(x_stack_loss,y_stack_loss)

stack_loss_betas_LQS, stack_loss_qth_resid_LQS = LeastQuantileSquares(x_stack_loss,y_stack_loss,2)

stack_loss_betas_LTS, stack_loss_outlier_indicators_LTS = LeastTrimmedSquares(4,x_stack_loss,y_stack_loss)

index = findall(x->x==1, stack_loss_outlier_indicators_LTS)

stack_loss_betas_LTSR, stack_loss_outlier_indicators_LTSR = LeastTrimmedSumResiduals(4,x_stack_loss,y_stack_loss)

findall(x->x==1, stack_loss_outlier_indicators_LTSR)

h_stack_loss, stack_loss_betas_BI, stack_loss_outlier_indicators_BI = BoundedInfluence(0.1,x_stack_loss,y_stack_loss)

findall(x->x==1, stack_loss_outlier_indicators_BI)

h_stack_loss

function standardize(df)
    n, p = size(df)
    for j in 1:p
        mean1 = mean(df[:,j])
        std1 = std(df[:,j])
        for i in 1:n
            df[i,j] = (df[i,j] - mean1)/std1
        end
    end
    return df
end

df1 = CSV.read("sds0.csv", header=false);

y1 = df1[:,1]
X1 = convert(Matrix, df1[:,2:end]);

size(X1)

t_OLS_df1 = @suppress @benchmark df1_betas_OLS = OrdinaryLeastSquares(X1,y1)

time(median(t_OLS_df1))/1000000000

t_LQS_df1 = @suppress @benchmark LeastQuantileSquares(X1,y1,2)

time(median(t_LQS_df1))/1000000000

t_LTS_df1 = @suppress @benchmark df1_LTS_betas, df1_LTS_outlier_indicators = LeastTrimmedSquares(0,X1,y1)

time(median(t_LTS_df1))/1000000000

t_LTSR_df1 = @suppress @benchmark df1_LTSR_betas, df1_LTSR_outlier_indicators = LeastTrimmedSumResiduals(0,X1,y1)

time(median(t_LTSR_df1))/1000000000

df1_betas_OLS = OrdinaryLeastSquares(X1,y1)

df1_LQS_betas, df1_LQS_qth_resid = LeastQuantileSquares(X1,y1,2)

df1_LTS_betas, df1_LTS_outlier_indicators = LeastTrimmedSquares(0,X1,y1)

df1_LTSR_betas, df1_LTSR_outlier_indicators = LeastTrimmedSumResiduals(0,X1,y1)

df1_BI_betas, df1_BI_outlier_indicators = BoundedInfluence(0.00001,X1,y1)

df2 = CSV.read("sds1.csv", header=false);

y2 = df2[:,1]
X2 = convert(Matrix, df2[:,2:end]);

size(X2)

sum(X2)

mean(X2)

t_OLS_df2 = @suppress @benchmark df2_betas_OLS = OrdinaryLeastSquares(X2,y2)

time(median(t_OLS_df2))/1000000000

t_LQS_df2 = @suppress @benchmark df2_LQS_betas, df2_LQS_qth_resid = LeastQuantileSquares(X2,y2,2)

time(median(t_LQS_df2))/1000000000

t_LTS_df2 = @suppress @benchmark df2_LTS_betas, df2_LTS_outlier_indicators = LeastTrimmedSquares(5,X2,y2)

time(median(t_LTS_df2))/1000000000

t_LTSR_df2 = @suppress @benchmark df2_LTSR_betas, df2_LTSR_outlier_indicators = LeastTrimmedSumResiduals(5,X2,y2)

time(median(t_LTSR_df2))/1000000000

t_BI_df2 = @suppress @benchmark df2_BI_betas, df2_BI_outlier_indicators = BoundedInfluence(0.4,X2,y2)

time(median(t_BI_df2))/1000000000

df2_betas_OLS = OrdinaryLeastSquares(X2,y2)

df2_LQS_betas, df2_LQS_qth_resid = LeastQuantileSquares(X2,y2,2)

df2_LTS_betas, df2_LTS_outlier_indicators = LeastTrimmedSquares(5,X2,y2)

findall(x->x==1, df2_LTS_outlier_indicators)

df2_LTSR_betas, df2_LTSR_outlier_indicators = LeastTrimmedSumResiduals(5,X2,y2)

findall(x->x==1, df2_LTSR_outlier_indicators)

df2_BI_betas, df2_BI_outlier_indicators = BoundedInfluence(0.4,X2,y2)

findall(x->x==1, df2_BI_outlier_indicators)

df3 = CSV.read("sds2.csv", header=false);

y3 = df3[:,1]
X3 = convert(Matrix, df3[:,2:end]);

t_OLS_df3 = @suppress @benchmark df3_betas_OLS = OrdinaryLeastSquares(X3,y3)

time(median(t_OLS_df3))/1000000000

t_LQS_df3 = @suppress @benchmark df3_LQS_betas, df3_LQS_qth_resid = LeastQuantileSquares(X3,y3,2)

time(median(t_LQS_df3))/1000000000

t_LTS_df3 = @suppress @benchmark df3_LTS_betas, df3_LTS_outlier_indicators = LeastTrimmedSquares(9,X3,y3)

time(median(t_LTS_df3))/1000000000

t_LTSR_df3 = @suppress @benchmark df3_LTSR_betas, df3_LTSR_outlier_indicators = LeastTrimmedSumResiduals(9,X3,y3)

time(median(t_LTSR_df3))/1000000000

t_BI_df3 = @suppress @benchmark df3_BI_betas, df3_BI_outlier_indicators = BoundedInfluence(0.82,X3,y3)

time(median(t_BI_df3))/1000000000

df3_betas_OLS = OrdinaryLeastSquares(X3,y3)

df3_LQS_betas, df3_LQS_qth_resid = LeastQuantileSquares(X3,y3,2)

df3_LTS_betas, df3_LTS_outlier_indicators = LeastTrimmedSquares(9,X3,y3)

findall(x->x==1, df3_LTS_outlier_indicators)

df3_LTSR_betas, df3_LTSR_outlier_indicators = LeastTrimmedSumResiduals(9,X3,y3)

findall(x->x==1, df3_LTSR_outlier_indicators)

df3_BI_betas, df3_BI_outlier_indicators = BoundedInfluence(0.82,X3,y3)

findall(x->x==1, df3_BI_outlier_indicators)

df4 = CSV.read("sds3.csv", header=false);

y4 = df4[:,1]
X4 = convert(Matrix, df4[:,2:end]);

size(X4)

t_OLS_df4 = @suppress @benchmark df4_betas_OLS = OrdinaryLeastSquares(X4,y4)

time(median(t_OLS_df4))/1000000000

t_LQS_df4 = @suppress @benchmark df4_LQS_betas, df4_LQS_qth_resid = LeastQuantileSquares(X4,y4,2)

time(median(t_LQS_df4))/1000000000

t_LTS_df4 = @suppress @benchmark df4_LTS_betas, df4_LTS_outlier_indicators = LeastTrimmedSquares(24,X4,y4)

time(median(t_LTS_df4))/1000000000

t_LTSR_df4 = @suppress @benchmark df4_LTSR_betas, df4_LTSR_outlier_indicators = LeastTrimmedSumResiduals(24,X4,y4)

time(median(t_LTSR_df4))/1000000000

t_BI_df4 = @suppress @benchmark df4_BI_betas, df4_BI_outlier_indicators = BoundedInfluence(0.85,X4,y4)

time(median(t_BI_df4))/1000000000

df4_betas_OLS = OrdinaryLeastSquares(X4,y4)

df4_LQS_betas, df4_LQS_qth_resid = LeastQuantileSquares(X4,y4,2)

df4_LTS_betas, df4_LTS_outlier_indicators = LeastTrimmedSquares(24,X4,y4)

findall(x->x==1, df4_LTS_outlier_indicators)

df4_LTSR_betas, df4_LTSR_outlier_indicators = LeastTrimmedSumResiduals(24,X4,y4)

findall(x->x==1, df4_LTSR_outlier_indicators)

df4_BI_betas, df4_BI_outlier_indicators = BoundedInfluence(0.89,X4,y4)

findall(x->x==1, df4_BI_outlier_indicators)

# df5 = CSV.read("./SyntheticData_&_Scripts/df5.csv", header=false)
# y5 = df5[:,1]
# X5 = convert(Matrix, df5[:,2:end]);

# size(X5)

# t_OLS_df5 = @suppress @benchmark df5_betas_OLS = OrdinaryLeastSquares(X5,y5)

# time(median(t_OLS_df5))/1000000000

# t_LQS_df5 = @suppress @benchmark df5_LQS_betas, df5_LQS_qth_resid = LeastQuantileSquares(X5,y5,2)

# time(median(t_LQS_df2))/1000000000

# t_LTS_df5 = @suppress @benchmark df5_LTS_betas, df5_LTS_outlier_indicators = LeastTrimmedSquares(4,X5,y5)

# time(median(t_LTS_df2))/1000000000

# t_ALTS_df5 = @suppress @benchmark df5_LTSR_betas, df5_LTSR_outlier_indicators = LeastTrimmedSumResiduals(4,X5,y5)

# time(median(t_LTSR_df2))/1000000000

# t_BI_df5 = @suppress @benchmark df5_BI_betas, df5_BI_outlier_indicators = BoundedInfluence(0.1,X5,y5)

# time(median(t_BI_df2))/1000000000

# df5_betas_OLS = OrdinaryLeastSquares(X5,y5)

# df5_LQS_betas, df5_LQS_qth_resid = LeastQuantileSquares(X5,y5,2)

# df5_LTS_betas, df5_LTS_outlier_indicators = LeastTrimmedSquares(5,X5,y5)

# findall(x->x==1, df5_LTS_outlier_indicators)

# df5_LTSR_betas, df5_LTSR_outlier_indicators = LeastTrimmedSumResiduals(5,X5,y5)

# findall(x->x==1, df5_LTSR_outlier_indicators)

# df5_BI_betas, df5_BI_outlier_indicators = BoundedInfluence(0.1,X5,y5)

# findall(x->x==1, df5_BI_outlier_indicators)

function standardize(df)
    n, p = size(df)
    for j in 1:p
        mean1 = mean(df[:,j])
        std1 = std(df[:,j])
        for i in 1:n
            df[i,j] = (df[i,j] - mean1)/std1
        end
    end
    return df
end

aq = CSV.read("Results/qsar_aquatic_toxicity.csv", header=false)
#Convert to matrices
Random.seed!(95)
aq = convert(Matrix, aq)
aq = aq[shuffle(1:end), :]

size(aq)

aq = standardize(aq);

# Split into 50%, 25%, and 25%
train_x = aq[1:40, 4:8]
train_y = aq[1:40, 9]

valid_x = aq[41:60, 4:8]
valid_y = aq[41:60, 9]

test_x = aq[61:80, 4:8]
test_y = aq[61:80, 9];

betas_OLS = OrdinaryLeastSquares(train_x,train_y)

betas_LQS, LQS_qth_resid = LeastQuantileSquares(train_x,train_y,2)

betas_LTS, LTS_outlier_indicators = LeastTrimmedSquares(4,train_x,train_y)

betas_LTSR, LTSR_outlier_indicators = LeastTrimmedSumResiduals(4,train_x,train_y)

betas_BI, BI_outlier_indicators = BoundedInfluence(0.4,train_x,train_y)

mean((valid_y .- betas_OLS[1] .- valid_x*betas_OLS[2:end]).^2)
mean((valid_y .- betas_LQS[1] .- valid_x*betas_LQS[2:end]).^2)
mean((valid_y .- betas_LTS[1] .- valid_x*betas_LTS[2:end]).^2)
mean((valid_y .- betas_LTSR[1] .- valid_x*betas_LTSR[2:end]).^2)
mean((valid_y .- betas_BI[1] .- valid_x*betas_BI[2:end]).^2)

full_x = vcat(train_x, valid_x)
full_y = vcat(train_y, valid_y);

time_OLS = @suppress @benchmark OrdinaryLeastSquares(full_x,full_y)
time(median(time_OLS))/1000000000

time_LQS = @suppress @benchmark betas_LQS, LQS_qth_resid = LeastQuantileSquares(full_x,full_y,2)
time(median(time_LQS))/1000000000

time_LTS = @suppress @benchmark  LeastTrimmedSquares(4,full_x,full_y)
time(median(time_LTS))/1000000000

time_LTSR = @suppress @benchmark LeastTrimmedSumResiduals(4,full_x,full_y)
time(median(time_LTSR))/1000000000

time_BI = @suppress @benchmark BoundedInfluence(0.4,full_x,full_y)
time(median(time_BI))/1000000000

max_stars_repo_path: DEC 9 - FINAL MIO Formulations - OLS, LMS, LTS, LTSR, BILS.jl


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## dummy stub for https://github.com/JuliaBinaryWrappers/p7zip_jll.jl
module p7zip_jll

const PATH_list = String[]
const LIBPATH_list = String[]

export p7zip

# These get calculated in __init__()
PATH = Ref("")
LIBPATH = Ref("")
artifact_dir = ""
p7zip_path = ""
if Sys.iswindows()
    const p7zip_exe = "7z.exe"
else
    const p7zip_exe = "7z"
end

# These functions look a little strange, but they're mimicking the JLLWrappers signature
p7zip(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true) = f(p7zip_path)
p7zip(; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true) = Cmd([p7zip_path])

function init_p7zip_path()
    # Prefer our own bundled p7zip, but if we don't have one, pick it up off of the PATH
    # If this is an in-tree build, `7z` will live in `bin`.  Otherwise, it'll be in `libexec`
    for bundled_p7zip_path in (joinpath(Sys.BINDIR, Base.LIBEXECDIR, p7zip_exe),
                               joinpath(Sys.BINDIR, p7zip_exe))
        if isfile(bundled_p7zip_path)
            global p7zip_path = abspath(bundled_p7zip_path)
            return
        end
    end
    global p7zip_path = Sys.which(p7zip_exe)
end

function __init__()
    global artifact_dir = dirname(Sys.BINDIR)
    global LIBPATH[] = joinpath(Sys.BINDIR, Base.LIBDIR, "julia")
    init_p7zip_path()
end

is_available() = true

end  # module p7zip_jll

max_stars_repo_name: saberdy/julia

max_issues_repo_name: saberdy/julia

max_forks_repo_name: saberdy/julia


text:
"""
$(TYPEDEF)
"""
abstract type AbstractDesign end

"""
$(TYPEDEF)
"""
abstract type AbstractScreeningDesign <: AbstractDesign end

"""
$(TYPEDEF)
"""
abstract type AbstractFactorialDesign <: AbstractDesign end

"""
$(TYPEDEF)
"""
abstract type AbstractOptimalDesign <: AbstractDesign end

"""
$(TYPEDEF)
"""
abstract type AbstractRandomDesign <: AbstractDesign end

"""
$(TYPEDEF)

Encapsulates  a  Plackett-Burman  screening  design  constructed  using  Paley's
method.  Factor  levels are  encoded as  `:high` and  `:low` symbols,  and extra
dummy variables, possibly  aliasing interactions between factors,  will be added
to pad a design.

$(TYPEDFIELDS)
"""
struct PlackettBurman <: AbstractScreeningDesign
    matrix::DataFrame
    factors::Tuple
    dummy_factors::Tuple
    formula::FormulaTerm
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> PlackettBurman(@formula(y ~ x1 + x2 + x3 + x4))
PlackettBurman
Dimension: (8, 7)
Factors: (:x1, :x2, :x3, :x4)
Dummy Factors: (:dummy1, :dummy2, :dummy3)
Formula: y ~ -1 + x1 + x2 + x3 + x4 + dummy1 + dummy2 + dummy3
Design Matrix:
8×7 DataFrame
│ Row │ x1    │ x2    │ x3    │ x4    │ dummy1 │ dummy2 │ dummy3 │
│     │ Int64 │ Int64 │ Int64 │ Int64 │ Int64  │ Int64  │ Int64  │
├─────┼───────┼───────┼───────┼───────┼────────┼────────┼────────┤
│ 1   │ 1     │ 1     │ 1     │ 1     │ 1      │ 1      │ 1      │
│ 2   │ -1    │ 1     │ -1    │ 1     │ 1      │ -1     │ -1     │
│ 3   │ 1     │ -1    │ 1     │ 1     │ -1     │ -1     │ -1     │
│ 4   │ -1    │ 1     │ 1     │ -1    │ -1     │ -1     │ 1      │
│ 5   │ 1     │ 1     │ -1    │ -1    │ -1     │ 1      │ -1     │
│ 6   │ 1     │ -1    │ -1    │ -1    │ 1      │ -1     │ 1      │
│ 7   │ -1    │ -1    │ -1    │ 1     │ -1     │ 1      │ 1      │
│ 8   │ -1    │ -1    │ 1     │ -1    │ 1      │ 1      │ -1     │

```
"""
function PlackettBurman(formula::FormulaTerm; symbol_encoding::Bool = false)
    symbol_factors = Tuple(r.sym for r in formula.rhs)

    initial_design = plackettburman(length(symbol_factors))
    categorical_design = similar(initial_design, Symbol)

    if symbol_encoding
        for i = 1:size(initial_design, 1), j = 1:size(initial_design, 2)
            categorical_design[i, j] = initial_design[i, j] == 1.0 ? :high : :low
        end

        design = DataFrame(categorical_design)
    else
        design = DataFrame(initial_design)
    end

    dummy_factors = Tuple(Symbol("dummy" * string(i)) for i = 1:(length(names(design)) -
                                                                 length(symbol_factors)))

    design_names = (symbol_factors..., dummy_factors...)
    rename!(design, collect(design_names))

    PlackettBurman(design,
                   symbol_factors,
                   dummy_factors,
                   term(formula.lhs) ~ term(-1) + sum(term.(design_names)))
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> fold!(PlackettBurman(2))
8×3 DataFrame
│ Row │ factor1 │ factor2 │ dummy1 │
│     │ Int64   │ Int64   │ Int64  │
├─────┼─────────┼─────────┼────────┤
│ 1   │ 1       │ 1       │ 1      │
│ 2   │ 1       │ -1      │ -1     │
│ 3   │ -1      │ -1      │ 1      │
│ 4   │ -1      │ 1       │ -1     │
│ 5   │ -1      │ -1      │ -1     │
│ 6   │ -1      │ 1       │ 1      │
│ 7   │ 1       │ 1       │ -1     │
│ 8   │ 1       │ -1      │ 1      │

```
"""
function fold!(design::AbstractScreeningDesign)
    mirror_matrix = deepcopy(design.matrix)
    if mirror_matrix[1,1] isa Symbol
        for col in names(mirror_matrix)
            replace!(mirror_matrix[col], :low => :high, :high => :low)
        end
    else
        for col in names(mirror_matrix)
            replace!(mirror_matrix[col], 1 => -1, -1 => 1)
        end
    end
    append!(design.matrix, mirror_matrix)
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> PlackettBurman(4)
PlackettBurman
Dimension: (8, 7)
Factors: (:factor1, :factor2, :factor3, :factor4)
Dummy Factors: (:dummy1, :dummy2, :dummy3)
Formula: 0 ~ -1 + factor1 + factor2 + factor3 + factor4 + dummy1 + dummy2 + dummy3
Design Matrix:
8×7 DataFrame
│ Row │ factor1 │ factor2 │ factor3 │ factor4 │ dummy1 │ dummy2 │ dummy3 │
│     │ Int64   │ Int64   │ Int64   │ Int64   │ Int64  │ Int64  │ Int64  │
├─────┼─────────┼─────────┼─────────┼─────────┼────────┼────────┼────────┤
│ 1   │ 1       │ 1       │ 1       │ 1       │ 1      │ 1      │ 1      │
│ 2   │ -1      │ 1       │ -1      │ 1       │ 1      │ -1     │ -1     │
│ 3   │ 1       │ -1      │ 1       │ 1       │ -1     │ -1     │ -1     │
│ 4   │ -1      │ 1       │ 1       │ -1      │ -1     │ -1     │ 1      │
│ 5   │ 1       │ 1       │ -1      │ -1      │ -1     │ 1      │ -1     │
│ 6   │ 1       │ -1      │ -1      │ -1      │ 1      │ -1     │ 1      │
│ 7   │ -1      │ -1      │ -1      │ 1       │ -1     │ 1      │ 1      │
│ 8   │ -1      │ -1      │ 1       │ -1      │ 1      │ 1      │ -1     │

```
"""
function PlackettBurman(factors::Int)
    PlackettBurman(ConstantTerm(0) ~ sum(term.(Symbol("factor" * string(i)) for i = 1:factors)))
end

"""
$(TYPEDEF)

Encapsulates  an *explicit*  full factorial  design, where  all experiments  are
completely computed  and stored.  Use [`IterableFullFactorial`](@ref)  to obtain
an *implicit* factorial  design with arbitrary access to the  elements of a full
factorial design.

$(TYPEDFIELDS)
"""
struct FullFactorial <: AbstractFactorialDesign
    matrix::DataFrame
    factors::NamedTuple
    formula::FormulaTerm
end

"""
$(TYPEDEF)

$(TYPEDFIELDS)
"""
struct IterableFullFactorial <: AbstractFactorialDesign
    matrix::DataFrame
    factors::NamedTuple
    formula::FormulaTerm
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> FullFactorial((A = [1, 2, 4], B = [:a, :b], C = [1.0, -1.0]), @formula(y ~ A + B +C))
FullFactorial
Dimension: (12, 3)
Factors: (A = [1, 2, 4], B = [:a, :b], C = [1.0, -1.0])
Formula: y ~ A + B + C
Design Matrix:
12×3 DataFrame
│ Row │ A   │ B   │ C    │
│     │ Any │ Any │ Any  │
├─────┼─────┼─────┼──────┤
│ 1   │ 1   │ a   │ 1.0  │
│ 2   │ 2   │ a   │ 1.0  │
│ 3   │ 4   │ a   │ 1.0  │
│ 4   │ 1   │ b   │ 1.0  │
│ 5   │ 2   │ b   │ 1.0  │
│ 6   │ 4   │ b   │ 1.0  │
│ 7   │ 1   │ a   │ -1.0 │
│ 8   │ 2   │ a   │ -1.0 │
│ 9   │ 4   │ a   │ -1.0 │
│ 10  │ 1   │ b   │ -1.0 │
│ 11  │ 2   │ b   │ -1.0 │
│ 12  │ 4   │ b   │ -1.0 │

```
"""
function FullFactorial(factors::NamedTuple, formula::FormulaTerm)
    iterator = fullfactorial(values(factors))
    matrix = DataFrame(explicit_fullfactorial(iterator))
    rename!(matrix, [r.sym for r in formula.rhs])
    FullFactorial(matrix, factors, formula)
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> FullFactorial((A = [1, 2, 4], B = [:a, :b], C = [1.0, -1.0]))
FullFactorial
Dimension: (12, 3)
Factors: (A = [1, 2, 4], B = [:a, :b], C = [1.0, -1.0])
Formula: 0 ~ A + B + C
Design Matrix:
12×3 DataFrame
│ Row │ A   │ B   │ C    │
│     │ Any │ Any │ Any  │
├─────┼─────┼─────┼──────┤
│ 1   │ 1   │ a   │ 1.0  │
│ 2   │ 2   │ a   │ 1.0  │
│ 3   │ 4   │ a   │ 1.0  │
│ 4   │ 1   │ b   │ 1.0  │
│ 5   │ 2   │ b   │ 1.0  │
│ 6   │ 4   │ b   │ 1.0  │
│ 7   │ 1   │ a   │ -1.0 │
│ 8   │ 2   │ a   │ -1.0 │
│ 9   │ 4   │ a   │ -1.0 │
│ 10  │ 1   │ b   │ -1.0 │
│ 11  │ 2   │ b   │ -1.0 │
│ 12  │ 4   │ b   │ -1.0 │

```
"""
function FullFactorial(factors::NamedTuple)
    FullFactorial(factors, ConstantTerm(0) ~ sum(term.(keys(factors))))
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> FullFactorial(([1, 2, 4], [:a, :b], [1.0, -1.0]))
FullFactorial
Dimension: (12, 3)
Factors: (factor1 = [1, 2, 4], factor2 = [:a, :b], factor3 = [1.0, -1.0])
Formula: 0 ~ factor1 + factor2 + factor3
Design Matrix:
12×3 DataFrame
│ Row │ factor1 │ factor2 │ factor3 │
│     │ Any     │ Any     │ Any     │
├─────┼─────────┼─────────┼─────────┤
│ 1   │ 1       │ a       │ 1.0     │
│ 2   │ 2       │ a       │ 1.0     │
│ 3   │ 4       │ a       │ 1.0     │
│ 4   │ 1       │ b       │ 1.0     │
│ 5   │ 2       │ b       │ 1.0     │
│ 6   │ 4       │ b       │ 1.0     │
│ 7   │ 1       │ a       │ -1.0    │
│ 8   │ 2       │ a       │ -1.0    │
│ 9   │ 4       │ a       │ -1.0    │
│ 10  │ 1       │ b       │ -1.0    │
│ 11  │ 2       │ b       │ -1.0    │
│ 12  │ 4       │ b       │ -1.0    │

```
"""
function FullFactorial(factors::Tuple)
    FullFactorial(NamedTuple{Tuple(Symbol("factor" * string(i)) for i = 1:length(factors))}(factors))
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> FullFactorial(fill([-1, 1], 3))
FullFactorial
Dimension: (8, 3)
Factors: (factor1 = [-1, 1], factor2 = [-1, 1], factor3 = [-1, 1])
Formula: 0 ~ factor1 + factor2 + factor3
Design Matrix:
8×3 DataFrame
│ Row │ factor1 │ factor2 │ factor3 │
│     │ Int64   │ Int64   │ Int64   │
├─────┼─────────┼─────────┼─────────┤
│ 1   │ -1      │ -1      │ -1      │
│ 2   │ 1       │ -1      │ -1      │
│ 3   │ -1      │ 1       │ -1      │
│ 4   │ 1       │ 1       │ -1      │
│ 5   │ -1      │ -1      │ 1       │
│ 6   │ 1       │ -1      │ 1       │
│ 7   │ -1      │ 1       │ 1       │
│ 8   │ 1       │ 1       │ 1       │

```
"""
function FullFactorial(factors::Array)
    FullFactorial(tuple(factors...))
end

"""
$(TYPEDEF)

$(TYPEDFIELDS)
"""
struct FractionalFactorial <: AbstractFactorialDesign
end

"""
$(TYPEDEF)

$(TYPEDFIELDS)
"""
struct RandomDesign <: AbstractRandomDesign
    matrix::DataFrame
    factors::NamedTuple
    formula::FormulaTerm
end

"""
$(TYPEDEF)

$(TYPEDFIELDS)

Encapsulates one or more `Distribution`s, generating random designs.  Receives a
`NamedTuple`  of  factor names  associated  with  `Distribution`s from  the  the
*Distributions* package.   After instantiating a `DesignDistribution`,  you must
request   samples   from   it    using   [`rand`](@ref),   which   generates   a
[`RandomDesign`](@ref)

"""
struct DesignDistribution
    factors::NamedTuple
    formula::FormulaTerm
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> DesignDistribution((f1 = Uniform(2, 3), f2 = DiscreteUniform(-1, 5), f3 = Uniform(5, 10)))
DesignDistribution
Formula: 0 ~ f1 + f2 + f3
Factor Distributions:
f1: Distributions.Uniform{Float64}(a=2.0, b=3.0)
f2: Distributions.DiscreteUniform(a=-1, b=5)
f3: Distributions.Uniform{Float64}(a=5.0, b=10.0)

```
"""
function DesignDistribution(factors::NamedTuple)
    DesignDistribution(factors, ConstantTerm(0) ~ sum(term.(keys(factors))))
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> DesignDistribution((Uniform(2, 3), DiscreteUniform(-1, 5), Uniform(5, 10)))
DesignDistribution
Formula: 0 ~ factor1 + factor2 + factor3
Factor Distributions:
factor1: Distributions.Uniform{Float64}(a=2.0, b=3.0)
factor2: Distributions.DiscreteUniform(a=-1, b=5)
factor3: Distributions.Uniform{Float64}(a=5.0, b=10.0)

```
"""
function DesignDistribution(factors::Tuple)
    DesignDistribution(NamedTuple{Tuple(Symbol("factor" * string(i)) for i = 1:length(factors))}(factors))
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> DesignDistribution([Uniform(2, 3), DiscreteUniform(-1, 5), Uniform(5, 10)])
DesignDistribution
Formula: 0 ~ factor1 + factor2 + factor3
Factor Distributions:
factor1: Distributions.Uniform{Float64}(a=2.0, b=3.0)
factor2: Distributions.DiscreteUniform(a=-1, b=5)
factor3: Distributions.Uniform{Float64}(a=5.0, b=10.0)

```
"""
function DesignDistribution(factors::Array)
    DesignDistribution(tuple(factors...))
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> DesignDistribution(DiscreteNonParametric([-1, 1], [0.5, 0.5]), 6)
DesignDistribution
Formula: 0 ~ factor1 + factor2 + factor3 + factor4 + factor5 + factor6
Factor Distributions:
factor1: Distributions.DiscreteNonParametric{Int64,Float64,Array{Int64,1},Array{Float64,1}}(support=[-1, 1], p=[0.5, 0.5])
factor2: Distributions.DiscreteNonParametric{Int64,Float64,Array{Int64,1},Array{Float64,1}}(support=[-1, 1], p=[0.5, 0.5])
factor3: Distributions.DiscreteNonParametric{Int64,Float64,Array{Int64,1},Array{Float64,1}}(support=[-1, 1], p=[0.5, 0.5])
factor4: Distributions.DiscreteNonParametric{Int64,Float64,Array{Int64,1},Array{Float64,1}}(support=[-1, 1], p=[0.5, 0.5])
factor5: Distributions.DiscreteNonParametric{Int64,Float64,Array{Int64,1},Array{Float64,1}}(support=[-1, 1], p=[0.5, 0.5])
factor6: Distributions.DiscreteNonParametric{Int64,Float64,Array{Int64,1},Array{Float64,1}}(support=[-1, 1], p=[0.5, 0.5])

```
"""
function DesignDistribution(distribution::D, n::Int) where D <: Distribution
    factors = tuple(fill(distribution, n)...)
    DesignDistribution(factors)
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> rand(DesignDistribution((f1 = Uniform(2, 3), f2 = DiscreteUniform(-1, 5), f3 = Uniform(5, 10))), 12)
ExperimentalDesign.RandomDesign
Dimension: (12, 3)
Factors: (f1 = Distributions.Uniform{Float64}(a=2.0, b=3.0), f2 = Distributions.DiscreteUniform(a=-1, b=5), f3 = Distributions.Uniform{Float64}(a=5.0, b=10.0))
Formula: 0 ~ f1 + f2 + f3
Design Matrix:
12×3 DataFrame
│ Row │ f1      │ f2      │ f3      │
│     │ Float64 │ Float64 │ Float64 │
├─────┼─────────┼─────────┼─────────┤
│ 1   │ 2.04922 │ 5.0     │ 5.72029 │
│ 2   │ 2.59117 │ 1.0     │ 8.48192 │
│ 3   │ 2.77148 │ 4.0     │ 5.42517 │
│ 4   │ 2.25659 │ 5.0     │ 6.75815 │
│ 5   │ 2.64968 │ 5.0     │ 8.94993 │
│ 6   │ 2.31523 │ 4.0     │ 5.89413 │
│ 7   │ 2.86526 │ 3.0     │ 9.8807  │
│ 8   │ 2.44753 │ -1.0    │ 9.24986 │
│ 9   │ 2.08284 │ 3.0     │ 9.99911 │
│ 10  │ 2.81201 │ 2.0     │ 7.28848 │
│ 11  │ 2.64232 │ 4.0     │ 5.34507 │
│ 12  │ 2.08189 │ 4.0     │ 8.303   │

```
"""
function rand(distribution::DesignDistribution, n::Int = 1)
    RandomDesign(DataFrame(random_design!(values(distribution.factors), n),
                           collect(keys(distribution.factors))),
                 distribution.factors,
                 distribution.formula)
end

"""
$(TYPEDEF)

$(TYPEDFIELDS)

Contains a  set of candidate experiments,  a target optimality criterion,  and a
model prior  formula.  A set  of experiments  that maximizes a  given optimality
criterion can selected from a candidate set using the [`kl_exchange`](@ref)
algorithm.

"""
struct OptimalDesign <: AbstractOptimalDesign
    matrix::DataFrame
    factors::NamedTuple
    formula::FormulaTerm
    selected_experiments::Array{Int}
    criteria::Dict{Symbol, Float64}
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> design_distribution = DesignDistribution((f1 = Uniform(2, 3), f2 = DiscreteUniform(-1, 5), f3 = Uniform(5, 10)))
DesignDistribution
Formula: 0 ~ f1 + f2 + f3
Factor Distributions:
f1: Distributions.Uniform{Float64}(a=2.0, b=3.0)
f2: Distributions.DiscreteUniform(a=-1, b=5)
f3: Distributions.Uniform{Float64}(a=5.0, b=10.0)

julia> design = rand(design_distribution, 400);

julia> f = @formula 0 ~ f1 + f2 + f3 + f2 ^ 2;

julia> OptimalDesign(design, f, 10)
OptimalDesign
Dimension: (10, 3)
Factors: (f1 = Distributions.Uniform{Float64}(a=2.0, b=3.0), f2 = Distributions.DiscreteUniform(a=-1, b=5), f3 = Distributions.Uniform{Float64}(a=5.0, b=10.0))
Formula: 0 ~ f1 + f2 + f3 + :(f2 ^ 2)
Selected Candidate Rows: [222, 103, 384, 140, 139, 156, 63, 184, 169, 54]
Optimality Criteria: Dict(:D => 2.418045723405036)
Design Matrix:
10×3 DataFrame
│ Row │ f1      │ f2      │ f3      │
│     │ Float64 │ Float64 │ Float64 │
├─────┼─────────┼─────────┼─────────┤
│ 1   │ 2.02407 │ -1.0    │ 5.07959 │
│ 2   │ 2.0203  │ 5.0     │ 9.50193 │
│ 3   │ 2.80676 │ -1.0    │ 9.87726 │
│ 4   │ 2.9996  │ 2.0     │ 5.03174 │
│ 5   │ 2.14846 │ 2.0     │ 9.72596 │
│ 6   │ 2.84392 │ 5.0     │ 5.01213 │
│ 7   │ 2.04046 │ 2.0     │ 5.00331 │
│ 8   │ 2.93432 │ 1.0     │ 9.71405 │
│ 9   │ 2.79858 │ -1.0    │ 5.35685 │
│ 10  │ 2.97635 │ 5.0     │ 8.10695 │

```
"""
function OptimalDesign(candidates::AbstractDesign,
                       formula::FormulaTerm,
                       experiments::Int;
                       tolerance::Float64 = 1e-9,
                       seed_design_size::Int = 2,
                       max_iterations::Int = 1000,
                       design_k::Int = experiments,
                       candidates_l::Int = size(candidates.matrix, 1) - design_k)
    selected = kl_exchange(formula,
                           candidates.matrix,
                           experiments = experiments,
                           tolerance = tolerance,
                           seed_design_size = seed_design_size,
                           max_iterations = max_iterations,
                           design_k = design_k,
                           candidates_l = candidates_l)

    optimality = d_criterion(selected)

    OptimalDesign(candidates.matrix[selected.indices[1], :],
                  candidates.factors,
                  formula,
                  selected.indices[1],
                  Dict(:D => optimality))
end


text:
using RiemannTheta
using Test, Random
using LinearAlgebra

######### some testing values

Ω1 = [ 1.690983006 + 0.9510565162im 1.5+0.363271264im ;
      1.5+0.363271264im 1.309016994+0.9510565162im ]
T1 = Matrix(cholesky(imag.(Ω1)).U)

Ω2 = -1/(2π * im) * [ 111.207 96.616 ; 96.616 83.943 ]
T2 = Matrix(cholesky(imag.(Ω2)).U)

Random.seed!(0)
tmp = 5*rand(10,10) .- 2.5
Ω3 = 5*rand(10,10) .- 2.5 + (tmp * tmp') * im
T3 = Matrix(cholesky(imag.(Ω3)).U)

derivs1 = [ rand(ComplexF64, 2) for i in 1:1 ]
derivs2 = [ rand(ComplexF64, 10) for i in 1:5 ]

ϵ1, ϵ2 = 1e-3, 1e-8

##########

include("internals.jl")
include("mainfuncs.jl")


text:
"""
    Apply a rotation CNN
"""
function transform_4corners_to_4cols(bb)
    top_left = bb[1:2, 1:2]
    top_right = [
            bb[1, 4] bb[2,4];
            bb[1, 3] bb[2,3]
        ]
    bottom_right =
        [
            bb[4, 4] bb[4,3];
            bb[3, 4] bb[3,3]
        ]
    bottom_left =
        [
            bb[4, 1] bb[3,1];
            bb[4, 2] bb[3,2]
        ]

    reduce(hcat, reshape.([top_left, top_right, bottom_left, bottom_right], :))
end


text:
module Stippler
using FileIO, ColorTypes, GZip

function extract_image_from_file(fileName::String)
	if !isfile(fileName)
		println("""The file \""""*fileName*"""\" is missing.""")
		exit()
	else
		image = RGB.(load(fileName))
		image = (Float64.(red.(image)) + Float64.(green.(image)) + Float64.(blue.(image))) ./ 3
		image = round.(UInt8, 255 .* image)
		return image::Array{UInt8, 2}
	end
end

function minimumish(image::Array{Float64, 2}, startInd::Tuple{Int, Int})

	function find(image::Array{Float64, 2}, ind::Tuple{Int, Int}, margin::Int)

		indMin = argmin(image[max(1, ind[1]-margin):min(end, ind[1]+margin), max(1, ind[2]-margin):min(end, ind[2]+margin)])
		indMin = (indMin[1]+max(1, ind[1]-margin)-1, indMin[2]+max(1, ind[2]-margin)-1)

		return indMin::Tuple{Int, Int}

	end

	indMinish = find(image, startInd, 2*HALO_MARGIN)
	if indMinish == startInd
		indMinish = Tuple(argmin(image))
		if indMinish == startInd
			global DONE = true
		end
	end

	return indMinish::Tuple{Int, Int}

end

function make_dot(image::Array{Float64, 2}, pImage::Array{Float64, 2}, ind::Tuple{Int, Int}, distance::Array{Float64, 2}, io::IOStream)

	if rand() < 0.5
		haloRadius = DOT_RADIUS + (HALO_MARGIN/2)*image[ind[1], ind[2]]
		dotRadius = DOT_RADIUS
	else
		haloRadius = 0.8*(DOT_RADIUS + (HALO_MARGIN/2)*image[ind[1], ind[2]])
		dotRadius = 0.8*DOT_RADIUS
	end

	halo = 4 .* 0.01 .* ((haloRadius ./ distance).^12 .- (haloRadius ./ distance).^6) #<-- slowest
	halo[HALO_MARGIN+1, HALO_MARGIN+1] = Inf

	pImage[ind[1]-HALO_MARGIN:ind[1]+HALO_MARGIN, ind[2]-HALO_MARGIN:ind[2]+HALO_MARGIN] .+= halo
	println(io, """<circle cx=\"""", ind[2]-HALO_MARGIN, """\" cy=\"""", ind[1]-HALO_MARGIN, """\" r=\"""", dotRadius, """\" />""")

	return pImage::Array{Float64, 2}

end

function stipple(fileName::String)

	data = Float64.(extract_image_from_file(fileName)) ./ 255
	image = Inf .* ones(size(data) .+ (2*HALO_MARGIN, 2*HALO_MARGIN))
	image[HALO_MARGIN+1:end-HALO_MARGIN, HALO_MARGIN+1:end-HALO_MARGIN] = data
	pImage = deepcopy(image)

	c1 = repeat(1:2*HALO_MARGIN+1, 1, 2*HALO_MARGIN+1)
	c2 = c1'
	distance = sqrt.((c1.-(HALO_MARGIN+1)).^2 + (c2.-(HALO_MARGIN+1)).^2)

	io = open("output/"*split(split(fileName, "/")[end], ".")[1]*"_stippled.svg", "w")

	println(io, """<?xml version="1.0" encoding="UTF-8" standalone="no"?>""")
	println(io, """<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">""")
	println(io, """<svg width=\"""", size(data)[2], """px" height=\"""", size(data)[1], """px" xmlns="http://www.w3.org/2000/svg">""")
	println(io, """<g fill="black">""")

	println("\nStippling \"", fileName, "\"...")

	indMinish = Tuple(argmin(pImage))
	while !DONE
		if pImage[indMinish[1], indMinish[2]] < 0.95
			pImage = make_dot(image, pImage, indMinish, distance, io)
		end
		indMinish = minimumish(pImage, indMinish)
	end

	println(io, """</g>""")
	println(io, """</svg>""")

	close(io)

	cio = GZip.open("output/"*split(split(fileName, "/")[end], ".")[1]*"_stippled.svgz", "w")
	write(cio, read("output/"*split(split(fileName, "/")[end], ".")[1]*"_stippled.svg", String))
	close(cio)

	println("...done.")

	return nothing

end

global DOT_RADIUS = 2.5
global HALO_MARGIN = 20

global DONE = false

cd("..")
println("\n------------------------------------------------")
println("\nCurrent path is \"", pwd(), "\".")
println("\nEnter input file path (without quotation):\n")
stipple(readline())
println("\n------------------------------------------------")

end # Stippler

max_stars_repo_name: JonPallbo/Stippler

max_issues_repo_name: JonPallbo/Stippler

max_forks_repo_name: JonPallbo/Stippler


text:
#--------------------------------------------------------------------
# DNSS.jl
# Soham 03-2022
# Spectral utitlies for changing the number of basis
# functions in each element.
#--------------------------------------------------------------------

function cheb(m::Int, x::Number)::Number
    if abs(x) <= 1
        return cos(m*acos(x))
    elseif x >= 1
        return cosh(m*acosh(x))
    else
        return ((-1)^m)*cosh(m*acosh(-x))
    end
end

function scaling(i::Int, N::Int)::Number
    (i == 1 || i == N+1) ? (return 1/2) : (return 1)
end

function linetransform(PS::Chebyshev{Tag, N, T}, modal::Array{T, 1}, nodal::Array{T,1})::Array{T, 1} where {Tag, N, T}
    PSC   = ChebyshevGL{Tag,N,T}(range(PS)...)
    for index in CartesianIndices(nodal)
        for ordr in 0:order(PSC)
            nodal[index] += scaling(ordr+1, order(PSC))*cheb(ordr, naturalcollocation(PSC, index.I[1]))*modal[ordr+1]
        end
    end
    return nodal
end

function linetransform(PS::ChebyshevGL{Tag, N, T}, nodal::Array{T, 1}, modal::Array{T,1})::Array{T, 1} where {Tag, N, T}
    for ordr in 0:order(PS)
        for index in CartesianIndices(modal)
            modal[ordr+1] += scaling(index.I[1], order(PS))*cheb(ordr, naturalcollocation(PS, index.I[1]))*nodal[index]
        end
    end
    return (2/order(PS))*modal
end

function embed!(u::Vector{T}, v::Vector{T})::Vector{T} where {T}
    a = length(u)
    b = length(v)
    for i in 1:min(a,b)
        v[i] = u[i]
    end
    return v
end

function embed!(u::Matrix{T}, v::Matrix{T})::Matrix{T} where {T}
    (a, b) = size(u)
    (c, d) = size(v)
    for i in 1:min(a,c), j in 1:min(b,d)
        v[i,j] = u[i,j]
    end
    return v
end

function basistransform(u::Field{ChebyshevGL{Tag, N, T}})::Field{Chebyshev{Tag, N, T}} where {Tag, N, T}
    nodal = u.value
    modal = zeros(T, size(nodal))
    return Field(Chebyshev{Tag, N, T}(range(u.space)...), linetransform(u.space, nodal, modal))
end

function basistransform(u::Field{Chebyshev{Tag, N, T}})::Field{ChebyshevGL{Tag, N, T}} where {Tag, N, T}
    modal = u.value
    nodal = zeros(T, size(modal))
    return Field(ChebyshevGL{Tag, N, T}(range(u.space)...), linetransform(u.space, modal, nodal))
end

function project(u::Field{Chebyshev{Tag, N1, T}}, S::Chebyshev{Tag, N2, T})::Field{Chebyshev{Tag, N2, T}} where {Tag, N1, N2, T}
    v = zeros(T, N2)
    return Field(Chebyshev{Tag, N2, T}(range(u.space)...), embed!(u.value, v))
end

function project(u::Field{ChebyshevGL{Tag, N1, T}}, S::ChebyshevGL{Tag, N2, T})::Field{ChebyshevGL{Tag, N2, T}} where {Tag, N1, N2, T}
    return basistransform(project(basistransform(u), basistransform(S)))
end

function Base. filter!(u::Field{ChebyshevGL{Tag, N, T}})::Field{ChebyshevGL{Tag, N, T}} where {Tag, N, T}
    SR = Chebyshev{Tag, Int(ceil(N * 2/3)), T}(range(u.space)...)
    S  = Chebyshev{Tag, N, T}(range(u.space)...)
    return basistransform(project(project(basistransform(u), SR), S))
end

function basistransform(u::Field{S}) where {S<:ProductSpace{ChebyshevGL{Tag1, N1, T},
                                                            ChebyshevGL{Tag2, N2, T}}} where {Tag1, Tag2, N1, N2, T}
    nodal   = u.value
    modalS1 = zeros(T, size(nodal))
    modalS2 = zeros(T, size(nodal))

    for index in 1:size(u.space.S2)
        modalS1[:, index] = linetransform(u.space.S1, nodal[:, index], modalS1[:, index])
    end

    for index in 1:size(u.space.S1)
        modalS2[index, :] = linetransform(u.space.S2, modalS1[index, :], modalS2[index, :])
    end

    return Field(ProductSpace(Chebyshev{Tag1, N1, T}(u.space.S1.min, u.space.S1.max),
                              Chebyshev{Tag2, N2, T}(u.space.S2.min, u.space.S2.max)), modalS2)
end

function basistransform(u::Field{S}) where {S<:ProductSpace{Chebyshev{Tag1, N1, T},
                                                            Chebyshev{Tag2, N2, T}}} where {Tag1, Tag2, N1, N2, T}
    modal   = u.value
    nodalS1 = zeros(T, size(modal))
    nodalS2 = zeros(T, size(modal))

    for index in 1:size(u.space.S1)
        nodalS1[index, :] = linetransform(u.space.S2, modal[index, :], nodalS1[index, :])
    end

    for index in 1:size(u.space.S2)
        nodalS2[:, index] = linetransform(u.space.S1, nodalS1[:, index], nodalS2[:, index])
    end

    return Field(ProductSpace(ChebyshevGL{Tag1, N1, T}(u.space.S1.min, u.space.S1.max),
                              ChebyshevGL{Tag2, N2, T}(u.space.S2.min, u.space.S2.max)), nodalS2)
end

function project(u::Field{ProductSpace{Chebyshev{Tag1, N1, T},
                                       Chebyshev{Tag2, N2, T}}}, PS::ProductSpace{Chebyshev{Tag1, N3, T},
                                                                                Chebyshev{Tag2, N4, T}})::Field{ProductSpace{Chebyshev{Tag1, N3, T},
                                                                                Chebyshev{Tag2, N4, T}}} where {Tag1, Tag2, N1, N2, N3, N4, T}
    @assert range(u.space) == range(PS)
    v = zeros(T, N3, N4)
    return Field(ProductSpace(Chebyshev{Tag1, N3, T}(range(u.space.S1)...),
                              Chebyshev{Tag2, N4, T}(range(u.space.S2)...)), embed!(u.value, v))
end

function project(u::Field{ProductSpace{ChebyshevGL{Tag1, N1, T},
                                       ChebyshevGL{Tag2, N2, T}}}, PS::ProductSpace{ChebyshevGL{Tag1, N3, T},
                                                                                ChebyshevGL{Tag2, N4, T}})::Field{ProductSpace{ChebyshevGL{Tag1, N3, T},
                                                                                ChebyshevGL{Tag2, N4, T}}} where {Tag1, Tag2, N1, N2, N3, N4, T}
    return basistransform(project(basistransform(u), basistransform(PS)))
end

function basistransform(PS::ProductSpace{ChebyshevGL{Tag1, N1, T}, ChebyshevGL{Tag2, N2, T}}) where {Tag1, Tag2, N1, N2, T}
    A = Chebyshev{Tag1, N1, T}(range(PS.S1)...)
    B = Chebyshev{Tag2, N2, T}(range(PS.S2)...)
    return ProductSpace(A,B)
end

function basistransform(PS::ProductSpace{Chebyshev{Tag1, N1, T}, Chebyshev{Tag2, N2, T}}) where {Tag1, Tag2, N1, N2, T}
    A = ChebyshevGL{Tag1, N1, T}(range(PS.S1)...)
    B = ChebyshevGL{Tag2, N2, T}(range(PS.S2)...)
    return ProductSpace(A,B)
end

function basistransform(S::ChebyshevGL{Tag, N, T}) where {Tag, N, T}
    return  Chebyshev{Tag, N, T}(range(S)...)
end

function basistransform(S::Chebyshev{Tag, N, T}) where {Tag, N, T}
    return  ChebyshevGL{Tag, N, T}(range(S)...)
end

function prolongate(S::ChebyshevGL{Tag, N, T}) where {Tag, N, T}
    return ChebyshevGL{Tag, 2*N, T}(range(S)...)
end

function restrict(S::ChebyshevGL{Tag, N, T}) where {Tag, N, T}
    N2 = iseven(N) ? div(N, 2) : div(N + 1, 2)
    return ChebyshevGL{Tag, N2, T}(range(S)...)
end

function prolongate(PS::ProductSpace{ChebyshevGL{Tag1, N1, T}, ChebyshevGL{Tag2, N2, T}}) where {Tag1, Tag2, N1, N2, T}
    A = ChebyshevGL{Tag1, 2*N1, T}(range(PS.S1)...)
    B = ChebyshevGL{Tag2, 2*N2, T}(range(PS.S2)...)
    return ProductSpace(A,B)
end

function restrict(PS::ProductSpace{ChebyshevGL{Tag1, N1, T}, ChebyshevGL{Tag2, N2, T}}) where {Tag1, Tag2, N1, N2, T}
    N3 = iseven(N1) ? div(N1, 2) : div(N1 + 1, 2)
    N4 = iseven(N1) ? div(N2, 2) : div(N2 + 1, 2)
    A = ChebyshevGL{Tag1, N3, T}(range(PS.S1)...)
    B = ChebyshevGL{Tag2, N4, T}(range(PS.S2)...)
    return ProductSpace(A,B)
end

function restrict(u::Field{S}) where {S}
    return project(u, restrict(u.space))
end

function prolongate(u::Field{S}) where {S}
    return project(u, prolongate(u.space))
end

function restrict(U::NTuple{N,Field{S}}) where {N,S}
    return map(restrict, U)
end

function prolongate(U::NTuple{N,Field{S}}) where {N,S}
    return map(prolongate, U)
end

function filtertophalf(u::Field{S}) where {S}
    return (prolongate ∘ restrict)(u)
end

function filtertophalf(U::NTuple{N,Field{S}}) where {N,S}
    return map(filtertophalf, U)
end


text:
type LongstaffSchwartzPathPricer{E <: EarlyExercisePathPricer, T} <: AbstractPathPricer
  calibrationPhase::Bool
  pathPricer::E
  coeff::Vector{Vector{Float64}}
  dF::Vector{Float64}
  paths::Vector{Path}
  len::Int
  exerciseProbability::NonWeightedStatistics
  v::Vector{T}
end

function LongstaffSchwartzPathPricer(tg::TimeGrid, ep::EarlyExercisePathPricer, yts::YieldTermStructure)
  v = basis_system(ep)
  len = length(tg.times)
  coeff = Vector{Vector{Float64}}(len - 2)
  dF = zeros(len - 1)
  paths = Vector{Path}()

  for i in eachindex(dF)
    dF[i] = discount(yts, tg.times[i+1]) / discount(yts, tg.times[i])
  end

  return LongstaffSchwartzPathPricer(true, ep, coeff, dF, paths, len, NonWeightedStatistics(), v)
end

type MCAmericanEngine{S <: AbstractBlackScholesProcess, P <: LsmBasisSystemPolynomType} <: MCLongstaffSchwartzEngine
  process::S
  timeSteps::Int
  timeStepsPerYear::Int
  requiredSamples::Int
  maxSamples::Int
  requiredTolerance::Float64
  brownianBridge::Bool
  seed::Int
  nCalibrationSamples::Int
  polynomOrder::Int
  polynomType::P
  mcSimulation::MCSimulation
  pathPricer::LongstaffSchwartzPathPricer

  MCAmericanEngine{S, P}(process::S, timeSteps::Int, timeStepsPerYear::Int, requiredSamples::Int, maxSamples::Int, requiredTolerance::Float64, brownianBridge::Bool,
                            seed::Int, nCalibrationSamples::Int, polynomOrder::Int, polynomType::P, mcSimulation::MCSimulation) =
                            new{S, P}(process, timeSteps, timeStepsPerYear, requiredSamples, maxSamples, requiredTolerance, brownianBridge, seed,
                                        nCalibrationSamples, polynomOrder, polynomType, mcSimulation)
end

function MCAmericanEngine{RSG <: AbstractRandomSequenceGenerator, S <: AbstractBlackScholesProcess}(process::S; timeSteps::Int = -1, timeStepsPerYear::Int = -1, brownianBridge::Bool = false,
                          antitheticVariate::Bool = false, requiredSamples::Int = -1, requiredTolerance::Float64 = -1.0, maxSamples::Int = typemax(Int), seed::Int = 0, rsg::RSG = InverseCumulativeRSG(seed),
                          nCalibrationSamples::Int = 2048, polynomOrder::Int = 2, polynomType::LsmBasisSystemPolynomType = Monomial())
  # build mc sim
  mcSim = MCSimulation{RSG, SingleVariate}(antitheticVariate, false, rsg, SingleVariate())

  return MCAmericanEngine{S, typeof(polynomType)}(process, timeSteps, timeStepsPerYear, requiredSamples, maxSamples, requiredTolerance, brownianBridge, seed,
                                                      nCalibrationSamples, polynomOrder, polynomType, mcSim)
end

type AmericanPathPricer{P <: StrikedTypePayoff, T, T1} <: EarlyExercisePathPricer
  scalingValue::Float64
  payoff::P
  v::Vector{Union{T, T1}}
end

function get_payoff(app::AmericanPathPricer)
  function _get_payoff(x::Float64)
    return app.payoff(x / app.scalingValue)
  end

  return _get_payoff
end

function AmericanPathPricer{P <: StrikedTypePayoff, L <: LsmBasisSystemPolynomType}(payoff::P, polynomOrder::Int, polynomType::L)
  T = get_type(polynomType)
  v = Vector{Union{T, P}}(polynomOrder + 1)
  path_basis_system!(polynomType, polynomOrder, v)

  scalingVal = 1.0 / payoff.strike

  app = AmericanPathPricer{P, T, Function}(scalingVal, payoff, v)

  push!(app.v, get_payoff(app))
  return app
end

basis_system(p::AmericanPathPricer) = p.v
get_state(p::AmericanPathPricer, path::Path, t::Int) = path[t] * p.scalingValue
get_payoff(p::AmericanPathPricer, st::Float64) = p.payoff(st / p.scalingValue)
call(p::AmericanPathPricer, path::Path, t::Int) = get_payoff(p, get_state(p, path, t))

function update_paths!(lpp::LongstaffSchwartzPathPricer, path::Path)
  push!(lpp.paths, copy(path))
  return lpp
end

function call(lpp::LongstaffSchwartzPathPricer, path::Path)
  if lpp.calibrationPhase
    update_paths!(lpp, path)
    return 0.0
  end
  price = lpp.pathPricer(path, lpp.len)
  # initialize with exercise on last date
  exercised = price != 0.0

  for i = lpp.len-1:-1:2
    price *= lpp.dF[i]
    exercise = lpp.pathPricer(path, i)
    if exercise > 0.0
      regValue = get_state(lpp.pathPricer, path, i)

      continuationValue = 0.0
      for l in eachindex(lpp.v)
        continuationValue += lpp.coeff[i-1][l] * lpp.v[l](regValue)
      end

      if continuationValue < exercise
        price = exercise
        exercised = true
      end
    end
  end
  add_sample!(lpp.exerciseProbability, exercised ? 1.0 : 0.0)

  return price * lpp.dF[1]
end

function calibrate!(lpp::LongstaffSchwartzPathPricer)
  n = length(lpp.paths)
  prices = zeros(n)
  exercise = zeros(n)
  len = length(lpp.paths[1])

  for i in eachindex(prices)
    prices[i] = lpp.pathPricer(lpp.paths[i], len)
  end

  for i = len-1:-1:2
    y = Float64[]
    x = Float64[]

    # rollback step
    for j in eachindex(exercise)
      exercise[j] = lpp.pathPricer(lpp.paths[j], i)

      if exercise[j] > 0.0
        push!(x, get_state(lpp.pathPricer, lpp.paths[j], i))
        push!(y, lpp.dF[i] * prices[j])
      end
    end

    if length(lpp.v) <= length(x)
      lpp.coeff[i-1] = get_coefficients(GeneralLinearLeastSquares(x, y, lpp.v))
    else
      # if number of itm paths is smaller than the number of calibration functions
      # then early exercise if earlyExercise > 0
      lpp.coeff[i-1] = zeros(length(lpp.v))
    end

    k = 1
    for j in eachindex(prices)
      prices[j] *= lpp.dF[i]
      if exercise[j] > 0.0
        continuationValue = 0.0
        for l in eachindex(lpp.v)
          continuationValue += lpp.coeff[i-1][l] * lpp.v[l](x[k])
        end
        if continuationValue < exercise[j]
          prices[j] = exercise[j]
        end
        k += 1
      end
    end
  end

  # remove calibration paths & release memory
  lpp.paths = Path[]
  # entering calculation phase
  lpp.calibrationPhase = false

  return lpp
end

function time_grid(pe::MCLongstaffSchwartzEngine, opt::VanillaOption)
  lastExerciseDate = opt.exercise.dates[end]
  t = get_time(pe.process, lastExerciseDate)
  if pe.timeSteps != -1
    return TimeGrid(t, pe.timeSteps)
  elseif pe.timeStepsPerYear != -1
    steps = round(Int, floor(pe.timeStepsPerYear * t))
    return TimeGrid(t, max(steps, 1))
  else
    error("time steps not specified")
  end
end

function lsm_path_pricer(pe::MCAmericanEngine, opt::AmericanOption)
  process = pe.process
  exercise = opt.exercise

  early_exercise_path_pricer = AmericanPathPricer(opt.payoff, pe.polynomOrder, pe.polynomType)

  return LongstaffSchwartzPathPricer(time_grid(pe, opt), early_exercise_path_pricer, process.riskFreeRate)
end

function path_generator(pe::MCLongstaffSchwartzEngine, opt::VanillaOption)
  dimensions = get_factors(pe.process)
  grid = time_grid(pe, opt)
  init_sequence_generator!(pe.mcSimulation.rsg, dimensions * (length(grid.times) - 1))

  return PathGenerator(pe.process, grid, pe.mcSimulation.rsg, pe.brownianBridge)
end

path_pricer(pe::MCLongstaffSchwartzEngine, ::VanillaOption) = pe.pathPricer

function _calculate!(pe::MCLongstaffSchwartzEngine, opt::VanillaOption)
  pe.pathPricer = lsm_path_pricer(pe, opt)
  pe.mcSimulation.mcModel = MonteCarloModel(path_generator(pe, opt), pe.pathPricer, gen_RiskStatistics(), pe.mcSimulation.antitheticVariate)
  add_samples!(pe.mcSimulation.mcModel, pe.nCalibrationSamples, 1)
  # calibration
  calibrate!(pe.pathPricer)
  _calculate!(pe.mcSimulation, pe, opt, pe.requiredTolerance, pe.requiredSamples, pe.maxSamples)

  opt.results.value = stats_mean(pe.mcSimulation.mcModel.sampleAccumulator)
  return pe, opt
end


text:
################## Binary Hamiltonian Monte Carlo  ##################

#################### Types ####################

mutable struct BHMCTune <: SamplerTune
  logf::Union{Function, Missing}
  traveltime::Float64
  position::Vector{Float64}
  velocity::Vector{Float64}
  wallhits::Int
  wallcrosses::Int

  BHMCTune() = new()
"""
    BHMCTune(x::Vector, traveltime::Real, logf::Union{Function, Missing})
Construct a `BHMCVariate` object that stores simulated values and tuning
 parameters for BHMC sampling.
Returns a BHMCVariate type object with fields set to the supplied
 x and tuning parameter values.
"""
  function BHMCTune(x::Vector, traveltime::Real, logf::Union{Function, Missing})
    n = length(x)
    new(logf, traveltime, randn(n), randn(n), 0, 0)
  end
end

BHMCTune(x::Vector, traveltime::Real) =
  BHMCTune(x, traveltime, missing)

const BHMCVariate = SamplerVariate{BHMCTune}

validate(v::BHMCVariate) = validatebinary(v)

#################### Sampler Constructor ####################
"""
    BHMC(params::ElementOrVector{Symbol}, traveltime::Real)
Construct a `Sampler` object for BHMC sampling. Parameters are assumed to
  have binary numerical values (0 or 1).

Returns a `Sampler{BHMCTune}` type object.

* `params`: stochastic node(s) to be updated with the sampler.

* `traveltime`: length of time over which particle paths are simulated. It is recommended that supplied values be of the form ``(n + \\frac{1}{2}) \\pi``, where optimal choices of n ``\\in \\mathbb{Z}^+`` are expected to grow with the parameter space dimensionality.

"""
function BHMC(params::ElementOrVector{Symbol}, traveltime::Real)
  samplerfx = function(model::Model, block::Integer)
    block = SamplingBlock(model, block)
    v = SamplerVariate(block, traveltime)
    sample!(v, x -> logpdf!(block, x))
    relist(block, v)
  end
  Sampler(params, samplerfx, BHMCTune())
end

#################### Sampling Functions ####################

sample!(v::BHMCVariate) = sample!(v, v.tune.logf)
"""
    sample!(v::BHMCVariate, logf::Function)
Draw one sample from a target distribution using the BHMC sampler.
Parameters are assumed to have binary numerical values (0 or 1).
Returns `v` updated with simulated values and associated tuning parameters.
"""
function sample!(v::BHMCVariate, logf::Function)
  tune = v.tune
  flag = false
  nearzero = 1e4 * eps()
  j = 0
  totaltime = 0.0                     ## time the particle already moved

  n = length(v)                       ## length of binary vector
  S = sign.(tune.position)

  while true
    a = tune.velocity[:]
    b = tune.position[:]
    phi = atan.(b, a)

    ## time to hit or cross wall
    walltime = -phi
    idx = findall(x-> x > 0.0, phi)
    walltime[idx] .= pi .- phi[idx]

    ## if there was a previous reflection (j > 0) and there is a potential
    ## reflection at the sample plane make sure that a new reflection at j
    ## is not found because of numerical error
    if j > 0
      if abs(walltime[j]) < nearzero || abs(walltime[j] - 2.0 * pi) < nearzero
        walltime[j] = Inf
      end
    end

    ## time till particle j hits or crosses wall
    movetime, j = findmin(walltime)
    if movetime == 0.0
      error("walking length zero")
    elseif movetime == Inf
      movetime = pi
    end

    totaltime += movetime
    if totaltime >= tune.traveltime
      movetime -= totaltime - tune.traveltime
      flag = true
    else
      tune.wallhits += 1
    end

    ## move the particle a time mt
    tune.velocity[:] = a * cos(movetime) - b * sin(movetime)
    tune.position[:] = a * sin(movetime) + b * cos(movetime)

    flag && break

    tune.position[j] = 0

    S1 = (S + ones(n)) / 2.0
    S1[j] = 0.0
    S2 = (S + ones(n)) / 2.0
    S2[j] = 1.0

    v2_new = tune.velocity[j]^2 +
             sign(tune.velocity[j]) * 2.0 * (logf(S2) - logf(S1))
    if v2_new > 0.0
      tune.velocity[j] = sqrt(v2_new) * sign(tune.velocity[j])
      S[j] *= -1.0
      tune.wallcrosses += 1
    else
      tune.velocity[j] *= -1.0
    end
  end

  ## convert from (-/+1) to (0/1)
  v[:] = (sign.(tune.position) + ones(n)) / 2.0
  v
end


text:
using DashBootstrapComponents

first_card = dbc_card(
    dbc_cardbody([
        html_h5("Card title", className = "card-title"),
        html_p("This card has some text content, but not much else"),
        dbc_button("Go somewhere", color = "primary"),
    ]),
);

second_card = dbc_card(
    dbc_cardbody([
        html_h5("Card title", className = "card-title"),
        html_p(
            "This card also has some text content and not much else, but " *
            "it is twice as wide as the first card.",
        ),
        dbc_button("Go somewhere", color = "primary"),
    ]),
);

cards = dbc_row([dbc_col(first_card, width = 4), dbc_col(second_card, width = 8)]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
using Sockets

bootstrap_script = raw"""
#!/bin/bash
if [ -z "${AWS_LAMBDA_RUNTIME_API}" ]; then
  LOCAL="127.0.0.1:9001"
  echo "AWS_LAMBDA_RUNTIME_API not found, starting AWS RIE on $LOCAL"
  exec ./aws-lambda-rie /usr/local/julia/bin/julia -e "using Jot; using $PKG_NAME; start_runtime(\\\"$LOCAL\\\", $FUNC_FULL_NAME, $FUNC_PARAM_TYPE)"
else
  echo "AWS_LAMBDA_RUNTIME_API = $AWS_LAMBDA_RUNTIME_API, running Julia"
  # touch /tmp/working
  # echo "CREATED TEMPORARY FILE"
  # exec /usr/local/julia/bin/julia -e "using Jot; println(pwd())"
  # echo "JULIA CREATED TEMP"
  exec /usr/local/julia/bin/julia -e "using Jot; using $PKG_NAME; start_runtime(\\\"$AWS_LAMBDA_RUNTIME_API\\\", $FUNC_FULL_NAME, $FUNC_PARAM_TYPE)"
fi
"""

function start_lambda_server(host::String, port::Int64)
  ROUTER = HTTP.Router()
  server = Sockets.listen(Sockets.InetAddr(parse(IPAddr, host), port))

  function get_respond(req::HTTP.Request)
    println("received get request")
    return HTTP.Response(
                         200,
                         ["Lambda-Runtime-Aws-Request-Id" => "dummy-request-id"];
                         body=JSON3.write(3),
                        )
  end

  function post_respond(req::HTTP.Request)
    answer = JSON3.read(IOBuffer(HTTP.payload(req)), Int64)
    println(answer)
    println("Closing server")
    close(server)
    return HTTP.Response(200, JSON3.write(answer))
  end

  HTTP.@register(ROUTER, "POST", "/*", post_respond)
  HTTP.@register(ROUTER, "GET", "/*", get_respond)
  @sync HTTP.serve(ROUTER, host, port, server=server)
end

function get_lambda_dummy_server_jl()::String
  """
  using HTTP
  using JSON3

  function get_respond(req::HTTP.Request)
    println("received get request")
    return HTTP.Response(
                         200,
                         ["Lambda-Runtime-Aws-Request-Id" => "dummy-request-id"];
                         body=JSON3.write(3),
                        )
  end

  function post_respond(req::HTTP.Request)
    answer = JSON3.read(IOBuffer(HTTP.payload(req)), Int64)
    println(answer)
    return HTTP.Response(200, JSON3.write(answer))
    exit(86)
  end

  const ROUTER = HTTP.Router()

  HTTP.@register(ROUTER, "POST", "/*", post_respond)
  HTTP.@register(ROUTER, "GET", "/*", get_respond)
  HTTP.serve(ROUTER, "127.0.0.1", 9001)
  """
end

function get_init_script(
    package_compile::Bool,
    cpu_target::String,
  )::String
  precomp = """
  import Pkg
  using Jot
  Pkg.precompile()
  """
  pc_script = """
  Pkg.add("PackageCompiler")
  using PackageCompiler
  @async Jot.start_lambda_server("127.0.0.1", 9001)
  create_sysimage(
                  :Jot,
                  precompile_execution_file="precompile.jl",
                  replace_default=true,
                  cpu_target="$cpu_target",
                 )
  """
  package_compile ? precomp * pc_script : precomp
end

function get_precompile_jl(
    package_name::String,
  )::String
  """
  using Jot
  using $package_name
  using Pkg

  Pkg.test("$package_name")

  rf(i::Int64) = i + 1

  Jot.start_runtime("127.0.0.1:9001", rf, Int64; single_shot=true)
  """
end

function get_lambda_function_tags_script(lambda_function::LambdaFunction)::String
  """
  aws lambda list-tags --resource $(lambda_function.FunctionArn)
  """
end

function get_images_in_ecr_repo_script(ecr_repo::ECRRepo)::String
  """
  aws ecr list-images --repository-name=$(ecr_repo.repositoryName)
  """
end

function get_delete_remote_image_script(remote_image::RemoteImage)::String
  """
  aws ecr batch-delete-image \\
    --repository-name=$(remote_image.ecr_repo.repositoryName) \\
    --image-ids imageDigest=$(remote_image.imageDigest)
  """
end

function get_ecr_login_script(aws_config::AWSConfig, image_suffix::String)
  image_full_name = get_image_full_name(aws_config, image_suffix)
  """
  aws ecr get-login-password --region $(aws_config.region) \\
    | docker login \\
    --username AWS \\
    --password-stdin \\
    $image_full_name
  """
end

get_docker_push_script(image_full_name_plus_tag::String) = """
docker push $image_full_name_plus_tag
"""

function get_create_ecr_repo_script(image_suffix::String, aws_region::String, labels::Labels)::String
  tags_json = to_json(labels)
  """
  aws ecr create-repository \\
    --repository-name $(image_suffix) \\
    --image-scanning-configuration scanOnPush=true \\
    --region $(aws_region) \\
    --tags '$tags_json'
  """
end

function get_delete_ecr_repo_script(repo_name)::String
  """
  aws ecr delete-repository \\
    --force \\
    --repository-name $(repo_name)
  """
end

function get_create_lambda_role_script(role_name)::String
  """
  TRUST_POLICY=\$(cat <<EOF
  {
    "Version": "2012-10-17",
    "Statement": [
      {
        "Effect": "Allow",
        "Principal": {
          "Service": "lambda.amazonaws.com"
        },
        "Action": "sts:AssumeRole"
      }
    ]
  }
  EOF
  )

  aws iam create-role \\
    --role-name $(role_name) \\
    --assume-role-policy-document "\$TRUST_POLICY"
  """
end

function get_delete_lambda_role_script(role_name)::String
  """
  aws iam delete-role --role-name $(role_name)
  """
end

function get_create_lambda_function_script(
    function_name::String,
    repo_uri::String,
    role_arn::String,
    timeout::Int64,
    memory_size::Int64;
    labels::Labels,
  )::String
  tags_shorthand = to_aws_shorthand(labels)
  """
  aws lambda create-function \\
    --function-name=$(function_name) \\
    --code ImageUri=$(repo_uri) \\
    --role $(role_arn) \\
    --package-type Image \\
    --timeout=$(timeout) \\
    --memory-size=$(memory_size) \\
    --tags $tags_shorthand
  """
end

function get_lambda_function_status(
    function_name::String,
  )::String
  """
  aws lambda get-function-configuration --function-name $function_name
  """
end

function get_delete_lambda_function_script(function_arn::String)::String
  """
  aws lambda delete-function \\
    --function-name=$(function_arn)
  """
end

function get_invoke_lambda_function_script(
    function_arn::String,
    request::String,
    outfile::String,
  )::String
  """
  aws lambda invoke \\
    --function-name=$(function_arn) \\
    --payload='$(request)' \\
    --cli-binary-format raw-in-base64-out \\
    $outfile
  """
end