text:
"""
spin16.jl - getting numerical accuracy for the long simulations
"""

using Dates
using DifferentialEquations
using HDF5
using LinearAlgebra
using Printf
import Plots
using Random
using StaticArrays
using Statistics

include(joinpath(ENV["ROBUST_QOC_PATH"], "spin.jl"))

EXPERIMENT_META = "spin"
EXPERIMENT_NAME = "spin16"
WDIR = ENV["ROBUST_QOC_PATH"]
META_PATH = joinpath(WDIR, "out", EXPERIMENT_META)
SAVE_PATH = joinpath(META_PATH, EXPERIMENT_NAME)
VE_SAVE_FILE_NAME = "ve_spin16"

# plotting configuration
ENV["GKSwstype"] = "nul"
Plots.gr()

PT_TO_STR = Dict(
    analytic => "Analytic",
    vanilla => "Vanilla QOC",
    t1_m1 => "T1 QOC"
)

# types
@enum PulseType begin
    analytic = 1
    vanilla = 2
    t1_m1 = 3
end

# data constants
NON_TRAJ_OPT_PULSE = -1
CIDX_KEY = 4

PULSE_DATA = Dict(
    analytic => Dict(
        zpiby2 => joinpath(META_PATH, "spin14", "00000_spin14.h5"),
        ypiby2 => joinpath(META_PATH, "spin14", "00001_spin14.h5"),
        xpiby2 => joinpath(META_PATH, "spin14", "00002_spin14.h5"),
        CIDX_KEY => NON_TRAJ_OPT_PULSE,
    ),
    vanilla => Dict(
        zpiby2 => joinpath(META_PATH, "spin15", "00013_spin15.h5"),
        ypiby2 => joinpath(META_PATH, "spin15", "00016_spin15.h5"),
        xpiby2 => joinpath(META_PATH, "spin15", "00012_spin15.h5"),
        CIDX_KEY => 10,
    ),
    t1_m1 => Dict(
        zpiby2 => joinpath(META_PATH, "spin15", "00023_spin15.h5"),
        ypiby2 => joinpath(META_PATH, "spin15", "00015_spin15.h5"),
        xpiby2 => joinpath(META_PATH, "spin15", "00024_spin15.h5"),
        CIDX_KEY => 10,
    ),
)

# other constants
MAXITERS = 1e10
DPI = 500

function grab_controls(gate_type, pulse_type)
    controls_file_path = PULSE_DATA[pulse_type][gate_type]
    controls_idx = PULSE_DATA[pulse_type][CIDX_KEY]
    (controls, gate_time) = h5open(controls_file_path, "r") do save_file
        if controls_idx == NON_TRAJ_OPT_PULSE
            controls = permutedims(read(save_file, "controls"), [2, 1])
            gate_time = read(save_file, "evolution_time")
        else
            controls = read(save_file, "states")[:, controls_idx]
            gate_time = read(save_file, "evolution_time")
        end
        controls = controls / (2 * pi)
        return (controls, gate_time)
    end

    data = (controls, controls_file_path, gate_time)
    return data
end

# Define integration.
DT_CONTROLS = 1e-2
DT_CONTROLS_INV = 1e2
DT = 1e-2
DT_INV = 1e2

function rk4_step(dynamics, x, u)
    k1 = dynamics(x, u) * DT
	k2 = dynamics(x + k1 / 2, u) * DT
	k3 = dynamics(x + k2 / 2, u) * DT
	k4 = dynamics(x + k3, u) * DT
	return x + (k1 + 2 * k2 + 2 * k3 + k4)/6
end

function dynamics_lindblad_hb(density, controls)
    gamma = (get_t1_poly(controls[1]))^(-1)
    neg_i_hamiltonian = OMEGA_NEG_I_H_S + controls[1] * NEG_I_2_PI_H_C1
    return (
        neg_i_hamiltonian * density - density * neg_i_hamiltonian
        + gamma * (G_E * density * E_G + NEG_E_E_BY2 * density + density * NEG_E_E_BY2)
        + gamma * (E_G * density * G_E + NEG_G_G_BY2 * density + density * NEG_G_G_BY2)
    )
end

function dynamics_lindblad_nodis_hb(density, controls)
    neg_i_hamiltonian = OMEGA_NEG_I_H_S + controls[1] * NEG_I_2_PI_H_C1
    return (
        neg_i_hamiltonian * density - density * neg_i_hamiltonian
    )
end

function run_sim_hb(dissipation_type, gate_count, gate_type, pulse_type;
                    save=true, seed=0, plot=false)
    # grab
    data = grab_controls(gate_type, pulse_type)
    (controls, controls_file_path, gate_time) = data
    gate_knot_count = Int(gate_time * DT_INV)
    total_knot_count = gate_knot_count * gate_count
    densities = zeros(gate_count + 1, STATE_SIZE_ISO, STATE_SIZE_ISO)
    if dissipation_type == dissipation
        dynamics = dynamics_lindblad_hb
    else
        dynamics = dynamics_lindblad_nodis_hb
    end

    # integrate
    densities[1, :, :] = density = initial_density = gen_rand_density_iso(seed)
    for i = 1:gate_count
        for j = 1:gate_knot_count
            density = rk4_step(dynamics, density, controls[j])
        end
        densities[i + 1, :, :] = density
    end
    densities[end, :, :] = density
    fidelities = zeros(gate_count + 1)
    g = GT_TO_GATE[gate_type]
    g2 = g^2
    g3 = g^3
    id0 = initial_density
    id1 = g * id0 * g'
    id2 = g2 * id0 * g2'
    id3 = g3 * id0 * g3'
    target_dag = id0_dag = id0'
    id1_dag = id1'
    id2_dag = id2'
    id3_dag = id3'
    target_fnorm = id0_fnorm = abs(tr(id0_dag * id0))
    id1_fnorm = abs(tr(id1_dag * id1))
    id2_fnorm = abs(tr(id2_dag * id2))
    id3_fnorm = abs(tr(id3_dag * id3))
    for i = 1:gate_count + 1
        i_eff = i - 1
        if i_eff % 4 == 0
            target = id0
            target_dag = id0_dag
            target_fnorm = id0_fnorm
        elseif i_eff % 4 == 1
            target = id1
            target_dag = id1_dag
            target_fnorm = id1_fnorm
        elseif i_eff % 4 == 2
            target = id2
            target_dag = id2_dag
            target_fnorm = id2_fnorm
        elseif i_eff % 4 == 3
            target = id3
            target_dag = id3_dag
            target_fnorm = id3_fnorm
        end
        fidelities[i] = abs(tr(target_dag * densities[i, :, :])) / target_fnorm
        # println("fidelity\n$(fidelities[i])")
        # println("density")
        # show_nice(density)
        # println("")
        # println("target")
        # show_nice(target)
        # println("")
    end

    # save data
    if save
        data_file_path = generate_save_file_path("h5", EXPERIMENT_NAME, SAVE_PATH)
        h5open(data_file_path, "cw") do data_file
            write(data_file, "controls_file_path", controls_file_path)
            write(data_file, "dissipation_type", Integer(dissipation_type))
            write(data_file, "gate_count", gate_count)
            write(data_file, "gate_time", gate_time)
            write(data_file, "gate_type", Integer(gate_type))
            write(data_file, "pulse_type", Integer(pulse_type))
            write(data_file, "seed", seed)
            write(data_file, "densities", densities)
            write(data_file, "fidelities", fidelities)
        end
        println("Saved to $(data_file_path)")
    end

    # plot
    if plot
        plot_file_path = generate_save_file_path("png", EXPERIMENT_NAME, SAVE_PATH)
        gate_axis = Array(0:1:gate_count)
        fig = Plots.plot(
            dpi=DPI, title="$(GT_TO_STR[gate_type]) $(PT_TO_STR[pulse_type])",
            ylims=(0, 1), yticks=(0:0.1:1),
        )
        Plots.plot!(gate_axis, fidelities, label=nothing)
        Plots.xlabel!("Gate Count")
        Plots.ylabel!("Fidelity")
        Plots.savefig(fig, plot_file_path)
        println("Plotted to $(plot_file_path)")
    end
end

function ve_hb(gate_count; seed=0)
    gate_knot_count = Int(ZPIBY2_GATE_TIME * DT_INV)
    total_knot_count = gate_knot_count * gate_count
    densities = zeros(gate_count + 1, STATE_SIZE_ISO, STATE_SIZE_ISO)
    densities[1, :, :] = density = initial_density = gen_rand_density_iso(seed)
    for i = 1:gate_count
        for j = 1:gate_knot_count
            density = rk4_step(dynamics_lindblad_hb, density, 0)
        end
        densities[i + 1, :, :] = density
    end
    densities[end, :, :] = density
    fidelities = zeros(gate_count + 1)
    z = ZPIBY2
    z2 = z^2
    z3 = z^3
    id0 = initial_density
    id1 = z * id0 * z'
    id2 = z2 * id0 * z2'
    id3 = z3 * id0 * z3'
    target_dag = id0_dag = id0'
    id1_dag = id1'
    id2_dag = id2'
    id3_dag = id3'
    target_fnorm = id0_fnorm = abs(tr(id0_dag * id0))
    id1_fnorm = abs(tr(id1_dag * id1))
    id2_fnorm = abs(tr(id2_dag * id2))
    id3_fnorm = abs(tr(id3_dag * id3))
    for i = 1:gate_count + 1
        density = densities[i, :, :]
        i_eff = i - 1
        if i_eff % 4 == 0
            target = id0
            target_dag = id0_dag
            target_fnorm = id0_fnorm
        elseif i_eff % 4 == 1
            target = id1
            target_dag = id1_dag
            target_fnorm = id1_fnorm
        elseif i_eff % 4 == 2
            target = id2
            target_dag = id2_dag
            target_fnorm = id2_fnorm
        elseif i_eff % 4 == 3
            target = id3
            target_dag = id3_dag
            target_fnorm = id3_fnorm
        end
        fidelities[i] = abs(tr(target_dag * density)) / target_fnorm
        # println("fidelity\n$(fidelities[i])")
        # println("density")
        # show_nice(density)
        # println("")
        # println("target")
        # show_nice(target)
        # println("")
    end

    # save data
    data_file_path = generate_save_file_path("h5", VE_SAVE_FILE_NAME, SAVE_PATH)
    h5open(data_file_path, "cw") do data_file
        write(data_file, "densities", densities)
        write(data_file, "fidelities", fidelities)
    end
    println("Saved to $(data_file_path)")

    # plot
    plot_file_path = generate_save_file_path("png", VE_SAVE_FILE_NAME, SAVE_PATH)
    gate_axis = Array(0:1:gate_count)
    fig = Plots.plot(
        dpi=DPI, title="Verify Empty",
        ylims=(0, 1), yticks=(0:0.1:1),
    )
    Plots.plot!(gate_axis, fidelities, label=nothing)
    Plots.xlabel!("Gate Count")
    Plots.ylabel!("Fidelity")
    Plots.savefig(fig, plot_file_path)
    println("Plotted to $(plot_file_path)")
end

function plot_fidelity_by_gate_count_single(fig, path; inds=nothing)
    (dissipation_type, fidelities, gate_type,
     pulse_type) = h5open(path) do data_file
         dissipation_type = DissipationType(read(data_file, "dissipation_type"))
         fidelities = read(data_file, "fidelities")
         gate_type = GateType(read(data_file, "gate_type"))
         pulse_type = PulseType(read(data_file, "pulse_type"))
        return (dissipation_type, fidelities, gate_type, pulse_type)
    end
    gate_count = size(fidelities)[1] - 1
    gate_count_axis = Array(0:1:gate_count)
    label = "$(GT_TO_STR[gate_type]) $(PT_TO_STR[pulse_type]) $(DT_TO_STR[dissipation_type])"
    if isnothing(inds)
        inds = 1:gate_count + 1
    end
    Plots.plot!(fig, gate_count_axis[inds], fidelities[inds], label=label, legend=:bottomleft)
end

function plot_fidelity_by_gate_count(paths; inds=nothing, title=nothing)
    plot_file_path = generate_save_file_path("png", EXPERIMENT_NAME, SAVE_PATH)
    fig = Plots.plot(dpi=DPI, ylims=(0, 1), yticks=(0:0.1:1))
    for path in paths
        plot_fidelity_by_gate_count_single(fig, path; inds=inds)
    end
    Plots.ylabel!("Fidelity")
    Plots.xlabel!("Gate Count")
    Plots.savefig(fig, plot_file_path)
    println("Saved plot to $(plot_file_path)")
end

function dynamics_lindblad_deqjl(density, (controls, gate_knot_count), t)
    knot_point = (Int(floor(t * DT_INV)) % gate_knot_count) + 1
    controls_ = controls[knot_point]
    gamma = (get_t1_poly(controls_[1]))^(-1)
    neg_i_hamiltonian = OMEGA_NEG_I_H_S + controls_[1] * NEG_I_2_PI_H_C1
    return (
        neg_i_hamiltonian * density - density * neg_i_hamiltonian
        + gamma * (G_E * density * E_G + NEG_E_E_BY2 * density + density * NEG_E_E_BY2)
        + gamma * (E_G * density * G_E + NEG_G_G_BY2 * density + density * NEG_G_G_BY2)
    )
end

function dynamics_lindblad_nodis_deqjl(density, (controls, gate_knot_count), t)
    knot_point = (Int(floor(t * DT_INV)) % gate_knot_count) + 1
    controls_ = controls[knot_point]
    neg_i_hamiltonian = OMEGA_NEG_I_H_S + controls_[1] * NEG_I_2_PI_H_C1
    return (
        neg_i_hamiltonian * density - density * neg_i_hamiltonian
    )
end

