text:
import Base.randn!

function advi_rand{T}(d::MvNormal{T})
    x = Array{T}(length(d))
    Distributions.add!(PDMats.unwhiten!(d.Σ, randn!(x)), d.μ)
end
function advi_rand{T}(d::MvNormal{T}, n::Int)
    x = Array{T}(length(d), n)
    Distributions.add!(PDMats.unwhiten!(d.Σ, randn!(x)), d.μ)
end
advi_rand(d::Normal) = d.μ + d.σ * randn()
function advi_rand{T}(d::Normal{T}, n::Int)
    x = Array(T, n)
    randn!(x)
    scale!(d.σ, x)
    x[:] += d.μ
end

function advi_rand{T<:Real}(x::Vector{T}, full=false)
    if full
        p = Int((-3 + sqrt(9 + 8 * length(x)))/2)
        L = LowerTriangular(x[p+1:endof(x)])
        out = L * randn(p)
    else
        p = Int(length(x)/2)
        out = exp(x[p+1:endof(x)]) .* randn(p)
    end
    for i in 1:p
        @inbounds out[i] += x[i]
    end
    out
end

function advi_rand{T<:Real}(x::Vector{T}, n::Int, full=false)
    if full
        p = Int((-3 + sqrt(9 + 8 * length(x)))/2)
        L = LowerTriangular(x[p+1:endof(x)])
        out = L * randn(p, n)
    else
        p = Int(length(x)/2)
        out = scale(exp(x[p+1:endof(x)]), randn(p, n))
    end
    for i in 1:p
        @inbounds out[i] += x[i]
    end
    out
end

function randn!{T<:ForwardDiff.Dual}(A::AbstractArray{T})
    for i in eachindex(A)
        @inbounds A[i] = randn()
    end
    A
end

"""
Calculate the entropy of a (multivariate) Normal distribution based on
a vector of unconstrained parameters for the mean and covariance.
"""
function H{T<:Real}(x::Vector{T}, full=false)
    if full
        # number of parameters: will throw InexactError if not an integer
        p = Int((-3 + sqrt(9 + 8 * length(x)))/2)

        # ½ logdet(Σ) = ∑ log L_ii, with L the Cholesky factor
        # but log L_ii is just the unconstrained parameter
        offset = p + 1  # first element of covariance part
        ldet = 0.
        skip = p
        while offset ≤ length(x)
            ldet += x[offset]
            offset += skip
            skip -= 1
        end
    else
        # return a multivariate normal with diagonal covariance
        p = Int(length(x)/2)
        ldet = sum(x[p+1:end])/2  # ½ ∑ log(σ)
    end
    p * (log2π + 1)/2 + ldet
end

"""
Initialize unconstrained and constrained arrays of advi variables.
x is the vector of unconstrained parameters being optimized over
offset is the initial offset within x to begin reading from
rv is the random variable (isa(T, RVType) == true)
ζ is the array of (sampled) unconstrained variables
v is the array of (sampled) constrained variables
Returns a tuple containing the total ELBO contribution due to
    entropy plus Jacobians and the updated offset
"""
function advi_initialize!(x, offset, rv, ζ, v)
    L = zero(eltype(x))
    npars = VinDsl.num_pars_advi(rv)

    ctr = offset
    for i in eachindex(v)
        L += H(x[ctr:ctr+npars-1])
        if isa(rv, RScalar)
            ζ[i] = VinDsl.advi_rand(x[ctr:ctr+npars-1])[1]
        else
            ζ[i] = VinDsl.advi_rand(x[ctr:ctr+npars-1])
        end
        v[i] = constrain(rv, ζ[i])
        L += logdetjac(rv, ζ[i])
        ctr += npars
    end
    (L, ctr)
end

"""
Add an ADVI random variable to the model. Add contributions to the ELBO
for entropy and the log determinant of the Jacobian for the transformation from
unconstrained to constrained variables. Return the ELBO contribution, an
offset in the parameter vector marking the next unused parameter, and
a sample value of the constrained variable.
"""
function advi_variable(x, offset, rv, dims)
    C = VinDsl.storage_type(rv, eltype(x))
    ζ = Array{C}(dims...)
    v = Array{C}(dims...)

    L, new_offset = advi_initialize!(x, offset, rv, ζ, v)

    (L, new_offset, v)
end

function advi_variable(x, offset, rv)
    L = zero(eltype(x))
    npars = VinDsl.num_pars_advi(rv)

    L += H(x[offset:offset+npars-1])
    if isa(rv, RScalar)
        ζ = VinDsl.advi_rand(x[offset:offset+npars-1])[1]
    else
        ζ = VinDsl.advi_rand(x[offset:offset+npars-1])
    end
    v = constrain(rv, ζ)
    L += logdetjac(rv, ζ)

    (L, offset + npars, v)
end

macro advi_declarations(adviparamlist)
    esc(_advi_declarations(adviparamlist))
    #_advi_declarations(adviparamlist)
end

function _advi_declarations(adviexpr::Expr, paramlist)
    if adviexpr.head == :(::)
        advivname = adviexpr.args[1]
        advitypearg = adviexpr.args[2]
        if isa(advitypearg, Expr)
            if advitypearg.head == :ref   # array of variables
                advitype = _convert_typename(advitypearg.args[1])
                paramdims = advitypearg.args[2:end]
                dimsexpr = :(
                begin
                    ΔLoweradvi, adviparamctr, $advivname = VinDsl.advi_variable(adviparamlist, adviparamctr, $advitype, tuple($(paramdims...)))
                    Loweradvi += ΔLoweradvi
                end
                    )
            else  # single variable
                advitype = _convert_typename(advitypearg)
                dimsexpr = :(
                begin
                    ΔLoweradvi, adviparamctr, $advivname = VinDsl.advi_variable(adviparamlist, adviparamctr, $advitype)
                    Loweradvi += ΔLoweradvi
                end
                    )
            end
            push!(paramlist, dimsexpr)
            #println(paramlist)
        end

    else
        for aadiviex in filter(adviparamlist -> isa(adviparamlist, Expr), adviexpr.args)
            _advi_declarations(aadiviex, paramlist)
        end
    end
    Expr(:block, paramlist...)
end

_advi_declarations(adviexpr::Expr) = _advi_declarations(adviexpr, [])

function _convert_typename(adviexpr)
    adviexprout = copy(adviexpr)
    adviexprout.args[1] = Symbol("R", adviexprout.args[1])
    adviexprout
end

macro advi_model(adviparamlist)
    esc(_advi_model(adviparamlist))
end

function _advi_model(adviexpr::Expr)
    # if we have a ~ expression, replace with increment of L
    if adviexpr.head == :macrocall && adviexpr.args[1] == Symbol("@~")
        adviexprout = :(Loweradvi += logpdf($(adviexpr.args[3]), $(adviexpr.args[2])))
    else  # recursively parse
        adviexprout = copy(adviexpr)
        for indxadiviex in eachindex(adviexprout.args)
            adviexprout.args[indxadiviex] = isa(adviexprout.args[indxadiviex], Expr) ? _advi_model(adviexprout.args[indxadiviex]) : adviexprout.args[indxadiviex]
        end
    end
    adviexprout
end

macro ELBO(adviexpr)
    adviexprout = quote
        function ELBO(adviparamlist::Vector)
            adviparamctr = 1
            Loweradvi = zero(eltype(adviparamlist))
            $adviexpr
            Loweradvi
        end
    end

    esc(adviexprout)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Mortar3D.jl/blob/master/LICENSE

module Mortar3D
include("calculate_normals.jl")
include("calculate_projections.jl")
include("calculate_polygon_clip.jl")
end

ext: jl


text:
using Primes

println(factor(600851475143).pe[end][1])

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
using Base.Test
using DataFrames

let
	i = Index()
	push!(i, "A")
	push!(i, "B")

	inds = {1,
			1.0,
			"A",
			:A,
			[true],
			trues(1),
			[1],
			[1.0],
			1:1,
			1.0:1.0,
			["A"],
			[:A],
			DataVector[true],
			DataVector[1],
			DataVector[1.0],
			DataVector["A"],
			DataVector[:A],
			PooledDataVector[true],
			PooledDataVector[1],
			PooledDataVector[1.0],
			PooledDataVector["A"],
			PooledDataVector[:A]}

	for ind in inds
		if isequal(ind, "A") || isequal(ind, :A) || ndims(ind) == 0
			@assert isequal(i[ind], 1)
		else
			@assert (i[ind] == [1])
		end
	end
end


text:
function get_intensities_over_sensors(
    problem::LinearOptics,
    β::AbstractArray,
    all_sensors::AbstractMatrix;
    scattering = :farField,
)
    @debug "start : get intensities over sensors - LinearOptics"

    n_sensors = size(all_sensors, 2)
    intensities = zeros(n_sensors)

    v_r = view(problem.atoms.r, :, :)
    if scattering == :nearField
        scattering_func = _scattering_nearField
    elseif scattering == :farField
        scattering_func = _scattering_farField
    else
        @error "Invalid scattering Value: should be `:nearField` or `:farField` "
    end
    Threads.@threads for i = 1:n_sensors
        oneSensor = view(all_sensors, :, i)
        intensities[i] = _get_intensity_over_sensor(
            problem.atoms.shape,
            problem.laser,
            v_r,
            oneSensor,
            β;
            scattering_func = scattering_func,
        )
    end

    @debug "end  : get intensities over sensors - LinearOptics"
    return intensities
end
#= WORKS ONLY FOR 3D CLOUD DISTRIBUTIONS =#
function _get_intensity_over_sensor(
    shape::T,
    laser::Laser,
    atoms::AbstractMatrix,
    sensor::AbstractArray,
    β::AbstractArray;
    scattering_func = _scattering_farField,
) where {T<:ThreeD}
    ## Laser Pump
    E_L = (im / Γ) * apply_laser_over_oneSensor(laser, sensor)

    ## Scattered
    E_scatt = scattering_func(sensor, atoms, β)

    return abs2(E_L + E_scatt)
end

@inline function _scattering_nearField(sensor::AbstractVector, atoms, β::Vector{ComplexF64})
    E_scatt = zero(eltype(β))

    j = 1
    @inbounds for atom in eachcol(atoms)
        d_SensorAtom = sqrt( (sensor[1] - atom[1])^2 + (sensor[2] - atom[2])^2 + (sensor[3] - atom[3])^2 )
        E_scatt += cis(k₀ * d_SensorAtom) * (β[j] / d_SensorAtom)
        j += 1
    end
    E_scatt = +(Γ / 2) * im * E_scatt
    return E_scatt
end
@inline function _scattering_farField(sensor::AbstractVector, atoms, β::Vector{ComplexF64})
    E_scatt = zero(ComplexF64)
    n̂ = sensor / norm(sensor)

    dot_n_r = zero(ComplexF64)
    j = 1
    @inbounds for atom in eachcol(atoms)
        dot_n_r = n̂[1] * atom[1] + n̂[2] * atom[2] + n̂[3] * atom[3]
        dot_n_r = cis(-k₀ * dot_n_r)
        E_scatt += dot_n_r * β[j]
        j += 1
    end

    ikr = im * k₀ * norm(sensor)
    E_scatt = -(Γ / 2) * E_scatt * exp(ikr) / ikr
    return E_scatt
end

function get_intensity_over_an_angle(
    problem::LinearOptics{Scalar},
    atoms_states::AbstractVector,
    θ::Number,
)
    N = problem.atoms.N
    β = view(atoms_states, 1:N)

    ϕ_range = range(0, 2π, length = 30)
    vr = view(problem.atoms.r, :, :)

    complex_intensity = zeros(ComplexF64, N)
    total_intensity = 0.0

    for k = 1:length(ϕ_range)
        ϕ = ϕ_range[k]
        Threads.@threads for j = 1:N
            complex_intensity[j] =
                cis(
                    -k₀ * (
                        vr[1, j] * sin(θ) * cos(ϕ) +
                        vr[2, j] * sin(θ) * sin(ϕ) +
                        vr[3, j] * cos(θ)
                    ),
                ) * β[j]
        end
        total_intensity += abs2(sum(complex_intensity))
    end
    return total_intensity / length(ϕ_range)