function run_sim_deqjl(dissipation_type, gate_count, gate_type, pulse_type;
                       save=true, seed=0)
    start_time = Dates.now()
    # grab
    data = grab_controls(gate_type, pulse_type)
    (controls, controls_file_path, gate_time) = data
    gate_knot_count = Int(gate_time * DT_CONTROLS_INV)
    gate_times = Array(0:1:gate_count) * gate_time

    # integrate
    if dissipation_type == dissipation
        f = dynamics_lindblad_deqjl
    else
        f = dynamics_lindblad_nodis_deqjl
    end
    initial_density = gen_rand_density_iso(seed)
    tspan = (0., gate_time * gate_count)
    p = (controls, gate_knot_count)
    prob = ODEProblem(f, initial_density, tspan, p)
    result = solve(prob, DifferentialEquations.Tsit5(), dt=DT, saveat=gate_times,
                   maxiters=MAXITERS, adaptive=false)

    # compute fidelity
    densities = zeros(gate_count + 1, STATE_SIZE_ISO, STATE_SIZE_ISO)
    fidelities = zeros(gate_count + 1)
    g = GT_TO_GATE[gate_type]
    g2 = g^2
    g3 = g^3
    id0 = initial_density
    id1 = g * id0 * g'
    id2 = g2 * id0 * g2'
    id3 = g3 * id0 * g3'
    target_dag = id0_dag = id0'
    id1_dag = id1'
    id2_dag = id2'
    id3_dag = id3'
    target_fnorm = id0_fnorm = abs(tr(id0_dag * id0))
    id1_fnorm = abs(tr(id1_dag * id1))
    id2_fnorm = abs(tr(id2_dag * id2))
    id3_fnorm = abs(tr(id3_dag * id3))
    for i = 1:gate_count + 1
        densities[i, :, :] = density = result.u[i]
        i_eff = i - 1
        if i_eff % 4 == 0
            target = id0
            target_dag = id0_dag
            target_fnorm = id0_fnorm
        elseif i_eff % 4 == 1
            target = id1
            target_dag = id1_dag
            target_fnorm = id1_fnorm
        elseif i_eff % 4 == 2
            target = id2
            target_dag = id2_dag
            target_fnorm = id2_fnorm
        elseif i_eff % 4 == 3
            target = id3
            target_dag = id3_dag
            target_fnorm = id3_fnorm
        end
        fidelities[i] = abs(tr(target_dag * density)) / target_fnorm
        # println("fidelity\n$(fidelities[i])")
        # println("density")
        # show_nice(density)
        # println("")
        # println("target")
        # show_nice(target)
        # println("")
    end
    end_time = Dates.now()
    run_time = end_time - start_time

    # save data
    if save
        data_file_path = generate_save_file_path("h5", EXPERIMENT_NAME, SAVE_PATH)
        h5open(data_file_path, "cw") do data_file
            write(data_file, "controls_file_path", controls_file_path)
            write(data_file, "dissipation_type", Integer(dissipation_type))
            write(data_file, "gate_count", gate_count)
            write(data_file, "gate_time", gate_time)
            write(data_file, "gate_type", Integer(gate_type))
            write(data_file, "pulse_type", Integer(pulse_type))
            write(data_file, "seed", seed)
            write(data_file, "densities", densities)
            write(data_file, "fidelities", fidelities)
            write(data_file, "run_time", string(run_time))
        end
        println("Saved to $(data_file_path)")
    end
end

function dynamics_ve_deqjl(u, p, t)
    return OMEGA_NEG_I_H_S * u - u * OMEGA_NEG_I_H_S
end

function ve_deqjl(gate_count; seed=0, data_file_path=nothing)
    gate_times = Array(0:1:gate_count) * ZPIBY2_GATE_TIME
    if isnothing(data_file_path)
        initial_density = u0 = gen_rand_density_iso(seed)
        tspan = (0., ZPIBY2_GATE_TIME * gate_count)
        prob = ODEProblem(dynamics_ve_deqjl, u0, tspan)
        result = solve(prob, DifferentialEquations.TsitPap8(), dt=1e-3, saveat=gate_times, adaptive=false)
    end
    densities = zeros(gate_count + 1, STATE_SIZE_ISO, STATE_SIZE_ISO)
    fidelities = zeros(gate_count + 1)
    z = ZPIBY2
    z2 = z^2
    z3 = z^3
    id0 = initial_density
    id1 = z * id0 * z'
    id2 = z2 * id0 * z2'
    id3 = z3 * id0 * z3'
    target_dag = id0_dag = id0'
    id1_dag = id1'
    id2_dag = id2'
    id3_dag = id3'
    target_fnorm = id0_fnorm = abs(tr(id0_dag * id0))
    id1_fnorm = abs(tr(id1_dag * id1))
    id2_fnorm = abs(tr(id2_dag * id2))
    id3_fnorm = abs(tr(id3_dag * id3))
    for i = 1:gate_count + 1
        density = result.u[i]
        densities[i, :, :] = Array(density)
        i_eff = i - 1
        if i_eff % 4 == 0
            target = id0
            target_dag = id0_dag
            target_fnorm = id0_fnorm
        elseif i_eff % 4 == 1
            target = id1
            target_dag = id1_dag
            target_fnorm = id1_fnorm
        elseif i_eff % 4 == 2
            target = id2
            target_dag = id2_dag
            target_fnorm = id2_fnorm
        elseif i_eff % 4 == 3
            target = id3
            target_dag = id3_dag
            target_fnorm = id3_fnorm
        end
        fidelities[i] = abs(tr(target_dag * density)) / target_fnorm
        # println("fidelity\n$(fidelities[i])")
        # println("density")
        # show_nice(density)
        # println("")
        # println("target")
        # show_nice(target)
        # println("")
    end

    # save data
    data_file_path = generate_save_file_path("h5", VE_SAVE_FILE_NAME, SAVE_PATH)
    h5open(data_file_path, "cw") do data_file
        write(data_file, "densities", densities)
        write(data_file, "fidelities", fidelities)
    end
    println("Saved to $(data_file_path)")

    # plot
    plot_file_path = generate_save_file_path("png", VE_SAVE_FILE_NAME, SAVE_PATH)
    gate_axis = Array(0:1:gate_count)
    fig = Plots.plot(
        dpi=DPI, title="Verify Empty",
        ylims=(0, 1), yticks=(0:0.1:1),
    )
    Plots.plot!(gate_axis, fidelities, label=nothing)
    Plots.xlabel!("Gate Count")
    Plots.ylabel!("Fidelity")
    Plots.savefig(fig, plot_file_path)
    println("Plotted to $(plot_file_path)")
end

max_stars_repo_name: SchusterLab/rbqoc

max_issues_repo_name: SchusterLab/rbqoc

max_forks_repo_name: SchusterLab/rbqoc


text:
"""
    NormalHistLoss

Create loss function from observed histogram and fitted probabilities from normal distribution.

# Notice
n interval boundaries, z(1), z(2), and z(3) for which the survey gives the following probabilities:
Pr(x < z(1)), Pr(z(1) <= x <z(2)), Pr(z(2) <= x< z(3)), and Pr(z(3) <= x),
so there are n+1 categories (and therefore probabilities)

Notice: both the lowest and the highest intervals are open

Paul Soderlind, February 2000, to Julia Dec 2016

"""
function NormalHistLoss(par,Probs,Bounds)

  n       = length(Bounds)                  #n bounds, n+1 categories
  (mu,s2) = (par[1],par[2]^2)               #square to avoid s<0

  xnorm        = (Bounds.-mu)./sqrt(s2)                  #(x-mu)/s
  TheoryProb_x = 0.5 .+ 0.5*erf.(xnorm/sqrt(2))          #theoretical Pr( z<=x[i] )

  TheoryProb_Interval = [ TheoryProb_x[1]                        ;
                          TheoryProb_x[2:n] - TheoryProb_x[1:n-1];
                          1.0 - TheoryProb_x[n]                  ]

  loss = 1.0 + 10000*sum( (TheoryProb_Interval-Probs).^2 )  #loss function

  return loss

end
#----------------------------------------------------------------------------

"""
    NormPdfPs

Pdf value of a normally distributed variable

"""
function NormPdfPs(x,mu=0,s2=1)

  s    = sqrt(s2)
  z    = (x - mu)/s
  pdfx = exp(-0.5*z^2)/(sqrt(2*pi)*s)

  return pdfx

end
#----------------------------------------------------------------------------

"""
  TriangularPdfPs

Pdf value of a triangularly distributed variable over (a,b) with mode c

"""
function TriangularPdfPs(x,a,b,c)

  if (x<a) || (x>b)     #outside (a,b)
    pdfx = 0
  elseif a <= x < c
    pdfx = 2*(x-a)/((b-a)*(c-a))
  elseif y == c
    pdfx = 2/(b-a)
  else
    pdfx = 2*(b-x)/((b-a)*(b-c))
  end

  return pdfx

end
#------------------------------------------------------------------------------

ext: jl

max_stars_repo_name: PaulSoderlind/PaulSoderlindCode

max_issues_repo_name: PaulSoderlind/PaulSoderlindCode

max_forks_repo_name: PaulSoderlind/PaulSoderlindCode


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

module RunCommandsClient

using Random
using Dates
using HTTP
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_VirtualMachineRunCommandsApi.jl")
include("api_VirtualMachinesApi.jl")

# export models
export convert, ApiError
export convert, ApiErrorBase
export convert, InnerError
export convert, OperationStatusResponse
export convert, RunCommandDocumentBase
export convert, RunCommandInput
export convert, RunCommandInputParameter
export convert, RunCommandListResult
export convert, RunCommandParameterDefinition
export convert, RunCommandResultProperties
export convert, RunCommandDocument
export convert, RunCommandResult

# export operations
export convert, VirtualMachineRunCommandsApi, VirtualMachinesApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
export minkowski_sum

"""
    minkowski_sum(P::LazySet{N}, Q::LazySet{N};
                  [backend]=nothing,
                  [algorithm]=nothing,
                  [prune]=true) where {N<:Real}

Concrete Minkowski sum for a pair of lazy sets using their constraint representation.

### Input

- `P`         -- lazy set
- `Q`         -- another lazy set
- `backend`   -- (optional, default: `nothing`) polyhedral computations backend
- `algorithm` -- (optional, default: `nothing`) algorithm to compute the elimination
                 of variables; available options are `Polyhedra.FourierMotzkin`,
                 `Polyhedra.BlockElimination`, and `Polyhedra.ProjectGenerators`
- `prune`     -- (optional, default: `true`) if `true`, apply a post-processing algorithm
                 to remove redundant constraints

### Output

An `HPolytope` that corresponds to the Minkowski sum of `P` and `Q` if both `P`
and `Q` are bounded; otherwise an `HPolyhedron`.

### Notes

This function requires that the list of constraints of both lazy sets `P` and
`Q` can be obtained. After obtaining the respective lists of constraints, the
`minkowski_sum` fucntion for polyhedral sets is used. For details see
[`minkowski_sum(::VPolytope, ::VPolytope)`](@ref).

This method requires `Polyhedra` and `CDDLib`, so you have to do:

```julia
julia> using LazySets, Polyhedra, CDDLib

julia> ...

julia> minkowski_sum(P, Q)
```
"""
function minkowski_sum(P::LazySet{N}, Q::LazySet{N};
                       backend=nothing,
                       algorithm=nothing,
                       prune=true) where {N<:Real}
    @assert applicable(constraints_list, P) &&
        applicable(constraints_list, Q) "this function requires that the " *
        "list of constraints is available for both arguments; try " *
        "overapproximating with an `HPolytope` first"

    res = _minkowski_sum_hrep_preprocess(P, Q, backend, algorithm, prune)
    if isbounded(P) && isbounded(Q)
        return convert(HPolytope, res)
    else
        return res
    end
end

"""
    minkowski_sum(P::AbstractPolyhedron{N}, Q::AbstractPolyhedron{N};
                  [backend]=nothing,
                  [algorithm]=nothing,
                  [prune]=true) where {N<:Real}

Compute the Minkowski sum between two polyhedra in constraint representation.

### Input

- `P`         -- polyhedron in constraint representation
- `Q`         -- another polyhedron in constraint representation
- `backend`   -- (optional, default: `nothing`) polyhedral computations backend
- `algorithm` -- (optional, default: `nothing`) algorithm to compute the elimination
                 of variables; available options are `Polyhedra.FourierMotzkin`,
                 `Polyhedra.BlockElimination`, and `Polyhedra.ProjectGenerators`
- `prune`     -- (optional, default: `true`) if `true`, apply a post-processing algorithm
                 to remove redundant constraints

### Output

A polyhedron in H-representation that corresponds to the Minkowski sum of `P` and `Q`.

### Notes

This method requires `Polyhedra` and `CDDLib`, so you have to do:

```julia
julia> using LazySets, Polyhedra, CDDLib

julia> ...

julia> minkowski_sum(P, Q)
```

### Algorithm

This function implements the concrete Minkowski sum by projection and variable
elimination as detailed in [1]. The idea is that if we write ``P`` and ``Q`` in
*simple H-representation*, that is, ``P = \\{x ∈ \\mathbb{R}^n : Ax ≤ b \\}``
and ``Q = \\{x ∈ \\mathbb{R}^n : Cx ≤ d \\}``, then their Minkowski sum can be
seen as the projection onto the first ``n``-dimensional coordinates of the polyhedron

```math
    \\begin{pmatrix} 0 & A \\ C & -C \\end{pmatrix} \\binom{x}{y} ≤ \binom{b}{d}
```
This is seen by noting that ``P ⊕ Q`` corresponds to the set of points
``x ∈ \\mathbb{R}^n`` such that ``x = y + z`` with ``Ay ≤ b`` and ``Cz ≤ d``;
hence it follows that ``Ay ≤ b`` and ``C(x-y) ≤ d``, and the inequality displayed
above follows by considering the ``2n``-dimensional space ``\\binom{x}{y}``.
The reduction from ``2n`` to ``n`` variables is performed using an elimination
algorithm as described next.

The elimination of variables depends on the concrete polyhedra library `Polyhedra`,
which itself uses `CDDLib` for variable elimination. The available algorithms are:

- `Polyhedra.FourierMotzkin`   -- computation of the projection by computing the
                                  H-representation and applying the Fourier-Motzkin
                                  elimination algorithm to it

- `Polyhedra.BlockElimination` -- computation of the projection by computing the
                                  H-representation and applying the block elimination
                                  algorithm to it

- `Polyhedra.ProjectGenerators` -- computation of the projection by computing the
                                   V-representation

[1] Kvasnica, Michal. "Minkowski addition of convex polytopes." (2005): 1-10.
"""
function minkowski_sum(P::AbstractPolyhedron{N}, Q::AbstractPolyhedron{N};
                       backend=nothing,
                       algorithm=nothing,
                       prune=true) where {N<:Real}
    return _minkowski_sum_hrep_preprocess(P, Q, backend, algorithm, prune)
end

function minkowski_sum(P::HPolytope{N}, Q::HPolytope{N};
                       backend=nothing,
                       algorithm=nothing,
                       prune=true) where {N<:Real}
    res = _minkowski_sum_hrep_preprocess(P, Q, backend, algorithm, prune)
    return convert(HPolytope, res)
end

# common code before calling _minkowski_sum_hrep
function _minkowski_sum_hrep_preprocess(P, Q, backend, algorithm, prune)
    require(:Polyhedra; fun_name="minkowski_sum")
    require(:CDDLib; fun_name="minkowski_sum")

    A, b = tosimplehrep(P)
    C, d = tosimplehrep(Q)
    return _minkowski_sum_hrep(A, b, C, d, backend=backend, algorithm=algorithm,
                               prune=prune)
end

# This function computes the concrete Minkowski sum between two polyhedra in
# simple H-representation,
# P = {x : Ax <= b} and Q = {x : Cx <= d}
# using the projection methods. See the documentation of `minkowski_sum` for details.
function _minkowski_sum_hrep(A::AbstractMatrix{N}, b::AbstractVector{N},
                             C::AbstractMatrix{N}, d::AbstractVector{N};
                             backend=nothing,
                             algorithm=nothing,
                             prune=true) where {N<:Real}

    if backend == nothing
        backend = default_cddlib_backend(N)
    end

    if algorithm == nothing
        algorithm = Polyhedra.FourierMotzkin()
    elseif !(algorithm <: EliminationAlgorithm)
        error("the algorithm $algorithm is not a valid elimination algorithm;
              choose among any of $(subtypes(Polyhedra.EliminationAlgorithm))")
    end

    mP, nP = size(A)
    mQ, nQ = size(C)
    E = [zeros(N, mP, nQ) A; C -C]
    f = [b; d]
    PQ = HPolyhedron(E, f)
    PQ_cdd = polyhedron(PQ, backend=backend)
    W = HPolyhedron(Polyhedra.eliminate(PQ_cdd, nP+1:2nP, algorithm))
    if prune
        success = remove_redundant_constraints!(W)
        if !success
            error("the constraints corresponding to the minkowski sum of the given " *
                  "sets are infeasible")
        end
    end
    return W
end

"""
    minkowski_sum(H1::AbstractHyperrectangle{N}, H2::AbstractHyperrectangle{N})
        where {N<:Real}

Concrete Minkowski sum of a pair of hyperrectangular sets.

### Input

- `H1` -- hyperrectangular set
- `H2` -- hyperrectangular set

### Output

A `Hyperrectangle` corresponding to the concrete Minkowski sum of `H1` and `H2`.

### Algorithm

The resulting hyperrectangle is obtained by summing up the centers and
radiuses of `H1` and `H2`.
"""
function minkowski_sum(H1::AbstractHyperrectangle{N},
                       H2::AbstractHyperrectangle{N}) where {N<:Real}
    return Hyperrectangle(center(H1) + center(H2),
                          radius_hyperrectangle(H1) + radius_hyperrectangle(H2))
end

"""
    minkowski_sum(Z1::AbstractZonotope{N}, Z2::AbstractZonotope{N})
        where {N<:Real}

Concrete Minkowski sum of a pair of zonotopic sets.

### Input

- `Z1` -- zonotopic set
- `Z2` -- zonotopic set

### Output

A `Zonotope` corresponding to the concrete Minkowski sum of `Z1` and `Z2`.

### Algorithm

The resulting zonotope is obtained by summing up the centers and concatenating
the generators of `Z1` and `Z2`.
"""
function minkowski_sum(Z1::AbstractZonotope{N},
                       Z2::AbstractZonotope{N}) where {N<:Real}
    return Zonotope(center(Z1) + center(Z2), [genmat(Z1) genmat(Z2)])
end


text:
module SeisProcessing
    using Interpolations
    using Requires,
    FFTW,
    LinearAlgebra,
    DSP
    import SeisMain

    include("Processing/Processing.jl")
    include("Modelling/Modelling.jl")
    include("Tools/Tools.jl")
end


text:
using MultipleViewGeometry
using Documenter

DocMeta.setdocmeta!(MultipleViewGeometry, :DocTestSetup, :(using MultipleViewGeometry); recursive=true)

makedocs(;
    modules=[MultipleViewGeometry],
    authors="sychen52 <sychen52@gmail.com> and contributors",
    repo="https://github.com/sychen52/MultipleViewGeometry.jl/blob/{commit}{path}#{line}",
    sitename="MultipleViewGeometry.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://sychen52.github.io/MultipleViewGeometry.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/sychen52/MultipleViewGeometry.jl",
    devbranch="main",
)


text:
using .Networks

using Dates
using JuMP
using Printf
using Parameters

@with_kw mutable struct RSOAssessmentConfigs
    BIG_M = 1e3 #max supposed overflow

    out_path = nothing
    problem_name = "RSOAssessment"
end

@with_kw struct RSOAssessment <: AbstractRunnable
    configs::RSOAssessmentConfigs = RSOAssessmentConfigs()
end

function run(runnable::RSOAssessment, ech::Dates.DateTime, TS::Vector{Dates.DateTime}, context::AbstractContext)
    problem_name_l = @sprintf("rso_assessment_%s", ech)
    runnable.configs.out_path = context.out_dir
    runnable.configs.problem_name = problem_name_l

    #FIXME : only create vars for started units

    return rso_assessment(get_network(context),
                        TS,
                        get_assessment_uncertainties(context),
                        get_tso_actions(context),
                        runnable.configs
                    )
end

function rso_assessment(network,
                        TS,
                        assessment_uncertainties,
                        tso_actions,
                        configs::RSOAssessmentConfigs
                    )
    model_container_l = formulate_rso_assessment(network, TS, assessment_uncertainties, tso_actions, configs)
    solve!(model_container_l, configs.problem_name, configs.out_path)

    for var in all_variables(model_container_l.model)
        println(name(var), " = ", value(var))
    end

    return model_container_l
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
@testset "device pointer" begin

# inner constructors

voidptr_a = CuPtr{Cvoid}(Int(0xDEADBEEF))
generic_voidptr_a = CUDA.DevicePtr{Cvoid,AS.Generic}(voidptr_a)
global_voidptr_a = CUDA.DevicePtr{Cvoid,AS.Global}(voidptr_a)
local_voidptr_a = CUDA.DevicePtr{Cvoid,AS.Local}(voidptr_a)

voidptr_b = CuPtr{Cvoid}(Int(0xCAFEBABE))
generic_voidptr_b = CUDA.DevicePtr{Cvoid,AS.Generic}(voidptr_b)
global_voidptr_b = CUDA.DevicePtr{Cvoid,AS.Global}(voidptr_b)
local_voidptr_b = CUDA.DevicePtr{Cvoid,AS.Local}(voidptr_b)

intptr_b = convert(CuPtr{Int}, voidptr_b)
generic_intptr_b = CUDA.DevicePtr{Int,AS.Generic}(intptr_b)
global_intptr_b = CUDA.DevicePtr{Int,AS.Global}(intptr_b)
local_intptr_b = CUDA.DevicePtr{Int,AS.Local}(intptr_b)

# outer constructors
@test CUDA.DevicePtr{Cvoid}(voidptr_a) == generic_voidptr_a
@test CUDA.DevicePtr(voidptr_a) == generic_voidptr_a

# getters
@test eltype(generic_voidptr_a) == Cvoid
@test eltype(global_intptr_b) == Int
@test addrspace(generic_voidptr_a) == AS.Generic
@test addrspace(global_voidptr_a) == AS.Global
@test addrspace(local_voidptr_a) == AS.Local

# comparisons
@test generic_voidptr_a != global_voidptr_a
@test generic_voidptr_a != generic_intptr_b

@testset "conversions" begin

# between host and device pointers

@test convert(CuPtr{Cvoid}, generic_voidptr_a) == voidptr_a
@test convert(CUDA.DevicePtr{Cvoid}, voidptr_a) == generic_voidptr_a
@test convert(CUDA.DevicePtr{Cvoid,AS.Global}, voidptr_a) == global_voidptr_a

# between device pointers

@test_throws ArgumentError convert(typeof(local_voidptr_a), global_voidptr_a)
@test convert(typeof(generic_voidptr_a), generic_voidptr_a) == generic_voidptr_a
@test convert(typeof(global_voidptr_a), global_voidptr_a) == global_voidptr_a
@test Base.unsafe_convert(typeof(local_voidptr_a), global_voidptr_a) == local_voidptr_a

@test convert(typeof(global_voidptr_a), global_intptr_b) == global_voidptr_b
@test convert(typeof(generic_voidptr_a), global_intptr_b) == generic_voidptr_b
@test convert(typeof(global_voidptr_a), generic_intptr_b) == global_voidptr_b

@test convert(CUDA.DevicePtr{Cvoid}, global_intptr_b) == global_voidptr_b

end

@testset "unsafe_load & unsafe_store!" begin

@testset for T in (Int8, UInt16, Int32, UInt32, Int64, UInt64, Int128, Float32, Float64),
             cached in (false, true)
    d_a = CuArray(ones(T))
    d_b = CuArray(zeros(T))

    ptr_a = CUDA.DevicePtr{T,AS.Global}(pointer(d_a))
    ptr_b = CUDA.DevicePtr{T,AS.Global}(pointer(d_b))
    @test Array(d_a) != Array(d_b)

    let ptr_a=ptr_a, ptr_b=ptr_b #JuliaLang/julia#15276
        if cached && capability(device()) >= v"3.2"
            @on_device unsafe_store!(ptr_b, unsafe_cached_load(ptr_a))
        else
            @on_device unsafe_store!(ptr_b, unsafe_load(ptr_a))
        end
    end
    @test Array(d_a) == Array(d_b)
end

@testset "indexing" begin
    function kernel(src, dst)
        unsafe_store!(dst, unsafe_cached_load(src, 4))
        return
    end

    T = Complex{Int8}
    # this also tests the fallback for unsafe_cached_load

    src = CuArray([T(1) T(9); T(3) T(4)])
    dst = CuArray([0])

    @cuda kernel(
        CUDA.DevicePtr{T,AS.Global}(pointer(src)),
        CUDA.DevicePtr{T,AS.Global}(pointer(dst))
    )

    @test Array(src)[4] == Array(dst)[1]
end

end

end


text:
"""
    bcvumps_env(model::MT, β, D; tol=1e-10, maxiter=20, verbose = false) where {MT <: HamiltonianModel}

return the bcvumps environment of the `model` as a function of the inverse
temperature `β` and the environment bonddimension `D` as calculated with
vumps. Save `env` in file `./data/model_β_D.jld2`. Requires that `model_tensor` are defined for `model`.
"""
function bcvumps_env(model::MT, β, D; tol=1e-10, maxiter=20, verbose = false) where {MT <: HamiltonianModel}
    M = model_tensor(model, β)
    mkpath("./data/")
    chkp_file = "./data/$(model)_β$(β)_D$(D).jld2"
    if isfile(chkp_file)
        rt = SquareBCVUMPSRuntime(M, chkp_file, D; verbose = verbose)
    else
        rt = SquareBCVUMPSRuntime(M, Val(:random), D; verbose = verbose)
    end
    env = bcvumps(rt; tol=tol, maxiter=maxiter, verbose = verbose)
    return env
end

"""
    Z(env::SquareBCVUMPSRuntime)

return the partition function of the `env`.
"""
function Z(env::SquareBCVUMPSRuntime)
    M,AL,C,FL,FR = env.M,env.AL,env.C,env.FL,env.FR
    Ni,Nj = size(M)
    ACij = [ein"asc,cb -> asb"(AL[i],C[i]) for i=1:Ni*Nj]
    AC = reshape(ACij,Ni,Nj)
    z_tol = 1
    for j = 1:Nj,i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        jr = j + 1 - Nj * (j==Nj)
        z = ein"αcβ,βsη,cpds,ηdγ,αpγ ->"(FL[i,j],AC[i,j],M[i,j],FR[i,j],conj(AC[ir,j]))[]
        λ = ein"αcβ,βη,ηcγ,αγ ->"(FL[i,jr],C[i,j],FR[i,j],conj(C[ir,j]))[]
        z_tol *= z/λ
    end
    return z_tol^(1/Ni/Nj)
end

"""
    magnetisation(env::SquareBCVUMPSRuntime, model::MT, β)

return the magnetisation of the `model`. Requires that `mag_tensor` are defined for `model`.
"""
function magnetisation(env::SquareBCVUMPSRuntime, model::MT, β) where {MT <: HamiltonianModel}
    M,AL,C,FL,FR = env.M,env.AL,env.C,env.FL,env.FR
    Ni,Nj = size(M)
    ACij = [ein"asc,cb -> asb"(AL[i],C[i]) for i=1:Ni*Nj]
    AC = reshape(ACij,Ni,Nj)
    Mag = mag_tensor(model, β)
    mag_tol = 0
    for j = 1:Nj,i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        mag = ein"αcβ,βsη,cpds,ηdγ,αpγ ->"(FL[i,j],AC[i,j],Mag[i,j],FR[i,j],conj(AC[ir,j]))[]
        λ = ein"αcβ,βsη,cpds,ηdγ,αpγ ->"(FL[i,j],AC[i,j],M[i,j],FR[i,j],conj(AC[ir,j]))[]
        mag_tol += mag/λ
    end
    return abs(mag_tol)/Ni/Nj
end

"""
    energy(env::SquareBCVUMPSRuntime, model::MT, β)

return the energy of the `model` as a function of the inverse
temperature `β` and the environment bonddimension `D` as calculated with
vumps. Requires that `model_tensor` are defined for `model`.
"""
function energy(env::SquareBCVUMPSRuntime, model::MT, β::Real) where {MT <: HamiltonianModel}
    M,AL,C,FL,FR = env.M,env.AL,env.C,env.FL,env.FR
    Ni,Nj = size(M)
    ACij = [ein"asc,cb -> asb"(AL[i],C[i]) for i=1:Ni*Nj]
    AC = reshape(ACij,Ni,Nj)
    Ene = energy_tensor(model, β)
    ene_tol = 0
    for j = 1:Nj,i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        ene = ein"αcβ,βsη,cpds,ηdγ,αpγ ->"(FL[i,j],AC[i,j],Ene[i,j],FR[i,j],conj(AC[ir,j]))[]
        λ = ein"αcβ,βsη,cpds,ηdγ,αpγ ->"(FL[i,j],AC[i,j],M[i,j],FR[i,j],conj(AC[ir,j]))[]
        ene_tol += ene/λ
    end
    return ene_tol/Ni/Nj
end

"""
    magofβ(::Ising,β)

return the analytical result for the magnetisation at inverse temperature
`β` for the 2d classical ising model.
"""
magofβ(::Ising, β) = β > isingβc ? (1-sinh(2*β)^-4)^(1/8) : 0.

"""
    eneofβ(::Ising,β)

return some the numerical integrations of analytical result for the energy at inverse temperature
`β` for the 2d classical ising model.
"""
function eneofβ(::Ising, β)
    if β == 0.2
        return -0.42822885693016843
    elseif β == 0.4
        return -1.1060792706185651
    elseif β == 0.6
        return -1.909085845408498
    elseif β == 0.8
        return -1.9848514445364174
    end
end

"""
    Zofβ(::Ising,β)

return some the numerical integrations of analytical result for the partition function at inverse temperature
`β` for the 2d classical ising model.
"""
function Zofβ(::Ising,β)
    if β == 0.2
        return 2.08450374046259
    elseif β == 0.4
        return 2.4093664345022363
    elseif β == 0.6
        return 3.3539286863974582
    elseif β == 0.8
        return 4.96201030069517
    end
end


text:
using statmech_tm_solver
using Documenter

DocMeta.setdocmeta!(statmech_tm_solver, :DocTestSetup, :(using statmech_tm_solver); recursive=true)

makedocs(;
    modules=[statmech_tm_solver],
    authors="Wei Tang <tangwei@smail.nju.edu.cn> and contributors",
    repo="https://github.com/wei.tang.nju@gmail.com/statmech_tm_solver.jl/blob/{commit}{path}#{line}",
    sitename="statmech_tm_solver.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://wei.tang.nju@gmail.com.github.io/statmech_tm_solver.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/wei.tang.nju@gmail.com/statmech_tm_solver.jl",
)


text:
function acronym(phrase)

end


text:
module Architectures

export AbstractArchitecture, AbstractCPUArchitecture, AbstractGPUArchitecture
export CPU, GPU
export device, architecture, array_type, arch_array

using CUDA
using KernelAbstractions
using Adapt
using OffsetArrays

"""
    AbstractArchitecture

Abstract supertype for architectures supported by Oceananigans.
"""
abstract type AbstractArchitecture end

"""
    AbstractCPUArchitecture

Abstract supertype for CPU architectures supported by Oceananigans.
"""
abstract type AbstractCPUArchitecture <: AbstractArchitecture end

"""
    AbstractGPUArchitecture

Abstract supertype for GPU architectures supported by Oceananigans.
"""
abstract type AbstractGPUArchitecture <: AbstractArchitecture end

"""
    CPU <: AbstractArchitecture

Run Oceananigans on one CPU node. Uses multiple threads if the environment
variable `JULIA_NUM_THREADS` is set.
"""
struct CPU <: AbstractCPUArchitecture end

"""
    GPU <: AbstractArchitecture

Run Oceananigans on a single NVIDIA CUDA GPU.
"""
struct GPU <: AbstractGPUArchitecture end

device(::AbstractCPUArchitecture) = KernelAbstractions.CPU()
device(::AbstractGPUArchitecture) = KernelAbstractions.CUDADevice()

architecture() = nothing
architecture(::Number) = nothing
architecture(::Array) = CPU()
architecture(::CuArray) = GPU()

array_type(::CPU) = Array
array_type(::GPU) = CuArray

arch_array(::AbstractCPUArchitecture, A::Array) = A
arch_array(::AbstractCPUArchitecture, A::CuArray) = Array(A)
arch_array(::AbstractGPUArchitecture, A::Array) = CuArray(A)
arch_array(::AbstractGPUArchitecture, A::CuArray) = A

const OffsetCPUArray = OffsetArray{FT, N, <:Array} where {FT, N}
const OffsetGPUArray = OffsetArray{FT, N, <:CuArray} where {FT, N}

Adapt.adapt_structure(::CPU, a::OffsetCPUArray) = a
Adapt.adapt_structure(::GPU, a::OffsetGPUArray) = a

Adapt.adapt_structure(::GPU, a::OffsetCPUArray) = OffsetArray(CuArray(a.parent), a.offsets...)
Adapt.adapt_structure(::CPU, a::OffsetGPUArray) = OffsetArray(Array(a.parent), a.offsets...)

device_event(arch) = Event(device(arch))

end


text:
function fibo_even_sum(N::Int)
    a = 0
    b = 1
    c = 0
    s = 0
    while(a+b<N)
        c = a + b
        if iseven(c)
            s+=c
        end
        a = b
        b = c
    end
    return s
end

println(fibo_even_sum(4000000))


text:
# # KissABC
# ## Usage guide
# The ingredients you need to use Approximate Bayesian Computation:
#
# 1. A simulation which depends on some parameters, able to generate datasets similar to your target dataset if parameters are tuned
# 1. A prior distribution over such parameters
# 1. A distance function to compare generated dataset to the true dataset
#
# We will start with a simple example, we have a dataset generated according to an Normal distribution whose parameters are unknown

tdata = randn(1000) .* 0.04 .+ 2;

# we are ofcourse able to simulate normal random numbers, so this constitutes our simulation

sim((μ, σ)) = randn(1000) .* σ .+ μ;

# The second ingredient is a prior over the parameters μ and σ

using KissABC
prior = Factored(Uniform(1, 3), Truncated(Normal(0, 0.1), 0, 100));

# we have chosen a uniform distribution over the interval [1,3] for μ and a normal distribution truncated over ℝ⁺ for σ.
#
# Now all that we need is a distance function to compare the true dataset to the simulated dataset, for this purpose comparing mean and std is optimal

function dist(x, y)
    d1 = mean(x) - mean(y)
    d2 = std(x) - std(y)
    hypot(d1, d2 * 50)
end

# Now we are all set, we can use `AIS` which is an Affine Invariant MC algorithm via the `sample` function, to simulate the posterior distribution for this model, inferring μ and σ
cost(x) = dist(tdata, sim(x))
approx_density = ApproxPosterior(prior, cost, 0.01)
res = sample(
    approx_density,
    AIS(50),
    500,
    discard_initial = 1000,
    ntransitions = 10,
    progress = true,
)

@show res

# You can also use Sequential Monte Carlo (SMC) to infer posterior parameters:

ressmc = smc(prior, cost, nparticles=500, epstol=0.01)

@show ressmc

# the parameters we chose are: a tolerance on distances equal to `0.1`, a number of samples equal to `2000`, the simulated posterior results are in `res`.
# We can now extract the inference results:

prsample = [rand(prior) for i = 1:2000] #some samples from the prior for comparison
μ_pr = getindex.(prsample, 1) # μ samples from the prior
σ_pr = getindex.(prsample, 2) # σ samples from the prior
μ_p = res[1].particles # μ samples from the posterior
σ_p = res[2].particles; # σ samples from the posterior

# and plotting prior and posterior side by side we get:

using Plots
a = stephist(
    μ_pr,
    xlims = (1, 3),
    xlabel = "μ prior",
    leg = false,
    lw = 2,
    normalize = true,
)
b = stephist(
    σ_pr,
    xlims = (0, 0.3),
    xlabel = "σ prior",
    leg = false,
    lw = 2,
    normalize = true,
)
ap = stephist(
    μ_p,
    xlims = (1, 3),
    xlabel = "μ posterior",
    leg = false,
    lw = 2,
    normalize = true,
)
bp = stephist(
    σ_p,
    xlims = (0, 0.3),
    xlabel = "σ posterior",
    leg = false,
    lw = 2,
    normalize = true,
)
plot(a, ap, b, bp)
savefig("inference.svg");
nothing; # hide

# ![inference_plot](inference.svg)
#
# we can see that the algorithm has correctly inferred both parameters, this exact recipe will work for much more complicated models and simulations, with some tuning.


text:
using ColorTypes, FixedPointNumbers
using Base.Test

if VERSION >= v"0.5.0-dev+2396"
    macro inferred5(ex)
        Expr(:macrocall, Symbol("@inferred"), esc(ex))
    end
else
    macro inferred5(ex)
        esc(ex)
    end
end

if VERSION >= v"0.5.0"
    if isdefined(Core, :UnionAll)
      #@test isempty(detect_ambiguities(ColorTypes, Base, Core, allow_bottom=false))
    else
      @test isempty(detect_ambiguities(ColorTypes, Base, Core))
    end
end

# Support pre- and post- julia #20288
if VERSION >= v"0.6.0-dev.2505"
    tformat(x...) = join(string.(x), ", ")
else
    tformat(x...) = join(map(string, x), ",")
end

@test ColorTypes.to_top(AGray32(.8)) == ColorTypes.Colorant{FixedPointNumbers.Normed{UInt8,8},2}
@test @inferred(eltype(Color{N0f8})) == N0f8
@test @inferred(eltype(RGB{Float32})) == Float32
@test @inferred(eltype(RGBA{Float64})) == Float64
# @test eltype(RGB) == TypeVar(:T, Fractional)
@inferred(eltype(RGB))      # just test that it doesn't error

@test length(RGB)    == 3
@test length(RGB1)   == 3
@test length(Gray)   == 1
@test length(ARGB)   == 4
@test length(RGB24)  == 3
@test length(ARGB32) == 4
@test length(AGray{Float32}) == 2

@test @inferred(color_type(RGB{N0f8})) == RGB{N0f8}
@test @inferred(color_type(RGB)) == RGB
@test @inferred(color_type(RGBA{Float32})) == RGB{Float32}
@test @inferred(color_type(GrayA{N0f8})) == Gray{N0f8}
@test @inferred(color_type(RGBA)) == RGB
@test @inferred(color_type(RGB24) ) == RGB24
@test @inferred(color_type(ARGB32)) == RGB24
@test @inferred(color_type(TransparentColor{RGB})) == RGB
@test @inferred(color_type(TransparentColor{RGB,Float64})) == RGB
@test @inferred(color_type(TransparentColor{RGB{Float64},Float64})) == RGB{Float64}
@test color_type(TransparentColor) <: Color
@test Color <: color_type(TransparentColor)
@test_throws MethodError color_type(Colorant{N0f8})

@test @inferred(base_color_type(RGBA{Float32})) == RGB
@test @inferred(base_color_type(ARGB{Float32})) == RGB
@test @inferred(base_color_type(BGR{N0f8})      ) == BGR
@test @inferred(base_color_type(HSV) ) == HSV
@test @inferred(base_color_type(HSVA)) == HSV
@test @inferred(base_color_type(TransparentColor{RGB{Float64},Float64})) == RGB

@test @inferred(base_colorant_type(RGBA{Float32})) == RGBA
@test @inferred(base_colorant_type(ARGB{Float32})) == ARGB
@test @inferred(base_colorant_type(BGR{N0f8})      ) == BGR
@test @inferred(base_colorant_type(HSV) ) == HSV
@test @inferred(base_colorant_type(HSVA)) == HSVA

@test @inferred(ccolor(Colorant{N0f8,3}, BGR{N0f8})) == BGR{N0f8}

@test @inferred(ccolor(RGB{Float32}, HSV{Float32})) == RGB{Float32}
@test @inferred(ccolor(RGB{N0f8},      HSV{Float32})) == RGB{N0f8}
@test @inferred(ccolor(RGB,          HSV{Float32})) == RGB{Float32}
@test @inferred(ccolor(ARGB{Float32}, HSV{Float32})) == ARGB{Float32}
@test @inferred(ccolor(ARGB{N0f8},      HSV{Float32})) == ARGB{N0f8}
@test @inferred(ccolor(ARGB,          HSV{Float32})) == ARGB{Float32}

@test @inferred(ccolor(Gray{N0f8}, Bool)) === Gray{N0f8}
@test @inferred(ccolor(Gray,     Bool)) === Gray{Bool}
@test @inferred(ccolor(Gray,     Int))  === Gray{N0f8}
# This tests the same thing as the last, but in a user-observable way
a = Array{Gray}(1)
a[1] = Gray(0)
a[1] = 1
@test a[1] === Gray(1)

@test @inferred(ccolor(RGB,  RGB))  == RGB
@test @inferred(ccolor(Gray, Gray)) == Gray

# Traits for instances (and their constructors)
@test @inferred(eltype(RGB{N0f8}(1,0,0))) == N0f8
@test @inferred(eltype(RGB(1.0,0,0))) == Float64
@test @inferred(eltype(ARGB(1.0,0.8,0.6,0.4))) == Float64
@test @inferred(eltype(RGBA{Float32}(1.0,0.8,0.6,0.4))) == Float32
@test @inferred(eltype(RGB(0x01,0x00,0x00))) == N0f8

@test length(ARGB(1.0,0.8,0.6,0.4)) == 4

@test @inferred(color_type(RGB{N0f8}(1,0,0))) == RGB{N0f8}
@test @inferred(color_type(ARGB(1.0,0.8,0.6,0.4))) == RGB{Float64}
@test @inferred(color_type(RGBA{Float32}(1.0,0.8,0.6,0.4))) == RGB{Float32}

@test @inferred(base_color_type(RGB{N0f8}(1,0,0))) == RGB
@test @inferred(base_color_type(ARGB(1.0,0.8,0.6,0.4))) == RGB
@test @inferred(base_color_type(RGBA{Float32}(1.0,0.8,0.6,0.4))) == RGB

@test @inferred(base_colorant_type(RGB{N0f8}(1,0,0))) == RGB
@test @inferred(base_colorant_type(ARGB(1.0,0.8,0.6,0.4))) == ARGB
@test @inferred(base_colorant_type(RGBA{Float32}(1.0,0.8,0.6,0.4))) == RGBA

@test N0f8 <: ColorTypes.eltypes_supported(RGB(1,0,0))

# Constructors
for val in (0.2, 0.2f0, N0f8(0.2), N4f12(0.2), N0f16(0.2),
            Gray{N0f8}(0.2), Gray{N4f12}(0.2), Gray24(0.2))
    @test isa(Gray(val), Gray)
    @test Gray{N0f8}(val) === Gray{N0f8}(0.2)
    @test Gray{N0f16}(val) === Gray{N0f16}(0.2)
    @test Gray24(val) === Gray24(0.2)
    @test AGray32(val) === AGray32(0.2, 1)
    @test AGray32(val, 0.8) === AGray32(0.2, 0.8)
end
for val in (1.2, 1.2f0, N4f12(1.2), Gray{N4f12}(1.2), 2)
    !isa(val, Int) && @test isa(Gray(val), Gray)
    @test_throws ArgumentError Gray{N0f8}(val)
    @test_throws ArgumentError Gray{N0f16}(val)
    @test_throws ArgumentError Gray24(val) == Gray24(0.2)
    @test_throws ArgumentError GrayA{N0f8}(val)
    @test_throws ArgumentError AGray{N0f8}(val)
    @test_throws ArgumentError AGray32(val)
    @test_throws ArgumentError AGray32(val, 0.8)
end
@test eltype(Gray()) == N0f8
@test Gray(Gray()) == Gray()  # no StackOverflowError
@test eltype(broadcast(Gray, rand(5))) == Gray{Float64}
@test eltype(broadcast(Gray, rand(Float32,5))) == Gray{Float32}

for C in ColorTypes.parametric3
    @test eltype(C{Float32}) == Float32
    et = (C <: AbstractRGB) ? N0f8 : Float32
    @test eltype(C(1,0,0)) == et
    @test color_type(C(1,0,0)) == C{et}
    @test color_type(C) == C
    @test color_type(C{Float32}) == C{Float32}
    @test eltype(C{Float32}(1,0,0)) == Float32
    @test C(C()) == C()  # no StackOverflowError
end

# Specifically test the AbstractRGB types
# This checks that the constructor order is the same, even if the
# storage order is not
for C in filter(T -> T <: AbstractRGB, ColorTypes.color3types)
    c = C(1, 0.5, 0)
    C == RGB24 && continue
    @test red(c)   == c.r == 1
    @test green(c) == c.g == 0.5
    @test blue(c)  == c.b == 0
end
# Check various input types and values
for C in filter(T -> T <: AbstractRGB, ColorTypes.color3types)
    for val1 in (0.2, 0.2f0, N0f8(0.2), N4f12(0.2), N0f16(0.2))
        for val2 in (0.2, 0.2f0, N0f8(0.2), N4f12(0.2), N0f16(0.2))
            c = C(val1,val2,val1)
            @test isa(c, C)
            @test isa(alphacolor(C)(val1,val2,val1), alphacolor(C))
            @test isa(alphacolor(C)(val1,val2,val1,0.2), alphacolor(C))
            @test alphacolor(C)(c) === alphacolor(C)(val1,val2,val1,convert(eltype(c), 1))
            @test alphacolor(C)(c, 0.2) === alphacolor(C)(val1,val2,val1,convert(eltype(c), 0.2))
            @test alphacolor(c, 0.2) === alphacolor(C)(val1,val2,val1,convert(eltype(c), 0.2))
            if C !== RGB24
                @test isa(coloralpha(C)(val1,val2,val1), coloralpha(C))
                @test isa(coloralpha(C)(val1,val2,val1,0.2), coloralpha(C))
                @test C{N0f8}(val1,val2,val1) === C{N0f8}(0.2,0.2,0.2)
                @test C{N0f16}(val1,val2,val1) === C{N0f16}(0.2,0.2,0.2)
                @test alphacolor(C){N0f8}(val1,val2,val1) === alphacolor(C){N0f8}(0.2,0.2,0.2,1)
                @test alphacolor(C){N0f8}(val1,val2,val1,0.2) === alphacolor(C){N0f8}(0.2,0.2,0.2,0.2)
                @test coloralpha(C){N0f8}(val1,val2,val1) === coloralpha(C){N0f8}(0.2,0.2,0.2,1)
                @test coloralpha(C){N0f8}(val1,val2,val1,0.2) === coloralpha(C){N0f8}(0.2,0.2,0.2,0.2)
                @test coloralpha(C)(c) === coloralpha(C){eltype(c)}(val1,val2,val1,1)
                @test coloralpha(C)(c, 0.2) === coloralpha(C)(val1,val2,val1,convert(eltype(c), 0.2))
                @test coloralpha(c, 0.2) === coloralpha(C)(val1,val2,val1,convert(eltype(c), 0.2))
            end
        end
    end
    @test isa(C(1,0,1), C)
    @test C() === C(0,0,0)
    @test isa(alphacolor(C)(1,0,1), alphacolor(C))
    @test alphacolor(C)() === alphacolor(C)(0,0,0,1)
    if C != RGB24
        @test C(1,0,1) === C{N0f8}(1,0,1)
        @test alphacolor(C)(1,0,1) === alphacolor(C){N0f8}(1,0,1,1)
        @test alphacolor(C)(1,0,1,0.2) === alphacolor(C){Float64}(1,0,1,0.2)
        @test coloralpha(C)(1,0,1) === coloralpha(C){N0f8}(1,0,1,1)
        @test coloralpha(C)(1,0,1,0.2) === coloralpha(C){Float64}(1,0,1,0.2)
        @test coloralpha(C)() === coloralpha(C)(0,0,0,1)
    end
    for val in (1.2, 1.2f0, N4f12(1.2), 2)
        if C !== RGB24
            @test_throws ArgumentError C{N0f8}(val,val,val)
            @test_throws ArgumentError C{N0f16}(val,val,val)
            @test_throws ArgumentError alphacolor(C){N0f8}(val,val,val)
            @test_throws ArgumentError coloralpha(C){N0f8}(val,val,val)
            if val != 2
                @test isa(C(val,val,val), C)
                c = C(val,val,val)
                @test_throws ArgumentError alphacolor(C){N0f8}(c)
                @test_throws ArgumentError coloralpha(C){N0f8}(c)
                if C != RGB1 && C != RGB4
                    @test_throws ArgumentError alphacolor(C){N0f8}(c, 0.2)
                    @test_throws ArgumentError coloralpha(C){N0f8}(c, 0.2)
                end
            end
        end
        if C == RGB24 || isa(val, Int)
            @test_throws ArgumentError C(val,val,val)
            @test_throws ArgumentError alphacolor(C)(val,val,val)
        end
    end
end
if VERSION >= v"0.5.0"
    ret = @test_throws ArgumentError RGB(255, 17, 48)
    @test contains(ret.value.msg, tformat(255,17,48))
    @test contains(ret.value.msg, "0-255")
    ret = @test_throws ArgumentError RGB(256, 17, 48)
    @test contains(ret.value.msg, tformat(256,17,48))
    @test !contains(ret.value.msg, "0-255")
end

c = Gray(0.8)
@test gray(c) == real(c) == 0.8
@test gray(0.8) == 0.8
c = convert(Gray, 0.8)
@test c === Gray{Float64}(0.8)

ac = convert(AGray, c)
@test ac === AGray{Float64}(0.8, 1.0)

c = AGray(0.8)
@test gray(c) == 0.8
@test color(c) == Gray(0.8)

c = convert(Gray, true)
@test c === Gray{Bool}(true)
@test gray(c) === true
@test gray(false) === false

# Transparency
@test alphacolor(Gray24(.2), .8) == AGray32(.2,.8)
@test alphacolor(RGB24(1,0,0), .8) == ARGB32(1,0,0,.8)
@test alphacolor(RGB(1,0,0), .8) == ARGB{N0f8}(1,0,0,.8)
@test coloralpha(RGB(1,0,0), .8) == RGBA{N0f8}(1,0,0,.8)
@test alphacolor(RGBA(1,0,0,.8)) == ARGB{Float64}(1,0,0,.8)
@test coloralpha(ARGB(1,0,0,.8)) == RGBA{Float64}(1,0,0,.8)
@test alphacolor(RGBA(1,0,0,.8)) == ARGB{Float64}(1,0,0,.8)
@test coloralpha(ARGB(1,0,0,.8)) == RGBA{Float64}(1,0,0,.8)
for C in setdiff(ColorTypes.parametric3, [RGB1,RGB4])
    for A in (alphacolor(C), coloralpha(C))
        @test eltype(A{Float32}) == Float32
        @test color_type(A) == C
        @test color_type(A{Float32}) == C{Float32}
        @test eltype(A(1,0.8,0.6,0.4)) == Float64
        c = A{Float64}(1,0.8,0.6,0.4)
        @test color_type(c) == C{Float64}
        cc = color(c)
        @test cc == C{Float64}(1,0.8,0.6)
        if VERSION >= v"0.4.0-dev"
            @test A(cc) == A{Float64}(1,0.8,0.6,1)
            @test A(cc, 0.4)  == c
            @test A(cc, 0x01) == A{Float64}(1,0.8,0.6,1)
            @test A{Float32}(cc, 0x01) == A{Float32}(1,0.8,0.6,1)
            @test C(c         ) == C{Float64}(1,0.8,0.6)
            @test C{Float32}(c) == C{Float32}(1,0.8,0.6)
        end
        @test convert(A, cc) == A{Float64}(1,0.8,0.6,1)
        @test convert(A, cc, 0.4)  == c
        @test convert(A, cc, 0x01) == A{Float64}(1,0.8,0.6,1)
        @test convert(A{Float32}, cc, 0x01) == A{Float32}(1,0.8,0.6,1)
        @test convert(C,          c) == C{Float64}(1,0.8,0.6)
        @test convert(C{Float32}, c) == C{Float32}(1,0.8,0.6)
    end
end
ac = reinterpret(ARGB32, rand(UInt32))
@test convert(ARGB32, ac) == ac
c = convert(RGB24, ac)
@test convert(RGB24, c) == c

crgb   = convert(RGB, c)
acargb = convert(ARGB, ac)
@test convert(Colorant,       crgb) === crgb
@test convert(Colorant{N0f8},   crgb) === crgb
@test convert(Colorant{Float32},       crgb) === convert(RGB{Float32}, crgb)
@test convert(Color,          crgb) === crgb
@test convert(Color{N0f8},      crgb) === crgb
@test convert(Color{Float32}, crgb) === convert(RGB{Float32}, crgb)
@test_throws ErrorException convert(TransparentColor, crgb)
@test_throws ErrorException convert(TransparentColor{RGB{N0f8}}, crgb)
@test_throws ErrorException convert(TransparentColor{RGB{N0f8},N0f8}, crgb)
@test_throws ErrorException convert(TransparentColor{RGB{N0f8},N0f8,2}, crgb)
@test convert(AlphaColor, crgb) === alphacolor(crgb)
@test convert(AlphaColor{RGB{N0f8}}, crgb) === alphacolor(crgb)
@test convert(AlphaColor{RGB{N0f8},N0f8}, crgb) === alphacolor(crgb)
@test convert(AlphaColor{RGB{N0f8},N0f8,4}, crgb) === alphacolor(crgb)
@test convert(ColorAlpha, crgb) === coloralpha(crgb)
@test convert(ColorAlpha{RGB{N0f8}}, crgb) === coloralpha(crgb)
@test convert(ColorAlpha{RGB{N0f8},N0f8}, crgb) === coloralpha(crgb)
@test convert(ColorAlpha{RGB{N0f8},N0f8,4}, crgb) === coloralpha(crgb)

@test convert(Colorant, acargb) === acargb
@test convert(Colorant{N0f8},   acargb) === acargb
@test convert(Colorant{N0f8,3}, acargb) === crgb
@test convert(TransparentColor, acargb) == acargb
@test convert(Color, acargb) == crgb

h = N0f8(0.5)
@test Gray24(0.5) == Gray24(h)
@test convert(Gray24, 0.5) == Gray24(h)
@test convert(AGray, Gray24(h)) === AGray{N0f8}(h, 1)
@test convert(AGray, Gray24(h), 0.8)  === AGray{N0f8}(h, 0.8)
@test convert(AGray, AGray32(h, 0.8)) === AGray{N0f8}(h, 0.8)
@test AGray32(0.5) == AGray32(h, 1)
@test convert(AGray32, 0.5, 0.8) == AGray32(h, N0f8(0.8))
@test RGB24(0.5) == RGB24(h, h, h)
@test convert(RGB24, 0.5) == RGB24(h, h, h)
@test ARGB32(0.5) == ARGB32(h, h, h, 1)

@test red(c)   == red(ac)
@test green(c) == green(ac)
@test blue(c)  == blue(ac)
ac2 = convert(ARGB32, c)
@test reinterpret(UInt32, ac2) == (c.color | 0xff000000)
@test color(c) == c
@test color(ac) == c
@test alpha(c) == N0f8(1)
@test alpha(ac) == N0f8(ac.color>>24, 0)
@test alpha(ac2) == N0f8(1)
@test reinterpret(UInt32, reinterpret(RGB24,  0xff020304)) == 0xff020304
@test reinterpret(UInt32, reinterpret(ARGB32, 0x01020304)) == 0x01020304
ac3 = convert(RGBA, ac)
@test convert(RGB24, ac3) == c
ac4 = AGray32(.2,.8)
@test alpha(ac4) == .8
@test gray(ac4) == .2

for C in filter(T -> T <: AbstractRGB, ColorTypes.color3types)
    rgb = convert(C, c)
    C == RGB24 && continue
    @test ccolor(Gray24, C) == Gray24
    @test ccolor(AGray32, C) == AGray32
    @test convert(AbstractRGB, c) == c
    @test convert(AbstractRGB{Float64}, rgb) === convert(C{Float64}, c)
    argb = convert(alphacolor(C), ac)
    rgba = convert(coloralpha(C), ac)
    @test rgb.r == red(c)
    @test rgb.g == green(c)
    @test rgb.b == blue(c)
    @test argb.alpha == alpha(ac)
    @test argb.r == red(ac)
    @test argb.g == green(ac)
    @test argb.b == blue(ac)
    @test rgba.alpha == alpha(ac)
    @test rgba.r == red(ac)
    @test rgba.g == green(ac)
    @test rgba.b == blue(ac)
end

@test_throws ErrorException convert(HSV, RGB(1,0,1))
@test_throws ErrorException convert(AHSV, RGB(1,0,1), 0.5)

@test convert(Float64, Gray(.3)) === .3
x = N0f8(0.3)
@test convert(N0f8, Gray24(0.3)) === x
@test convert(GrayA{N0f8}, .2) == GrayA{N0f8}(.2)
@test convert(AGray{N0f8}, .2) == AGray{N0f8}(.2)
@test Gray{N0f8}(0.37).val           == N0f8(0.37)
@test convert(Gray{N0f8}, 0.37).val  == N0f8(0.37)
@test reinterpret(UInt32, Gray24(reinterpret(N0f8, 0x0d)))           == 0x000d0d0d
@test reinterpret(UInt32, convert(Gray24, reinterpret(N0f8, 0x0d)))  == 0x000d0d0d
@test reinterpret(UInt32, AGray32(reinterpret(N0f8, 0x0d)))          == 0xff0d0d0d
@test reinterpret(UInt32, convert(AGray32, reinterpret(N0f8, 0x0d))) == 0xff0d0d0d
@test reinterpret(UInt32, AGray32(reinterpret(N0f8, 0x0d), reinterpret(N0f8, 0x80))) == 0x800d0d0d
@test convert(Gray{N0f16}, Gray24(reinterpret(N0f8, 0x0d))) == Gray{N0f16}(0.05098)

@test promote(Gray{N0f8}(0.2), Gray24(0.3)) === (Gray{N0f8}(0.2), Gray{N0f8}(0.3))
@test promote(Gray(0.2f0), Gray24(0.3)) === (Gray{Float32}(0.2), Gray{Float32}(N0f8(0.3)))
@test promote(RGB{N0f8}(0.2,0.3,0.4), RGB24(0.3,0.8,0.1)) === (RGB{N0f8}(0.2,0.3,0.4), RGB{N0f8}(0.3,0.8,0.1))
@test promote(RGB{Float32}(0.2,0.3,0.4), RGB24(0.3,0.8,0.1)) === (RGB{Float32}(0.2,0.3,0.4), RGB{Float32}(N0f8(0.3),N0f8(0.8),N0f8(0.1)))

iob = IOBuffer()
cf = RGB{Float32}(0.32218,0.14983,0.87819)
c  = convert(RGB{N0f8}, cf)
show(iob, c)
@test String(take!(iob)) == "RGB{N0f8}(0.322,0.149,0.878)"
c = RGB{N0f16}(0.32218,0.14983,0.87819)
show(iob, c)
@test String(take!(iob)) == "RGB{N0f16}(0.32218,0.14983,0.87819)"
c = RGBA{N0f8}(0.32218,0.14983,0.87819,0.99241)
show(iob, c)
@test String(take!(iob)) == "RGBA{N0f8}(0.322,0.149,0.878,0.992)"
showcompact(iob, c)
@test String(take!(iob)) == "RGBA{N0f8}(0.322,0.149,0.878,0.992)"
show(iob, cf)
@test String(take!(iob)) == "RGB{Float32}(0.32218f0,0.14983f0,0.87819f0)"
showcompact(iob, cf)
@test String(take!(iob)) == "RGB{Float32}(0.32218,0.14983,0.87819)"

@test one(Gray{N0f8}) == Gray{N0f8}(1)
@test zero(Gray{N0f8}) == Gray{N0f8}(0)

c = Gray(0.8)
@test c == 0.8
show(iob, c)
@test String(take!(iob)) == "Gray{Float64}(0.8)"
show(iob, AGray(0.8))
@test String(take!(iob)) == "AGray{Float64}(0.8,1.0)"
show(iob, AGray32(0.8))
@test String(take!(iob)) == "AGray32{N0f8}(0.8,1.0)"

# if the test below fails, please extend the list of types at the call to
# make_alpha in types.jl (this is the price of making that list explicit)
@test Set(ColorTypes.ctypes) ==
Set([DIN99d, DIN99o, DIN99, HSI, HSL, HSV, LCHab, LCHuv,
     LMS, Lab, Luv, XYZ, YCbCr, YIQ, xyY, BGR, RGB, Gray])

## operations
for T in (Gray{N0f8}, Gray{N2f6}, Gray{N0f16}, Gray{N2f14}, Gray{N0f32}, Gray{N2f30},
          Gray{Float16}, Gray{Float32}, Gray{Float64},
          RGB{N0f8}, RGB{N2f6}, RGB{N0f16}, RGB{N2f14}, RGB{N0f32}, RGB{N2f30},
          RGB{Float16}, RGB{Float32}, RGB{Float64},
          AGray{Float32}, GrayA{Float64},
          RGBA{Float32}, ARGB{N0f16}, RGB1{N0f8}, RGB4{Float64},
          BGR{Float16}, ABGR{N0f32}, BGRA{N2f14},
          Gray, AGray, GrayA,
          RGB, ARGB, RGBA, BGR, ABGR, BGRA, RGB1, RGB4,
          HSV, HSL, Lab, LCHab, YIQ)
    a = rand(T)
    @test all(x->x[2]<=getfield(a,x[1])<=x[3],
                zip(ColorTypes.colorfields(T),gamutmin(T),gamutmax(T)))
    @test isa(a, T)
    a = rand(T, (3, 5))
    for el in a
        @test all(x->x[2]<=getfield(el,x[1])<=x[3],
                    zip(ColorTypes.colorfields(T),gamutmin(T),gamutmax(T)))
    end
    @test eltype(a) <: T
    @test size(a) == (3,5)
end
for T in (Gray24, AGray32)
    a = rand(T)
    b = a.color
    @test b&0xff == (b>>8)&0xff == (b>>16)&0xff
    @test isa(a, T)
    a = rand(T,3,5)
    for el in a
        b = el.color
        @test b&0xff == (b>>8)&0xff == (b>>16)&0xff
    end
    @test eltype(a) <: T
    @test size(a) == (3,5)
end

a = [BGR(1,0,0)]
@test eltype(broadcast(RGB, a)) == RGB{N0f8}

# colorfields
@test ColorTypes.colorfields(AGray32(.2)) == (:color,:alpha)
@test ColorTypes.colorfields(Gray) == (:val,)
@test ColorTypes.colorfields(RGB1) == (:r, :g, :b)
@test ColorTypes.colorfields(RGB4) == (:r, :g, :b)
@test ColorTypes.colorfields(BGR) == (:r, :g, :b)

# UInt32 comparison
@test reinterpret(UInt32, Gray24()) == 0x00000000
@test reinterpret(UInt32, Gray24(.2)) == 0x00333333
@test reinterpret(UInt32, Gray24(reinterpret(N0f8, 0x23))) == 0x00232323
@test reinterpret(UInt32, AGray32()) == 0xff000000
@test reinterpret(UInt32, AGray32(.2)) == 0xff333333
@test reinterpret(UInt32, convert(AGray32, .2, 0.)) == 0x00333333
@test reinterpret(UInt32, AGray32(reinterpret(N0f8, 0x23))) == 0xff232323
@test reinterpret(UInt32, RGB24()) == 0x00000000
@test reinterpret(UInt32, ARGB32()) == 0xff000000
@test reinterpret(UInt32, ARGB32(1,.2,.3)) == 0xffff334c

@test @inferred5(mapc(sqrt, Gray{N0f8}(0.04))) == Gray(sqrt(N0f8(0.04)))
@test @inferred5(mapc(sqrt, AGray{N0f8}(0.04, 0.4))) == AGray(sqrt(N0f8(0.04)), sqrt(N0f8(0.4)))
@test @inferred5(mapc(sqrt, GrayA{N0f8}(0.04, 0.4))) == GrayA(sqrt(N0f8(0.04)), sqrt(N0f8(0.4)))
@test @inferred5(mapc(x->2x, RGB{N0f8}(0.04,0.2,0.3))) == RGB(map(x->2*N0f8(x), (0.04,0.2,0.3))...)
@test @inferred5(mapc(sqrt, RGBA{N0f8}(0.04,0.2,0.3,0.7))) == RGBA(map(x->sqrt(N0f8(x)), (0.04,0.2,0.3,0.7))...)
@test @inferred5(mapc(x->1.5f0x, RGBA{N0f8}(0.04,0.2,0.3,0.4))) == RGBA(map(x->1.5f0*N0f8(x), (0.04,0.2,0.3,0.4))...)

@test @inferred5(mapc(max, Gray{N0f8}(0.2), Gray{N0f8}(0.3))) == Gray{N0f8}(0.3)
@test @inferred5(mapc(-, AGray{Float32}(0.3), AGray{Float32}(0.2))) == AGray{Float32}(0.3f0-0.2f0,0.0)
@test @inferred5(mapc(min, RGB{N0f8}(0.2,0.8,0.7), RGB{N0f8}(0.5,0.2,0.99))) == RGB{N0f8}(0.2,0.2,0.7)
@test @inferred5(mapc(+, RGBA{N0f8}(0.2,0.8,0.7,0.3), RGBA{Float32}(0.5,0.2,0.99,0.5))) == RGBA(0.5f0+N0f8(0.2),0.2f0+N0f8(0.8),0.99f0+N0f8(0.7),0.5f0+N0f8(0.3))
@test @inferred5(mapc(+, HSVA(0.1,0.8,0.3,0.5), HSVA(0.5,0.5,0.5,0.3))) == HSVA(0.1+0.5,0.8+0.5,0.3+0.5,0.5+0.3)
@test_throws ArgumentError mapc(min, RGB{N0f8}(0.2,0.8,0.7), BGR{N0f8}(0.5,0.2,0.99))
@test @inferred(mapc(abs, -2)) == 2

@test @inferred(reducec(+, 0.0, Gray(0.3))) === 0.3
@test @inferred(reducec(+, 1.0, Gray(0.3))) === 1.3
@test @inferred(reducec(+, 0, Gray(0.3))) === 0.3
@test @inferred(reducec(+, 0.0, AGray(0.3, 0.8))) === 0.3 + 0.8
@test @inferred(reducec(+, 0.0, RGB(0.3, 0.8, 0.5))) === (0.3 + 0.8) + 0.5
@test @inferred(reducec(+, 0.0, RGBA(0.3, 0.8, 0.5, 0.7))) === ((0.3 + 0.8) + 0.5) + 0.7
@test @inferred(reducec(&, true, Gray(true)))
@test !(@inferred(reducec(&, false, Gray(true))))
@test !(@inferred(reducec(&, true, Gray(false))))

@test @inferred(reducec(+, 0.0, 0.3)) === 0.3
@test @inferred(reducec(+, 0, 0.3)) === 0.3
@test @inferred(reducec(&, true, true))
@test !(@inferred(reducec(&, false, true)))
@test !(@inferred(reducec(&, true, false)))

@test @inferred(mapreducec(x->x^2, +, 0.0, Gray(0.3))) === 0.3^2
@test @inferred(mapreducec(x->x^2, +, 1.0, Gray(0.3))) === 1 + 0.3^2
@test @inferred(mapreducec(x->x^2, +, 0, Gray(0.3))) === 0.3^2
@test @inferred(mapreducec(x->x^2, +, 0.0, AGray(0.3, 0.8))) === 0.3^2 + 0.8^2
@test @inferred(mapreducec(x->x^2, +, 0.0, RGB(0.3, 0.8, 0.5))) === (0.3^2 + 0.8^2) + 0.5^2
@test @inferred(mapreducec(x->x^2, +, 0.0, RGBA(0.3, 0.8, 0.5, 0.7))) === ((0.3^2 + 0.8^2) + 0.5^2) + 0.7^2
@test !(@inferred(mapreducec(x->!x, &, true, Gray(true))))
@test !(@inferred(mapreducec(x->!x, &, false, Gray(true))))
@test @inferred(mapreducec(x->!x, &, true, Gray(false)))
@test !@inferred(mapreducec(x->!x, &, false, Gray(false)))

@test @inferred(mapreducec(x->x^2, +, 0.0, 0.3)) === 0.3^2
@test @inferred(mapreducec(x->x^2, +, 1.0, 0.3)) === 1 + 0.3^2
@test @inferred(mapreducec(x->x^2, +, 0, 0.3)) === 0.3^2
@test !(@inferred(mapreducec(x->!x, &, true, true)))
@test !(@inferred(mapreducec(x->!x, &, false, true)))
@test @inferred(mapreducec(x->!x, &, true, false))
@test !@inferred(mapreducec(x->!x, &, false, false))

# issue #52
@test AGray{BigFloat}(0.5,0.25) == AGray{BigFloat}(0.5,0.25)
@test RGBA{BigFloat}(0.5, 0.25, 0.5, 0.5) == RGBA{BigFloat}(0.5, 0.25, 0.5, 0.5)

for (a, b) in ((Gray(1.0), Gray(1)),
               (GrayA(0.8, 0.6), AGray(0.8, 0.6)),
               (RGB(1, 0.5, 0), BGR(1, 0.5, 0)),
               (RGBA(1, 0.5, 0, 0.8), ABGR(1, 0.5, 0, 0.8)))
    @test a == b
    @test hash(a) == hash(b)
end
for (a, b) in ((RGB(1, 0.5, 0), RGBA(1, 0.5, 0, 0.9)),)
    @test a != b
    @test hash(a) != hash(b)
end
# It's not obvious whether we want these to compare as equal, but
# whatever happens, you want hashing and equality-testing to yield the
# same result
for (a, b) in ((RGB(1, 0.5, 0), RGBA(1, 0.5, 0, 1)),)
    @test (a == b) == (hash(a) == hash(b))
end

### Test deprecations
mktemp() do tmpfile, io
    redirect_stderr(io) do
        @test convert(UInt32, RGB24(1,1,1))    == 0x00ffffff
        @test convert(UInt32, ARGB32(1,1,1,0)) == 0x00ffffff
        @test convert(UInt32, Gray24(1))    == 0x00ffffff
        @test convert(UInt32, AGray32(1,0)) == 0x00ffffff
        @test RGB24(1,1,1)    == 0x00ffffff
        @test ARGB32(1,1,1,0) == 0x00ffffff
        @test Gray24(1)    == 0x00ffffff
        @test AGray32(1,0) == 0x00ffffff
        @test RGB24(0x00ffffff)  === RGB24(1,1,1)
        @test ARGB32(0x00ffffff) === ARGB32(1,1,1,0)
        @test Gray24(0x00ffffff)  === Gray24(1)
        @test AGray32(0x00ffffff) === AGray32(1,0)
        @test RGB24[0x00000000] == [RGB24(0)]
        @test RGB24[0x00000000,0x00808080] == [RGB24(0), RGB24(0.5)]
        @test RGB24[0x00000000,0x00808080,0x00ffffff] == [RGB24(0), RGB24(0.5), RGB24(1)]
        @test RGB24[0x00000000,0x00808080,0x00ffffff,0x000000ff] == [RGB24(0), RGB24(0.5), RGB24(1), RGB24(0,0,1)]
    end
    close(io)
    Base.JLOptions().depwarn==1 && @test sum(x->contains(x, "WARNING"), readlines(tmpfile)) == 16
end

nothing


text:
@testset "toggle docs" begin
    # basic
    let doc = Atom.docs("push!")
        @test !doc[:error]
    end

    # context module
    let doc = Atom.docs("handle", "Atom")
        @test !doc[:error]
    end

    # keyword
    let doc = Atom.docs("begin", "Atom")
        @test !doc[:error]
    end
end


text:
# Join works on collections of collections (e.g. a table is a collection of
# rows).

# This seems to be a generically useful primitive on indexable collections,
# and is a useful definition for "natural" joins
function overlaps(a, b)
    for i in keys(a)
        if haskey(b, i)
            if a[i] != b[i]
                return false
            end
        end
    end
    return true
end

# TODO simplified signatures

innerjoin(left, right) = innerjoin(identity, identity, left, right)
innerjoin(lkey, rkey, left, right) = innerjoin(lkey, rkey, merge, left, right)
innerjoin(lkey, rkey, f, left, right) = innerjoin(lkey, rkey, f, overlaps, left, right)

const ⨝ = innerjoin

"""
    innerjoin(lkey, rkey, f, comparison, left, right)

Performs a relational-style join operation between iterables `left` and `right`, returning
a collection of elements `f(l, r)` for which `comparison(lkey(l), rkey(r))` is `true` where
`l ∈ left`, `r ∈ right.`

# Example

```jldoctest
julia> innerjoin(iseven, iseven, tuple, ==, [1,2,3,4], [0,1,2])
6-element Array{Tuple{Int64,Int64},1}:
 (1, 1)
 (2, 0)
 (2, 2)
 (3, 1)
 (4, 0)
 (4, 2)
```
"""
function innerjoin(lkey, rkey, f, comparison, left, right)
    # The O(length(left)*length(right)) generic method when nothing about `comparison` is known

    # TODO Do this inference-free, like comprehensions...
    T = promote_op(f, eltype(left), eltype(right))
    out = T[]
    for a ∈ left
        for b ∈ right
            if comparison(lkey(a), rkey(b))
                push!(out, f(a, b))
            end
        end
    end
    return out
end

function innerjoin(lkey, rkey, f, ::typeof(isequal), left, right)
    # isequal heralds a hash-based approach, roughly O(length(left) * log(length(right)))

    # TODO Do this inference-free, like comprehensions...
    T = promote_op(f, eltype(left), eltype(right))
    K = promote_op(rkey, eltype(right))
    V = eltype(right)
    dict = Dict{K,Vector{V}}() # maybe a different stategy if right is unique
    for b ∈ right
        key = rkey(b)
        push!(get!(()->Vector{V}(), dict, key), b)
    end

    out = T[]
    for a ∈ left
        key = lkey(a)
        dict_index = Base.ht_keyindex(dict, key)
        if dict_index > 0 # -1 if key not found
            for b ∈ dict.vals[dict_index]
                push!(out, f(a, b))
            end
        end
    end
    return out
end

# TODO more specialized methods for comparisons: ==, <, isless, etc - via sorting strategies

# TODO perhaps a better version would be:
# function innerjoin(left, right; lkey = identity, rkey = identity, f = tuple, comparison = isequal)
#     ...
# end


text:
# ---
# title: 40. Combination Sum II
# id: problem40
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Backtracking
# link: <https://leetcode.com/problems/combination-sum-ii/description/>
# hidden: true
# ---
#
# Given a collection of candidate numbers (`candidates`) and a target number
# (`target`), find all unique combinations in `candidates` where the candidate
# numbers sum to `target`.
#
# Each number in `candidates` may only be used **once** in the combination.
#
# **Note:**  The solution set must not contain duplicate combinations.
#
#
#
# **Example 1:**
#
#
#
#     Input: candidates = [10,1,2,7,6,1,5], target = 8
#     Output:
#     [
#     [1,1,6],
#     [1,2,5],
#     [1,7],
#     [2,6]
#     ]
#
#
# **Example 2:**
#
#
#
#     Input: candidates = [2,5,2,1,2], target = 5
#     Output:
#     [
#     [1,2,2],
#     [5]
#     ]
#
#
#
#
# **Constraints:**
#
#   * `1 <= candidates.length <= 100`
#   * `1 <= candidates[i] <= 50`
#   * `1 <= target <= 30`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#########################
# Sampler I/O Interface #
#########################

##########
# Sample #
##########

mutable struct Sample
  weight :: Float64     # particle weight
  value :: Dict{Symbol,Any}
end

Base.getindex(s::Sample, v::Symbol) = getjuliatype(s, v)

getjuliatype(s::Sample, v::Symbol, cached_syms=nothing) = begin
  # NOTE: cached_syms is used to cache the filter entiries in svalue. This is helpful when the dimension of model is huge.
  if cached_syms == nothing
    # Get all keys associated with the given symbol
    syms = collect(Iterators.filter(k -> occursin(string(v)*"[", string(k)), keys(s.value)))
  else
    syms = collect((Iterators.filter(k -> occursin(string(v), string(k)), cached_syms)))
  end

  # Map to the corresponding indices part
  idx_str = map(sym -> replace(string(sym), string(v) => ""), syms)
  # Get the indexing component
  idx_comp = map(idx -> collect(Iterators.filter(str -> str != "", split(string(idx), [']','[']))), idx_str)

  # Deal with v is really a symbol, e.g. :x
  if isempty(idx_comp)
    @assert haskey(s.value, v)
    return Base.getindex(s.value, v)
  end

  # Construct container for the frist nesting layer
  dim = length(split(idx_comp[1][1], ','))
  if dim == 1
    sample = Vector(undef, length(unique(map(c -> c[1], idx_comp))))
  else
    d = max(map(c -> eval(parse(c[1])), idx_comp)...)
    sample = Array{Any, length(d)}(undef, d)
  end

  # Fill sample
  for i = 1:length(syms)
    # Get indexing
    idx = Main.eval(parse(idx_comp[i][1]))
    # Determine if nesting
    nested_dim = length(idx_comp[1]) # how many nested layers?
    if nested_dim == 1
      setindex!(sample, getindex(s.value, syms[i]), idx...)
    else  # nested case, iteratively evaluation
      v_indexed = Symbol("$v[$(idx_comp[i][1])]")
      setindex!(sample, getjuliatype(s, v_indexed, syms), idx...)
    end
  end
  sample
end

#########
# Chain #
#########

"""
    Chain(weight::Float64, value::Array{Sample})

A wrapper of output trajactory of samplers.

Example:

```julia
# Define a model
@model xxx begin
  ...
  return(mu,sigma)
end

# Run the inference engine
chain = sample(xxx, SMC(1000))

chain[:logevidence]   # show the log model evidence
chain[:mu]            # show the weighted trajactory for :mu
chain[:sigma]         # show the weighted trajactory for :sigma
mean(chain[:mu])      # find the mean of :mu
mean(chain[:sigma])   # find the mean of :sigma
```
"""
mutable struct Chain{R<:AbstractRange{Int}} <: AbstractChains
  weight  ::  Float64                 # log model evidence
  value2  ::  Array{Sample}
  value   ::  Array{Float64, 3}
  range   ::  R # TODO: Perhaps change to UnitRange?
  names   ::  Vector{String}
  chains  ::  Vector{Int}
  info    ::  Dict{Symbol,Any}
end

Chain() = Chain{AbstractRange{Int}}(0.0, Vector{Sample}(), Array{Float64, 3}(undef, 0, 0, 0), 0:0,
                Vector{String}(), Vector{Int}(), Dict{Symbol,Any}())

Chain(w::Real, s::Array{Sample}) = begin
  chn = Chain()
  chn.weight = w
  chn.value2 = deepcopy(s)

  chn = flatten!(chn)
end

flatten!(chn::Chain) = begin
  ## Flatten samples into Mamba's chain type.
  local names = Vector{Array{AbstractString}}()
  local vals  = Vector{Array}()
  for s in chn.value2
    v, n = flatten(s)
    push!(vals, v)
    push!(names, n)
  end

  # Assuming that names[i] == names[j] for all (i,j)
  vals2 = [v[i] for v in vals, i=1:length(names[1])]
  vals2 = reshape(vals2, length(vals), length(names[1]), 1)
  c = Chains(vals2, names = names[1])
  chn.value = c.value
  chn.range = c.range
  chn.names = c.names
  chn.chains = c.chains
  chn
end

# ind2sub is deprecated in Julia 1.0
ind2sub(v, i) = Tuple(CartesianIndices(v)[i])

flatten(s::Sample) = begin
  vals  = Vector{Float64}()
  names = Vector{AbstractString}()
  for (k, v) in s.value
    flatten(names, vals, string(k), v)
  end
  return vals, names
end
flatten(names, value :: Array{Float64}, k :: String, v) = begin
    if isa(v, Number)
      name = k
      push!(value, v)
      push!(names, name)
    elseif isa(v, Array)
      for i = eachindex(v)
        if isa(v[i], Number)
          name = k * string(ind2sub(size(v), i))
          name = replace(name, "(" => "[");
          name = replace(name, ",)" => "]");
          name = replace(name, ")" => "]");
          isa(v[i], Nothing) && println(v, i, v[i])
          push!(value, Float64(v[i]))
          push!(names, name)
        elseif isa(v[i], Array)
          name = k * string(ind2sub(size(v), i))
          flatten(names, value, name, v[i])
        else
          error("Unknown var type: typeof($v[i])=$(typeof(v[i]))")
        end
      end
  else
    error("Unknown var type: typeof($v)=$(typeof(v))")
  end
end

function Base.getindex(c::Chain, v::Symbol)
  # This strange implementation is mostly to keep backward compatability.
  #  Needs some refactoring a better format for storing results is available.
  if v == :logevidence
    log(c.weight)
  elseif v==:samples
    c.value2
  elseif v==:logweights
    c[:lp]
  else
    map((s)->Base.getindex(s, v), c.value2)
  end
end

Base.getindex(c::Chain, expr::Expr) = begin
  str = replace(string(expr), r"\(|\)" => "")
  @assert match(r"^\w+(\[(\d\,?)*\])+$", str) != nothing "[Turing.jl] $expr invalid for getindex(::Chain, ::Expr)"
  c[Symbol(str)]
end

function Base.vcat(c1::Chain, args::Chain...)

  names = c1.names
  all(c -> c.names == names, args) ||
    throw(ArgumentError("chain names differ"))

  chains = c1.chains
  all(c -> c.chains == chains, args) ||
    throw(ArgumentError("sets of chains differ"))

  value2 = cat(c1.value2, map(c -> c.value2, args)..., dims=1)
  Chain(0, value2)
end

save!(c::Chain, spl::Sampler, model, vi) = begin
  c.info[:spl] = spl
  c.info[:model] = model
  c.info[:vi] = deepcopy(vi)
end

resume(c::Chain, n_iter::Int) = begin
  @assert !isempty(c.info) "[Turing] cannot resume from a chain without state info"
  sample(c.info[:model],
         c.info[:spl].alg;    # this is actually not used
         resume_from=c,
         reuse_spl_n=n_iter
        )
end


text:
module Bickley

# The Bickley jet
Ψ(y, L) = - tanh(y) + 2y / L
U(y, L) = sech(y)^2 + 2 / L

# A sinusoidal tracer
C(y, L) = sin(2π * y / L)

# Slightly off-center vortical perturbations
ψ̃(x, y, ℓ, k) = exp(-(y + ℓ/10)^2 / 2ℓ^2) * cos(k * x) * cos(k * y)

# Vortical velocity fields (ũ, ṽ) = (-∂_y, +∂_x) ψ̃
ũ(x, y, ℓ, k) = + ψ̃(x, y, ℓ, k) * (k * tan(k * y) + y / ℓ^2)
ṽ(x, y, ℓ, k) = - ψ̃(x, y, ℓ, k) * k * tan(k * x)

end # module


text:
include("solution.jl")

lines = readlines("./data/input.txt")

rows = length(lines)
cols = length(lines[1])
input = zeros(Int, rows, cols)

for row in 1:rows
  for col in 1:cols
    input[row, col] = parse(Int, lines[row][col])
  end
end

print("PART 1: ")
println(power_consumption(input))

print("PART 2: ")
println(life_support_rating(input))


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

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Random, StatsBase, DataFrames, Plots
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.13"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
dists = [
    DiscreteUniform(10,20),
    Binomial(10,0.5),
    Geometric(0.5),
    NegativeBinomial(10,0.5),
    Hypergeometric(30, 40, 10),
    Poisson(5.5)]

# ╔═╡ b3515100-2155-11eb-3ae3-55afcaa336ba
begin
	df = DataFrame()
	df.distributions = dists
	df.params = params.(dists)
	df.support = ( (d)->( minimum(d),maximum(d) ) ).(dists)
	df
end

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.13"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═b3515100-2155-11eb-3ae3-55afcaa336ba
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
abstract type AbstractCommand end

struct SetObject{O <: AbstractObject} <: AbstractCommand
    object::O
    path::Path
end

struct SetTransform{T <: Transformation} <: AbstractCommand
    tform::T
    path::Path
end

struct Delete <: AbstractCommand
    path::Path
end

struct SetProperty{T} <: AbstractCommand
    path::Path
    property::String
    value::T
end

abstract type AbstractControl end

struct Button <: AbstractControl
    observer::Observable
    name::String
end

struct NumericControl{T} <: AbstractControl
    observer::Observable
    name::String
    value::T
    min::T
    max::T
end

struct SetControl <: AbstractCommand
    control::AbstractControl
end

struct SetAnimation{A <: Animation} <: AbstractCommand
    animation::A
    play::Bool
    repetitions::Int
end

SetAnimation(anim::Animation; play=true, repetitions=1) = SetAnimation(anim, play, repetitions)


text:
using DFControl, Test

DFControl.removedefault_pseudos(:test)

@test DFControl.getdefault_pseudo(:Si, :test) == nothing

setdefault_server(prevdefault)


text:
"""
Functionality for handling the defautl conent of the various Genie files (migrations, models, controllers, etc).
"""
module FileTemplates

import Inflector

"""
    newtask(module_name::String) :: String

Default content for a new Genie Toolbox task.
"""
function newtask(module_name::String) :: String
  """
  module $module_name

  \"\"\"
  Description of the task here
  \"\"\"
  function runtask()
    # Build something great
  end

  end
  """
end

"""
    newcontroller(controller_name::String) :: String

Default content for a new Genie controller.
"""
function newcontroller(controller_name::String) :: String
  """
  module $(controller_name)Controller
    # Build something great
  end
  """
end

"""
    newtest(plural_name::String, singular_name::String) :: String

Default content for a new test file.
"""
function newtest(plural_name::String, singular_name::String) :: String
  """
  using Genie, App.$(plural_name)

  ### Your tests here
  @test 1 == 1
  """
end

"""
    appmodule(path::String)

Generates a custom app module when a new app is bootstrapped.
"""
function appmodule(path::String)
  path = replace(path, '-'=>'_') |> strip
  appname = split(path, '/', keepempty = false)[end] |> String |> Inflector.from_underscores

  content = """
  module $appname

  using Genie, Logging, LoggingExtras

  function main()
    Core.eval(Main, :(const UserApp = \$(@__MODULE__)))

    Genie.genie(; context = @__MODULE__)

    Core.eval(Main, :(const Genie = UserApp.Genie))
    Core.eval(Main, :(using Genie))
  end

  end
  """

  (appname, content)
end

"""
    dockerfile(; user::String = "genie", supervisor::Bool = false, nginx::Bool = false, env::String = "dev",
                      filename::String = "Dockerfile", port::Int = 8000, dockerport::Int = 80, host::String = "0.0.0.0",
                      websockets_port::Int = port, websockets_dockerport::Int = dockerport)

Generates dockerfile for the Genie app.
"""
function dockerfile(; user::String = "genie", supervisor::Bool = false, nginx::Bool = false, env::String = "dev",
                      filename::String = "Dockerfile", port::Int = Genie.config.server_port, dockerport::Int = 80,
                      host::String = "0.0.0.0", websockets_port::Int = port, platform::String = "",
                      websockets_dockerport::Int = dockerport, earlybind::Bool = true)
  appdir = "/home/$user/app"

  string(
  """
  # pull latest julia image
  FROM $(isempty(platform) ? "" : "--platform=$platform") julia:latest

  # create dedicated user
  RUN useradd --create-home --shell /bin/bash $user

  # set up the app
  RUN mkdir $appdir
  COPY . $appdir
  WORKDIR $appdir

  # configure permissions
  RUN chown $user:$user -R *

  RUN chmod +x bin/repl
  RUN chmod +x bin/server
  RUN chmod +x bin/runtask

  # switch user
  USER $user

  # instantiate Julia packages
  RUN julia -e "using Pkg; Pkg.activate(\\".\\"); Pkg.instantiate(); Pkg.precompile(); "

  # ports
  EXPOSE $port
  EXPOSE $dockerport
  """,

  (websockets_port != port ?
  """

  # websockets ports
  EXPOSE $websockets_port
  EXPOSE $websockets_dockerport
  """ : ""),

  """

  # set up app environment
  ENV JULIA_DEPOT_PATH "/home/$user/.julia"
  ENV GENIE_ENV "$env"
  ENV HOST "$host"
  ENV PORT "$port"
  ENV WSPORT "$websockets_port"
  ENV EARLYBIND "$earlybind"
  """,

  """

  # run app
  CMD ["bin/server"]

  # or maybe include a Julia file
  # CMD julia -e 'using Pkg; Pkg.activate("."); include("IrisClustering.jl"); '
  """)
end

end


text:
using PowerSystems
using JuMP
base_dir = dirname(dirname(pathof(PowerSystems)))
include(joinpath(base_dir,"data/data_5bus_uc.jl"))
sys5 = PowerSystem(nodes5, generators5, loads5_DA, branches5, nothing,  100.0);
using PowerSimulations
const PS = PowerSimulations

simple_reserve = PowerSystems.StaticReserve("test_reserve",vcat(sys5.generators.thermal,sys5.generators.renewable[2]),60.0,[gen.tech for gen in sys5.generators.thermal])
#simple_reserve = PowerSystems.StaticReserve("test_reserve",sys5.generators.thermal,60.0,[sys5.generators.thermal[1].tech])

@test try
    Net = PS.CopperPlatePowerModel
    m = Model();
    netinjection = PS.instantiate_network(Net, sys5);
    PS.constructdevice!(m, netinjection, ThermalGen, PS.ThermalDispatch, Net, sys5);
    PS.constructdevice!(m, netinjection, RenewableGen, PS.RenewableCurtail, Net, sys5);
    PS.constructdevice!(m, netinjection, ElectricLoad, PS.InterruptibleLoad, Net, sys5);
    PS.constructnetwork!(m, [(device=Branch, formulation=PS.PiLine)], netinjection, Net, sys5)
    PS.constructservice!(m, simple_reserve, PS.RampLimitedReserve, [(device = ThermalGen, formulation =PS.ThermalDispatch),
                                                              (device = RenewableGen, formulation = PS.RenewableCurtail)],
                                                              sys5)
    m.obj_dict
true finally end

@test try
    Net = PS.CopperPlatePowerModel
    m = Model();
    netinjection = PS.instantiate_network(Net, sys5);
    PS.constructdevice!(m, netinjection, ThermalGen, PS.StandardThermalCommitment, Net, sys5);
    PS.constructdevice!(m, netinjection, RenewableGen, PS.RenewableCurtail, Net, sys5);
    PS.constructdevice!(m, netinjection, ElectricLoad, PS.InterruptibleLoad, Net, sys5);
    PS.constructnetwork!(m, [(device=Branch, formulation=PS.PiLine)], netinjection, Net, sys5)
    PS.constructservice!(m, simple_reserve, PS.RampLimitedReserve, [(device = ThermalGen, formulation =PS.StandardThermalCommitment),
                                                              (device = RenewableGen, formulation = PS.RenewableCurtail)],
                                                              sys5)
    m.obj_dict
true finally end


text:
"""
Pretty printing of everything metabolic-modelish.
"""
function Base.show(io::IO, ::MIME"text/plain", m::MetabolicModel)
    _pretty_print_keyvals(io, "", "Metabolic model of type $(typeof(m))")
    if n_reactions(m) <= _constants.default_stoich_show_size
        println(io, stoichiometry(m))
    else # too big to display nicely
        println(io, "S = [...]")
    end
    _pretty_print_keyvals(io, "Number of reactions: ", string(n_reactions(m)))
    _pretty_print_keyvals(io, "Number of metabolites: ", string(n_metabolites(m)))
end


text:
export iauTpxes
"""
In the tangent plane projection, given celestial spherical
coordinates for a star and the tangent point, solve for the star's
rectangular coordinates in the tangent plane.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   a,b       double  star's spherical coordinates
   a0,b0     double  tangent point's spherical coordinates

Returned:
   *xi,*eta  double  rectangular coordinates of star image (Note 2)

Returned (function value):
            int     status:  0 = OK
                              1 = star too far from axis
                              2 = antistar on tangent plane
                              3 = antistar too far from axis

Notes:

   1. The tangent plane projection is also called the "gnomonic
   projection" and the "central projection".

   2. The eta axis points due north in the adopted coordinate system.
   If the spherical coordinates are observed (RA,Dec), the tangent
   plane coordinates (xi,eta) are conventionally called the
   "standard coordinates".  For right-handed spherical coordinates,
   (xi,eta) are also right-handed.  The units of (xi,eta) are,
   effectively, radians at the tangent point.

   3. All angular arguments are in radians.

   4. This function is a member of the following set:

      spherical      vector         solve for

      > iauTpxes <    iauTpxev         xi,eta
      iauTpsts      iauTpstv          star
      iauTpors      iauTporv         origin

References:

   Calabretta M.R. & Greisen, E.W., 2002, "Representations of
   celestial coordinates in FITS", Astron.Astrophys. 395, 1077

   Green, R.M., "Spherical Astronomy", Cambridge University Press,
   1987, Chapter 13.

This revision:   2018 January 2

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauTpxes(double a, double b, double a0, double b0,
#              double *xi, double *eta)

function iauTpxes(a::Real, b::Real,
                  a0::Real, b0::Real)
   ref_xi = Ref{Float64}(0.0)
   ref_eta = Ref{Float64}(0.0)

   status = ccall((:iauTpxes, libsofa_c), Cint,
       (Cdouble, Cdouble, Cdouble, Cdouble,
       Ref{Cdouble}, Ref{Cdouble}),
       convert(Float64, a),
       convert(Float64, b),
       convert(Float64, a0),
       convert(Float64, b0),
       ref_xi, ref_eta)

   return status, ref_xi[], ref_eta[]
end


text:
include("C:\\Users\\darien.shannon\\Documents\\Code\\Julia\\JuliaSE\\src\\2DBeam.jl")

COORD = [Node(0, 0)
         Node(6, 0)
         Node(10, 0)
         Node(20, 0)]

MSUP = [Support(1, 1, 1)
        Support(2, 0, 0)
        Support(3, 1, 0)
        Support(4, 1, 0)]

EM = Dict( 1 => Material(28e6))

CP = Dict( 1 => Section(5.8e9/(1000^4)),
           2 => Section(1.5*5.8e9/(1000^4)))

MPRP = [Element(1, 2, 1, 2)
        Element(2, 3, 1, 1)
        Element(3, 4, 1, 1)]

memberLoads = [(1, DLoad(30, 0, 0, 6))
               (3, CLoad(150, 5))]

nodeLoads = [(2, NLoad((-200, 0)))
             (3, NLoad((0, -90)))]

forces, reactions = run_analysis(COORD, MSUP, EM, CP, MPRP, memberLoads, nodeLoads)
println(forces)
println(reactions)

max_stars_repo_name: Julia-SE/matrix-analysis

max_issues_repo_name: Julia-SE/matrix-analysis

max_forks_repo_name: Julia-SE/matrix-analysis


text:
include("../src/FluxFederated.jl")
using .FluxFederated
using Sockets
using Flux # Note: this takes awhile to load

"Adds a small number to model parameters"
function fake_update(model)
  W, restructure = Flux.destructure(model)
  num = 0.01 * rand()
  return restructure(W .+ Float32(num))
end

# Connect to localhost port 3000
con = connect(3000)
println("Connected to localhost port 3000")

# Establish shared AES key via RSA
shared_key = handshake_client(con)

# Receive a Flux AI model via socket connection
model_bytes = secure_receive_bytes(con, shared_key)
model = model_from_bytes(model_bytes) # --> Chain(Dense(10,5,relu),Dense(5,2),softmax)

# Update model
# ...
# TODO: Perform training here
model = fake_update(model)

# Send back updated model
model_bytes = model_to_bytes(model)
secure_send_bytes(con, model_bytes, shared_key)

println("Received model: $model\n")

println("Model parameters:\n", params(model))

close(con)