end

function get_intensity_over_an_angle_shared(
    problem::LinearOptics{Scalar},
    β::AbstractVector,
    θ::Number,
    r_shared,
)
    N = problem.atoms.N

    ϕ_range = range(0, 2π, length = 30)

    complex_intensity = zeros(ComplexF64, N)
    total_intensity = 0.0

    for k = 1:length(ϕ_range)
        ϕ = ϕ_range[k]
        for j = 1:N
            complex_intensity[j] =
                cis(
                    -k₀ * (
                        r_shared[1, j] * sin(θ) * cos(ϕ) +
                        r_shared[2, j] * sin(θ) * sin(ϕ) +
                        r_shared[3, j] * cos(θ)
                    ),
                ) * β[j]
        end
        total_intensity += abs2(sum(complex_intensity))
    end
    return total_intensity / length(ϕ_range)
end

function get_intensity_over_an_angle(
    problem::LinearOptics{Scalar},
    atoms_states::Matrix,
    θ::Number,
)
    timeSteps = size(atoms_states, 2)
    intensities = zeros(timeSteps)

    r_shared = SharedArray{Float64,2}(3, problem.atoms.N)
    r_shared .= problem.atoms.r

    Threads.@threads for i = 1:timeSteps
        oneState = view(atoms_states, :, i)
        intensities[i] = get_intensity_over_an_angle_shared(problem, oneState, θ, r_shared)
    end

    return intensities
end

# ### --------------- MEAN FIELD ---------------
"""
	get_intensities_over_sensors(problem::NonLinearOptics{MeanField}, atoms_states::AbstractArray, sensors::AbstractArray

atoms_states: vcat(β, z)
sensor: [ [sensor1], [sensor2], ... ,[sensorN]  ]

**important**: check benchmark folders to understand how this code was produced
"""
function get_intensities_over_sensors(
    problem::NonLinearOptics{MeanField},
    atoms_states::AbstractArray,
    sensors::AbstractArray,
)
    @debug "start : get intensities over sensors - NonLinearOptics"

    N, r = problem.atoms.N, problem.atoms.r
    β, z = view(atoms_states, 1:N), view(atoms_states, (N+1):2N)

    number_configurations = ((N^2) ÷ 2 - N ÷ 2)

    βₙₘ = Array{ComplexF64}(undef, number_configurations)
    cont = 1
    for n = 1:N
        for m = (n+1):N
            βₙₘ[cont] = conj(β[n]) * β[m]
            cont += 1
        end
    end
    vβₙₘ = view(βₙₘ, :)

    rₙₘ = Array{ComplexF64}(undef, 3, number_configurations)
    cont = 1
    for n = 1:N
        r_n = @view(r[:, n])
        for m = (n+1):N
            rₙₘ[1, cont] = r_n[1] - r[1, m]
            rₙₘ[2, cont] = r_n[2] - r[2, m]
            rₙₘ[3, cont] = r_n[3] - r[3, m]
            cont += 1
        end
    end
    vrₙₘ = view(rₙₘ, :, :)

    n_sensors = size(sensors, 2)
    intensities = Float64[]

    if n_sensors == 1
        n_hat = sensors ./ norm(sensors)
        intensity =
            _oneSensor_MeanField_scattering(n_hat, vβₙₘ, vrₙₘ, number_configurations)
        intensity += sum((1 .+ z) ./ 4)

        push!(intensities, 2real(intensity))
    else
        const_sum = sum((1 .+ z) ./ 4)
        for oneSensor in eachcol(sensors)
            n_hat = oneSensor ./ norm(oneSensor)

            intensity =
                _manySensors_MeanField_scattering(n_hat, vβₙₘ, vrₙₘ, number_configurations)
            intensity += const_sum

            push!(intensities, 2real(intensity))
        end
    end

    @debug "end  : get intensities over sensors - NonLinearOptics"
    return intensities
end
function _oneSensor_MeanField_scattering(n_hat, vβₙₘ, vrₙₘ, number_configurations; k₀ = 1)
    intensity = ComplexF64(0)
    for cont = 1:number_configurations
        dot_n_r =
            n_hat[1] * vrₙₘ[1, cont] + n_hat[2] * vrₙₘ[2, cont] + n_hat[3] * vrₙₘ[3, cont]
        intensity += vβₙₘ[cont] * cis(k₀ * dot_n_r)
    end
    return intensity
end

function _manySensors_MeanField_scattering(n_hat, vβₙₘ, vrₙₘ, number_configurations; k₀ = 1)
    intensity = ComplexF64(0)
    intensity = Folds.mapreduce(+, 1:number_configurations) do k
        (
            begin
                @inbounds vβₙₘ[k] * cis(
                    k₀ * (
                        n_hat[1] * vrₙₘ[1, k] +
                        n_hat[2] * vrₙₘ[2, k] +
                        n_hat[3] * vrₙₘ[3, k]
                    ),
                )
            end
        )
    end
    return intensity
end

function get_intensity_over_an_angle(
    problem::NonLinearOptics{MeanField},
    atoms_states::Vector,
    θ::Float64,
)
    @debug "start : get intensity over an angle - NonLinearOptics{MeanField}"

    if is_integration_const_term_available(problem)
        Gₙₘ = problem.data[:Gₙₘ]
    else
        Gₙₘ = _get_meanfield_constant_term(problem.atoms, θ)
        problem.data[:Gₙₘ] = Gₙₘ
    end

    N = problem.atoms.N
    β = view(atoms_states, 1:N)
    z = view(atoms_states, (N+1):2N)

    βₙₘ = transpose(β * β') # I have to do "transpose" and NOT "adjoint = complex+tranpose"
    βₙₘ[diagind(βₙₘ)] .= (1 .+ z) ./ 2

    # IMPORTANT FOR NEXT LINE: you should use ELEMENT WISE multiplication.
    # Also, you can reduce memory allocation with inplace multiplication
    βₙₘ .*= Gₙₘ
    intensity = real(sum(βₙₘ))

    @debug "end  : get intensity over an angle - NonLinearOptics{MeanField}"
    return intensity
end

function is_integration_const_term_available(problem)
    if haskey(problem.data, :Gₙₘ)
        return true
    else
        return false
    end
end

function _get_meanfield_constant_term(atoms, Θ)
    N, r = atoms.N, atoms.r

    xₙₘ, yₙₘ, zₙₘ = get_xyz_distances(r)
    k₀sinΘ = k₀ * sin(Θ)
    cos_Θ = cos(Θ)

    Gₙₘ_shared = SharedArray{ComplexF64,2}(N, N)
    @sync for n = 1:N
        Threads.@spawn for m = 1:N # we had to compute all terms, and not the upper part
            @inbounds Gₙₘ_shared[n, m] =
                _constant_term_core_computation(xₙₘ, yₙₘ, zₙₘ, n, m, cos_Θ, k₀sinΘ)
        end
    end
    Gₙₘ = Array(Gₙₘ_shared)
    #=
        IF I want to compute only the upper part,
        I have to multiply all terms by π/2:   Gₙₘ = (π/2)Array(Gₙₘ_shared)

        I decided to don't make this, to don't appear with factors
        not mentioned on theory.
    =#
    #=
        Before returning, we HAVE to do some memory cleaning,
        EVEN losing some performance.

        Without this cleaning, the garbage collector gets lost
        outside this function - when many simulation occurs at the same time.
    =#
    Gₙₘ_shared = xₙₘ = yₙₘ = zₙₘ = 1
    GC.gc()  # DO NOT DELETE
    return Gₙₘ
end
function _constant_term_core_computation(xₙₘ, yₙₘ, zₙₘ, n, m, cos_Θ, k₀sinΘ)
    a = zero(ComplexF64)
    a = cis(k₀ * zₙₘ[n, m] * cos_Θ) * besselj(0, k₀sinΘ * sqrt(xₙₘ[n, m]^2 + yₙₘ[n, m]^2))
    return a
end
function get_xyz_distances(r) # @memoize  --> creating warning, let's ignore it right now
    dimensions = size(r, 1)
    N = size(r, 2)

    xₙₘ = SharedArray{Float64,2}(N, N)
    yₙₘ = SharedArray{Float64,2}(N, N)
    zₙₘ = SharedArray{Float64,2}(N, N)

    r_shared = SharedArray{Float64,2}(dimensions, N)
    r_shared .= r
    @sync for n = 1:N
        r_n = view(r_shared, :, n)
        Threads.@spawn for m = 1:N
            xₙₘ[n, m] = r_n[1] - r_shared[1, m]
            yₙₘ[n, m] = r_n[2] - r_shared[2, m]
            zₙₘ[n, m] = r_n[3] - r_shared[3, m]
        end
    end

    return xₙₘ, yₙₘ, zₙₘ
end


text:
println("This is a test on the subreferencing of Associative Array")

A = ReadCSV("A.csv")

A1r = A["a,b,",:];               # Get rows a and b.
#TODO A2r = A("a *,",1:3)             # Get rows containing a and columns 1 thru 3.
#TODO A3r = A("a,:,b,",:)             # Get rows a to b.
#A4r = A(StartsWith("a,c,"),:)   # Get rows starting with a or c.
A1c = A[:,"a,b,"]               # Get cols a and b.
#TODO A2c = A(1:3,"a *,")             # Get rows 1 thru 3 and cols containing a.
#TODO A3c = A(:,"a,:,b,")             # Get cols a to b.
#TODO A4c = A(:,StartsWith("a,c,"))   # Get cols starting with a or b.
#TODO A1v = (A < "b,")                # Get all values less than b.

#print(A1v)


text:
using Stheno, Plots, Distributions, ColorTypes, Random, FixedPointNumbers
using Stheno: @model, Titsias

rng = MersenneTwister(12345);

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

# Specify generative model.
@model function gp(σ²)
    f = 1.5 * GP(EQ())
    y = f + GP(Noise(σ²))
    return f, y
end

# Construct generative model.
σ² = 1.0;
f, y = gp(σ²);

# Sample prior at random locations.
Xy = rand(Uniform(-5, 5), 100);
y_ = rand(y(Xy));

# Approximately condition processes. Play around with different lengths of Z.
Z = range(-5, stop=5, length=10);
f′ = f | Titsias(f(Xy)←y_, f(Z), sqrt(σ²));
f′_exact = f | (y(Xy)←y_);

#######################  Do posterior inference give a few samples  #######################

Np, S = 500, 10;
Xp = range(-6; stop=6, length=Np);
f′Xp = rand(rng, f′(Xp), S);

# Get posterior marginals for approximate and exact conditioning.
μ′f, σ′f = marginals(f′(Xp));
μf′_exact, σf′_exact = marginals(f′_exact(Xp));

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

plotly();

posterior_plot = plot(
    legend=:topright,
    legendfont=Plots.Font(
        "sans-serif",
        10,
        :hcenter,
        :vcenter,
        0.0,
        RGB{Normed{UInt8, 8}}(0.0,0.0,0.0)
    ),
    background_color_legend=RGBA(1, 1, 1, 0),
    foreground_color_legend=RGBA(1, 1, 1, 0),
);

plot!(posterior_plot, Xp, μf′_exact;
    linewidth=2.0,
    linecolor=:blue,
    label="");
plot!(posterior_plot, Xp, μf′_exact;
    linewidth=0.0,
    linecolor=:blue,
    linealpha=0.2,
    fillrange=[μf′_exact.- 3  .* σf′_exact],
    fillalpha=0.2,
    fillcolor=:blue,
    label=["exact posterior, f", ""]);
plot!(posterior_plot, Xp, μf′_exact;
    linewidth=0.0,
    linecolor=:blue,
    linealpha=0.2,
    fillrange=[μf′_exact .+ 3 .* σf′_exact],
    fillalpha=0.2,
    fillcolor=:blue,
    label="");

plot!(posterior_plot, Xp, μ′f;
    linewidth=2.0,
    linecolor=:red,
    label="");
plot!(posterior_plot, Xp, μ′f;
    linewidth=0.0,
    linecolor=:red,
    linealpha=0.2,
    fillrange=[μ′f.- 3  .* σ′f],
    fillalpha=0.2,
    fillcolor=:red,
    label=["posterior, f", ""]);
plot!(posterior_plot, Xp, μ′f;
    linewidth=0.0,
    linecolor=:red,
    linealpha=0.2,
    fillrange=[μ′f .+ 3 .* σ′f],
    fillalpha=0.2,
    fillcolor=:red,
    label="");

scatter!(posterior_plot, Xy, y_;
    markercolor=:red,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=6,
    markeralpha=0.3,
    label="y_");

display(posterior_plot);


text:
# ******************************************************************************************
# Notices:
#
# Copyright © 2022 United States Government as represented by the Administrator of the
# National Aeronautics and Space Administration.  All Rights Reserved.
#
# Disclaimers
#
# No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND,
# EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY
# THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY
# WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION,
# IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER,
# CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS,
# RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM
# USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND
# LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND
# DISTRIBUTES IT "AS IS."
#
# Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED
# STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.
# IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES,
# EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON,
# OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND
# HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL
# AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY
# SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
# ******************************************************************************************

using AdaStress
using NBInclude
using Test

@testset "Interface" begin
    @test AdaStress.GrayBox <: Interface.AbstractSimulation
end

@testset "Solvers" begin
    @test begin
        AdaStress.enable("SoftActorCritic")
        using AdaStress.SoftActorCritic
        ac = SoftActorCritic.MLPActorCritic(1, 1, [-1.0], [1.0])
        ac isa AdaStress.GlobalResult
    end
end

@testset "Analysis" begin
    @test begin
        AdaStress.enable("PolicyValueVerification")
        using AdaStress.PolicyValueVerification
        ac = SoftActorCritic.MLPActorCritic(1, 1, [-1.0], [1.0])
        SoftActorCritic.to_cpu!(ac)
        nnet = PolicyValueVerification.policy_network(ac; act_mins=[-1.0], act_maxs=[1.0])
        nnet isa PolicyValueVerification.ExtendedNetwork
    end
end

# Example notebooks to be tested (cells beginning with `# autoskip` are not executed)
envs = ["walk1d", "walk2d", "pvv"]

for env in envs
    @testset "Example: $env" begin
        dir = joinpath(@__DIR__, "..", "examples", env)
        cd(dir) do
            @test begin
                @nbinclude(joinpath(dir, "$env.ipynb"); regex=r"^(?!# autoskip).*")
                true
            end
        end
    end
end


text:
using Ju
using BenchmarkTools

state_size = (80, 80, 1)
actions = 1:10

buffer = CircularSARDBuffer(10^5; state_type=Array{Float64, 3}, state_size=state_size)

push!(buffer, rand(state_size...), rand(actions))

println("\n", repeat('=', 50))
println("\n push! buffer\n")
display(@benchmark push!($buffer, 1.0, false, $(rand(state_size...)), $(rand(actions))))

batch_size = 32
println("\n", repeat('=', 50))
println("\n sample buffer (batch_size=$batch_size)\n")
display(@benchmark sample($buffer, $batch_size))

batch_size = 64
println("\n", repeat('=', 50))
println("\n sample buffer (batch_size=$batch_size)\n")
display(@benchmark sample($buffer, $batch_size))

function sample_N_batches(buffer, N, batch_size)
    for _ in 1:N
        sample(buffer, batch_size)
    end
end

N, batch_size = 10, 32
println("\n", repeat('=', 50))
println("\n sample buffer (batch_size=$batch_size) $N times\n")
display(@benchmark sample_N_batches($buffer, $N, $batch_size))

N, batch_size = 10, 64
println("\n", repeat('=', 50))
println("\n sample buffer (batch_size=$batch_size) $N times\n")
display(@benchmark sample_N_batches($buffer, $N, $batch_size))

# ==================================================

#  push! buffer

# BenchmarkTools.Trial:
#   memory estimate:  0 bytes
#   allocs estimate:  0
#   --------------
#   minimum time:     4.657 μs (0.00% GC)
#   median time:      5.114 μs (0.00% GC)
#   mean time:        5.138 μs (0.00% GC)
#   maximum time:     9.500 μs (0.00% GC)
#   --------------
#   samples:          10000
#   evals/sample:     7

# ==================================================

#  sample buffer (batch_size=32)

# BenchmarkTools.Trial:
#   memory estimate:  3.13 KiB
#   allocs estimate:  34
#   --------------
#   minimum time:     958.163 ns (0.00% GC)
#   median time:      1.081 μs (0.00% GC)
#   mean time:        1.239 μs (11.31% GC)
#   maximum time:     25.488 μs (94.94% GC)
#   --------------
#   samples:          10000
#   evals/sample:     43

# ==================================================

#  sample buffer (batch_size=64)

# BenchmarkTools.Trial:
#   memory estimate:  5.09 KiB
#   allocs estimate:  34
#   --------------
#   minimum time:     1.270 μs (0.00% GC)
#   median time:      1.460 μs (0.00% GC)
#   mean time:        1.709 μs (9.25% GC)
#   maximum time:     80.030 μs (95.69% GC)
#   --------------
#   samples:          10000
#   evals/sample:     10

# ==================================================

#  sample buffer (batch_size=32) 10 times

# BenchmarkTools.Trial:
#   memory estimate:  30.47 KiB
#   allocs estimate:  320
#   --------------
#   minimum time:     8.699 μs (0.00% GC)
#   median time:      10.000 μs (0.00% GC)
#   mean time:        13.832 μs (9.60% GC)
#   maximum time:     1.497 ms (98.13% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1

# ==================================================

#  sample buffer (batch_size=64) 10 times

# BenchmarkTools.Trial:
#   memory estimate:  50.16 KiB
#   allocs estimate:  320
#   --------------
#   minimum time:     12.599 μs (0.00% GC)
#   median time:      14.301 μs (0.00% GC)
#   mean time:        16.655 μs (9.19% GC)
#   maximum time:     741.201 μs (96.80% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1


text:
using Plots, Test

@testset "Shapes" begin
    @testset "Type" begin
        square = Shape([(0, 0.0), (1, 0.0), (1, 1.0), (0, 1.0)])
        @test isa(square, Shape{Int64,Float64})
        @test coords(square) isa Tuple{Vector{S},Vector{T}} where {T,S}
    end

    @testset "Copy" begin
        square = Shape([(0, 0), (1, 0), (1, 1), (0, 1)])
        square2 = Shape(square)
        @test square2.x == square.x
        @test square2.y == square.y
    end

    @testset "Center" begin
        square = Shape([(0, 0), (1, 0), (1, 1), (0, 1)])
        @test Plots.center(square) == (0.5, 0.5)
    end

    @testset "Translate" begin
        square = Shape([(0, 0), (1, 0), (1, 1), (0, 1)])
        squareUp = Shape([(0, 1), (1, 1), (1, 2), (0, 2)])
        squareUpRight = Shape([(1, 1), (2, 1), (2, 2), (1, 2)])

        @test Plots.translate(square, 0, 1).x == squareUp.x
        @test Plots.translate(square, 0, 1).y == squareUp.y

        @test Plots.center(translate!(square, 1)) == (1.5, 1.5)
    end

    @testset "Rotate" begin
        # 2 radians rotation matrix
        R2 = [cos(2) sin(2); -sin(2) cos(2)]
        coords = [0 0; 1 0; 1 1; 0 1]'
        coordsRotated2 = R2 * (coords .- 0.5) .+ 0.5

        square = Shape([(0, 0), (1, 0), (1, 1), (0, 1)])

        # make a new, rotated square
        square2 = Plots.rotate(square, -2)
        @test square2.x ≈ coordsRotated2[1, :]
        @test square2.y ≈ coordsRotated2[2, :]

        # unrotate the new square in place
        rotate!(square2, 2)
        @test square2.x ≈ coords[1, :]
        @test square2.y ≈ coords[2, :]
    end

    @testset "Plot" begin
        ang = range(0, 2π, length = 60)
        ellipse(x, y, w, h) = Shape(w * sin.(ang) .+ x, h * cos.(ang) .+ y)
        myshapes = [ellipse(x, rand(), rand(), rand()) for x = 1:4]
        @test coords(myshapes) isa Tuple{Vector{Vector{S}},Vector{Vector{T}}} where {T,S}
        local p
        @test_nowarn p = plot(myshapes)
        @test p[1][1][:seriestype] == :shape
    end
end

@testset "Brush" begin
    @testset "Colors" begin
        baseline = brush(1, RGB(0, 0, 0))
        @test brush(:black) == baseline
        @test brush("black") == baseline
    end
    @testset "Weight" begin
        @test brush(10).size == 10
        @test brush(0.1).size == 1
    end
    @testset "Alpha" begin
        @test brush(0.4).alpha == 0.4
        @test brush(20).alpha == nothing
    end
    @testset "Bad Argument" begin
        # using test_logs because test_warn seems to not work anymore
        @test_logs (:warn, "Unused brush arg: nothing (Nothing)") begin
            brush(nothing)
        end
    end
end

@testset "Fonts" begin
    @testset "Scaling" begin
        sizesToCheck = [
            :titlefontsize,
            :legendfontsize,
            :legendtitlefontsize,
            :xtickfontsize,
            :ytickfontsize,
            :ztickfontsize,
            :xguidefontsize,
            :yguidefontsize,
            :zguidefontsize,
        ]
        # get inital font sizes
        initialSizes = [Plots.default(s) for s in sizesToCheck]

        #scale up font sizes
        scalefontsizes(2)

        # get inital font sizes
        doubledSizes = [Plots.default(s) for s in sizesToCheck]

        @test doubledSizes == initialSizes * 2

        # reset font sizes
        resetfontsizes()

        finalSizes = [Plots.default(s) for s in sizesToCheck]

        @test finalSizes == initialSizes
    end
end

@testset "Series Annotations" begin
    square = Shape([(0, 0), (1, 0), (1, 1), (0, 1)])
    @test_logs (:warn, "Unused SeriesAnnotations arg: triangle (Symbol)") begin
        p = plot(
            [1, 2, 3],
            series_annotations = (
                ["a"],
                2,    # pass a scale factor
                (1, 4), # pass two scale factors (overwrites first one)
                square, # pass a shape
                font(:courier), # pass an annotation font
                :triangle,  # pass an incorrect argument
            ),
        )
        sa = p.series_list[1].plotattributes[:series_annotations]
        @test sa.strs == ["a"]
        @test sa.font.family == "courier"
        @test sa.baseshape == square
        @test sa.scalefactor == (1, 4)
    end
    spl = scatter(
        4.53 .* [1/1 1/2 1/3 1/4 1/5],
        [0 0 0 0 0],
        layout = (5, 1),
        ylims = (-1.1, 1.1),
        xlims = (0, 5),
        series_annotations = permutedims([["1/1"],["1/2"],["1/3"],["1/4"],["1/5"]]),
    )
    @test spl.series_list[1].plotattributes[:series_annotations].strs == ["1/1"]
    @test spl.series_list[2].plotattributes[:series_annotations].strs == ["1/2"]
    @test spl.series_list[3].plotattributes[:series_annotations].strs == ["1/3"]
    @test spl.series_list[4].plotattributes[:series_annotations].strs == ["1/4"]
    @test spl.series_list[5].plotattributes[:series_annotations].strs == ["1/5"]
end


text:
module FCG

include("../Misc.jl")
using .Misc

using Distributions, SpecialFunctions, Einsum, Combinatorics
import Base.Iterators: product

export generate, log_marginal

function generate(T::Ƶ, dims::Vararg{Ƶ}; γ::ℜ=0.1) where {ℜ<:Real, Ƶ<:Int}
    len = prod(dims)
    θ_flat = rand(Dirichlet(γ,len))

    X::Array{ℜ} = reshape(rand(Multinomial(T,θ_flat)),dims)

    return X, reshape(θ_flat,dims)
end

function log_marginal(X::Array{ℜ,N}; μ::ℜ=nanmean(X), γ::ℜ=0.1) where {ℜ<:Real, Ƶ<:Int, N}
    dims, len = size(X), length(X)
    X₊ = sum(X)

    a::ℜ, b::ℜ = len*γ, γ/μ

    log_PX = a*log(b) - (a + X₊)*log(b + 1.0) + lgamma(a + X₊) - lgamma(a) - sum(lgamma.(X .+ 1.0))
    log_PX += sum(lgamma.(γ .+ X)) - lgamma(a + X₊) + lgamma(a) - len*lgamma(γ)

    return log_PX
end

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
@testset "Enmap I/O" begin
    for trim in (true, false)
        imap = read_map("data/test.fits"; trim=trim)
        @test size(imap) == (100, 100, 3)
        @test imap.wcs.naxis == 2
        @test collect(imap.wcs.cdelt) == [-1, 1]
        @test collect(imap.wcs.crval) == [0.5, 0.0]
        @test sum(imap) ≈ 14967.2985
        # read with sel
        imap = read_map("data/test.fits", sel=(11:20,21:40,1:2); trim=trim)
        @test size(imap) == (10,20,2)
    end
    # todo: add tests on IAU conversion w/ and w/o sel
end


text:
using OptimalApplication
using Test
using Random

# Number of markets to generate for each test
const n_markets = 3

@testset verbose = true "OptimalApplication.jl" begin

    @testset verbose = true "Same app. costs" begin
        @testset verbose = true "t int" begin
            m = 20

            for _ in 1:n_markets
                mkt = SameCostsMarket(m)

                X_enum, vX_enum = optimalportfolio_enumerate(mkt)
                sort!(X_enum)
                X_heap, vX_heap = applicationorder_heap(mkt)
                X_list, vX_list = applicationorder_list(mkt)
                @test X_enum == sort(X_heap)
                @test vX_enum ≈ last(vX_heap)
                @test X_enum == sort(X_list)
                @test vX_enum ≈ last(vX_list)
            end
        end

        @testset verbose = true "t float" begin
            m = 20

            for _ in 1:n_markets
                mkt = SameCostsMarket(rand(16), 12 * rand(16), m ÷ 4)

                X_enum, vX_enum = optimalportfolio_enumerate(mkt)
                sort!(X_enum)
                X_heap, vX_heap = applicationorder_heap(mkt)
                X_list, vX_list = applicationorder_list(mkt)
                @test X_enum == sort(X_heap)
                @test vX_enum ≈ last(vX_heap)
                @test X_enum == sort(X_list)
                @test vX_enum ≈ last(vX_list)
            end
        end
    end

    @testset verbose = true "Varied app. costs" begin
        m = 16

        @testset "Exact algorithms" begin
            for _ in 1:n_markets
                mkt = VariedCostsMarket(m)

                X_enum, vX_enum = optimalportfolio_enumerate(mkt)
                sort!(X_enum)
                X_dp, vX_dp = optimalportfolio_dynamicprogram(mkt)
                X_bnb, vX_bnb = optimalportfolio_branchbound(mkt)

                @test X_enum == sort(X_dp)
                @test vX_enum ≈ vX_dp
                @test X_enum == sort(X_bnb)
                @test vX_enum ≈ vX_bnb
            end
        end

        m = 50
        ε = 0.1

        @testset "FPTAS" begin
            for _ in 1:n_markets
                mkt = VariedCostsMarket(m)

                W, vW = optimalportfolio_fptas(mkt, ε)
                Y, vY = optimalportfolio_dynamicprogram(mkt)

                @test vW / vY ≥ 1 - ε
            end
        end
    end

    # Verify that permuting the problem data doesn't change the optimum
    @testset verbose = true "Permutation" begin
        @testset verbose = true "Same costs" begin
            m = 6
            f = rand(m)
            t = rand(m)
            h = m ÷ 2

            perm = randperm(m)
            sp = sortperm(t)

            # Randomly sorted market
            mkt1 = SameCostsMarket(f, t, h)
            # A random permutation thereof
            mkt2 = SameCostsMarket(f[perm], t[perm], h)
            # Now sorted by t
            mkt3 = SameCostsMarket(f[sp], t[sp], h)

            # A random portfolio
            x = rand(Bool, m)
            @test valuation((1:m)[x], mkt1) ≈ valuation((1:m)[x[perm]], mkt2)
            @test valuation((1:m)[x], mkt1) ≈ valuation((1:m)[x[sp]], mkt3)

            # Optimal valuations all agree
            x, v = optimalportfolio_enumerate(mkt3)
            @test v ≈ valuation(applicationorder_list(mkt1)[1], mkt1)
            @test v ≈ valuation(applicationorder_list(mkt2)[1], mkt2)
            @test v ≈ valuation(applicationorder_list(mkt3)[1], mkt3)
            @test v ≈ valuation(applicationorder_heap(mkt1)[1], mkt1)
            @test v ≈ valuation(applicationorder_heap(mkt2)[1], mkt2)
            @test v ≈ valuation(applicationorder_heap(mkt3)[1], mkt3)
            @test v ≈ valuation(optimalportfolio_enumerate(mkt1)[1], mkt1)
            @test v ≈ valuation(optimalportfolio_enumerate(mkt2)[1], mkt2)
        end

        @testset verbose = true "Het. costs" begin
            m = 6
            f = rand(m)
            t = rand(50:60, m)
            g = rand(5:10, m)
            H = sum(g) ÷ 2

            perm = randperm(m)
            sp = sortperm(t)

            # Randomly sorted market
            mkt1 = VariedCostsMarket(f, t, g, H)
            # A random permutation thereof
            mkt2 = VariedCostsMarket(f[perm], t[perm], g[perm], H)
            # Now sorted by t
            mkt3 = VariedCostsMarket(f[sp], t[sp], g[sp], H)

            # A random portfolio
            x = rand(Bool, m)
            @test valuation((1:m)[x], mkt1) ≈ valuation((1:m)[x[perm]], mkt2)
            @test valuation((1:m)[x], mkt1) ≈ valuation((1:m)[x[sp]], mkt3)

            # Optimal valuations all agree
            x, v = optimalportfolio_enumerate(mkt3)
            @test v ≈ valuation(optimalportfolio_branchbound(mkt1)[1], mkt1)
            @test v ≈ valuation(optimalportfolio_branchbound(mkt2)[1], mkt2)
            @test v ≈ valuation(optimalportfolio_branchbound(mkt3)[1], mkt3)
            @test v ≈ valuation(optimalportfolio_dynamicprogram(mkt1)[1], mkt1)
            @test v ≈ valuation(optimalportfolio_dynamicprogram(mkt2)[1], mkt2)
            @test v ≈ valuation(optimalportfolio_dynamicprogram(mkt3)[1], mkt3)
            @test v ≈ valuation(optimalportfolio_enumerate(mkt1)[1], mkt1)
            @test v ≈ valuation(optimalportfolio_enumerate(mkt2)[1], mkt2)

            @test v ≤ 2 * valuation(optimalportfolio_fptas(mkt1, 0.5)[1], mkt1)
            @test v ≤ 2 * valuation(optimalportfolio_fptas(mkt2, 0.5)[1], mkt2)
            @test v ≤ 2 * valuation(optimalportfolio_fptas(mkt3, 0.5)[1], mkt3)
        end
    end

    @testset verbose = true "Edge cases" begin
        @testset verbose = true "Same costs" begin
            # h > m should have all schools
            mkt = SameCostsMarket(rand(4), rand(4), 6)

            @test 1:4 == sort(applicationorder_heap(mkt)[1])
            @test 1:4 == sort(applicationorder_list(mkt)[1])
            @test 1:4 == sort(optimalportfolio_enumerate(mkt)[1])
        end

        @testset verbose = true "Het. costs" begin
            # H ≥ sum(g) should have all schools
            mkt = VariedCostsMarket(fill(0.5, 4), fill(2, 4), fill(4, 4), 16)

            @test 1:4 == sort(optimalportfolio_branchbound(mkt)[1])
            @test 1:4 == sort(optimalportfolio_dynamicprogram(mkt)[1])
            @test 1:4 == sort(optimalportfolio_fptas(mkt, 0.05)[1])
            @test 1:4 == sort(optimalportfolio_enumerate(mkt)[1])
            @test 1:4 == sort(optimalportfolio_greedy(mkt)[1])
            @test 1:4 == sort(optimalportfolio_simulatedannealing(mkt; nit=4)[1])

            # g[j] > H should be excluded
            mkt = VariedCostsMarket([0.5, 0.5, 1.0], [1, 2, 100], [3, 3, 5], 4)

            @test [2] == optimalportfolio_branchbound(mkt)[1]
            @test [2] == optimalportfolio_dynamicprogram(mkt)[1]
            @test [2] == optimalportfolio_fptas(mkt, 0.25)[1]
            @test [2] == optimalportfolio_enumerate(mkt)[1]
            @test [2] == sort(optimalportfolio_greedy(mkt)[1])
            @test [2] == sort(optimalportfolio_simulatedannealing(mkt; nit=4)[1])
        end
    end

    @testset verbose = true "Heuristic algorithms" begin
        mkt = VariedCostsMarket([0.34, 0.8, 0.1, 0.5], [9, 2, 2, 4], [3, 4, 4, 3], 7)
        X, v = optimalportfolio_greedy(mkt)
        @test v == valuation(X, mkt)
        @test sort(X) == [1, 4]
        Y, w = optimalportfolio_simulatedannealing(mkt; nit=10)
        @test w == valuation(Y, mkt)
        @test w ≥ v

        mkt = VariedCostsMarket([0.131, 0.136, 0.09, 0.016], [2, 3, 9, 54], [9, 5, 5, 9], 14)
        X, v = optimalportfolio_greedy(mkt)
        @test v == valuation(X, mkt)
        @test sort(X) == [3, 4]
        Y, w = optimalportfolio_simulatedannealing(mkt; nit=10)
        @test w == valuation(Y, mkt)
        @test w ≥ v

        m = 20
        for _ in 1:n_markets
            mkt = VariedCostsMarket(m)

            X, v = optimalportfolio_greedy(mkt)
            @test v == valuation(X, mkt)
            Y, w = optimalportfolio_simulatedannealing(mkt; nit=50)
            @test w == valuation(Y, mkt)
            @test w ≥ v
        end
    end

    @testset verbose = true "Large problems" begin
        mkt = SameCostsMarket(5000)
        X, V = applicationorder_list(mkt)
        @test !isempty(X)

        mkt = VariedCostsMarket(500)
        W, vW = optimalportfolio_fptas(mkt, 0.5)
        Y, vY = optimalportfolio_dynamicprogram(mkt)
        @test vW / vY ≥ 0.5
        @test !isempty(optimalportfolio_greedy(mkt)[1])
        @test !isempty(optimalportfolio_simulatedannealing(mkt; nit=10)[1])
    end

    @testset verbose = true "Bad markets" begin
        # Dim mismatch
        f = [0.1, 0.1]
        t = [1, 2, 4]
        g = [2, 2]
        H = 3

        @test_throws DimensionMismatch Market(f, t, 1)
        @test_throws DimensionMismatch Market(f, t, g, H)

        # f not in (0, 1]
        f = [5.0, 1.0]
        t = [4, 7]

        @test_throws DomainError Market(f, t, 1)
        @test_throws DomainError Market(f, t, g, H)

        # t negative
        f = [5.0, 1.0]
        t = [4, -7]

        @test_throws DomainError Market(f, t, 1)
        @test_throws DomainError Market(f, t, g, H)
    end
end

max_stars_repo_name: maxkapur/OptimalApplication

max_issues_repo_name: maxkapur/OptimalApplication

max_forks_repo_name: maxkapur/OptimalApplication


text:
using OffsetArrays

#============================================================================================
function copyq1(meqn::Int,mbc::Int,mx::Int,q1::OffsetArray{Float64},q2::OffsetArray{Float64})
#============================================================================================

     # copy the contents of q1 into q2

#      implicit double precision (a-h,o-z)
#      dimension q1(meqn,1-mbc:mx+mbc)
#      dimension q2(meqn,1-mbc:mx+mbc)

    for i = 1-mbc : mx+mbc
        for m=1:meqn
            q2[m,i] = q1[m,i]
        end
    end

end


text:
## ---

### Rate calculation used in mingpt
function cosRate( tokens, lr=6e-4, warmuptoken=512*20, len_train_dataset=1115394, blocksize=128)
    finalTok=2*len_train_dataset*blocksize
    progress = (tokens - warmuptoken) / max(1, finalTok - warmuptoken)
    lr_mult = max(0.1, 0.5 * (1.0 + cos(pi * progress)))
    lr_mult*lr
end

@show cosRate(60*512*128)
@show cosRate(125*512*128)
@show cosRate(235*512*128)
@show cosRate(350*512*128)

max_stars_repo_name: daryoush/CudaAndFluxExperiments

max_issues_repo_name: daryoush/CudaAndFluxExperiments

max_forks_repo_name: daryoush/CudaAndFluxExperiments


text:
julia_version = v"1.3.1"

include("../common.jl")


text:
immutable OrBind
  pat1
  pat2
end

or_(a, b) = OrBind(a, b)
or_(p...) = foldl(or_, p)

function match_inner(pat::OrBind, ex, env)
  env′ = trymatch(pat.pat1, ex)
  env′ == nothing ? match(pat.pat2, ex, env) : merge!(env, env′)
end

isor(ex) = isexpr(ex, :call) && ex.args[1] in (:or_, :|)

function ornew(ex)
  isor(ex) || return ex
  or_(ex.args[2:end]...)
end

subor(s) = s
subor(s::Symbol) = s
subor(s::Expr) = isor(s) ? subor(ornew(s)) : Expr(s.head, map(subor, s.args)...)
subor(s::OrBind) = OrBind(subor(s.pat1), subor(s.pat2))


text:
function TronDenseMatrix(I::VI, J::VI, V::oneArray, n) where {VI}
    @assert n >= 1
    @assert length(I) == length(J) == length(V)

    A = TronDenseMatrix{oneArray{Float64, 2}}(n, n, tron_zeros(oneArray{eltype(V)}, (n, n)))
    for i=1:length(I)
        @assert 1 <= I[i] <= n && 1 <= J[i] <= n && I[i] >= J[i]
        @inbounds A.vals[I[i], J[i]] += V[i]
    end

    return A
end

@inline function ExaTronKernels.Base.fill!(w::oneDeviceArray{Float64,1}, val::Float64)
    tx = get_local_id()
    bx = get_group_id()

    if tx <= length(w) && bx == 1
        @inbounds w[tx] = val
    end
    barrier()

    return
end


text:
using Clang

# Currently https://github.com/JaneliaSciComp/ImarisWriter/tree/disablePyramidOption
const IMARISWRITER_INCLUDE = raw"C:\Users\kittisopikulm\source\repos\ImarisWriter\interfaceC"
# Might need to focus this on bpImageConverterInterfaceC.h or include <stdbool.h> before loading bpConverterTypesC.h
const IMARISWRITER_HEADERS = [joinpath(IMARISWRITER_INCLUDE, header) for header in readdir(IMARISWRITER_INCLUDE) if endswith(header, ".h")]

wc = init(; headers = IMARISWRITER_HEADERS,
            output_file = joinpath(@__DIR__, "imariswriter_api.jl"),
            common_file = joinpath(@__DIR__, "imariswriter_common.jl"),
            clang_includes = vcat(IMARISWRITER_INCLUDE, CLANG_INCLUDE),
            clang_args = ["-I", joinpath(IMARISWRITER_INCLUDE, "..")],
            header_wrapped = (root, current)->root == current,
            header_library = x->"bpImarisWriter96",
            clang_diagnostics = true,
            )

run(wc)


text:
heaviside(x::T) where {T<:Number} = ifelse(x <= 0, zero(x), ifelse(x > 0, one(x), oftype(x,0)))
heaviside0(x::T) where {T<:Number} = ifelse(x < 0, zero(x), ifelse(x >= 0, one(x), oftype(x, 0)))

func_name1 = Symbol(string("heaviside"))

# generate heaviside functions for variably sized arrays
for t in 1:5
    @eval function $(func_name1)(v::AbstractArray{T, $t}) where T <: Number
        h = similar(v)
        for i in eachindex(v)
            h[i] = heaviside(v[i])
        end
        return(h)
    end
end

@eval function $(func_name1)(v::Vector{T}) where T <: Number
    h = similar(v)
    for i in eachindex(v)
        h[i] = heaviside(v[i])
    end
    return(h)
end

@eval function $(func_name1)(v::Adjoint{T}) where T <: Number
    h = similar(v)
    for i in eachindex(v)
        h[i] = heaviside(v[i])
    end
    return(h)
end

func_name2 = Symbol(string("heaviside0"))

# generate heaviside0 functions for variably sized arrays
for t in 1:5
    @eval function $(func_name2)(v::AbstractArray{T, $t}) where T <: Number
        h = similar(v)
        for i in eachindex(v)
            h[i] = heaviside0(v[i])
        end
        return(h)
    end
end

@eval function $(func_name2)(v::Vector{T}) where T <: Number
    h = similar(v)
    for i in eachindex(v)
        h[i] = heaviside0(v[i])
    end
    return(h)
end

@eval function $(func_name2)(v::Adjoint{T}) where T <: Number
    h = similar(v)
    for i in eachindex(v)
        h[i] = heaviside0(v[i])
    end
    return(h)
end

export heaviside, heaviside0


text:
# save symbols
compress_then_write(b::Vector{Symbol}, io) = begin
    string_meta = compress_then_write(String.(b), io)
    (type = Symbol, len = string_meta.len, metadata = string_meta)
end

# load a Symbol column
column_loader(::Type{Symbol}, io, metadata) = begin
    strs = column_loader(String, io, metadata.metadata)
    Symbol.(strs)
end

column_loader!(_, ::Type{Symbol}, io, metadata) = begin
    column_loader(Symbol, io, metadata)
end

some_elm(::Type{Symbol}) = :JDF


text:
######################################################################
# batch.jl
# writing my own batch sim code
######################################################################

using FEBOL

export my_batch
#export my_parsim

# for running batches of simulations on a vector of sim units
function my_batch(m::SearchDomain, suv, n_sims::Int)

    theta_start = m.theta
    x0 = suv[1].x.x
    y0 = suv[1].x.y
    h0 = suv[1].x.heading

    costs = zeros(n_sims, length(suv))
    nmacs = zeros(n_sims, length(suv))
    d2_errors = zeros(n_sims, length(suv))
    for sim_ind = 1:n_sims

        # reset jammer to random location

        print("Simulation ", sim_ind, ": ")
        for (su_ind, su) in enumerate(suv)
            m.theta = theta_start
            print(su_ind, ",")

            # reset the filter, vehicle, and policy
            #reset!(m, su)
            reset!(su.f)

            # reset vehicle
            reset!(m, su.x)
            su.x.x = 100.0
            su.x.y = 10.0
            su.x.heading = 0.0
            #su.x.x = x0
            #su.x.y = y0
            #su.x.heading = h0

            # does policy need to be reset?

            c1, c2, c3 = my_simulate(m, su)
            costs[sim_ind, su_ind] = c1
            nmacs[sim_ind, su_ind] = c2
            d2_errors[sim_ind, su_ind] = c3
        end
        println("complete.")
    end
    return costs, nmacs, d2_errors
end

function my_simulate(m::SearchDomain, su::SimUnit)

    # Don't reset the filter, vehicle, and policy
    # I think I assume the SimUnit comes in clean and ready to go

    # What was the cost to getting this first observation?
    #cost_sum = get_cost(su, m)
	cost_sum = 0.0			# really just sum of entropy
    nmac_cost = 0.0
    error_cost = 0.0
	total_cost = 0.0		# sum of entropy and collision costs, including lambda

    # before doing anything else, we observe
    #  and update filter once
    o = observe(m, su.x)
    update!(su, o)

    # This was our first step; steps count number of observations
    step_count = 1

    while !is_complete(su.f, su.tc, step_count)
        # act
        a = FEBOL.action(m, su, o)
        act!(m, su.x, a)

        move_target!(m)

        # is there an nmac?
        dx = su.x.x - m.theta[1]
        dy = su.x.y - m.theta[2]
        if sqrt(dx*dx + dy*dy) < 15.0
            nmac_cost += 1.0
			total_cost += su.p.bmdp.lambda
        end

        # observe and update
        o = observe(m, su.x)
        update!(su, o)

        # get cost and update step count
		ec = get_cost(su, m, a)
        cost_sum += ec
		total_cost += ec
        step_count += 1

        if step_count > 20
            c = centroid(su.f)
            dx = c[1] - m.theta[1]
            dy = c[2] - m.theta[2]
            error_cost += sqrt(dx*dx + dy*dy)
        end
    end

    return cost_sum, nmac_cost, error_cost, total_cost
end
# like simulation, but special for parallel simulations
# main change is deep-copying of arguments so as to not disturb them
function psim(idx::Int, m::SearchDomain, vsu::Vector{SimUnit}, theta_start::NTuple{4, Float64})

    # copy these so they don't get messed up by other cores
    m = deepcopy(m)
    vsu = deepcopy(vsu)

    # change the target location
    theta!(m)

    costs = zeros(length(vsu))
    nmacs = zeros(length(vsu))
    d2s = zeros(length(vsu))
    totals = zeros(length(vsu))

    print("Simulation ", idx, ": ")
    for (su_ind, su) in enumerate(vsu)
        m.theta = theta_start
        print(su_ind, ",")

        # reset the filter, vehicle, and policy
        reset!(m, su)
        reset!(su.f)

        # reset vehicle
        reset!(m, su.x)
        su.x.x = 100.0
        su.x.y = 10.0
        su.x.heading = 0.0

        c1, c2, c3, c4 = my_simulate(m, su)
        costs[su_ind] = c1
        nmacs[su_ind] = c2
        d2s[su_ind] = c3
        totals[su_ind] = c4
    end
    println("complete.")
    return costs, nmacs, d2s, totals
end

function my_parsim(m::SearchDomain, vsu::Vector{SimUnit}, n_sims::Int)
    np = nprocs()       # number of processes available
    i = 1
    nextidx() = (idx=i; i+=1; idx)
    #costs = zeros(n_sims, length(vsu))

    costs = zeros(n_sims, length(vsu))
    nmacs = zeros(n_sims, length(vsu))
    d2_errors = zeros(n_sims, length(vsu))
    totals = zeros(n_sims, length(vsu))

    theta_start = m.theta
    @sync begin
        for p=1:np
            if p != myid() || np == 1
                @async begin
                    while true
                        idx = nextidx()
                        idx > n_sims && break
                        c1,c2,c3,c4=remotecall_fetch(psim, p, idx, m, vsu, theta_start)
                        costs[idx,:] = c1
                        nmacs[idx,:] = c2
                        d2_errors[idx,:] = c3
                        totals[idx,:] = c4
                    end
                end
            end
        end
    end
    return costs, nmacs, d2_errors, totals
end
function my_parsim(m::SearchDomain, su::SimUnit, n_sims::Int)
    return vec( my_parsim(m, [su], n_sims) )
end


text:
@testset "Record Options" begin
    # helper to get debug as string
    io = IOBuffer()
    M = Euclidean(2)
    x = RnPoint([4.,2.])
    o = GradientDescentOptions(x,stopAfterIteration(20), ConstantStepsize(1.))
    f = y -> distance(M,y,x).^2
    ∇f = y -> -2*log(M,y,x)
    p = GradientProblem(M,f,∇f)
    a = RecordIteration()
    # constructors
    rO = RecordOptions(o,a)
    @test getOptions(o) == o
    @test getOptions(rO) == o
    @test_throws ErrorException getOptions(p)
    #
    @test getInitialStepsize(p,rO) == 1.
    @test getStepsize!(p,rO,1) == 1.
    @test getLastStepsize(p,rO,1) == 1.
    #
    @test rO.recordDictionary[:All] == a
    @test RecordOptions(o,[a]).recordDictionary[:All].group[1] == a
    @test RecordOptions(o,Dict(:A => a)).recordDictionary[:A] == a
    @test isa( RecordOptions(o, [:Iteration]).recordDictionary[:All].group[1], RecordIteration)
    @test !hasRecord(o)
    @test !hasRecord(DebugOptions(o,[]))
    @test hasRecord(rO)
    @test_throws ErrorException getRecord(o)
    @test length( getRecord(rO,:All) ) == 0
    @test length( getRecord(rO) ) == 0
    @test length(  getRecord(DebugOptions( rO,[]) )  ) == 0
    @test_throws ErrorException getRecord(RecordOptions(o,Dict{Symbol,RecordAction}()))
    @test_throws ErrorException getRecord(o)
    @test getRecord(rO) == Array{Int64,1}()
    # RecordIteration
    @test a(p,o,0) == nothing # inactive
    @test a(p,o,1) == [1]
    @test a(p,o,2) == [1,2]
    @test a(p,o,9) == [1,2,9]
    @test a(p,o,-1) == []
    # RecordGroup
    @test length(RecordGroup().group) == 0
    b = RecordGroup([ RecordIteration(), RecordIteration() ])
    b(p,o,1)
    b(p,o,2)
    @test b.group[1].recordedValues == [1,2]
    @test b.group[2].recordedValues == [1,2]
    @test getRecord(b) == [ (1,1), (2,2) ]
    #RecordEvery
    c = RecordEvery(a,10,true)
    @test c(p,o,0) == nothing
    @test c(p,o,1) == nothing
    @test c(p,o,10) == [10]
    @test c(p,o,20) == [10,20]
    @test c(p,o,-1) == []
    @test getRecord(c) == []
    # RecordChange
    c = RecordChange()
    c(p,o,1)
    @test c.recordedValues == [0.] # no x0 -> assume x0 is the first iterate
    o.x = RnPoint([3.,2.])
    c(p,o,2)
    @test c.recordedValues == [0.,1.] # no x0 -> assume x0 is the first iterate
    c = RecordChange(RnPoint([4.,2.]))
    c(p,o,1)
    @test c.recordedValues == [1.] # no x0 -> assume x0 is the first iterate
    # RecordEntry
    o.x = x
    d = RecordEntry(x,:x)
    d2 = RecordEntry(RnPoint{Float64},:x)
    @test d.field == d2.field
    d(p,o,1)
    @test d.recordedValues == [x]
    # RecordEntryChange
    o.x = x
    e = RecordEntryChange(:x, (p,o,x,y) -> distance(p.M,x,y))
    @test updateStorage!(e.storage,o) == (:x,)
    e2 = RecordEntryChange(x,:x,(p,o,x,y) -> distance(p.M,x,y))
    @test e.field == e2.field
    e(p,o,1)
    @test e.recordedValues == [0.]
    o.x = RnPoint([3.0,2.0])
    e(p,o,2)
    @test e.recordedValues == [0.,1.]
    # RecordIterate
    o.x = x
    f = RecordIterate(x)
    f2 = RecordIterate(RnPoint{Float64})
    @test_throws ErrorException RecordIterate()
    @test f.recordedValues == f2.recordedValues
    f(p,o,1)
    @test f.recordedValues == [x]
    # RecordCost
    g = RecordCost()
    g(p,o,1)
    @test g.recordedValues == [0.]
    o.x = RnPoint([3.,2.])
    g(p,o,2)
    @test g.recordedValues == [0.,1.]
    #RecordFactory
    o.∇ = RnTVector([0.,0.])
    rf = RecordFactory(o,[:Cost,:∇])
    @test isa(rf[:All], RecordGroup)
    @test isa(rf[:All].group[1], RecordCost)
    @test isa(rf[:All].group[2], RecordEntry)
    @test isa( RecordFactory(o,[2])[:All], RecordEvery)
    @test rf[:All].group[2].field == :∇
    @test length(rf[:All].group) == 2
    @test all(  isa.( RecordFactory(o,[:Cost, :Iteration, :Change, :Iterate])[:All].group,
         [RecordCost, RecordIteration, RecordChange, RecordIterate]) )
    @test RecordActionFactory(o,g) == g
end


text:
using LinearAlgebra

A = [1 1 1 1;
     2 3 5 4;
     4 9 25 16;
     8 27 125 64]

function vandermonde(A)
    det = 1
    for j in 2:4
        for i in 1:j-1
            det *= (A[2,j]-A[2,i])
        end
    end

    det
end

println("vandermonde: ", vandermonde(A))
println("det(A): ", det(A))


text:
# Copyright (c) 2018 Maarten Pronk, MIT license
# `https://github.com/evetion/GeoArrays.jl/blob/master/src/crs.jl`
# https://github.com/evetion/GeoArrays.jl/blob/master/src/geoutils.jl

# import GeoArrays: epsg2wkt, proj2wkt, wkt2wkt, str2wkt;
# import GeoArrays: get_affine_map, geotransform_to_affine, affine_to_geotransform,
#     is_rotated, unitrange_to_affine;

st_affine(ga::AbstractGeoArray) = ga.f

"""Check wether the AffineMap of a AbstractGeoArray contains rotations."""
function is_rotated(ga::AbstractGeoArray)
    ga.f.linear[2] != 0. || ga.f.linear[3] != 0.
end

"""
    st_crs!(ga::AbstractGeoArray, crs::WellKnownText{GeoFormatTypes.CRS,<:AbstractString})
    st_crs!(ga::AbstractGeoArray, crs::GFT.CoordinateReferenceSystemFormat)
    st_crs!(ga::AbstractGeoArray, epsgcode::Int)

Set CRS on AbstractGeoArray by epsgcode or string
"""
function st_crs!(ga::AbstractGeoArray, crs::WellKnownText)
    ga.crs = crs
    ga
end

function st_crs!(ga::AbstractGeoArray, crs::GFT.CoordinateReferenceSystemFormat)
    ga.crs = convert(GFT.WellKnownText, GFT.CRS(), crs)
    ga
end

st_crs!(ga::AbstractGeoArray, epsgcode::Int) = st_crs!(ga, EPSG(epsgcode))

st_crs(ga::AbstractGeoArray) = ga.crs

epsg!(ga::AbstractGeoArray, epsgcode::Int) = st_crs!(ga, EPSG(epsgcode))
epsg!(ga::AbstractGeoArray, epsgstring::AbstractString) = st_crs!(ga, EPSG(epsgstring))

crs! = st_crs!

export st_crs, st_crs!, st_affine


text:
module NaiveBayes
    using Pluto

    run_notebook() = Pluto.run(notebook=joinpath(@__DIR__, "../notebook/naive_bayes.jl"))
end


text:
# NOTE, user variables and manifolds will require the same definitions, TODO perhaps add into `@defVariable`
# unusual definitions, but all they need to do is pack and unpack as one-to-one
# this step actually occurs separate from the actual variables or factors (with their own manifolds)
# relies on later use of getManifold to give back the same <:AbstractManifold
# NOTE added to DFG.@defVariable
getVariableType(M::Euclidean{Tuple{N}}) where N = ContinuousEuclid(N)
getVariableType(M::TranslationGroup{Tuple{N}}) where N = ContinuousEuclid(N)

# getVariableType(M::RealCircleGroup) = Circular()
# getVariableType(M::Circle) = error("Circle manifold is deprecated use RealCircleGroup, will come back when we generalize to non-group Riemannian")

# Type converters for MKD
Base.convert(::Type{<:SamplableBelief}, ::Type{<:PackedManifoldKernelDensity}) = ManifoldKernelDensity
Base.convert(::Type{<:PackedSamplableBelief}, ::Type{<:ManifoldKernelDensity}) = PackedManifoldKernelDensity

"""
    $SIGNATURES

Parching refers to drying/hollowing out the object to reduce its memory size to the bare minimum.

Notes
- Likely to be used in combination with [stashing](@ref section_stash_unstash) where large data blobs are independently stored.
- For example, a point cloud stored as a MKD is very large and likely to duplicate the already stored point cloud object values,
  - When storing the MKD object, it might make sense to parch the MKD first and just persisting the context of the data.
  - Reconstituting a full MKD object would then require the inverse, where the parched shell is sized and filled from a separate large data blob.
"""
function parchDistribution(mkd::ManifoldKernelDensity)
  pts = getPoints(mkd)
  bw = getBW(mkd)[:,1]
  manikde!(mkd.manifold, pts[1:1], mkd._u0; bw, partial=mkd._partial, infoPerCoord=mkd.infoPerCoord)
end

# Data converters for MKD
function packDistribution( mkd::ManifoldKernelDensity )
  #
  pts = getPoints(mkd)

  PackedManifoldKernelDensity(
    "IncrementalInference.PackedManifoldKernelDensity",
    # piggy back on InferenceVariable serialization rather than try serialize anything Manifolds.jl
    DFG.typeModuleName(getVariableType(mkd.manifold)),
    [AMP.makeCoordsFromPoint(mkd.manifold, pt) for pt in pts],
    getBW(mkd.belief)[:,1],
    mkd._partial isa Nothing ? collect(1:manifold_dimension(mkd.manifold)) : mkd._partial ,
    mkd.infoPerCoord
  )
end

function unpackDistribution(dtr::PackedManifoldKernelDensity)
  # find InferenceVariable type from string (anything Manifolds.jl?)
  M = DFG.getTypeFromSerializationModule(dtr.varType) |> getManifold
  vecP = [AMP.makePointFromCoords(M, pt) for pt in dtr.pts]
  bw = length(dtr.bw) === 0 ? nothing : dtr.bw
  partial = length(dtr.partial) == manifold_dimension(M) || length(dtr.partial) === 0 ? nothing : dtr.partial

  return manikde!( M, vecP; bw, partial, infoPerCoord=dtr.infoPerCoord )
end

function Base.convert(::Type{String},
                      mkd::ManifoldKernelDensity )
  #
  packedMKD = packDistribution(mkd)
  JSON2.write(packedMKD)
end

# Use general dispatch
# Base.convert(::Type{<:PackedSamplableBelief}, mkd::ManifoldKernelDensity) = convert(String, mkd)

# make module specific
# good references:
#  https://discourse.julialang.org/t/converting-string-to-datatype-with-meta-parse/33024/2
#  https://discourse.julialang.org/t/is-there-a-way-to-import-modules-with-a-string/15723/6
function Base.convert(::Type{<:ManifoldKernelDensity}, str::AbstractString)
  dtr = JSON2.read(str, PackedManifoldKernelDensity)
  unpackDistribution(dtr)
end

#


text:
type BiMapInternal{Bool,L,R}
    left::Dict{L,R}
    right::Dict{R,L}
end

haskey(bimap::BiMapInternal{true},  lkey) = in(lkey, keys(bimap.left))
haskey(bimap::BiMapInternal{false}, rkey) = in(rkey, keys(bimap.right))

isempty(bimap::BiMapInternal) = isempty(bimap.left)
length(bimap::BiMapInternal) = length(bimap.left)

getindex(bimap::BiMapInternal{true},  lkey) = bimap.left[lkey]
getindex(bimap::BiMapInternal{false}, rkey) = bimap.right[rkey]

function setindex!(bimap::BiMapInternal{true}, rkey, lkey)
    # TODO: atomic
    if haskey(bimap.left, lkey)
        delete!(bimap.right, bimap.left[lkey])
    end
    if haskey(bimap.right, rkey)
        delete!(bimap.left, bimap.right[rkey])
    end
    bimap.left[lkey] = rkey
    bimap.right[rkey] = lkey
    rkey
end

function setindex!(bimap::BiMapInternal{false}, lkey, rkey)
    # TODO: atomic
    if haskey(bimap.left, lkey)
        delete!(bimap.right, bimap.left[lkey])
    end
    if haskey(bimap.right, rkey)
        delete!(bimap.left, bimap.right[rkey])
    end
    bimap.left[lkey] = rkey
    bimap.right[rkey] = lkey
    lkey
end

start(bimap::BiMapInternal{true}) = start(bimap.left)
done(bimap::BiMapInternal{true}, s) = done(bimap.left, s)
next(bimap::BiMapInternal{true}, s) = next(bimap.left, s)


text:
@testset "Simulator: simulate" begin
    ## hopper
    q1 = nominal_configuration(hopper)
    v1 = zeros(hopper.nq)

    h = 0.01
    T = 100

    # simulator
    s = Simulator(hopper, T, h=h, diff_sol=true, sim_opts=RoboDojo.SimulatorOptions(record=true))

    # step
    q2 = step!(s, q1, v1, zeros(hopper.nu), 1)
    @test sum(q2) != 0.0

    # simulate
    status = simulate!(s, q1, v1, reset_traj=true)
    @test status

    ## quadruped
    q1 = nominal_configuration(quadruped)
    v1 = zeros(quadruped.nq)

    h = 0.01
    T = 100

    # simulator
    s = Simulator(quadruped, T, h=h, diff_sol=true, sim_opts=RoboDojo.SimulatorOptions(record=true))

    # step
    q2 = step!(s, q1, v1, zeros(quadruped.nu), 1)
    @test sum(q2) != 0.0

    # simulate
    status = simulate!(s, q1, v1, reset_traj=true)
    @test status

    ## halfquadruped
    q1 = nominal_configuration(halfquadruped)
    v1 = zeros(halfquadruped.nq)

    h = 0.01
    T = 100

    # simulator
    s = Simulator(halfquadruped, T, h=h, diff_sol=true, sim_opts=RoboDojo.SimulatorOptions(record=true))

    # step
    q2 = step!(s, q1, v1, zeros(halfquadruped.nu), 1)
    @test sum(q2) != 0.0

    # simulate
    status = simulate!(s, q1, v1, reset_traj=true)
    @test status

    ## halfcheetah
    q1 = nominal_configuration(halfcheetah)
    v1 = zeros(halfcheetah.nq)

    h = 0.01
    T = 100

    # simulator
    s = Simulator(halfcheetah, T, h=h, diff_sol=true, sim_opts=RoboDojo.SimulatorOptions(record=true))

    # step
    q2 = step!(s, q1, v1, zeros(halfcheetah.nu), 1)
    @test sum(q2) != 0.0

    # simulate
    status = simulate!(s, q1, v1, reset_traj=true)
    @test status

    ## biped
    q1 = nominal_configuration(biped)
    v1 = zeros(biped.nq)

    h = 0.01
    T = 100

    # simulator
    s = Simulator(biped, T, h=h, diff_sol=true, sim_opts=RoboDojo.SimulatorOptions(record=true))

    # step
    q2 = step!(s, q1, v1, zeros(biped.nu), 1)
    @test sum(q2) != 0.0

    # simulate
    status = simulate!(s, q1, v1, reset_traj=true)
    @test status

    # process timing results
    RoboDojo.process!(s.stats, 1)
    @test sum(s.stats.policy_time) > 0.0
    @test s.stats.policy_mean[1] > 0.0

    @test sum(s.stats.sim_time) > 0.0
    @test s.stats.sim_mean[1] > 0.0

    ## box
    q1 = nominal_configuration(box)
    v1 = zeros(box.nq)

    h = 0.01
    T = 100

    # simulator
    s = Simulator(box, T, h=h, diff_sol=true, sim_opts=RoboDojo.SimulatorOptions(record=true))

    # step
    q2 = step!(s, q1, v1, zeros(box.nu), 1)
    @test sum(q2) != 0.0

    # simulate
    status = simulate!(s, q1, v1, reset_traj=true)
    @test status

    # process timing results
    RoboDojo.process!(s.stats, 1)
    @test sum(s.stats.policy_time) > 0.0
    @test s.stats.policy_mean[1] > 0.0

    @test sum(s.stats.sim_time) > 0.0
    @test s.stats.sim_mean[1] > 0.0

    ## particle
    q1 = nominal_configuration(particle)
    v1 = zeros(particle.nq)

    h = 0.01
    T = 100

    # simulator
    s = Simulator(particle, T, h=h, diff_sol=true, sim_opts=RoboDojo.SimulatorOptions(record=true))

    # step
    q2 = step!(s, q1, v1, zeros(particle.nu), 1)
    @test sum(q2) != 0.0

    # simulate
    status = simulate!(s, q1, v1, reset_traj=true)
    @test status

    # process timing results
    RoboDojo.process!(s.stats, 1)
    @test sum(s.stats.policy_time) > 0.0
    @test s.stats.policy_mean[1] > 0.0

    @test sum(s.stats.sim_time) > 0.0
    @test s.stats.sim_mean[1] > 0.0

    ## particle (time budget)
    q1 = nominal_configuration(particle)
    v1 = zeros(particle.nq)

    h = 0.01
    T = 100

    # simulator
    s = Simulator(particle, T, h=h, diff_sol=true,
        sim_opts=RoboDojo.SimulatorOptions(record=true),
        solver_opts=solver_opts=RoboDojo.InteriorPointOptions(
            max_time=0.49,
            undercut=Inf,
            γ_reg=0.1,
            r_tol=1e-8,
            κ_tol=1e-8,
            max_ls=25,
            ϵ_min=0.25,
            diff_sol=false,
            verbose=false),)

    # step
    q2 = step!(s, q1, v1, zeros(particle.nu), 1)
    @test sum(q2) != 0.0

    # simulate
    status = simulate!(s, q1, v1, reset_traj=true)
    @test status

     # simulator
     s = Simulator(particle, T, h=h, diff_sol=true,
     sim_opts=RoboDojo.SimulatorOptions(record=true),
     solver_opts=solver_opts=RoboDojo.InteriorPointOptions(
         max_time=1.0e-6,
         undercut=Inf,
         γ_reg=0.1,
         r_tol=1e-8,
         κ_tol=1e-8,
         max_ls=25,
         ϵ_min=0.25,
         diff_sol=false,
         verbose=false),)

    # simulate
    status = simulate!(s, q1, v1, reset_traj=true)
    @test !status

    ## centroidal quadruped
    q1 = nominal_configuration(centroidal_quadruped)
    v1 = zeros(centroidal_quadruped.nq)

    h = 0.01
    T = 100

    # simulator
    s = Simulator(centroidal_quadruped, T, h=h, diff_sol=true, sim_opts=RoboDojo.SimulatorOptions(record=true))

    # step
    q2 = step!(s, q1, v1, zeros(centroidal_quadruped.nu), 1)
    @test sum(q2) != 0.0

    # simulate
    status = simulate!(s, q1, v1, reset_traj=true)
    @test status

    # process timing results
    RoboDojo.process!(s.stats, 1)
    @test sum(s.stats.policy_time) > 0.0
    @test s.stats.policy_mean[1] > 0.0

    @test sum(s.stats.sim_time) > 0.0
    @test s.stats.sim_mean[1] > 0.0
end


text:
# https://docs.julialang.org/en/stable/manual/integers-and-floating-point-numbers/
println(1 |> UInt8)
println(1 |> Int128)
println(1. |> Float16)
println(typeof(1))
println(Sys.WORD_SIZE)
println(typemin(UInt8), " ", typemax(UInt8))

println(bitstring(0.0))
println(eps(Float64))
println(eps(BigFloat))
println(nextfloat(0.0))

println(factorial(BigInt(40)))
println(parse(BigInt, "123123123123"))

x = 3
println(2x^2 -3x + 1)
println(2(x-1)^2)
println((x-1)x)

zero(Float64)
one(BigFloat)

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
using COSMO
using Distributions
using Plots
using PortfolioOpt
using PortfolioOpt.TestUtils:
    backtest_po, get_test_data, mean_variance,
    percentchange, timestamp, rename!

############ Read Prices #############
prices = get_test_data()
numD, numA = size(prices) # A: Assets    D: Days

############# Calculating returns #####################
returns_series = percentchange(prices)

# risk free asset
# rf =  fill(3.2e-4, numD)

############# backtest Parameters #####################
DEFAULT_SOLVER = optimizer_with_attributes(
    COSMO.Optimizer, "verbose" => false, "max_iter" => 900000
)

start_date = timestamp(returns_series)[100]

## Strategies
# Look back days
k_back = 60

############# backtest with limit return strategies #####################

wealth_markowitz_limit_R, strategy_returns =
    backtest_po(returns_series; start_date=start_date) do past_returns, current_wealth, rf
        # Prep
        numD, numA = size(past_returns)
        returns = values(past_returns)
        Σ, r̄ = mean_variance(returns[(end - k_back):end, :])

        # Parameters
        R = 0.0015 / current_wealth
        formulation = MeanVariance(;
            predicted_mean = r̄,
            predicted_covariance = Σ,
        )

        # model
        model = po_min_variance_limit_return(formulation, R; rf = rf)
        x = compute_solution(model, DEFAULT_SOLVER)
        return x * current_wealth
    end;
rename!(wealth_markowitz_limit_R, :markowitz_limit_R);

wealth_markowitz_infeasible_limit_R, strategy_returns =
    backtest_po(returns_series; start_date=start_date) do past_returns, current_wealth, rf
        # Prep
        numD, numA = size(past_returns)
        returns = values(past_returns)
        Σ, r̄ = mean_variance(returns[(end - k_back):end, :])

        # Parameters
        R = 0.005 / current_wealth
        formulation = MeanVariance(;
            predicted_mean = r̄,
            predicted_covariance = Σ,
        )

        # model
        model = po_min_variance_limit_return(formulation, R; rf = rf)
        x = compute_solution(model, DEFAULT_SOLVER)
        return x * current_wealth
    end;
rename!(wealth_markowitz_infeasible_limit_R, :markowitz_infeasible_limit_R);

wealth_soyster_limit_R, strategy_returns =
    backtest_po(returns_series; start_date=start_date) do past_returns, current_wealth, rf
        # Prep
        numD, numA = size(past_returns)
        returns = values(past_returns)
        Σ, r̄ = mean_variance(returns[(end - k_back):end, :])

        # Parameters
        R = 0.0015 / current_wealth
        formulation = RobustBertsimas(;
            predicted_mean = r̄,
            predicted_covariance = Σ,
            uncertainty_delta = std(returns[(end - k_back):end, :]; dims=1)'[:,1] / 5,
            bertsimas_budget = Float64(numA),
        )

        # model
        model = po_min_variance_limit_return(formulation, R; rf = rf)
        x = compute_solution(model, DEFAULT_SOLVER)
        return x * current_wealth
    end;
rename!(wealth_soyster_limit_R, :soyster_limit_R);

wealth_bertsimas_limit_R, strategy_returns =
    backtest_po(returns_series; start_date=start_date) do past_returns, current_wealth, rf
        # Prep
        numD, numA = size(past_returns)
        returns = values(past_returns)
        Σ, r̄ = mean_variance(returns[(end - k_back):end, :])

        # Parameters
        R = 0.0015 / current_wealth
        formulation = RobustBertsimas(;
            predicted_mean = r̄,
            predicted_covariance = Σ,
            uncertainty_delta = std(returns[(end - k_back):end, :]; dims=1)'[:,1] / 5,
            bertsimas_budget = 3.0,
        )

        # model
        model = po_min_variance_limit_return(formulation, R; rf = rf)
        x = compute_solution(model, DEFAULT_SOLVER)
        return x * current_wealth
    end;
rename!(wealth_bertsimas_limit_R, :bertsimas_limit_R);

wealth_bental_limit_R, strategy_returns =
backtest_po(returns_series; start_date=start_date) do past_returns, current_wealth, rf
        # Prep
        numD, numA = size(past_returns)
        returns = values(past_returns)
        Σ, r̄ = mean_variance(returns[(end - k_back):end, :])

        # Parameters
        R = 0.0015 / current_wealth
        formulation = RobustBenTal(;
            predicted_mean = r̄,
            predicted_covariance = Σ,
            uncertainty_delta = 0.025
        )

        # model
        model = po_min_variance_limit_return(formulation, R; rf = rf)
        x = compute_solution(model, DEFAULT_SOLVER)
        return x * current_wealth
    end;
rename!(wealth_bental_limit_R, :bental_limit_R);

############# backtest with limit variance strategies #####################

wealth_markowitz_limit_var, returns_markowitz_limit_var =
backtest_po(returns_series; start_date=start_date) do past_returns, current_wealth, rf
        # Prep
        numD, numA = size(past_returns)
        returns = values(past_returns)
        Σ_s, r̄_s = mean_variance(returns[(end - 14):end, :])
        Σ, r̄ = mean_variance(returns[(end - k_back):end, :])

        # Parameters
        max_risk = 0.8
        formulation = MeanVariance(;
            predicted_mean = r̄_s,
            predicted_covariance = Σ,
        )

        # model
        model = po_max_return_limit_variance(formulation, max_risk; rf = rf)
        x = compute_solution(model, DEFAULT_SOLVER)
        return x * current_wealth
    end;
rename!(wealth_markowitz_limit_var, :markowitz_limit_var);

wealth_soyster_limit_var, returns_soyster_limit_var =
    backtest_po(returns_series; start_date=start_date) do past_returns, current_wealth, rf
        # Prep
        numD, numA = size(past_returns)
        returns = values(past_returns)
        Σ_s, r̄_s = mean_variance(returns[(end - 14):end, :])
        Σ, r̄ = mean_variance(returns[(end - k_back):end, :])

        # Parameters
        max_risk = 0.8
        formulation = RobustBertsimas(;
            predicted_mean = r̄_s,
            predicted_covariance = Σ,
            uncertainty_delta = std(returns[(end - k_back):end, :]; dims=1)'[:,1] / 3,
            bertsimas_budget = Float64(numA),
        )

        # model
        model = po_max_return_limit_variance(formulation, max_risk; rf = rf)
        x = compute_solution(model, DEFAULT_SOLVER)
        return x * current_wealth
    end;
rename!(wealth_soyster_limit_var, :soyster_limit_var);

wealth_bertsimas_2_limit_var, returns_bertsimas_2_limit_var =
    backtest_po(returns_series; start_date=start_date) do past_returns, current_wealth, rf
        # Prep
        numD, numA = size(past_returns)
        returns = values(past_returns)
        Σ_s, r̄_s = mean_variance(returns[(end - 14):end, :])
        Σ, r̄ = mean_variance(returns[(end - k_back):end, :])

        # Parameters
        max_risk = 0.8
        formulation = RobustBertsimas(;
            predicted_mean = r̄_s,
            predicted_covariance = Σ,
            uncertainty_delta = std(returns[(end - k_back):end, :]; dims=1)'[:,1] / 3,
            bertsimas_budget = 2.0,
        )

        # model
        model = po_max_return_limit_variance(formulation, max_risk; rf = rf)
        x = compute_solution(model, DEFAULT_SOLVER)
        return x * current_wealth
    end;
rename!(wealth_bertsimas_2_limit_var, :bertsimas_2_limit_var);

wealth_bertsimas_4_limit_var, returns_bertsimas_4_limit_var =
    backtest_po(returns_series; start_date=start_date) do past_returns, current_wealth, rf
        # Prep
        numD, numA = size(past_returns)
        returns = values(past_returns)
        Σ_s, r̄_s = mean_variance(returns[(end - 14):end, :])
        Σ, r̄ = mean_variance(returns[(end - k_back):end, :])

        # Parameters
        max_risk = 0.8
        formulation = RobustBertsimas(;
            predicted_mean = r̄_s,
            predicted_covariance = Σ,
            uncertainty_delta = std(returns[(end - k_back):end, :]; dims=1)'[:,1] / 3,
            bertsimas_budget = 4.0,
        )

        # model
        model = po_max_return_limit_variance(formulation, max_risk; rf = rf)
        x = compute_solution(model, DEFAULT_SOLVER)
        return x * current_wealth
    end;
rename!(wealth_bertsimas_4_limit_var, :bertsimas_4_limit_var);

wealth_bental_limit_var, returns_bental_limit_var =
    backtest_po(returns_series; start_date=start_date) do past_returns, current_wealth, rf
        # Prep
        numD, numA = size(past_returns)
        returns = values(past_returns)
        Σ_s, r̄_s = mean_variance(returns[(end - 14):end, :])
        Σ, r̄ = mean_variance(returns[(end - k_back):end, :])

        # Parameters
        max_risk = 0.8
        formulation = RobustBenTal(;
            predicted_mean = r̄_s,
            predicted_covariance = Σ,
            uncertainty_delta = 0.028
        )

        # model
        model = po_max_return_limit_variance(formulation, max_risk; rf = rf)
        x = compute_solution(model, DEFAULT_SOLVER)
        return x * current_wealth
    end;
rename!(wealth_bental_limit_var, :bental_limit_var);

############# backtest with Data-Driven RO strategies #####################
wealth_betina, returns_betina =
    backtest_po(returns_series; start_date=start_date) do past_returns, current_wealth, rf
        # Prep
        numD, numA = size(past_returns)
        returns = values(past_returns)
        # Parameters
        j_robust = 25
        R = -0.015 / current_wealth
        num_train = 20
        kst = 14
        klt = 30
        kmom = 20
        Q = 2

        formulation = RobustBetina(;
            sampled_returns = returns[(end - j_robust):end, :]
        )

        # Predict return
        r_bar_t = zeros(numA)
        for i in 1:numA
            r_bar_t[i] = mixed_signals_predict_return(
                returns[:, i], num_train, kst, klt, kmom, Q
            )
        end

        # model
        model = po_max_predicted_return_limit_return(formulation, R; predicted_mean = r_bar_t)
        x = compute_solution(model, DEFAULT_SOLVER)
        return x * current_wealth
    end;
rename!(wealth_betina, :betina);

############# backtest with DRO strategies #####################

wealth_delague, returns_delague_limit_var =
    backtest_po(returns_series; start_date=start_date) do past_returns, current_wealth, rf
        # Prep
        numD, numA = size(past_returns)
        returns = values(past_returns)
        Σ_s, r̄_s = mean_variance(returns[(end - 14):end, :])
        Σ, r̄ = mean_variance(returns[(end - k_back):end, :])
        # Parameters
        # a =  [0.020646446303211906; 0.01456965970660176; 0.011268181574938444; 0.009189943838602721; 0.007760228688512961; 0.006716054241470283; 0.005919840492702679; 0.005292563196132596; 0.004785576102825234; 0.004367285293323039; 0.004367285293323039]
        # b = [12351.417663268794; 11949.933298064285; 11797.83904955914; 11743.66244416854; 11734.986173090612; 11749.533056761096; 11776.549758689594; 11810.379744602318; 11847.862045043921; 11887.152715141867; 11887.152715141867]

        formulation = RobustDelague(;
            predicted_mean = r̄_s,
            predicted_covariance = Σ,
            γ1 = 0.0065,
            γ2 = 5.5,
            utility_coeficients = [1.0],
            utility_intercepts = [0.0],
        )

        model = po_max_utility_return(formulation)
        x = compute_solution(model, DEFAULT_SOLVER)
        return x * current_wealth
    end;
rename!(wealth_delague, :delague);

############# extra strategies #####################

wealth_sharpe, returns_sharpe =
    backtest_po(returns_series; start_date=start_date) do past_returns, current_wealth, rf
        # Prep
        numD, numA = size(past_returns)
        returns = values(past_returns)
        Σ, r̄ = mean_variance(returns[(end - k_back):end, :])
        # Parameters
        x = max_sharpe(Σ, r̄, rf)
        return x * current_wealth
    end;
rename!(wealth_sharpe, :sharpe);

wealth_equal_weights, returns_equal =
    backtest_po(returns_series; start_date=start_date) do past_returns, current_wealth, rf
        # Prep
        numD, numA = size(past_returns)
        # Parameters
        x = equal_weights(numA)
        return x * current_wealth
    end;
rename!(wealth_equal_weights, :equal_weights);

############# plot results  #####################
## Wealth
plt = plot(
    wealth_sharpe;
    title="Culmulative Wealth",
    xlabel="Time",
    ylabel="Wealth",
    legend=:outertopright,
);
plot!(plt, wealth_equal_weights);
# plot!(plt, wealth_markowitz_limit_R);
# plot!(plt, wealth_soyster_limit_R);
# plot!(plt, wealth_bertsimas_limit_R);
# plot!(plt, wealth_bental_limit_R);
# plot!(plt, wealth_markowitz_infeasible);
plot!(plt, wealth_markowitz_limit_var);
plot!(plt, wealth_soyster_limit_var);
plot!(plt, wealth_bertsimas_2_limit_var);
plot!(plt, wealth_bertsimas_4_limit_var);
plot!(plt, wealth_bental_limit_var);
plot!(plt, wealth_delague);
plot!(plt, wealth_betina);
plt

## Stats
plt = scatter(
    [std(returns_sharpe)],
    [mean(returns_sharpe)];
    label="Sharpe",
    title="Mean Vs Std",
    xlabel="σ",
    ylabel="E[r]",
    legend=:outertopright,
);
scatter!(plt, [std(returns_equal)], [mean(returns_equal)]; label="Equal");
scatter!(
    plt,
    [std(returns_markowitz_limit_var)],
    [mean(returns_markowitz_limit_var)];
    label="Markowitz",
);
scatter!(
    plt,
    [std(returns_soyster_limit_var)],
    [mean(returns_soyster_limit_var)];
    label="Soyster",
);
scatter!(
    plt,
    [std(returns_bertsimas_2_limit_var)],
    [mean(returns_bertsimas_2_limit_var)];
    label="Bertsimas 2",
);
scatter!(
    plt,
    [std(returns_bertsimas_4_limit_var)],
    [mean(returns_bertsimas_4_limit_var)];
    label="Bertsimas 4",
);
scatter!(
    plt, [std(returns_bental_limit_var)], [mean(returns_bental_limit_var)]; label="Bental"
);
scatter!(
    plt,
    [std(returns_delague_limit_var)],
    [mean(returns_delague_limit_var)];
    label="Delague",
);
scatter!(plt, [std(returns_betina)], [mean(returns_betina)]; label="Betina");
plt