text:
function evaluate_tc(obj::IERSConventions, arg0::AbsoluteDate)
    return jcall(obj, "evaluateTC", jdouble, (AbsoluteDate,), arg0)
end

function evaluate_tc(obj::IERSConventions, arg0::AbsoluteDate, arg1::TimeScales)
    return jcall(obj, "evaluateTC", jdouble, (AbsoluteDate, TimeScales), arg0, arg1)
end

function evaluate_tc(obj::IERSConventions, arg0::FieldAbsoluteDate)
    return jcall(obj, "evaluateTC", RealFieldElement, (FieldAbsoluteDate,), arg0)
end

function evaluate_tc(obj::IERSConventions, arg0::FieldAbsoluteDate, arg1::TimeScales)
    return jcall(obj, "evaluateTC", RealFieldElement, (FieldAbsoluteDate, TimeScales), arg0, arg1)
end

function get_earth_orientation_angle_function(obj::IERSConventions, arg0::TimeScale)
    return jcall(obj, "getEarthOrientationAngleFunction", TimeScalarFunction, (TimeScale,), arg0)
end

function get_earth_orientation_angle_function(obj::IERSConventions, arg0::TimeScale, arg1::TimeScale)
    return jcall(obj, "getEarthOrientationAngleFunction", TimeScalarFunction, (TimeScale, TimeScale), arg0, arg1)
end

function get_eop_tidal_correction(obj::IERSConventions)
    return jcall(obj, "getEOPTidalCorrection", TimeVectorFunction, ())
end

function get_eop_tidal_correction(obj::IERSConventions, arg0::TimeScales)
    return jcall(obj, "getEOPTidalCorrection", TimeVectorFunction, (TimeScales,), arg0)
end

function get_gast_function(obj::IERSConventions, arg0::TimeScale, arg1::EOPHistory)
    return jcall(obj, "getGASTFunction", TimeScalarFunction, (TimeScale, EOPHistory), arg0, arg1)
end

function get_gast_function(obj::IERSConventions, arg0::TimeScale, arg1::EOPHistory, arg2::TimeScales)
    return jcall(obj, "getGASTFunction", TimeScalarFunction, (TimeScale, EOPHistory, TimeScales), arg0, arg1, arg2)
end

function get_gmst_function(obj::IERSConventions, arg0::TimeScale)
    return jcall(obj, "getGMSTFunction", TimeScalarFunction, (TimeScale,), arg0)
end

function get_gmst_function(obj::IERSConventions, arg0::TimeScale, arg1::TimeScales)
    return jcall(obj, "getGMSTFunction", TimeScalarFunction, (TimeScale, TimeScales), arg0, arg1)
end

function get_gmst_rate_function(obj::IERSConventions, arg0::TimeScale)
    return jcall(obj, "getGMSTRateFunction", TimeScalarFunction, (TimeScale,), arg0)
end

function get_gmst_rate_function(obj::IERSConventions, arg0::TimeScale, arg1::TimeScales)
    return jcall(obj, "getGMSTRateFunction", TimeScalarFunction, (TimeScale, TimeScales), arg0, arg1)
end

function get_love_numbers(obj::IERSConventions)
    return jcall(obj, "getLoveNumbers", LoveNumbers, ())
end

function get_mean_obliquity_function(obj::IERSConventions)
    return jcall(obj, "getMeanObliquityFunction", TimeScalarFunction, ())
end

function get_mean_obliquity_function(obj::IERSConventions, arg0::TimeScales)
    return jcall(obj, "getMeanObliquityFunction", TimeScalarFunction, (TimeScales,), arg0)
end

function get_nominal_tidal_displacement(obj::IERSConventions)
    return jcall(obj, "getNominalTidalDisplacement", Vector{jdouble}, ())
end

function get_nutation_arguments(obj::IERSConventions, arg0::TimeScale)
    return jcall(obj, "getNutationArguments", FundamentalNutationArguments, (TimeScale,), arg0)
end

function get_nutation_arguments(obj::IERSConventions, arg0::TimeScale, arg1::TimeScales)
    return jcall(obj, "getNutationArguments", FundamentalNutationArguments, (TimeScale, TimeScales), arg0, arg1)
end

function get_nutation_correction_converter(obj::IERSConventions)
    return jcall(obj, "getNutationCorrectionConverter", IERSConventions_NutationCorrectionConverter, ())
end

function get_nutation_correction_converter(obj::IERSConventions, arg0::TimeScales)
    return jcall(obj, "getNutationCorrectionConverter", IERSConventions_NutationCorrectionConverter, (TimeScales,), arg0)
end

function get_nutation_function(obj::IERSConventions)
    return jcall(obj, "getNutationFunction", TimeVectorFunction, ())
end

function get_nutation_function(obj::IERSConventions, arg0::TimeScales)
    return jcall(obj, "getNutationFunction", TimeVectorFunction, (TimeScales,), arg0)
end

function get_nutation_reference_epoch(obj::IERSConventions)
    return jcall(obj, "getNutationReferenceEpoch", AbsoluteDate, ())
end

function get_nutation_reference_epoch(obj::IERSConventions, arg0::TimeScales)
    return jcall(obj, "getNutationReferenceEpoch", AbsoluteDate, (TimeScales,), arg0)
end

function get_ocean_pole_tide(obj::IERSConventions, arg0::EOPHistory)
    return jcall(obj, "getOceanPoleTide", TimeVectorFunction, (EOPHistory,), arg0)
end

function get_permanent_tide(obj::IERSConventions)
    return jcall(obj, "getPermanentTide", jdouble, ())
end

function get_precession_function(obj::IERSConventions)
    return jcall(obj, "getPrecessionFunction", TimeVectorFunction, ())
end

function get_precession_function(obj::IERSConventions, arg0::TimeScales)
    return jcall(obj, "getPrecessionFunction", TimeVectorFunction, (TimeScales,), arg0)
end

function get_solid_pole_tide(obj::IERSConventions, arg0::EOPHistory)
    return jcall(obj, "getSolidPoleTide", TimeVectorFunction, (EOPHistory,), arg0)
end

function get_tidal_displacement_frequency_correction_diurnal(obj::IERSConventions)
    return jcall(obj, "getTidalDisplacementFrequencyCorrectionDiurnal", PoissonSeries_CompiledSeries, ())
end

function get_tidal_displacement_frequency_correction_zonal(obj::IERSConventions)
    return jcall(obj, "getTidalDisplacementFrequencyCorrectionZonal", PoissonSeries_CompiledSeries, ())
end

function get_tide_frequency_dependence_function(obj::IERSConventions, arg0::TimeScale)
    return jcall(obj, "getTideFrequencyDependenceFunction", TimeVectorFunction, (TimeScale,), arg0)
end

function get_tide_frequency_dependence_function(obj::IERSConventions, arg0::TimeScale, arg1::TimeScales)
    return jcall(obj, "getTideFrequencyDependenceFunction", TimeVectorFunction, (TimeScale, TimeScales), arg0, arg1)
end

function get_xy_sp_xy2_function(obj::IERSConventions)
    return jcall(obj, "getXYSpXY2Function", TimeVectorFunction, ())
end

function get_xy_sp_xy2_function(obj::IERSConventions, arg0::TimeScales)
    return jcall(obj, "getXYSpXY2Function", TimeVectorFunction, (TimeScales,), arg0)
end

function value_of(::Type{IERSConventions}, arg0::JString)
    return jcall(IERSConventions, "valueOf", IERSConventions, (JString,), arg0)
end

function values(::Type{IERSConventions})
    return jcall(IERSConventions, "values", Vector{IERSConventions}, ())
end


text:
## Definition od different type of networks

# Network is a weighted directed graph and use the Arc type
type Network
  source::Int
  target::Int # Also called sink
  tails::Vector{Int} # Used to store the start of the arcs incident to each node
  links::Vector{Arc}
end

function findArc(g::Network,t::Int,h::Int)
  for (id,e) in enumerate(g.links[g.tails[t]:(g.tails[t+1]-1)])
    if e.head == h
      return g.tails[t] + id - 1
    end
  end
  return 0
end

function minmaxCapacity(g::Network)
  cmin, cmax = typemax(Float64), 0.
  for a in g.links
    if a.cap != 0.
      cmin = min(cmin,a.cap)
    end
    cmax = max(cmax,a.cap)
  end
  return cmin, cmax
end

ext: jl


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

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

	v_r = view(problem.atoms.r, :, :)

	Threads.@threads for i = 1:n_sensors
		one_sensors = view(all_sensors, :, i)
		intensities[i] = _get_intensity_over_sensor(problem.atoms.shape, problem.laser, v_r, one_sensors, β)
    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) where T <: ThreeD
	## Laser Pump
	E_L = (im/Γ)*apply_laser_over_oneAtom(laser, sensor)

    # ## Scattered
	E_scatt = zero(ComplexF64)
	n̂ = sensor/norm(sensor)

    dot_n_r = zero(ComplexF64)
    j = 1
	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 = E_scatt*exp(ikr)/ikr
	return abs2(E_L + E_scatt)
end

# function get_intensity_over_angle(problem::SimulationScalar, atoms_states, θ::Number; Gₙₘ=nothing)
#     if isnothing(Gₙₘ)
#         Gₙₘ = get_geometric_factor(problem.atoms, θ)
#     end

#     β = atoms_states

#     βₙₘ = transpose(β*β')
#     βₙₘ[diagind(βₙₘ)] .= abs2.(β)

#     intensity = real(sum(βₙₘ.*Gₙₘ))
# 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_scattered_intensity(problem::T, atoms_states, θ::Number; Gₙₘ=nothing) where {T<:MeanFieldProblem}
#     if isnothing(Gₙₘ)
#         Gₙₘ = view(get_geometric_factor(problem.atoms, θ),:,:)
#     end

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

#     βₙₘ = transpose(β*β')
#     βₙₘ[diagind(βₙₘ)] .= (1 .+ z)./2

#     intensity = real(sum(βₙₘ.*Gₙₘ))
#     return intensity
# end

# function get_geometric_factor(atoms, Θ)
#     N = atoms.N
#     r = atoms.r

#     Gₙₘ = zeros(ComplexF64, N,N)
#     for n=1:N
#         for m=n:N
#             rx = r[n][1] - r[m][1]
#             ry = r[n][2] - r[m][2]
#             rz = r[n][3] - r[m][3]

#             argument_bessel = sqrt( Complex( -k₀^2*(sin(Θ)^2)*(rx^2 +ry^2)) )
#             Gₙₘ[n,m] = cis(k₀*rz*cos(Θ))*besseli(0, argument_bessel )
#         end
#     end
#     return Hermitian(Gₙₘ)
# end


text:
isunitary(op) = op' * op ≈ IMatrix(size(op, 1))
isunitary(op::Number) = op' * op ≈ one(op)

isreflexive(op) = op * op ≈ IMatrix(size(op, 1))
isreflexive(op::Number) = op * op ≈ one(op)

"""
    ishermitian(op) -> Bool

check if this operator is hermitian.
"""
LinearAlgebra.ishermitian(op) = op' ≈ op

function iscommute(ops...)
    n = length(ops)
    for i = 1:n
        for j = i+1:n
            iscommute(ops[i], ops[j]) || return false
        end
    end
    true
end

iscommute(op1, op2) = op1 * op2 ≈ op2 * op1


text:
function get_profit_for_single_day(u::String="*", eval_date::String=string(today()))
    yr, m, d = split(eval_date, "-")
    tmp = glob(joinpath(Pinguu.base_dir, "users/" * u * "/*/closed_orders/$yr/$m.json"));
    profit = 0
    for file in tmp
        month = JSON.parsefile(file)
        for (day, info) in month
            if day == eval_date
                for (id_str, trade) in info["trades"]
                    profit += haskey(trade, "profit_EUR") ? trade["profit_EUR"] : trade["profit"]
                end
            end
        end
    end
    grid_bot_paths = glob(joinpath(Pinguu.base_dir, "grid_bots/user/" * u * "/*/config.json"))
    for grid_bot_path in grid_bot_paths
        tmp_config = JSON.parsefile(grid_bot_path)
        if tmp_config["testing"] == "0"
            tmp = glob(joinpath(dirname(grid_bot_path), "history/$yr/$m.json"))
            for file in tmp
                month = JSON.parsefile(file)
                for (day, info) in month
                    if day == eval_date
                        for trade in info
                            profit += trade["profit"]
                        end
                    end
                end
            end
        end
    end
    return profit
end
export get_profit_for_single_day

function get_profit_for_past_days(u::String="*", days_num=30)
    past_days = []
    profits = []
    for num in 0:1:days_num-1
        eval_date = string(today() - Day(num))
        push!(past_days, eval_date)
        push!(profits, Pinguu.get_profit_for_single_day(u, eval_date))
    end
    return past_days, profits
end
export get_profit_for_past_days

function get_profit_for_past_months(u::String="*", months_num=7)
    past_months = []
    profits = []

    for num in 0:1:months_num-1
        eval_date = string(today()-Month(num))
        yr, m, d = split(eval_date, "-")
        #
        ## Pinguus
        tmp = glob(joinpath(Pinguu.base_dir, "users/" * u * "/*/closed_orders/$yr/$m.json"));
        profit = 0
        for file in tmp
            month = JSON.parsefile(file)
            for (day, info) in month
                if haskey(info, "trades")
                    for (id_str, trade) in info["trades"]
                        profit += haskey(trade, "profit_EUR") ? trade["profit_EUR"] : trade["profit"]
                    end
                end
            end
        end
        #
        ## Grid bots
        grid_bot_paths = glob(joinpath(Pinguu.base_dir, "grid_bots/user/" * u * "/*/config.json"))
        for grid_bot_path in grid_bot_paths
            tmp_config = JSON.parsefile(grid_bot_path)
            if tmp_config["testing"] == "0"
                tmp = glob(joinpath(dirname(grid_bot_path), "history/$yr/$m.json"))
                for file in tmp
                    month = JSON.parsefile(file)
                    for (day, info) in month
                        for trade in info
                            profit += trade["profit"]
                        end
                    end
                end
            end
        end

        push!(past_months, yr * "-" * m)
        push!(profits, profit)
    end
    return past_months, profits
end
export get_profit_for_past_months


text:
# Autogenerated wrapper script for ADIOS2_jll for x86_64-w64-mingw32-libgfortran4-cxx03
export bpls, libadios2_c, libadios2_c_mpi, libadios2_core, libadios2_core_mpi, libadios2_cxx11, libadios2_cxx11_mpi

using Blosc_jll
using Bzip2_jll
using CompilerSupportLibraries_jll
using MicrosoftMPI_jll
using ZeroMQ_jll
using libpng_jll
using zfp_jll
JLLWrappers.@generate_wrapper_header("ADIOS2")
JLLWrappers.@declare_library_product(libadios2_c, "libadios2_c.dll")
JLLWrappers.@declare_library_product(libadios2_c_mpi, "libadios2_c_mpi.dll")
JLLWrappers.@declare_library_product(libadios2_core, "libadios2_core.dll")
JLLWrappers.@declare_library_product(libadios2_core_mpi, "libadios2_core_mpi.dll")
JLLWrappers.@declare_library_product(libadios2_cxx11, "libadios2_cxx11.dll")
JLLWrappers.@declare_library_product(libadios2_cxx11_mpi, "libadios2_cxx11_mpi.dll")
JLLWrappers.@declare_executable_product(bpls)
function __init__()
    JLLWrappers.@generate_init_header(Blosc_jll, Bzip2_jll, CompilerSupportLibraries_jll, MicrosoftMPI_jll, ZeroMQ_jll, libpng_jll, zfp_jll)
    JLLWrappers.@init_library_product(
        libadios2_c,
        "bin\\libadios2_c.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libadios2_c_mpi,
        "bin\\libadios2_c_mpi.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libadios2_core,
        "bin\\libadios2_core.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libadios2_core_mpi,
        "bin\\libadios2_core_mpi.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libadios2_cxx11,
        "bin\\libadios2_cxx11.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libadios2_cxx11_mpi,
        "bin\\libadios2_cxx11_mpi.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        bpls,
        "bin\\bpls.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Documenter, Revise

makedocs(
    modules = [Revise],
    clean = false,
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "Revise.jl",
    authors = "Tim Holy",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "index.md",
        "config.md",
        "limitations.md",
        "debugging.md",
        "internals.md",
        "user_reference.md",
        "dev_reference.md",
    ],
)

deploydocs(
    repo = "github.com/timholy/Revise.jl.git",
)


text:
import YAML
config = YAML.load(open("config.yaml"))

using DiskJockey.constants
using DiskJockey.model
using DiskJockey.visibilities
using DiskJockey.image
using DiskJockey.gridding
using HDF5

import PyPlot.plt
using LaTeXStrings

dv = DataVis(config["data_file"], true)[5]

# Conjugation is necessary for the SMA and ALMA
visibilities.conj!(dv)

qq = get_qq(dv)

# Find the indices that sort these in increasing order
ind = sortperm(qq)

qq = qq[ind]
VV = dv.VV[ind]

fig, ax = plt[:subplots](nrows=3, figsize=(10, 9), sharex=true)

ax[1][:plot](qq, abs(VV), "k.")
ax[2][:plot](qq, real(VV), "k.")
ax[3][:plot](qq, imag(VV), "k.")

ax[3][:set_xlabel](L"$qq\quad k\lambda$")

ax[1][:set_ylabel]("Mod")
ax[2][:set_ylabel]("Real")
ax[3][:set_ylabel]("Imag")

plt[:savefig]("sorted_visibilities.png", dpi=600)

max_stars_repo_name: iancze/DiskEmulator

max_issues_repo_name: iancze/DiskEmulator

max_forks_repo_name: iancze/DiskEmulator


text:
# Coulomb's Law
# F = k * ((q1 * q2) / r²)
function coulombs(particle1, particle2)
  k = 0.0003
  force = k * ((charge(particle1) * charge(particle2)) / max(magnitude(particle1.position - particle2.position), 0.0001))
  q = (charge(particle1) * charge(particle2))
  ma = max(magnitude(particle1.position - particle2.position), 0.0001)
  positive = Vector3(
    force * (particle1.position.x - particle2.position.x),
    force * (particle1.position.y - particle2.position.y),
    force * (particle1.position.z - particle2.position.z)
  )
  negative = Vector3(
    -force * (particle1.position.x - particle2.position.x),
    -force * (particle1.position.y - particle2.position.y),
    -force * (particle1.position.z - particle2.position.z)
  )
  return (positive, negative)
end

## Velocity
# Δx = v.x * Δt
# Δy = v.y * Δt
# Δz = v.z * Δt
function applyVelocity(particle, deltaTime)
  particle.position = particle.position + Vector3(
    particle.velocity.x * deltaTime,
    particle.velocity.y * deltaTime,
    particle.velocity.z * deltaTime
  )
end

## Acceleration
# Δv.x = a.x * Δt
# Δv.y = a.y * Δt
# Δv.z = a.z * Δt
function applyAcceleration(particle, acceleration, deltaTime)
  particle.velocity = particle.velocity + Vector3(
    acceleration.x * deltaTime,
    acceleration.y * deltaTime,
    acceleration.z * deltaTime
  )
end

## F = ma
## a = F/m
function applyForce(particle, force, deltaTime)
  mass = particleMass(particle)
  acceleration = Vector3(
    force.x / mass,
    force.y / mass,
    force.z / mass
  )
  return acceleration
end

max_stars_repo_name: jonathanbuchanan/PhysicsEngine

max_issues_repo_name: jonathanbuchanan/PhysicsEngine

max_forks_repo_name: jonathanbuchanan/PhysicsEngine


text:
module GetindexArrays

export GetindexArray

struct GetindexArray{T,N,F,Axs<:NTuple{N,AbstractUnitRange{Int}}} <: AbstractArray{T,N}
    getindex::F
    axes::Axs
end
GetindexArray{T}(getindex::F, axes) where {T,F} =
    GetindexArray{T, length(axes), F, typeof(axes)}(getindex, axes)
GetindexArray(getindex::F, axes) where F =
    GetindexArray{Base._return_type(getindex, NTuple{length(axes), Int})}(getindex, axes)

Base.@propagate_inbounds function Base.getindex(A::GetindexArray{T,N}, i::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A, i...)
    convert(T, A.getindex(i...))::T
end
Base.size(A::GetindexArray) = map(length, A.axes)
Base.axes(A::GetindexArray) = A.axes

end


text:
print("\tdenseBlockAngular.jl")
Base.BLAS.set_num_threads(1)

srand(0)
m = 24
n = 32
R = 16
u = [rand(m, n) for _ in 1:R]

# Constructor test
A = Tulip.Cholesky.DenseBlockAngular(u)
@test m == A.m
@test R*n == A.n
@test R == A.R
for r in 1:R
    @test u[r] == A.cols[:, A.colptr[r]:(A.colptr[r+1]-1)]
end

# Base interface tests
@test size(A) == (m+R, n*R)
@test A[1, 1] == 1.0
@test A[1, end] == (R == 1 ? 1.0 : 0.0)
@test A[end, end] == u[end][end, end]
@test A[R+1, 1] == u[1][1, 1]

# Matrix-Vector multiplication tests
x = rand(A.n)
# these two should not throw any error
y = A * x
Base.LinAlg.A_mul_B!(y, A, x)
y_ = hcat(u...) * x
@test y[(R+1):end] == y_

# factorization tests
F = Tulip.Cholesky.cholesky(A, ones(A.n))
@test m == F.m
@test R == F.R
@test n*R == F.n
@test F.colptr[end] == (F.n+1)
# factor update
θ = rand(A.n)
Tulip.Cholesky.cholesky!(A, θ, F)

# Left division tests
A_ = sparse(A)  # sparse representation of A
b = rand(m+R)

y = F \ b
err = maximum(abs.(A_ * (θ .* (A_' * y)) - b))
@test err < 10.0^-10

Base.LinAlg.A_ldiv_B!(y, F, b)
err = maximum(abs.(A_ * (θ .* (A_' * y)) - b))
@test err < 10.0^-10

y = copy(b)
Base.LinAlg.A_ldiv_B!(F, y)
err = maximum(abs.(A_ * (θ .* (A_' * y)) - b))
@test err < 10.0^-10

# Tulip tests
# create and solve model
m, n, R = 2, 2, 4
u = 1.0 - 2.0 * rand(m, n*R)
for r in 1:R
    u[:, 1+n*(r-1)] = 0.0
end
colptr = [1+n*(r-1) for r in 1:(R+1)]
A = Tulip.Cholesky.DenseBlockAngular(m, n*R, R, colptr, u)
b = vcat(ones(R), zeros(m))
c = rand(n*R)
colub_ind = collect(1:(n*R))
colub_val = 10.0 * ones(n*R)

# solve model
model = Tulip.Model(A, b, c, colub_ind, colub_val)
Tulip.solve!(model, verbose=0, tol=10.0^-8)

println("\tPassed.")


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

function bsplinesurface()
    vknots = KnotVector{Float64}([0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 1.0, 1.0, 1.0])
    uknots = KnotVector{Float64}([0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0])
    points = map(
        x -> collect(x),
        [
            (0.0, 0.0, 0.0) (0.33, 0.0, 0.0) (0.66, 0.0, 0.0) (1.0, 0.0, 0.0) (1.33, 0.0, 0.0)
            (0.0, 0.33, 0.0) (0.33, 0.33, 1.0) (0.66, 0.33, 1.0) (1.0, 0.33, 0.5) (1.33, 0.33, 0.0)
            (0.0, 0.66, 0.0) (0.33, 0.66, 1.0) (0.66, 0.66, 1.0) (1.0, 0.66, 0.5) (1.33, 0.66, 0.0)
            (0.0, 1.0, 0.0) (0.33, 1.0, 0.0) (0.66, 1.0, 0.0) (1.0, 1.0, 0.0) (1.33, 1.0, 0.0)
        ],
    )
    return BSplineSurface{OpticSim.Euclidean,Float64,3,3}(uknots, vknots, points)
end

function beziersurface()
    # only simple way I could think of to write this array literal as 2D array of 1D arrays.
    # Julia doesn't parse a 2D array of 1D arrays properly. Maybe should use tuples instead though.
    points = map(
        x -> collect(x),
        [
            (0.0, 0.0, 0.0) (0.0, 0.33, 0.0) (0.0, 0.66, 0.0) (0.0, 1.0, 0.0)
            (0.33, 0.0, 0.0) (0.33, 0.33, 1.0) (0.33, 0.66, 1.0) (0.33, 1.0, 0.0)
            (0.66, 0.0, 0.0) (0.66, 0.33, 1.0) (0.66, 0.66, 1.0) (0.66, 1.0, 0.0)
            (1.0, 0.0, 0.0) (1.0, 0.33, 0.0) (1.0, 0.66, 0.0) (1.0, 1.0, 0.0)
        ],
    )
    return BezierSurface{OpticSim.Euclidean,Float64,3,3}(points)
end

function upsidedownbeziersurface()
    points = map(
        x -> collect(x),
        [
            (0.0, 0.0, 0.0) (0.33, 0.0, 0.0) (0.66, 0.0, 0.0) (1.0, 0.0, 0.0)
            (0.0, 0.33, 0.0) (0.33, 0.33, -1.0) (0.66, 0.33, -1.0) (1.0, 0.33, -.5)
            (0.0, 0.66, 0.0) (0.33, 0.66, -1.0) (0.66, 0.66, -1.0) (1.0, 0.66, -.5)
            (0.0, 1.0, 0.0) (0.33, 1.0, 0.0) (0.66, 1.0, 0.0) (1.0, 1.0, 0.0)
        ],
    )
    return BezierSurface{OpticSim.Euclidean,Float64,3,3}(points)
end

function wavybeziersurface()
    points = map(
        x -> collect(x),
        [
            (0.0, 0.0, -0.3) (0.0, 0.33, 1.0) (0.0, 0.66, -1.0) (0.0, 1.0, 0.3)
            (0.33, 0.0, -0.3) (0.33, 0.33, 1.0) (0.33, 0.66, -1.0) (0.33, 1.0, 0.3)
            (0.66, 0.0, -0.3) (0.66, 0.33, 1.0) (0.66, 0.66, -1.0) (0.66, 1.0, 0.3)
            (1.0, 0.0, -0.3) (1.0, 0.33, 1.0) (1.0, 0.66, -1.0) (1.0, 1.0, 0.3)
        ],
    )
    return BezierSurface{OpticSim.Euclidean,Float64,3,3}(points)
end

function verywavybeziersurface()
    points = map(
        x -> collect(x),
        [
            (0.0, 0.0, 0.5) (0.0, 0.2, 1.0) (0.0, 0.4, -3.0) (0.0, 0.6, 3.0) (0.0, 0.8, -1.0) (0.0, 1.0, 0.5)
            (0.2, 0.0, 0.5) (0.2, 0.2, 1.0) (0.2, 0.4, -3.0) (0.2, 0.6, 3.0) (0.2, 0.8, -1.0) (0.2, 1.0, 0.5)
            (0.4, 0.0, 0.0) (0.4, 0.2, 1.0) (0.4, 0.4, -3.0) (0.4, 0.6, 3.0) (0.4, 0.8, -1.0) (0.4, 1.0, 0.5)
            (0.6, 0.0, 0.0) (0.6, 0.2, 1.0) (0.6, 0.4, -3.0) (0.6, 0.6, 3.0) (0.6, 0.8, -1.0) (0.6, 1.0, 0.5)
            (0.8, 0.0, -0.5) (0.8, 0.2, 1.0) (0.8, 0.4, -3.0) (0.8, 0.6, 3.0) (0.8, 0.8, -1.0) (0.8, 1.0, 0.5)
            (1.0, 0.0, -0.5) (1.0, 0.2, 1.0) (1.0, 0.4, -3.0) (1.0, 0.6, 3.0) (1.0, 0.8, -1.0) (1.0, 1.0, 0.5)
        ],
    )
    return BezierSurface{OpticSim.Euclidean,Float64,3,5}(points)
end

qtypesurface1() = QTypeSurface(1.5, radius = 5.0, conic = 1.0, αcoeffs = [(0, 1, -0.3), (5, 4, 0.1)], βcoeffs = [(3, 7, 0.1)], normradius = 1.6)

zernikesurface1() = ZernikeSurface(1.5, radius = 5.0, conic = 1.0, zcoeff = [(1, 0.0), (4, -0.1), (7, 0.05), (33, 0.03)], aspherics = [(2, 0.0), (10, 0.01)])

zernikesurface1a() = ZernikeSurface(1.5, radius = 5.0, conic = 1.0, zcoeff = [(1, 0.0), (4, -0.1), (7, 0.05), (33, 0.03)], aspherics = [(2, 0.0), (10, 0.01)], normradius = 1.6)

zernikesurface2() = ZernikeSurface(1.5, radius = 5.0, conic = 1.0, zcoeff = [(7, 0.2), (9, 0.2)], aspherics = [(10, -0.01)])

zernikesurface3() = ZernikeSurface(9.5, zcoeff = [(1, 0.0), (4, -1.0), (7, 0.5)])

conicsurface() = AcceleratedParametricSurface(ZernikeSurface(9.5, radius = -15.0, conic = -5.0))

function simplesaggrid()
    points = map(
        x -> collect(x),
        [
            (0.1, 0.0, 0.0, 0.0) (0.3, 0.0, 0.0, 0.0) (0.0, 0.0, 0.0, 0.0) (0.3, 0.0, 0.0, 0.0)
            (0.4, 0.0, 0.0, 0.0) (0.1, 0.0, 0.0, 0.0) (0.0, 0.0, 0.0, 0.0) (0.1, 0.0, 0.0, 0.0)
            (-0.3, 0.0, 0.0, 0.0) (0.1, 0.0, 0.0, 0.0) (0.2, 0.0, 0.0, 0.0) (0.3, 0.0, 0.0, 0.0)
            (0.2, 0.0, 0.0, 0.0) (0.3, 0.0, 0.0, 0.0) (-0.1, 0.0, 0.0, 0.0) (-0.1, 0.0, 0.0, 0.0)
            (0.1, 0.0, 0.0, 0.0) (0.3, 0.0, 0.0, 0.0) (0.0, 0.0, 0.0, 0.0) (0.3, 0.0, 0.0, 0.0)
            (0.4, 0.0, 0.0, 0.0) (0.1, 0.0, 0.0, 0.0) (0.0, 0.0, 0.0, 0.0) (0.1, 0.0, 0.0, 0.0)
        ],
    )
    return points
end

gridsagsurfacelinear() = GridSagSurface(ZernikeSurface(1.5, radius = 5.0, conic = 1.0, zcoeff = [(1, 0.0), (4, -0.1), (7, 0.05), (33, 0.03)], aspherics = [(2, 0.0), (10, 0.01)]), simplesaggrid(), interpolation = GridSagLinear)

gridsagsurfacebicubic() = GridSagSurface(ZernikeSurface(1.5, radius = 25.0, conic = 1.0, zcoeff = [(1, 0.0), (4, -0.1), (7, 0.05), (33, 0.03)], aspherics = [(2, 0.0), (10, 0.01)]), simplesaggrid(), interpolation = GridSagBicubic)

gridsagsurfacebicubiccheby() = GridSagSurface(ChebyshevSurface(1.5, 2.0, radius = 25.0, conic = 0.2, [(0, 1, 0.03), (2, 1, -0.01)]), simplesaggrid(), interpolation = GridSagBicubic)

chebyshevsurface1() = ChebyshevSurface(2.0, 2.0, [(1, 1, 0.1), (4, 5, -0.3), (2, 1, 0.0), (22, 23, 0.01)], radius = 25.0, conic = 0.2)

chebyshevsurface2() = ChebyshevSurface(2.0, 2.0, [(1, 1, 0.1), (0, 2, 0.3), (4, 2, -0.1)], radius = 10.0, conic = 0.2)


text:
using Stingray
using Test

@testset "Stingray.jl" begin
    @test typeof(2) <: Real
end


text:
using MedicalImagingUtils
using Test
using TestSetExtensions


text:
module DEACreject
using Random
const functionNames = ["smallerFit",
                       "smallChange1",
                       "keepAll"]

function smallerFit(rng::MersenneTwister,rInd::BitArray{1},cFitness::Array{Float64,1},fitnessP::Array{Float64,1})
    @inbounds for i = 1:size(rInd,1)
        rInd[i] = cFitness[i] < fitnessP[i];
    end
    nothing
end

function smallChange1(rng::MersenneTwister,rInd::BitArray{1},cFitness::Array{Float64,1},fitnessP::Array{Float64,1})
    @inbounds for i = 1:size(rInd,1)
        rInd[i] = rand(rng) < exp(-abs(1.0-cFitness[i]/fitnessP[i]));
    end
    nothing
end

function keepAll(rng::MersenneTwister,rInd::BitArray{1},cFitness::Array{Float64,1},fitnessP::Array{Float64,1})
    @inbounds for i = 1:size(rInd,1)
        rInd[i] = true;
    end
    nothing
end

end


text:
# (Stoecker, 1971; modified)
# solution is: x1=6.2934300, x[2]=3.8218391, x[3]=201.1593341
# objective = 201.1593341

function nlp2(gm::Bool = false)
    m = JuMP.Model()
    @variable(m, 0 <= x[1:3])
    for i = 1:3
        JuMP.set_upper_bound(x[i], [9.422, 5.9023, 267.417085245][i])
    end
    @objective(m, Min, x[3])
    if !gm
        @NLconstraint(m, e1, 250 + 30*x[1] - 6*x[1]^2 - x[3] == 0)
        @NLconstraint(m, e2, 300 + 20*x[2] - 12*x[2]^2 - x[3] == 0)
        @NLconstraint(m, e3, 150 + 0.5*(x[1]+x[2])^2 - x[3] == 0)
        set_optimizer(m, BARON_SILENT)
        return m
    else
        gm = GlobalModel(model = m, name = "nlp2")
        set_optimizer(gm, CPLEX_SILENT)
        add_nonlinear_constraint(gm, :(x -> 250 + 30*x[1] - 6*x[1]^2 - x[3]), vars = [x[1], x[3]], name = "e1", equality = true)
        add_nonlinear_constraint(gm, :(x -> 300 + 20*x[2] - 12*x[2]^2 - x[3]), vars = [x[2], x[3]], name = "e2", equality = true)
        add_nonlinear_constraint(gm, :(x ->  150 + 0.5*(x[1]+x[2])^2 - x[3]), name = "e3", equality = true)
        return gm
    end
end


text:
function PoleCorrection(arg0::jdouble, arg1::jdouble)
    return PoleCorrection((jdouble, jdouble), arg0, arg1)
end

function get_xp(obj::PoleCorrection)
    return jcall(obj, "getXp", jdouble, ())
end

function get_yp(obj::PoleCorrection)
    return jcall(obj, "getYp", jdouble, ())
end


text:
using Pkg
Pkg.activate(@__DIR__)
cd(@__DIR__)
#using Revise
using ZigZagBoomerang
using DataStructures
using LinearAlgebra
using Random
using SparseArrays
using Test
using FileIO
using Statistics
using Makie, AbstractPlotting

function gridlaplacian(T, m, n)
    S = sparse(T(0.0)I, n*m, n*m)
    linear = LinearIndices((1:m, 1:n))
    for i in 1:m
        for j in 1:n
            for (i2, j2) in ((i + 1, j), (i, j + 1))
                if i2 <= m && j2 <= n
                    S[linear[i, j], linear[i2, j2]] -= 1.
                    S[linear[i2, j2], linear[i, j]] -= 1.

                    S[linear[i, j], linear[i, j]] += 1.
                    S[linear[i2, j2], linear[i2, j2]] += 1.
                end
            end
        end
    end
    S
end

Random.seed!(1)

# Define precision operator of a Gaussian random field (sparse matrix operating on `vec`s of `n*n` matrices)
#n = 100
n = 100
const σ2 = 0.5
Γ0 = 2gridlaplacian(Float64, n, n)
Γ = 0.1I + Γ0
mat(x) = reshape(x, (n, n)) # vector to matrix
function mat0(y)
    mat(y  .- 0.1)
end
# Γ is very sparse
@show nnz(Γ)/length(Γ) # 0.000496

# Corresponding Gaussian potential
# ϕ(x, Γ, y) = 0.5*x'*Γ*x  + dot(x - y, x - y)/(2*σ2) # not used by the program

# Define ∇ϕ(x, i, Γ) giving the partial derivative of ϕ(x) with respect to x[i]
∇ϕ(x, i, Γ, y) = ZigZagBoomerang.idot(Γ, i, x)  + (x[i]-y[i])/σ2 # more efficient that dot(Γ[:, i], x)

# Random initial values
t0 = 0.0
h(x, y) = x^2+(5y/4-sqrt(abs(x)))^2
heart = [ max.(1 - h(x, y), 0) for x in range(-1.5,1.5, length=n), y   in range(-1.1,1.9, length=n)]
image(heart)
μ0 = 5.0*vec(heart)
y = μ = μ0 + randn(n*n)
μpost = yhat = (I + Γ)\y
Γpost = (Γ + I)/σ2

x0 = μpost
θ0 = rand([-1.0,1.0], n*n)

# Rejection bounds
c = [norm(Γpost[:, i], 2) for i in 1:n*n]

# Define ZigZag
Z = ZigZag(Γpost, μpost)
# or try the FactBoomerang
#Z = FactBoomerang(Γ, x0*0, 0.1)

κ = 0.4*ones(length(x0))
# Run sparse ZigZag for T time units and collect trajectory
T = 800.0
@time trace0, (tT, xT, θT), (acc, num) = spdmp(∇ϕ, t0, x0, θ0, T, c, Z, Γ, μ; structured = true, adapt = false)

@time traj0 = collect(discretize(trace0, 0.2))

su = false
adapt = false
trace, (t, x, θ), (acc, num), c = @time sspdmp(∇ϕ, t0, x0, θ0, T, c, Z, κ, Γ, μ; structured = true,
                                                strong_upperbounds = su ,
                                                adapt = adapt)

@time traj = collect(discretize(trace, 0.2))

#display(scene)

#scene1, _ = heatmap(abs.([mat(mean(last.(traj0[end÷2:end]))) mat(mean(last.(traj[end÷2:end]))) ]))
scene2a, _ = heatmap(abs.(mat(μ0 - mean(last.(traj0[end÷2:end])))))
scene2b, _ = heatmap(abs.(mat(μ0 - mean(last.(traj[end÷2:end])))))

#error("stop")
if false

# Prepare surface plot
    M = Node(mat0(traj[end].second))
    scene = Scene()
    surface!(scene, M; shading=false, show_axis=false, colormap = :oleron, colorrange = (-3,3))
    Makie.scale!(scene, 1.0, 1.0, 2.0)
    zlims!(scene, -5.0, 5.0)
    # Movie frames
    mkpath(joinpath(@__DIR__, "output"))
    for i in eachindex(traj)
        M[] = mat0(traj[i].second)
        FileIO.save(joinpath(@__DIR__, "output", "surfs$i.png"), scene)
    end

    # Make video
    dir = joinpath(@__DIR__, "output")
    run(`ffmpeg -y -r 40 -f image2 -i $dir/surfs%d.png -vcodec libx264 -crf 25  -pix_fmt yuv420p $(typeof(Z).name)field.mp4`)
end

# Save hearts
save("hearttrue.png", image(mat(μ0)))
save("heart.png", image(mat(y)))
save("hearthat.png", image(mat(yhat)))
save("heartpostmeana.png", image(mat(mean(last.(traj0[end÷2:end])))))
save("heartpostmeanb.png", image(mat(mean(last.(traj[end÷2:end])))))
save("hearterrora.png", scene2a)
save("hearterrorb.png", scene2b)

mean((mat(μ0 - mean(last.(traj0[end÷2:end])))).^2)
mean((mat(μ0 - mean(last.(traj[end÷2:end])))).^2)

@show mean(abs.(mat(μ0 - mean(last.(traj0[end÷2:end])))))
@show mean(abs.(mat(μ0 - mean(last.(traj[end÷2:end])))))
@show extrema(μ0)

@show  extrema(y)

@show extrema(yhat)

@show extrema((mat(mean(last.(traj0[end÷2:end])))))
@show extrema((mat(mean(last.(traj[end÷2:end])))))

@show extrema(abs.(mat(μ0 - mean(last.(traj0[end÷2:end])))))

@show extrema(abs.(mat(μ0 - mean(last.(traj[end÷2:end])))))


text:
using DynamicalSystemsBase

"""
    example_uncertain_indexvalue_datasets(system::DynamicalSystems.DiscreteDynamicalSystem, n::Int, vars; Ttr = 1000,
        d_xval = Uniform(0.01, 0.4), d_yval = Uniform(0.01, 0.5),
        d_xind = Uniform(0.5, 1.5), d_yind = Uniform(0.5, 1.5))

Generate a pair of `UncertainIndexValueDataset`s from a discrete dynamical `system`, generated
by iterating the system `n` time after a transient run of `Ttr` steps, then gathering the
columns at positions `vars` (should be two column indices) as separate time series.

Each of the time series, call them `x` and `y`, are then converted to uncertain values.
Specifically, replace `x[i]` and `y[i]` with `UncertainValue(Normal, x[i], rand(d_xval)`
and `UncertainValue(Normal, y[i], rand(d_xval)`. Because the time series don't have
explicit time indices associated with them, we'll create some time indices as the range
`1:tstep:length(x)*tstep`, call them `x_inds` and `y_inds`. The time indices for `x` and `y`
are also normally distributed, such that `x_inds[i] = UncertainValue(Normal, i, rand(d_xind)`,
and the same for `y_inds`.

Returns a tuple of `UncertainIndexValueDataset` instances, one for `x` and one for `y`.
"""
function example_uncertain_indexvalue_datasets(system::DiscreteDynamicalSystem, n::Int, vars; Ttr = 1000,
        tstep = 3,
        d_xval = Uniform(0.01, 0.4), d_yval = Uniform(0.01, 0.5),
        d_xind = Uniform(0.5, 1.5), d_yind = Uniform(0.5, 1.5))

    orbit = trajectory(system, n - 1, Ttr = 1000);
    x_vals = UncertainValueDataset([UncertainValue(Normal, x, rand(d_xval)) for x in orbit[:, vars[1]]])
    y_vals = UncertainValueDataset([UncertainValue(Normal, y, rand(d_yval)) for y in orbit[:, vars[2]]])

    x_inds = UncertainIndexDataset([UncertainValue(Normal, i, rand(d_xind)) for i = 1:tstep:length(x_vals)*tstep]);
    y_inds = UncertainIndexDataset([UncertainValue(Normal, i, rand(d_yind)) for i = 1:tstep:length(y_vals)*tstep]);

    X = UncertainIndexValueDataset(x_inds, x_vals)
    Y = UncertainIndexValueDataset(y_inds, y_vals);

    return X, Y
end

export example_uncertain_indexvalue_datasets


text:
gc()

#Tests des packages
include("../ALAS/include.jl")

precision=1e-15

using NLPModels
#package pour plot
using PyPlot
#package pour plot
#using Gadfly

# test : 0 pour tout tester, sinon sélectionner entre 1 et 7.
test=0
done=false

	(test==0 || test==1) ? include("run_test_Thetafunc.jl") : print("")

	(test==0 || test==2) ? include("run_test_Relaxation.jl") : print("")

	# Un exemple simple pour tester :
	 f(x)=x[1]-x[2]
	 G(x)=x[1]
	 H(x)=x[2]
	 c(x)=[1-x[2]]
	 lcon=zeros(1)

	done && (test==0 || test==3) ? include("run_test_ActifMPCCmod.jl") : print("")

	done && (test==0 || test==4) ? include("run_test_UnconstrainedMPCCActif.jl") : print("")

	done && (test==0 || test==5) ? include("run_test_ALASMPCCmod.jl") : print("")

	done && (test==0 || test==6) ? include("run_test_MPCCmod.jl") : print("")

	done && (test==0 || test==7) ? include("run_test_MPCCsolve.jl") : print("")
#FIN des tests


text:
# solve lyapunov equation (dH/dt) = B(t) H(t) + H(t) B(t)' - a(t), backward, conditional on Hend
function lyapunovpsdbackward_step(t, y, dt, P)
    B = Bridge.B(t - dt/2, P)
    ϕ = (I + 1/2*dt*B)\(I - 1/2*dt*B)
    ϕ *(y + 1/2*dt*a(t - dt, P))* ϕ' + 1/2*dt*Bridge.a(t, P)
end

lyapunovpsdbackward(t, P, Hend) = lyapunovpsdbackward!(samplepath(t, zero(Hend)), P, Hend)
function lyapunovpsdbackward!(H, P, Hend) # backward solve on grid tt, with final value Hend
    H.yy[end] = Hend
    t = H.tt
    for i in length(t)-1:-1:1
         H.yy[i] = lyapunovpsdbackward_step(t[i+1], H.yy[i+1], t[i+1] - t[i], P)
    end
    H
end


text:
# Michelle Medina
# DATE

using PyPlot

						# 1D TESTS
# ---------------------------------------------------------
p=Params(101,1.0)
dX(p)

#101 points and identity map gives the coordinates
x = functionOnGrid(x->x,Params(100,2));

sinx=functionOnGrid(x->sin(pi*x),Params(100,2));
plot(x,sinx)
show()

					# 1D Convergence Test
H = []
E = []

for N in 50:200
    grid = Params(N,4)
    u = functionOnGrid(x->sin(x),grid)
	#known derivative of sin(x)
    dUdx = functionOnGrid(x->cos(x),grid)

    #L infinity norm: difference between the numerical derivative and
    #the analytical one, while taking the maximum of it's abs. value
    error = maximum(abs.(Deriv1D(u,grid) - dUdx))

    h = dX(grid)
    #append h to the array H and error to the array E
    append!(H,h)
    append!(E,error)

end

#plotting error vs h
loglog(H,E,"r-")

#first order convergence
loglog(H,H,"g--",label="1")

#second order convergence
loglog(H,H.^2,"b--",label="2")
xlabel("h")
ylabel("L Infinity Norm")
legend(loc="best",frameon=false)
show()

# USING L2 NORM
H = []
E = []

for N in 50:200
    grid = Params(N,4)
    u = functionOnGrid(x->sin(x),grid)
    dUdx = functionOnGrid(x->cos(x),grid) #known derivative of sin(x)
    derivative = Deriv1D(u,grid)

    #L 2 norm: sum of all the differences between the numerical derivative and
    #the analytical one squared and then take the square root
    error = sqrt(sum((derivative[2:end-1] - dUdx[2:end-1]).^2))

    h = dX(grid)
    #append h to the array H and error to the array E
    append!(H,h)
    append!(E,error)

end

#plotting error vs h
loglog(H,E,"r-")

#first order convergence
loglog(H,H,"g--",label="1")

#second order convergence
loglog(H,H.^2,"b--",label="2")
xlabel("h")
ylabel("L 2 Norm")
legend(loc="best",frameon=false)
show()

						# 2D TESTS
# ---------------------------------------------------------

sinxcosy = functionOnGrid((x,y)->sin(x)*cos(y), Params2D(100,200,10.0,10.0));

contourf(sinxcosy)
colorbar()

x2y2= functionOnGrid((x,y)->x^2+y^2,Params2D(10,10,1.0,1.0));

contourf(x2y2)
colorbar()

cosxcosy = DerivX_2D(sinxcosy,Params2D(100,100,1.0,1.0));

					# 2D Convergence Test
# Norm for X ; f(x) = sin(pi*x)*cos(pi*y)
Hx = []
Ex = []
Hy = []
Ey = []

for Nx in 50:100, Ny in 50:100
	grid = Params2D(Nx,Ny,4,4)
	u = functionOnGrid((x,y)->sin(pi*x)*cos(pi*y),grid)

	#analytical derivatives
	dUdx = functionOnGrid((x,y)->pi*cos(pi*x)*cos(pi*y),grid)
	dUdy = functionOnGrid((x,y)->sin(pi*x)*(-pi)*sin(pi*y),grid)
	errorx = maximum(abs.(DerivX_2D(u,grid)-dUdx))
	errory = maximum(abs.(DerivY_2D(u,grid)-dUdy))
	hx=dX_2D(grid)
	hy=dY_2D(grid)
	append!(Hx,hx)
	append!(Ex,errorx)
	append!(Hy,hy)
	append!(Ey,errory)
end

loglog(Hx,Ex,"r-")
loglog(Hx,Hx.^2,"--", label="2")
loglog(Hx,Hx,"--",label="1")
xlabel("hx")
ylabel("L Infinity Norm")
legend(loc="best", frameon=false)
show()

# Norm for Y ; f(x) = sin(pi*x)*cos(pi*y)
Hx = []
Ex = []
Hy = []
Ey = []

for Nx in 50:100, Ny in 50:100
    grid = Params2D(Nx,Ny,4,4)
    u = functionOnGrid((x,y)->sin(pi*x)*cos(pi*y),grid)

	#analytical derivatives
    dUdx = functionOnGrid((x,y)->pi*cos(pi*x)*cos(pi*y),grid)
    dUdy = functionOnGrid((x,y)->sin(pi*x)*(-pi)*sin(pi*y),grid)
    errorx = maximum(abs.(DerivX_2D(u,grid)-dUdx))
    errory = maximum(abs.(DerivY_2D(u,grid)-dUdy))
    hx=dX_2D(grid)
    hy=dY_2D(grid)
    append!(Hx,hx)
    append!(Ex,errorx)
    append!(Hy,hy)
    append!(Ey,errory)
end

loglog(Hy,Ey,"r-")
loglog(Hy,Hy.^2,"--", label="2")
loglog(Hy,Hy,"--",label="1")
xlabel("hy")
ylabel("L Infinity Norm")
legend(loc="best", frameon=false)
show()

# Norm for X ; f(x) = x^2+y^2
Hx = []
Ex = []

for Nx in 50:100, Ny in 50:100
    grid = Params2D(Nx,Ny,4,4)
    u = functionOnGrid((x,y)->x^2+y^2,grid)
    #analytical derivatives
    dUdx = functionOnGrid((x,y)->2x,grid)

    errorx = maximum(abs.(DerivX_2D(u,grid)-dUdx))

    hx=dX_2D(grid)

    append!(Hx,hx)
    append!(Ex,errorx)

end

loglog(Hx,Ex,"r-")
loglog(Hx,Hx.^2,"--", label="2")
loglog(Hx,Hx,"--",label="1")
xlabel("hx")
ylabel("L Infinity Norm")
legend(loc="best", frameon=false)
show()

# Norm for Y ; f(x) = x^2+y^2
Hy = []
Ey = []

for Nx in 50:100, Ny in 50:100
    grid = Params2D(Nx,Ny,4,4)
    u = functionOnGrid((x,y)->x^2+y^2,grid)
    #analytical derivatives
    dUdy = functionOnGrid((x,y)->2y,grid)

    errory = maximum(abs.(DerivY_2D(u,grid)-dUdy))

    hy=dY_2D(grid)

    append!(Hy,hy)
    append!(Ey,errory)

end

loglog(Hy,Ey,"r-")
loglog(Hy,Hy.^2,"--", label="2")
loglog(Hy,Hy,"--",label="1")
xlabel("hy")
ylabel("L Infinity Norm")
legend(loc="best", frameon=false)
show()


text:
"""
Interface between `Iris` and `Arpack`
"""
module ArpackInterface

using Arpack
using ..Common

################################################################################
# LEP

using LinearAlgebra
using SparseArrays
import ..AbstractLinearEigenproblem
import ..HelmholtzLEP
import ..MaxwellLEP

"""
    eigs(lep, ω, [ωs, ψs; nev=6, ncv=max(20,2*nev+1), which=:LM, tol=0.0, maxiter=300, sigma=ω², ritzvec=true, v0=zeros((0,))]) -> (ω, [ψ,], nconv, niter, nmult, resid)

`lep` is an `AbstractLinearEigenproblem`, constructed from `HelmholtzLEP` or `MaxwellLEP`.
If frequencies `ωs::Vector` and fields `ψs::VectorField` are provided, susceptability is *saturated*.
Linear eigenproblem evaluates any boundary matching and dispersive domains at `ω`.
"""
function Arpack.eigs(
			lep::AbstractLinearEigenproblem,
			ω::Number,
			args...;
			kwargs...)

	# get A, B matrices from lep
    A, B = lep(ω, args...)
	# solve B\A for eigenpair
    ω²::Vector{ComplexF64}, ψ::Matrix{ComplexF64}, nconv, niter, nmult, resid = eigs(spdiagm(0=>1 ./diag(B))*A; kwargs..., sigma=ω^2)
	# wrap ψ as a VectorField, depending on Helmhotz vs Maxwell
	if typeof(lep) <: HelmholtzLEP
		Ψ = ScalarField(lep,ψ)
	elseif typeof(lep) <: MaxwellLEP
		Ψ = ElectricField(lep,ψ)
	end
	# return ω (eval is ω²)
    return sqrt.(ω²), Ψ, nconv, niter, nmult, resid
end

################################################################################
# CF

using LinearAlgebra
using SparseArrays
import ..AbstractCFEigenproblem
import ..HelmholtzCF
import ..MaxwellCF

"""
    eigs(cf, ω, [ωs, ψs; η=0, nev=6, ncv=max(20,2*nev+1), which=:LM, tol=0.0, maxiter=300, sigma=η, ritzvec=true, v0=zeros((0,))]) -> (η, [u,], nconv, niter, nmult, resid)

`cf` is an `AbstractCFEigenproblem`, constructed from `HelmholtzCF` or `MaxwellCF`.
If frequencies `ωs::Vector` and fields `ψs::ElectricField` are provided, susceptability is *saturated*.
CF eigenproblem evaluates any boundary matching and dispersive domains at `ω`.
"""
function Arpack.eigs(
			cf::AbstractCFEigenproblem,
			ω::Number,
			args...;
			η::Number=0,
			kwargs...)

	# get A, B from cf
    A, B = cf(ω, args...)
	# note Arpack can't handle B which is not positive semidefinite
		lmul!(-1,A)
		lmul!(-1,B)
		n = length(cf.simulation)
		sgn = sparse(I*1,n,n)
		rows = rowvals(B)
		vals = nonzeros(B)
		@inbounds for col ∈ 1:n
			@inbounds for j ∈ nzrange(B, col)
				row = rows[j]
				sgn[row,col] = flipsign(sgn[row,col],B.nzval[j])
				B.nzval[j] = abs(B.nzval[j])
			end
		end
	# now do actual solve
    ηs::Vector{ComplexF64}, u::Matrix{ComplexF64}, nconv, niter, nmult, resid = eigs(sgn*A, B; kwargs..., sigma=η)
	# wrap u as a VectorField, depending on Helmhotz vs Maxwell
	if typeof(cf) <: HelmholtzCF
		U = ScalarField(cf, u)
	elseif typeof(cf) <: MaxwellCF
		U = ElectricField(cf, u)
	end
    return ηs, U, nconv, niter, nmult, resid
end

################################################################################
# maxwelleigen and helmholtzeigen, LEP

using LinearAlgebra
import ..orthogonalize!
import ..INDEX_OFFSET

# applies when Arpack is default interface
function iris_eigen_arpack(
            lep::AbstractLinearEigenproblem,
            ω::Number,
			args...;
            verbose::Bool=false,
            kwargs...)

	# use previously defined eigs
	ωs, ψs, nconv, niter, nmult, resid = eigs(lep, ω, args...; kwargs...)
    length(ωs) ≤ nconv || @warn "$(length(ωs) - nconv) evecs did not converge"
	# Normalize according to (ψ₁,ψ₂)=δ₁₂
    normalize!(lep.simulation, ψs.values, lep.αεpFχ, size(ψs,1)÷2+INDEX_OFFSET)
	# Orthogonalize degenerate evecs
	if typeof(lep) <: HelmholtzLEP
		orthogonalize!(ψs,lep.simulation, ωs, lep.αεpFχ, 0, 0)
	elseif typeof(lep) <: MaxwellLEP
		orthogonalize!(ψs,lep.simulation, ωs, lep.αεpFχ, lep.ky, lep.kz)
		# for Maxwell ky=kz=0, there is inherent two-fold polarization degeneracy, this is cheap way to treat
		if all(iszero,(lep.ky,lep.kz)) normalize!(lep.simulation, ψs.values, lep.αεpFχ, size(ψs,1)÷2+INDEX_OFFSET) end
	end
	return ωs, ψs
end

import ..DEFAULT_LINEAR_EIGENSOLVER
import ..DOC_LEP_H
import ..DOC_LEP_M
import ..helmholtzeigen
import ..maxwelleigen

if DEFAULT_LINEAR_EIGENSOLVER == :Arpack
	@doc """
	$DOC_LEP_H
	`nev` Number of eigenvalues (`6`) |
	`v0` Starting vector (`zeros((0,))`) |
	`maxiter` Maximum iterations (`300`) |
	`ncv` Number of Krylov vectors (`max(20,2*nev+1)`) |
	`tol` Tolerance is max of ε and `tol` (`0.0`)
	"""
	helmholtzeigen(lep::HelmholtzLEP, args...;kwargs...) = iris_eigen_arpack(lep,args...;kwargs...)

	@doc """
	$DOC_LEP_M
	`nev` Number of eigenvalues (`6`) |
	`v0` Starting vector (`zeros((0,))`) |
	`maxiter` Maximum iterations (`300`) |
	`ncv` Number of Krylov vectors (`max(20,2*nev+1)`) |
	`tol` Tolerance is max of ε and `tol` (`0.0`)
	"""
	maxwelleigen(lep::MaxwellLEP, args...;kwargs...) = iris_eigen_arpack(lep,args...;kwargs...)
end

################################################################################
# maxwelleigen and helmholtzeigen, CF

function iris_eigen_arpack(
            cf::AbstractCFEigenproblem,
            ω::Number,
			args...;
			η::Number=0,
            verbose::Bool=false,
            kwargs...)

	# use previously defined eigs
	ηs, us, nconv, niter, nmult, resid = eigs(cf, ω, args...; η=η, kwargs...)
    length(ηs) ≤ nconv || @warn "$(length(ηs) - nconv) evecs did not converge"
	# Normalize according to (ψ₁,ψ₂)=δ₁₂ via CF.F, which may be saturated
    normalize!(cf.simulation, us, cf.F)
	# Orthogonalize degenerate evecs
	orthogonalize!(us,cf.simulation, ηs, cf.F, cf.ky, cf.kz)
	# for Maxwell ky=kz=0, there is inherent two-fold polarization degeneracy, this is cheap way to treat
	if all(iszero,(cf.ky,cf.kz)) normalize!(cf.simulation, us, cf.F) end
	return ηs, us
end

import ..DOC_CF_H
import ..DOC_CF_M

if DEFAULT_LINEAR_EIGENSOLVER == :Arpack
	@doc """
	$DOC_CF_H
	`nev` Number of eigenvalues (`6`) |
	`v0` Starting vector (`zeros((0,))`) |
	`maxiter` Maximum iterations (`300`) |
	`ncv` Number of Krylov vectors (`max(20,2*nev+1)`) |
	`tol` Tolerance is max of ε and `tol` (`0.0`)
	""" ->
	helmholtzeigen(cf::HelmholtzCF, args...;kwargs...) = iris_eigen_arpack(cf,args...;kwargs...)

	@doc """
	$DOC_CF_M
	`nev` Number of eigenvalues (`6`) |
	`v0` Starting vector (`zeros((0,))`) |
	`maxiter` Maximum iterations (`300`) |
	`ncv` Number of Krylov vectors (`max(20,2*nev+1)`) |
	`tol` Tolerance is max of ε and `tol` (`0.0`)
	""" ->
	maxwelleigen(cf::MaxwellCF, args...;kwargs...) = iris_eigen_arpack(cf,args...;kwargs...)
end

end #module

using .ArpackInterface


text:
using POMDPs

POMDPs.add("BasicPOMCP")
POMDPs.add("QMDP")
POMDPs.add("POMCPOW")
POMDPs.add("DESPOT")
POMDPs.add("DiscreteValueIteration")
POMDPs.add("ARDESPOT")

try Pkg.clone("https://github.com/slundberg/PmapProgressMeter.jl.git") end

try Pkg.clone("https://github.com/zsunberg/LightDarkPOMDPs.jl") end
try Pkg.clone("https://github.com/zsunberg/VDPTag.jl.git") end
try Pkg.clone("https://github.com/zsunberg/VDPTag2.jl.git") end
try Pkg.clone("https://github.com/zsunberg/LaserTag.jl.git") end
try Pkg.clone("https://github.com/zsunberg/SubHunt.jl.git") end


text:
@testset "preprocess for training sqth_th" begin
    n = 6
    file = joinpath(SqState.training_data_path(), "test.jld2")

    points, 𝛒s, args, σs = gen_data_sqth_th(n_data=n)
    jldsave(file; points, 𝛒s, args, σs)

    loader = SqState.preprocess_q2σs(file, batch_size=2)
    x, y = first(loader)
    @test size(x) == (2, 4096, 2)
    @test size(y) == (4096, 2)

    loader = SqState.preprocess_q2ρ(file, batch_size=2)
    x, y = first(loader)
    @test size(x) == (2, 4096, 2)
    @test size(y) == (2, 100, 100, 2)

    loader = SqState.preprocess_q2args(file, batch_size=2)
    x, y = first(loader)
    @test size(x) == (2, 4096, 2)
    @test size(y) == (6, 2)

    rm(file)
end

@testset "preprocess for training sqth" begin
    n = 6
    file = joinpath(SqState.training_data_path(), "test.jld2")

    points, 𝛒s, args, σs = gen_data_sqth(n_data=n)
    jldsave(file; points, 𝛒s, args, σs)

    loader = SqState.preprocess_q2σs(file, batch_size=2)
    x, y = first(loader)
    @test size(x) == (2, 4096, 2)
    @test size(y) == (4096, 2)

    loader = SqState.preprocess_q2ρ(file, batch_size=2)
    x, y = first(loader)
    @test size(x) == (2, 4096, 2)
    @test size(y) == (2, 100, 100, 2)

    loader = SqState.preprocess_q2args(file, batch_size=2)
    x, y = first(loader)
    @test size(x) == (2, 4096, 2)
    @test size(y) == (3, 2)

    rm(file)
end


text:
using Base.Test
using MathProgBase

function mixintprogtest(solver=MathProgBase.defaultMIPsolver)
    println("Testing mixintprog with solver ", string(typeof(solver)))

    # integer knapsack problem
    sol = mixintprog(-[5,3,2,7,4],Float64[2 8 4 2 5],'<',10,:Int,0,1,solver)
    @test sol.status == :Optimal
    @test_approx_eq_eps sol.objval -16.0 1e-6
    @test_approx_eq_eps norm(sol.sol - [1.0, 0.0, 0.0, 1.0, 1.0]) 0.0 1e-4

    println("Done")
end


text:
using Test
using SequentialRaytrace

@testset "Raytracing" begin
    X = 0.0
    Y = 0.0

    ray = Ray(0.0, 20.0, 0.0, X, Y)
    @testset "paraxial surface x" begin
        lens = Lens(
            "",
            Object(Air, 200.0),
            [OpticalComponent(:base, Nothing, :meta, [Paraxial(150.5, nothing, Air, 0.0)])],
        )
        result = gen_result(Vector{typeof(ray)}, lens)
        @test trace!(result, lens, ray, 1.0)[end] ≈ Ray(0.0, 20.0, 0.0, 0.0, -0.1317322700) atol =
            0.00000000001

        ray2 = ray_from_NA(0.0, 0.1)
        r = trace!(result, lens, ray2, 1.0)
        @test r[end] ≈ Ray(0.0, 20.100756305, 0.0, 0.0, -0.0330380156) atol = 0.000000001
        y = r[end].y
        cy = r[end].cy
        cz = r[end].cz
        sl = cy / cz
        @test -y / sl ≈ 608.081 atol = 1e-3
    end

    ray = Ray(5.0, 20.0, 0.0, X, Y)
    @testset "without raytracing errors" begin
        lens = Lens(
            "",
            Object(Air, 200.0),
            [
                OpticalComponent(
                    :c1,
                    Nothing,
                    :meta1,
                    [
                        Sphere(1 / 50, nothing, Silica, 15.0)
                        EvenAsphere(-1 / 50, 0.0, 1e-7, 1e-9, 0.0, nothing, Air, 35.0)
                    ],
                ),
                OpticalComponent(
                    :c2,
                    Nothing,
                    "meta2",
                    [
                        Plano(nothing, Air, 10.0)
                        Paraxial(-120.3, nothing, Air, 30.0)
                    ],
                ),
            ],
        )
        result = gen_result(Vector{typeof(ray)}, lens)
        @test trace!(result, lens, ray, 1.0)[end] ≈
              Ray(-4.0353445, -16.1413780, 0.0, -0.103684, -0.4147364) atol = 0.000001
    end

    X = 0.1
    Y = -0.1
    ray = Ray(0.0, 20.0, 0.0, X, Y)
    @testset "with raytracing errors" begin
        lens = Lens(
            "",
            Object(Air, 200.0),
            [OpticalComponent(
                :base,
                Nothing,
                nothing,
                [
                    Sphere(1 / 50, ClearDiameter(2 * sqrt(40.0^2 + 40.0^2)), Silica, 15.0)
                    EvenAsphere(
                        -1 / 50,
                        0.0,
                        1e-7,
                        1e-9,
                        0.0,
                        ClearDiameter(2 * sqrt(20.0^2 + (-1.0)^2)),
                        Air,
                        35.0,
                    )
                    Plano(nothing, Air, 15.0)
                    Plano(ClearDiameter(2 * sqrt(10.0^2 + 10.0^2)), Air, 15.0)
                ],
            )],
        )
        result = gen_result(Vector{typeof(ray)}, lens)
        @test iserror(trace!(result, lens, ray, 1.0))
        @test trace!(result, lens, ray, 1.0).data[3] ≈
              Ray(20.189725, -1.065485, -4.184919, -0.313750, -0.082502) atol = 0.000001
        @test trace!(result, lens, ray, 1.0).global_index == 3
        @test trace!(result, lens, ray, 1.0).component_id == :base
        @test trace!(result, lens, ray, 1.0).local_index == 2
        @test iserror(trace!(result, lens, Ray(0.0, 0.0, 0.0, X, Y), 1.0))
        ray = Ray(0.0, 20.0, 0.0, X, Y)
        @test error_type(trace!(result, lens, Ray(0.0, 0.0, 0.0, X, Y), 1.0)) === :vignetted
    end

end


text:
function get_date(obj::FieldTurnSpan)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_turn_duration(obj::FieldTurnSpan)
    return jcall(obj, "getTurnDuration", RealFieldElement, ())
end

function get_turn_end_date(obj::FieldTurnSpan)
    return jcall(obj, "getTurnEndDate", FieldAbsoluteDate, ())
end

function get_turn_start_date(obj::FieldTurnSpan)
    return jcall(obj, "getTurnStartDate", FieldAbsoluteDate, ())
end

function in_turn_time_range(obj::FieldTurnSpan, arg0::AbsoluteDate)
    return jcall(obj, "inTurnTimeRange", jboolean, (AbsoluteDate,), arg0)
end

function update_end(obj::FieldTurnSpan, arg0::FieldAbsoluteDate, arg1::AbsoluteDate)
    return jcall(obj, "updateEnd", void, (FieldAbsoluteDate, AbsoluteDate), arg0, arg1)
end

function update_start(obj::FieldTurnSpan, arg0::FieldAbsoluteDate, arg1::AbsoluteDate)
    return jcall(obj, "updateStart", void, (FieldAbsoluteDate, AbsoluteDate), arg0, arg1)
end


text:
using MolecularDynamicsIntegrators: move!, RattleIntegrator, FixedDistanceConstraint

@testset "RATTLE" begin
    N = 4
    d = 1.0
    r = [
        Vector3(0, 0, 0),
        Vector3(0, 1, 0),
        Vector3(0, 2, 0),
        Vector3(0, 3, 0)
    ]
    v = [
        Vector3(0, 0, 0),
        Vector3(0.5, 0, 0.8),
        Vector3(0.9, 0, -0.4),
        Vector3(0.6, 0, -0.6)
    ]
    dt = 0.001
    mass = i -> 1.0
    force = lj_potential_forces
    steps = 2000
    tol = 0.000001
    max_iter = 10000
    K = 3
    gs = [
        FixedDistanceConstraint(1, 2, d),
        FixedDistanceConstraint(2, 3, d),
        FixedDistanceConstraint(3, 4, d)
    ]

    setup = RattleIntegrator(r, v, dt, mass, force, tol, max_iter, gs)

    function test_constraints(r, v, gs, tol)
        K = length(gs)
        @inbounds for ki = 1:K
            i, j, d = gs[ki].i, gs[ki].j, gs[ki].distance
            @test abs(norm_sqr(r[i] - r[j]) - d ^ 2) ≤ 2 * tol * d^2
            @test abs((r[i] - r[j]) ⋅ (v[i] - v[j])) ≤ tol * d^2
        end
    end

    test_constraints(setup.positions, setup.velocities, gs, tol)

    for i = 1:steps
        move!(setup)
        test_constraints(setup.positions, setup.velocities, gs, tol)
    end

    @test vector_equal(setup.positions, [
        [0.09737654288485667, 0.7887420889546933, 0.5214751803156499],
        [0.9743787172580586, 1.2551728024457112, 0.4061033839318231],
        [1.5416316823754233, 1.815657545994333, -0.19728780411825142],
        [1.3866130574816684, 2.140427562605264, -1.1302907601292107]
    ])
    @test vector_equal(setup.velocities, [
        [0.22170214165089389, 0.48444541940694, 0.5318306598452084],
        [0.35308534346411097, 0.022864444526390852, -0.3355580167319755],
        [0.5589660298465181, 0.13590564532722463, -0.03700640122784894],
        [0.8662464850384826, -0.6432155092605506, -0.35926624188538064]
    ])
end


text:
export PPOPolicy, PPOTrajectory, MaskedPPOTrajectory

const PPOTrajectory = Trajectory{
    <:NamedTuple{
        (:action_log_prob, SART...),
        <:Tuple{
            <:CircularArrayBuffer,
            <:CircularArrayBuffer,
            <:CircularArrayBuffer,
            <:CircularArrayBuffer,
            <:CircularArrayBuffer,
        },
    },
}

function PPOTrajectory(; capacity, action_log_prob, kwargs...)
    merge(
        CircularArrayTrajectory(;
            capacity = capacity + 1,
            action_log_prob = action_log_prob,
        ),
        CircularArraySARTTrajectory(; capacity = capacity, kwargs...),
    )
end

const MaskedPPOTrajectory = Trajectory{
    <:NamedTuple{
        (:action_log_prob, SLART...),
        <:Tuple{
            <:CircularArrayBuffer,
            <:CircularArrayBuffer,
            <:CircularArrayBuffer,
            <:CircularArrayBuffer,
            <:CircularArrayBuffer,
            <:CircularArrayBuffer,
        },
    },
}

function MaskedPPOTrajectory(; capacity, action_log_prob, kwargs...)
    merge(
        CircularArrayTrajectory(;
            capacity = capacity + 1,
            action_log_prob = action_log_prob,
        ),
        CircularArraySLARTTrajectory(; capacity = capacity, kwargs...),
    )
end

function Base.length(t::Union{PPOTrajectory,MaskedPPOTrajectory})
    x = t[:terminal]
    size(x, ndims(x))
end

"""
    PPOPolicy(;kwargs)

# Keyword arguments

- `approximator`,
- `γ = 0.99f0`,
- `λ = 0.95f0`,
- `clip_range = 0.2f0`,
- `max_grad_norm = 0.5f0`,
- `n_microbatches = 4`,
- `n_epochs = 4`,
- `actor_loss_weight = 1.0f0`,
- `critic_loss_weight = 0.5f0`,
- `entropy_loss_weight = 0.01f0`,
- `dist = Categorical`,
- `rng = Random.GLOBAL_RNG`,

By default, `dist` is set to `Categorical`, which means it will only works
on environments of discrete actions. To work with environments of continuous
actions `dist` should be set to `Normal` and the `actor` in the `approximator`
should be a `GaussianNetwork`. Using it with a `GaussianNetwork` supports
multi-dimensional action spaces, though it only supports it under the assumption
that the dimensions are independent since the `GaussianNetwork` outputs a single
`μ` and `σ` for each dimension which is used to simplify the calculations.
"""
mutable struct PPOPolicy{A<:ActorCritic,D,R} <: AbstractPolicy
    approximator::A
    γ::Float32
    λ::Float32
    clip_range::Float32
    max_grad_norm::Float32
    n_microbatches::Int
    n_epochs::Int
    actor_loss_weight::Float32
    critic_loss_weight::Float32
    entropy_loss_weight::Float32
    rng::R
    n_random_start::Int
    update_freq::Int
    update_step::Int
    # for logging
    norm::Matrix{Float32}
    actor_loss::Matrix{Float32}
    critic_loss::Matrix{Float32}
    entropy_loss::Matrix{Float32}
    loss::Matrix{Float32}
end

function PPOPolicy(;
    approximator,
    update_freq,
    n_random_start = 0,
    update_step = 0,
    γ = 0.99f0,
    λ = 0.95f0,
    clip_range = 0.2f0,
    max_grad_norm = 0.5f0,
    n_microbatches = 4,
    n_epochs = 4,
    actor_loss_weight = 1.0f0,
    critic_loss_weight = 0.5f0,
    entropy_loss_weight = 0.01f0,
    dist = Categorical,
    rng = Random.GLOBAL_RNG,
)
    PPOPolicy{typeof(approximator),dist,typeof(rng)}(
        approximator,
        γ,
        λ,
        clip_range,
        max_grad_norm,
        n_microbatches,
        n_epochs,
        actor_loss_weight,
        critic_loss_weight,
        entropy_loss_weight,
        rng,
        n_random_start,
        update_freq,
        update_step,
        zeros(Float32, n_microbatches, n_epochs),
        zeros(Float32, n_microbatches, n_epochs),
        zeros(Float32, n_microbatches, n_epochs),
        zeros(Float32, n_microbatches, n_epochs),
        zeros(Float32, n_microbatches, n_epochs),
    )
end

function RLBase.prob(
    p::PPOPolicy{<:ActorCritic{<:GaussianNetwork},Normal},
    state::AbstractArray,
    mask,
)
    if p.update_step < p.n_random_start
        @error "todo"
    else
        μ, logσ = p.approximator.actor(send_to_device(device(p.approximator), state)) |> send_to_host
        StructArray{Normal}((μ, exp.(logσ)))
    end
end

function RLBase.prob(p::PPOPolicy{<:ActorCritic,Categorical}, state::AbstractArray, mask)
    logits = p.approximator.actor(send_to_device(device(p.approximator), state))
    if !isnothing(mask)
        logits .+= ifelse.(mask, 0f0, typemin(Float32))
    end
    logits = logits |> softmax |> send_to_host
    if p.update_step < p.n_random_start
        [
            Categorical(fill(1 / length(x), length(x)); check_args = false) for
            x in eachcol(logits)
        ]
    else
        [Categorical(x; check_args = false) for x in eachcol(logits)]
    end
end

function RLBase.prob(p::PPOPolicy, env::MultiThreadEnv)
    mask =  ActionStyle(env) === FULL_ACTION_SET ? legal_action_space_mask(env) : nothing
    prob(p, state(env), mask)
end

function RLBase.prob(p::PPOPolicy, env::AbstractEnv)
    s = state(env)
    s = Flux.unsqueeze(s, ndims(s) + 1)
    mask =  ActionStyle(env) === FULL_ACTION_SET ? legal_action_space_mask(env) : nothing
    prob(p, s, mask)
end

(p::PPOPolicy)(env::MultiThreadEnv) = rand.(p.rng, prob(p, env))

# !!! https://github.com/JuliaReinforcementLearning/ReinforcementLearning.jl/pull/533/files#r728920324
(p::PPOPolicy)(env::AbstractEnv) = rand.(p.rng, prob(p, env))

function (agent::Agent{<:PPOPolicy})(env::MultiThreadEnv)
    dist = prob(agent.policy, env)
    action = rand.(agent.policy.rng, dist)
    if ndims(action) == 2
        action_log_prob = sum(logpdf.(dist, action), dims = 1)
    else
        action_log_prob = logpdf.(dist, action)
    end
    EnrichedAction(action; action_log_prob = vec(action_log_prob))
end

function RLBase.update!(
    p::PPOPolicy,
    t::Union{PPOTrajectory,MaskedPPOTrajectory},
    ::AbstractEnv,
    ::PreActStage,
)
    length(t) == 0 && return  # in the first update, only state & action are inserted into trajectory
    p.update_step += 1
    if p.update_step % p.update_freq == 0
        _update!(p, t)
    end
end

function _update!(p::PPOPolicy, t::AbstractTrajectory)
    rng = p.rng
    AC = p.approximator
    γ = p.γ
    λ = p.λ
    n_epochs = p.n_epochs
    n_microbatches = p.n_microbatches
    clip_range = p.clip_range
    w₁ = p.actor_loss_weight
    w₂ = p.critic_loss_weight
    w₃ = p.entropy_loss_weight
    D = device(AC)

    n_envs, n_rollout = size(t[:terminal])
    @assert n_envs * n_rollout % n_microbatches == 0 "size mismatch"
    microbatch_size = n_envs * n_rollout ÷ n_microbatches

    n = length(t)
    states_plus = send_to_device(D, t[:state])

    states_flatten = flatten_batch(select_last_dim(states_plus, 1:n))
    states_plus_values =
        reshape(send_to_host(AC.critic(flatten_batch(states_plus))), n_envs, :)
    advantages = generalized_advantage_estimation(
        t[:reward],
        states_plus_values,
        γ,
        λ;
        dims = 2,
        terminal = t[:terminal],
    )
    returns = advantages .+ select_last_dim(states_plus_values, 1:n_rollout)

    actions_flatten = flatten_batch(select_last_dim(t[:action], 1:n))
    action_log_probs = select_last_dim(t[:action_log_prob], 1:n)

    # TODO: normalize advantage
    for epoch in 1:n_epochs
        rand_inds = shuffle!(rng, Vector(1:n_envs*n_rollout))
        for i in 1:n_microbatches
            inds = rand_inds[(i-1)*microbatch_size+1:i*microbatch_size]
            if t isa MaskedPPOTrajectory
                lam = send_to_device(
                    D,
                    select_last_dim(
                        flatten_batch(select_last_dim(t[:legal_actions_mask], 2:n+1)),
                        inds,
                    ),
                )
            else
                lam = nothing
            end
            s = send_to_device(D, select_last_dim(states_flatten, inds))  # !!! performance critical
            a = send_to_device(D, select_last_dim(actions_flatten, inds))

            if eltype(a) === Int
                a = CartesianIndex.(a, 1:length(a))
            end

            r = send_to_device(D, vec(returns)[inds])
            log_p = send_to_device(D, vec(action_log_probs)[inds])
            adv = send_to_device(D, vec(advantages)[inds])

            ps = Flux.params(AC)
            gs = gradient(ps) do
                v′ = AC.critic(s) |> vec
                if AC.actor isa GaussianNetwork
                    μ, logσ = AC.actor(s)
                    if ndims(a) == 2
                        log_p′ₐ = vec(sum(normlogpdf(μ, exp.(logσ), a), dims = 1))
                    else
                        log_p′ₐ = normlogpdf(μ, exp.(logσ), a)
                    end
                    entropy_loss = mean(size(logσ, 1) * (log(2.0f0π) + 1) .+ sum(logσ; dims = 1)) / 2
                else
                    # actor is assumed to return discrete logits
                    raw_logit′ = AC.actor(s)
                    if isnothing(lam)
                        logit′ = raw_logit′
                    else
                        logit′ = raw_logit′ .+ ifelse.(lam, 0.0f0, typemin(Float32))
                    end
                    p′ = softmax(logit′)
                    log_p′ = logsoftmax(logit′)
                    log_p′ₐ = log_p′[a]
                    entropy_loss = -sum(p′ .* log_p′) * 1 // size(p′, 2)
                end
                ratio = exp.(log_p′ₐ .- log_p)
                surr1 = ratio .* adv
                surr2 = clamp.(ratio, 1.0f0 - clip_range, 1.0f0 + clip_range) .* adv

                actor_loss = -mean(min.(surr1, surr2))
                critic_loss = mean((r .- v′) .^ 2)
                loss = w₁ * actor_loss + w₂ * critic_loss - w₃ * entropy_loss

                ignore() do
                    p.actor_loss[i, epoch] = actor_loss
                    p.critic_loss[i, epoch] = critic_loss
                    p.entropy_loss[i, epoch] = entropy_loss
                    p.loss[i, epoch] = loss
                end

                loss
            end

            p.norm[i, epoch] = clip_by_global_norm!(gs, ps, p.max_grad_norm)
            update!(AC, gs)
        end
    end
end

function RLBase.update!(
    trajectory::Union{PPOTrajectory,MaskedPPOTrajectory},
    ::PPOPolicy,
    env::MultiThreadEnv,
    ::PreActStage,
    action::EnrichedAction,
)
    push!(
        trajectory;
        state = state(env),
        action = action.action,
        action_log_prob = action.meta.action_log_prob,
    )

    if trajectory isa MaskedPPOTrajectory
        push!(trajectory; legal_actions_mask = legal_action_space_mask(env))
    end
end


text:
using Multilane
using GenerativeModels
using JLD

N = 5000

# NOTE: assumes 4 lanes, desired lane is lane 4
desired_lane_reward = 10.

lambdas = Float64[0.1, 1., 10., 100., 1000.]

nb_lanes = 4 # XXX assumption

rmodels = Multilane.NoCrashRewardModel[
                Multilane.NoCrashRewardModel(desired_lane_reward*lambda,
                                             desired_lane_reward,2.0,nb_lanes)
                for lambda in lambdas]

nb_lanes = 4
pp = PhysicalParam(nb_lanes,lane_length=100.) #2.=>col_length=8\n",
_discount = 0.95
nb_cars=10
dmodel = NoCrashIDMMOBILModel(nb_cars, pp, lane_terminate=false)

curve_problems = []
for i in 1:length(lambdas)
    rmodel = rmodels[i]
    mdp = NoCrashMDP(dmodel, rmodel, _discount);
    push!(curve_problems, mdp)
end

isrng = MersenneTwister(123)
initial_states = MLState[initial_state(curve_problems[1], isrng) for i in 1:N]

initials = assign_keys(curve_problems, initial_states, rng=isrng)

filename = string("initials_", Dates.format(Dates.now(),"u_d_HH_MM"), ".jld")
save(filename, initials)
println("problems and initial conditions saved to $filename")


text:
# https://github.com/altercation/solarized
const _solarized_palette = (:red, :yellow, :blue, :green, :orange, :magenta, :violet, :cyan)

const _solarized_colors = Dict(
    :base03 => RGB255(0, 43, 54),
    :base02 => RGB255(7, 54, 66),
    :base01 => RGB255(88, 110, 117),
    :base00 => RGB255(101, 123, 131),
    :base0  => RGB255(131, 148, 150),
    :base1  => RGB255(147, 161, 161),
    :base2  => RGB255(238, 232, 213),
    :base3  => RGB255(253, 246, 227),
    :blue   => RGB255(38, 139, 210),
    :orange => RGB255(203, 75, 22),
    :red    => RGB255(220, 50, 47),
    :green  => RGB255(133, 153, 0),
    :yellow => RGB255(181, 137, 0),
    :magenta=> RGB255(211, 54, 130),
    :violet => RGB255(108, 113, 196),
    :cyan   => RGB255(42, 161, 152)
)

_themes[:solarized] = PlotTheme(
    bg = _solarized_colors[:base03],
    bginside = _solarized_colors[:base02],
    fg = _solarized_colors[:base00],
    fgtext = _solarized_colors[:base01],
    fgguide = _solarized_colors[:base01],
    fglegend = _solarized_colors[:base01],
    palette = expand_palette(_solarized_colors[:base03], [_solarized_colors[c] for c in _solarized_palette]),
    colorgradient = cgrad(:YlOrRd_r).colors
)

_themes[:solarized_light] = PlotTheme(
    bg = _solarized_colors[:base3],
    bginside = _solarized_colors[:base2],
    fg = _solarized_colors[:base0],
    fgtext = _solarized_colors[:base1],
    fgguide = _solarized_colors[:base1],
    fglegend = _solarized_colors[:base1],
    palette = expand_palette(_solarized_colors[:base3], [_solarized_colors[c] for c in _solarized_palette]),
    colorgradient = cgrad(:YlOrRd_r).colors
)


text:
struct Affine{M,V}
    A::M
    u::V
    """
        Affine(A, b)

    Represents an (inverse) affine transform of the form `x -> A * (x - u)`.
    """
    function Affine(A::AbstractMatrix, u::AbstractVector)
        @assert size(A, 2) == length(u)
        return new{typeof(A), typeof(u)}(A, u)
    end
end

function Affine(A::AbstractMatrix)
    u = zeros(eltype(A), size(A, 2))
    return Affine(A, u)
end

function Affine(u::AbstractVector)
    A = Diagonal(ones(eltype(u), length(u)))
    return Affine(A, u)
end

Affine(::Type{T}, n::Int) where {T} = Affine(zeros(T, n))
Affine(n::Int) = Affine(Float64, n)

"""
    one(t::Affine)

Returns the identity affine transformation, which maps all points to themselves.
"""
Base.one(t::Affine) = Affine(one(t.A), zero(t.u))

"""
    zero(t::Affine)

Returns the zero affine transformation, which maps all points to zero.
"""
Base.zero(t::Affine) = Affine(zero(t.A), zero(t.u))

function Base.copy!(dst::Affine, src::Affine)
    dst.A .= src.A
    dst.u .= src.u
    return dst
end

Base.:*(t::Affine, x::AbstractVecOrMat) = t.A * (x .- t.u)
Base.:\(t::Affine, y::AbstractVecOrMat) = t.A \ y .+ t.u
Base.:*(s::Affine, t::Affine) = Affine(s.A * t.A, t \ s.u)
Base.:\(s::Affine, t::Affine) = Affine(s.A \ t.A, t.u - t.A \ s.A * s.u)

Base.:*(t::Affine, c::Number) = Affine(t.A * c, t.u / c)
Base.:*(c::Number, t::Affine) = Affine(c * t.A, t.u)
Base.:-(t::Affine) = Affine(-t.A, t.u)

Base.:+(s::Affine, t::Affine) = Affine(s.A + t.A, (s.A + t.A) \ (s.A * s.u + t.A * t.u))
Base.:-(s::Affine, t::Affine) = Affine(s.A - t.A, (s.A - t.A) \ (s.A * s.u - t.A * t.u))

Base.inv(t::Affine) = Affine(inv(t.A), -t.A * t.u)

"""
    whitening_transform(μ, C)

Returns the `Affine` transform that whitens data with mean `μ` and covariance `C`.
"""
whitening_transform(u::AbstractVector, C::AbstractMatrix) = Affine(inv(cholesky(C).L), u)
whitening_transform(C::AbstractMatrix) = Affine(inv(cholesky(C).L))
whitening_transform(u::AbstractVector) = Affine(u)

function whitening_transform!(affine::Affine, μ::AbstractVector, C::AbstractMatrix)
    affine.A .= inv(cholesky(C).L)
    affine.u .= μ
    return affine
end

function whitening_transform!(affine::Affine, μ::AbstractVector)
    copyto!(affine.A, I)
    affine.u .= μ
    return affine
end

function whitening_transform!(affine::Affine, C::AbstractMatrix)
    affine.A .= inv(cholesky(C).L)
    affine.u .= 0
    return affine
end


text:
using RunTests
using Base.Test

@testmodule MultinomialBayesTest begin
  using MLLibrary

  function test_getCounts()
    A=[1. 2. 1. 5.;
       4. 2. 1. 10.;
       3. 3. 2. 6.;]

    @test MLLibrary.LearningAlgorithms.MultinomialBayes.getCounts(A)==[8. 7. 4. 21.];
  end

  function test_getDictionarySize()
    A=[1. 2. 1. 5.;
       4. 2. 1. 10.;
       3. 3. 2. 6.;]

    @test MLLibrary.LearningAlgorithms.MultinomialBayes.getDictionarySize(A)==4
  end

  function test_getClassiferCount()
    A=[1;2;2;1;0;2;3;1;5]
    @test MLLibrary.LearningAlgorithms.MultinomialBayes.getClassifierCount(A)==Dict(0=>1,1=>3,2=>3,3=>1,5=>1)

    A=[1;0;1;0;1;1;0;0;0;1;1]
    @test MLLibrary.LearningAlgorithms.MultinomialBayes.getClassifierCount(A)==Dict(0=>5,1=>6)
  end

  function test_getDictionary()
    A=[1. 2. 1. 5. 1.;
       4. 2. 1. 10. 0.;
       3. 3. 2. 6. 1.;]

    dict,dictProb = MLLibrary.LearningAlgorithms.MultinomialBayes.getDictionary(A[:,1:end-1],A[:,end])
    @test dict==Dict(1=>[4. 5. 3. 11.],0=>[4. 2. 1. 10.])
    for key in keys(dictProb)
      dictProb[key] = round(dictProb[key],2)
    end
    @test dictProb ==Dict(0=>[.24 .14 .10 .52],1=>[.19 .22 .15 .44])
  end

  function test_classifyDocuments()
    A=[1. 2. 1. 5. 1.;
       4. 2. 1. 10. 0.;
       3. 3. 2. 6. 1.;]

    B=[1. 1. 1. 3. 1;
       2. 2. 1. 4. 0.;
       4. 2. 1. 1. 0.;]
    dict,dictProb = MLLibrary.LearningAlgorithms.MultinomialBayes.getDictionary(A[:,1:end-1],A[:,end])
    @test MLLibrary.LearningAlgorithms.MultinomialBayes.classifyDocuments(B[:,1:end-1],dictProb) == [1.;1.;1.]''
  end

end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
function voronoi_metric(x, support)
    a = colwise(PeriodicEuclidean(ones(2)), support, x)
    b = partialsort!(a, 1:2)
    abs(b[2] - b[1])
end

function circle_metric(x, support, r)
    a = abs.(colwise(PeriodicEuclidean(ones(2)), support, x) .- r)
    partialsort!(a, 1)
end


text:
@testset "Execution Control Intrinsics" begin

function completion_signal_kernel(X)
    X[1] = AMDGPU._completion_signal()
    nothing
end

RA = ROCArray(rand(UInt64, 1))

ev = @roc completion_signal_kernel(RA)
wait(ev)
@test collect(RA)[1] == ev.event.signal.signal[].handle

function exec_ctl_kernel(X)
    AMDGPU.sendmsg(5)
    AMDGPU.sendmsghalt(6)
    AMDGPU.endpgm()
end

iob = IOBuffer()
AMDGPU.code_native(iob, exec_ctl_kernel, Tuple{Int64}; kernel=true)
str = String(take!(iob))
@test occursin("s_sendmsg ", str)
@test occursin("s_sendmsghalt ", str)
@test occursin("s_endpgm", str) # TODO: Test for multiple occurences

end


text:
# load library
using GeometricIntegrators

# define vector field
f = (t, x, ẋ) -> ẋ[1] = x[1]

# create ODE for x = f(t,x) with initial condition x(0) = 1.0
ode = ODE(f, [1.0])

# create integrator for explicit Euler method with time step 0.1
int = Integrator(ode, getTableauExplicitEuler(), 0.1)

# compute 10 time steps and retrieve solution
sol = integrate(int, 10)

# plot and compare with the exact solution
using Plots
plot(xlim=[0,1], xlab="t", ylab="x(t)", legend=:bottomright)
plot!(sol.t.t, sol.q.d[1,:], label="numeric")
plot!(sol.t.t, exp.(sol.t.t), label="exact")

savefig("exponential_growth.pdf")


text:
##
using ACE, Printf, Test, LinearAlgebra
using ACE: evaluate, evaluate_d, read_dict, write_dict
using ACE.Testing

verbose = false
maxdeg = 10

##

@info("Testing Transforms and TransformedPolys")
for p in 2:4
   @info("p = $p, random transform")
   trans = PolyTransform(1+rand(), 1+rand())
   @info("      test (de-)dictionisation")
   @test all(ACE.Testing.test_fio(trans))
   B1 = transformed_jacobi(maxdeg, trans, 3.0; pcut = p)
   B2 = transformed_jacobi(maxdeg, trans, 3.0, 0.5, pin = p, pcut = p)
   for B in [B1, B2]
      B == B1 && @info("basis = 1s")
      B == B2 && @info("basis = 2s")
      for r in [3 * rand(10); [3.0]]
         P = evaluate(B, r)
         dP = evaluate_d(B, r)
         errs = []
         verbose && @printf("     h    |     error  \n")
         for p = 2:10
            local h = 0.1^p
            dPh = (evaluate(B, r+h) - P) / h
            push!(errs, norm(dPh - dP, Inf))
            verbose && @printf(" %.2e | %2e \n", h, errs[end])
         end
         print_tf(@test (/(extrema(errs)...) < 1e-3) || (minimum(errs) < 1e-10) )
      end
      println()
   end
end

##

@info("Testing PolyTransforms")
for p = 2:4
   r0 = 1+rand()
   trans = PolyTransform(p, r0)
   ACE.Testing.test_transform(trans, [r0/2, 3*r0])
end
println()

##

@info("Testing Morse Transform")
for lam = 1.0:3.0
   r0 = 1+rand()
   trans = ACE.Transforms.MorseTransform(lam, r0)
   ACE.Testing.test_transform(trans, [r0/2, 3*r0])
end
println()

##

@info("Testing Agnesi Transform")
for p = 2:4
   r0 = 1+rand()
   trans = ACE.Transforms.AgnesiTransform(r0, p)
   ACE.Testing.test_transform(trans, [r0/2, 3*r0])
end
println()

##

# TODO: This could be moved to some Tools package
#       to visualize the good transforms

# using Plots
# r0 = 1.0
# rr = range(0.0, 3*r0, length=200)
# plot(; size = (500, 300))
# for p = 2:4
#    tpoly = PolyTransform(p, r0)
#    tagnesi = ACE.Transforms.AgnesiTransform(r0, p)
#    plot!(rr, tagnesi.(rr), lw=2, c=p-1, label = "p = $p")
#    plot!(rr, tpoly.(rr), lw=2, c=p-1, ls = :dash, label = "")
# end
# xlabel!("r")
# ylabel!("x")
# title!("solid = Agnesi, dashed = Poly")
# vline!([1.0], lw=2, c=:black, label = "r0")
# ylims!(0.0, 2.0)
# #---
#
# function visualize_transform(T, rrange; nrays = 50, c = 1,
#                              inverse = false, straightr0 = false)
#    rnn = T.r0
#    xnn = T(rnn)
#    r0, r1 = extrema(rrange)
#    x0, x1 = T(r0), T(r1)
#    # plot r -> x
#    if !inverse
#       rr = range(r0, r1, length=nrays)
#       xx = T.(rr)
#    else
#       xx = range(x0, x1, length=nrays)
#       rr = ACE.Transforms.inv_transform.(Ref(T), xx)
#    end
#    rnn = (rnn - rr[1]) / (rr[end] - rr[1])
#    xnn = (xnn - xx[1]) / (xx[end] - xx[1])
#    rr = (rr .- rr[1]) / (rr[end] - rr[1])
#    xx = (xx .- xx[1]) / (xx[end] - xx[1])
#    if straightr0
#       rr = rr .- rnn; rnn = 0.0
#       xx = xx .- xnn; xnn = 0.0
#    end
#    plot(; size = (500, 200), ylims = (-1, 2), xticks = [])
#    for (r, x) in zip(rr, xx)
#       plot!([r, x], [0, 1], lw=1, label = "", c = c, )
#    end
#    plot!([rnn, xnn], [0, 1], lw=3, c = :red, label = "")
#    plot!([rr[1], rr[end]], [0, 0], lw=3, label = "", c = :black)
#    plot!([xx[1], xx[end]], [1, 1], lw=3, label = "", c = :black)
#    yticks!([0, 1], ["r", "x"])
# end
#
#
#
# tagnesi = AgnesiTransform(1.0, 3)
# P1 = visualize_transform(tagnesi, (0.3*r0, 3*r0), straightr0=true)
# title!(P1, "Agnesi-3")
#
# tagnesi4 = AgnesiTransform(1.0, 4)
# P3 = visualize_transform(tagnesi4, (0.3*r0, 3*r0), straightr0=true)
# title!(P3, "Agnesi-4")
#
# tpoly = PolyTransform(2, 1.0)
# P2 = visualize_transform(tpoly, (0.3*r0, 3*r0), straightr0=true)
# title!(P2, "Poly-2")
# plot(P2, P1, P3, layout = (3,1), size = (500, 600))
#

#---


text:
# Utils

"Get best Bid offer from the buy order book (order with highest price)"
function get_best(buy_book::Array{BidLimitOrder, 1})
    index_for_best = argmax(map(x -> x.price, buy_book))
    return buy_book[index_for_best]
end

"Get best Ask offer from the buy order book (order with lowest price)"
function get_best(sell_book::Array{AskLimitOrder, 1})
    index_for_best = argmin(map(x -> x.price, sell_book))
    return sell_book[index_for_best]
end

function minprice(orders::Array{<:LimitOrder, 1})
    if isempty(orders)
        return NaN
    end
    prices = map(order -> order.price, orders)
    minimum(prices)
end

function maxprice(orders::Array{<:LimitOrder, 1})
    if isempty(orders)
        return NaN
    end
    prices = map(order -> order.price, orders)
    maximum(prices)
end

best_order_index(book::Array{BidLimitOrder, 1}) = argmax(map(x -> x.price, book))
best_order_index(book::Array{AskLimitOrder, 1}) = argmin(map(x -> x.price, book))

"Clear empty orders (orders with quantity of 0)"
function clear_empty!(book::Array{<:LimitOrder, 1})
    empty_orders = findall(x -> x.quantity <= 0, book)
    deleteat!(book, empty_orders)
end

function get_prices(orders::Array{<:LimitOrder, 1})
    return [order.price for order in orders]
end

function get_equlibrium(bids::Array{BidLimitOrder, 1}, asks::Array{AskLimitOrder, 1})
    clearable_quantities = Dict{Int64, Int64}()

    prices = Set(vcat(get_prices(bids), get_prices(asks)))

    for price_level in prices
        #println(price_level)
        bids_clearable = [bid.quantity for bid in bids if bid.price >= price_level]

        asks_clearable = [ask.quantity for ask in asks if ask.price <= price_level]
        clearable_quantity = minimum([sum(bids_clearable), sum(asks_clearable)])

        #clearable_quantities = vcat(clearable_quantities, [price_level clearable_quantity])
        clearable_quantities[price_level] = clearable_quantity
    end

    #eq_quantity = findmax(clearable_quantities, dims=2)
    val, ind = findmax(clearable_quantities)
    eq_prices = [key for (key, val) in clearable_quantities if val == eq_quantity]

    eq_quantity = val
    eq_price = round(mean([minimum(eq_prices), maximum(eq_prices)]))
    return eq_quantity, eq_price
end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((2^(1/2)+cos(z)+sin(z))^-1, z) == :(-1*(-1*sin(z)+cos(z))^-1*(1+-1*2^(1/2)*sin(z)))
@test integrate(((1+x)^(1/2)+(1+-1x)^(1/2))^-2, x) == :((1/2)*arcsin(x)+-1//2*x^-1+(1/2)*x^-1*(1+-1*x^2)^(1/2))
@test integrate((1+cos(x))^-2, x) == :((3+3*cos(x))^-1*sin(x)+1//3*(1+cos(x))^-2*sin(x))
@test integrate((1+x)^-1//2*sin(x), x) == :(2^(1/2)*pi^(1/2)*cos(1)*fresnels(2^(1/2)*pi^-1//2*(1+x)^(1/2))+-1*2^(1/2)*pi^(1/2)*fresnelc(2^(1/2)*pi^-1//2*(1+x)^(1/2))*sin(1))
@test integrate((cos(x)+sin(x))^-6, x) == :(2*(15*cos(x)+15*sin(x))^-1*sin(x)+-1//10*(cos(x)+sin(x))^-5*(-1*sin(x)+cos(x))+-1//15*(cos(x)+sin(x))^-3*(-1*sin(x)+cos(x)))
@test integrate(log(x^-4+x^4), x) == :(-4x+x*log(x^-4+x^4)+(2+2^(1/2))^(1/2)*arctan((2+2^(1/2))^-1//2*((2+-1*2^(1/2))^(1/2)+2x))+(2+-1*2^(1/2))^(1/2)*arctan((2+-1*2^(1/2))^-1//2*((2+2^(1/2))^(1/2)+2x))+(1/2)*(2+2^(1/2))^(1/2)*log(1+x^2+x*(2+2^(1/2))^(1/2))+(1/2)*(2+-1*2^(1/2))^(1/2)*log(1+x^2+x*(2+-1*2^(1/2))^(1/2))+-1*(2+2^(1/2))^(1/2)*arctan((2+2^(1/2))^-1//2*((2+-1*2^(1/2))^(1/2)+-2x))+-1*(2+-1*2^(1/2))^(1/2)*arctan((2+-1*2^(1/2))^-1//2*((2+2^(1/2))^(1/2)+-2x))+-1//2*(2+2^(1/2))^(1/2)*log(1+x^2+-1*x*(2+2^(1/2))^(1/2))+-1//2*(2+-1*2^(1/2))^(1/2)*log(1+x^2+-1*x*(2+-1*2^(1/2))^(1/2)))
@test integrate(x^-1*(1+(1+x)^(1/2))^-1//2*log(1+x), x) == :(-8*arctanh((1+(1+x)^(1/2))^(1/2))+2^(1/2)*Polylogarithms.polylog(2,2^(1/2)*(2+2^(1/2))^-1*(1+(1+(1+x)^(1/2))^(1/2)))+2^(1/2)*Polylogarithms.polylog(2,-1*2^(1/2)*(2+-1*2^(1/2))^-1*(1+-1*(1+(1+x)^(1/2))^(1/2)))+-1*2^(1/2)*Polylogarithms.polylog(2,2^(1/2)*(2+2^(1/2))^-1*(1+-1*(1+(1+x)^(1/2))^(1/2)))+-1*2^(1/2)*Polylogarithms.polylog(2,-1*2^(1/2)*(2+-1*2^(1/2))^-1*(1+(1+(1+x)^(1/2))^(1/2)))+-2*(1+(1+x)^(1/2))^-1//2*log(1+x)+-1*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+(1+x)^(1/2))^(1/2))*log(1+x)+-2*2^(1/2)*arctanh((1/2)*2^(1/2))*log(1+(1+(1+x)^(1/2))^(1/2))+2*2^(1/2)*arctanh((1/2)*2^(1/2))*log(1+-1*(1+(1+x)^(1/2))^(1/2)))
@test integrate(x^-1*(1+(1+x)^(1/2))^(1/2)*log(1+x), x) == :(-16*(1+(1+x)^(1/2))^(1/2)+16*arctanh((1+(1+x)^(1/2))^(1/2))+-2*2^(1/2)*Polylogarithms.polylog(2,2^(1/2)*(2+2^(1/2))^-1*(1+-1*(1+(1+x)^(1/2))^(1/2)))+-2*2^(1/2)*Polylogarithms.polylog(2,-1*2^(1/2)*(2+-1*2^(1/2))^-1*(1+(1+(1+x)^(1/2))^(1/2)))+2*2^(1/2)*Polylogarithms.polylog(2,2^(1/2)*(2+2^(1/2))^-1*(1+(1+(1+x)^(1/2))^(1/2)))+2*2^(1/2)*Polylogarithms.polylog(2,-1*2^(1/2)*(2+-1*2^(1/2))^-1*(1+-1*(1+(1+x)^(1/2))^(1/2)))+4*(1+(1+x)^(1/2))^(1/2)*log(1+x)+-4*2^(1/2)*arctanh((1/2)*2^(1/2))*log(1+(1+(1+x)^(1/2))^(1/2))+-2*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+(1+x)^(1/2))^(1/2))*log(1+x)+4*2^(1/2)*arctanh((1/2)*2^(1/2))*log(1+-1*(1+(1+x)^(1/2))^(1/2)))
@test integrate((1+(x+(1+x^2)^(1/2))^(1/2))^-1, x) == :((x+(1+x^2)^(1/2))^(1/2)+(x+(1+x^2)^(1/2))^-1//2+(1/2)*log(x+(1+x^2)^(1/2))+-1*(2x+2*(1+x^2)^(1/2))^-1+-2*log(1+(x+(1+x^2)^(1/2))^(1/2)))
@test integrate((1+x)^(1/2)*(x+(1+(1+x)^(1/2))^(1/2))^-1, x) == :(2*(1+x)^(1/2)+8//5*5^(1/2)*arctanh(1//5*5^(1/2)*(1+2*(1+(1+x)^(1/2))^(1/2))))
@test integrate((x+-1*(1+(1+x)^(1/2))^(1/2))^-1, x) == :((2+-2//5*5^(1/2))*log(1+5^(1/2)+-2*(1+(1+x)^(1/2))^(1/2))+(2+2//5*5^(1/2))*log(1+-1*5^(1/2)+-2*(1+(1+x)^(1/2))^(1/2)))
@test integrate(x*(x+(1+-1*(1+x)^(1/2))^(1/2))^-1, x) == :((1+-1*(1+x)^(1/2))^2+-4*(1+-1*(1+x)^(1/2))^(1/2)+2*(1+x)^(1/2)+8//5*5^(1/2)*arctanh(1//5*5^(1/2)*(1+2*(1+-1*(1+x)^(1/2))^(1/2))))
@test integrate((1+x)^-1//2*(1+x^2)^-1*(x+(1+x)^(1/2))^(1/2), x) == :((1/2)*im*(-1*(im+-1*(1+im)^(1/2))^-1*(1+im))^-1//2*arctan((1/2)*(x+(1+x)^(1/2))^-1//2*((1+im)^(1/2)+-1im)^-1//2*(2+(1+im)^(1/2)+-1*(1+x)^(1/2)*(1+-2*(1+im)^(1/2))))+(1/2)*im*(-1*(im+-1*(1+-1im)^(1/2))^-1*(1+-1im))^-1//2*arctanh((1/2)*(x+(1+x)^(1/2))^-1//2*((1+-1im)^(1/2)+-1im)^-1//2*(2+-1*(1+-1im)^(1/2)+-1*(1+x)^(1/2)*(1+2*(1+-1im)^(1/2))))+-1//2*im*((im+(1+im)^(1/2))^-1*(1+im))^-1//2*arctanh((1/2)*(im+(1+im)^(1/2))^-1//2*(x+(1+x)^(1/2))^-1//2*(2+-1*(1+im)^(1/2)+-1*(1+x)^(1/2)*(1+2*(1+im)^(1/2))))+-1//2*im*((im+(1+-1im)^(1/2))^-1*(1+-1im))^-1//2*arctan((1/2)*(im+(1+-1im)^(1/2))^-1//2*(x+(1+x)^(1/2))^-1//2*(2+(1+-1im)^(1/2)+-1*(1+x)^(1/2)*(1+-2*(1+-1im)^(1/2)))))
@test integrate((1+x^2)^-1*(x+(1+x)^(1/2))^(1/2), x) == :((1/2)*im*(im+(1+-1im)^(1/2))^(1/2)*arctan((1/2)*(im+(1+-1im)^(1/2))^-1//2*(x+(1+x)^(1/2))^-1//2*(2+(1+-1im)^(1/2)+-1*(1+x)^(1/2)*(1+-2*(1+-1im)^(1/2))))+(1/2)*im*((1+-1im)^(1/2)+-1im)^(1/2)*arctanh((1/2)*(x+(1+x)^(1/2))^-1//2*((1+-1im)^(1/2)+-1im)^-1//2*(2+-1*(1+-1im)^(1/2)+-1*(1+x)^(1/2)*(1+2*(1+-1im)^(1/2))))+-1//2*im*(im+(1+im)^(1/2))^(1/2)*arctanh((1/2)*(im+(1+im)^(1/2))^-1//2*(x+(1+x)^(1/2))^-1//2*(2+-1*(1+im)^(1/2)+-1*(1+x)^(1/2)*(1+2*(1+im)^(1/2))))+-1//2*im*((1+im)^(1/2)+-1im)^(1/2)*arctan((1/2)*(x+(1+x)^(1/2))^-1//2*((1+im)^(1/2)+-1im)^-1//2*(2+(1+im)^(1/2)+-1*(1+x)^(1/2)*(1+-2*(1+im)^(1/2)))))
@test integrate((1+x^(1/2)+(1+2x+2*x^(1/2))^(1/2))^(1/2), x) == :(2//15*x^-1//2*(1+x^(1/2)+(1+2x+2*x^(1/2))^(1/2))^(1/2)*(2+x^(1/2)+6*x^3//2+-1*(1+2x+2*x^(1/2))^(1/2)*(2+-1*x^(1/2))))
@test integrate((2^(1/2)+x^(1/2)+(2+2x+2*2^(1/2)*x^(1/2))^(1/2))^(1/2), x) == :(2//15*2^(1/2)*x^-1//2*(2^(1/2)+x^(1/2)+2^(1/2)*(1+x+2^(1/2)*x^(1/2))^(1/2))^(1/2)*(4+2^(1/2)*x^(1/2)+3*2^(1/2)*x^3//2+-1*2^(1/2)*(1+x+2^(1/2)*x^(1/2))^(1/2)*(-1*x^(1/2)+2*2^(1/2))))
@test integrate(x^-2*(x+(1+x)^(1/2))^(1/2), x) == :(-1//4*arctan((1/2)*(x+(1+x)^(1/2))^-1//2*(3+(1+x)^(1/2)))+3//4*arctanh((1/2)*(x+(1+x)^(1/2))^-1//2*(1+-3*(1+x)^(1/2)))+-1*x^-1*(x+(1+x)^(1/2))^(1/2))
@test integrate((x^-1+(1+x^-1)^(1/2))^(1/2), x) == :(-3//4*arctanh((1/2)*(x^-1+(1+x^-1)^(1/2))^-1//2*(1+-3*(1+x^-1)^(1/2)))+1//4*arctan((1/2)*(x^-1+(1+x^-1)^(1/2))^-1//2*(3+(1+x^-1)^(1/2)))+x*(x^-1+(1+x^-1)^(1/2))^(1/2))
@test integrate((1+exp(-1x))^(1/2)*(-1*exp(-1x)+exp(x))^-1, x) == :(-1*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+exp(-1x))^(1/2)))
@test integrate((1+exp(-1x))^(1/2)*sinh(x)^-1, x) == :(-2*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+exp(-1x))^(1/2)))
@test integrate((cos(x)+cos(3x))^-5, x) == :(-523//256*arctanh(sin(x))+-437*(512+-1024*sin(x)^2)^-1*sin(x)+-43//256*sec(x)*tan(x)+-17//192*(1+-2*sin(x)^2)^-3*sin(x)+-1//128*sec(x)^3*tan(x)+1//32*(1+-2*sin(x)^2)^-4*sin(x)+203//768*(1+-2*sin(x)^2)^-2*sin(x)+1483//1024*2^(1/2)*arctanh(2^(1/2)*sin(x)))
@test integrate((1+cos(x)+sin(x))^-2, x) == :(-1*log(1+tan((1/2)*x))+-1*(1+cos(x)+sin(x))^-1*(-1*sin(x)+cos(x)))
@test integrate((1+tanh(4x))^(1/2), x) == :(1//4*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+tanh(4x))^(1/2)))
@test integrate((exp(x)+exp(2x))^-1//2*tanh(x), x) == :((1+-1im)^-1//2*arctan((1/2)*(1+-1im)^-1//2*(exp(x)+exp(2x))^-1//2*(im+(1+2im)*exp(x)))+-1*(1+im)^-1//2*arctan((1/2)*(1+im)^-1//2*(exp(x)+exp(2x))^-1//2*(im+-1*(1+-2im)*exp(x)))+2*(exp(x)+exp(2x))^(1/2)*exp(-1x))
@test integrate(log(x^2+(1+-1*x^2)^(1/2)), x) == :(-1*arcsin(x)+-2x+x*log(x^2+(1+-1*x^2)^(1/2))+(1/2+(1/2)*5^(1/2))^(1/2)*arctan(x*2^(1/2)*(1+5^(1/2))^-1//2)+(1/2+(1/2)*5^(1/2))^(1/2)*arctan(x*(1+-1*x^2)^-1//2*(1/2+(1/2)*5^(1/2))^(1/2))+(-1//2+(1/2)*5^(1/2))^(1/2)*arctanh(x*2^(1/2)*(-1+5^(1/2))^-1//2)+-1*(-1//2+(1/2)*5^(1/2))^(1/2)*arctanh(x*(1+-1*x^2)^-1//2*(-1//2+(1/2)*5^(1/2))^(1/2)))
@test integrate((1+exp(2x))^-1*log(1+exp(x)), x) == :(-1*Polylogarithms.polylog(2,-1*exp(x))+-1//2*Polylogarithms.polylog(2,(1+exp(x))*(1/2+(1/2)*im))+-1//2*Polylogarithms.polylog(2,(1+exp(x))*(1/2+-1//2*im))+-1//2*log((1/2+-1//2*im)*(im+-1*exp(x)))*log(1+exp(x))+-1//2*log((-1//2+-1//2*im)*(im+exp(x)))*log(1+exp(x)))
@test integrate(log(1+cosh(x)^2)^2*cosh(x), x) == :(8*sinh(x)+log(2+sinh(x)^2)^2*sinh(x)+-8*2^(1/2)*arctan((1/2)*2^(1/2)*sinh(x))+-4*log(2+sinh(x)^2)*sinh(x)+4*im*2^(1/2)*arctan((1/2)*2^(1/2)*sinh(x))^2+4*im*2^(1/2)*Polylogarithms.polylog(2,1+-2*2^(1/2)*(2^(1/2)+im*sinh(x))^-1)+4*2^(1/2)*arctan((1/2)*2^(1/2)*sinh(x))*log(2+sinh(x)^2)+8*2^(1/2)*arctan((1/2)*2^(1/2)*sinh(x))*log(2*2^(1/2)*(2^(1/2)+im*sinh(x))^-1))
@test integrate(log(cosh(x)^2+sinh(x))^2*cosh(x), x) == :(-2*log(1+sinh(x)^2+sinh(x))+8*sinh(x)+log(1+sinh(x)^2+sinh(x))^2*sinh(x)+-1*log(1+2*sinh(x)+im*3^(1/2))^2*(1/2+(1/2)*im*3^(1/2))+-1*log(1+2*sinh(x)+-1*im*3^(1/2))^2*(1/2+-1//2*im*3^(1/2))+-1*(1+im*3^(1/2))*Polylogarithms.polylog(2,-1//6*3^(1/2)*(im+-1*3^(1/2)+2*im*sinh(x)))+-1*(1+-1*im*3^(1/2))*Polylogarithms.polylog(2,1//6*3^(1/2)*(im+3^(1/2)+2*im*sinh(x)))+-4*3^(1/2)*arctan(1//3*3^(1/2)*(1+2*sinh(x)))+-4*log(1+sinh(x)^2+sinh(x))*sinh(x)+(1+im*3^(1/2))*log(1+sinh(x)^2+sinh(x))*log(1+2*sinh(x)+im*3^(1/2))+(1+-1*im*3^(1/2))*log(1+sinh(x)^2+sinh(x))*log(1+2*sinh(x)+-1*im*3^(1/2))+-1*(1+im*3^(1/2))*log(1//6*im*3^(1/2)*(1+2*sinh(x)+-1*im*3^(1/2)))*log(1+2*sinh(x)+im*3^(1/2))+-1*(1+-1*im*3^(1/2))*log(-1//6*im*3^(1/2)*(1+2*sinh(x)+im*3^(1/2)))*log(1+2*sinh(x)+-1*im*3^(1/2)))
@test integrate((1+x^2)^-1*log(x+(1+x)^(1/2)), x) == :((1/2)*im*Polylogarithms.polylog(2,(1+5^(1/2)+2*(1+im)^(1/2))^-1*(-2*(1+x)^(1/2)+2*(1+im)^(1/2)))+(1/2)*im*Polylogarithms.polylog(2,(1+-1*5^(1/2)+2*(1+im)^(1/2))^-1*(-2*(1+x)^(1/2)+2*(1+im)^(1/2)))+(1/2)*im*Polylogarithms.polylog(2,-1*(1+5^(1/2)+-2*(1+im)^(1/2))^-1*(2*(1+im)^(1/2)+2*(1+x)^(1/2)))+(1/2)*im*Polylogarithms.polylog(2,-1*(1+-1*5^(1/2)+-2*(1+im)^(1/2))^-1*(2*(1+im)^(1/2)+2*(1+x)^(1/2)))+-1//2*im*Polylogarithms.polylog(2,(1+5^(1/2)+2*(1+-1im)^(1/2))^-1*(-2*(1+x)^(1/2)+2*(1+-1im)^(1/2)))+-1//2*im*Polylogarithms.polylog(2,(1+-1*5^(1/2)+2*(1+-1im)^(1/2))^-1*(-2*(1+x)^(1/2)+2*(1+-1im)^(1/2)))+-1//2*im*Polylogarithms.polylog(2,-1*(1+5^(1/2)+-2*(1+-1im)^(1/2))^-1*(2*(1+x)^(1/2)+2*(1+-1im)^(1/2)))+-1//2*im*Polylogarithms.polylog(2,-1*(1+-1*5^(1/2)+-2*(1+-1im)^(1/2))^-1*(2*(1+x)^(1/2)+2*(1+-1im)^(1/2)))+(1/2)*im*log((1+5^(1/2)+-2*(1+im)^(1/2))^-1*(1+5^(1/2)+2*(1+x)^(1/2)))*log((1+im)^(1/2)+(1+x)^(1/2))+(1/2)*im*log((1+5^(1/2)+2*(1+im)^(1/2))^-1*(1+5^(1/2)+2*(1+x)^(1/2)))*log((1+im)^(1/2)+-1*(1+x)^(1/2))+(1/2)*im*log((1+-1*5^(1/2)+-2*(1+im)^(1/2))^-1*(1+-1*5^(1/2)+2*(1+x)^(1/2)))*log((1+im)^(1/2)+(1+x)^(1/2))+(1/2)*im*log((1+-1*5^(1/2)+2*(1+im)^(1/2))^-1*(1+-1*5^(1/2)+2*(1+x)^(1/2)))*log((1+im)^(1/2)+-1*(1+x)^(1/2))+(1/2)*im*log(x+(1+x)^(1/2))*log((1+x)^(1/2)+(1+-1im)^(1/2))+(1/2)*im*log(x+(1+x)^(1/2))*log((1+-1im)^(1/2)+-1*(1+x)^(1/2))+-1//2*im*log((1+5^(1/2)+-2*(1+-1im)^(1/2))^-1*(1+5^(1/2)+2*(1+x)^(1/2)))*log((1+x)^(1/2)+(1+-1im)^(1/2))+-1//2*im*log((1+5^(1/2)+2*(1+-1im)^(1/2))^-1*(1+5^(1/2)+2*(1+x)^(1/2)))*log((1+-1im)^(1/2)+-1*(1+x)^(1/2))+-1//2*im*log((1+-1*5^(1/2)+-2*(1+-1im)^(1/2))^-1*(1+-1*5^(1/2)+2*(1+x)^(1/2)))*log((1+x)^(1/2)+(1+-1im)^(1/2))+-1//2*im*log((1+-1*5^(1/2)+2*(1+-1im)^(1/2))^-1*(1+-1*5^(1/2)+2*(1+x)^(1/2)))*log((1+-1im)^(1/2)+-1*(1+x)^(1/2))+-1//2*im*log(x+(1+x)^(1/2))*log((1+im)^(1/2)+(1+x)^(1/2))+-1//2*im*log(x+(1+x)^(1/2))*log((1+im)^(1/2)+-1*(1+x)^(1/2)))
@test integrate((1+x)^-2*log(x+(1+x)^(1/2))^2, x) == :(-6*Polylogarithms.polylog(2,1+2*(1+x)^(1/2)*(1+-1*5^(1/2))^-1)+6*Polylogarithms.polylog(2,-2*(1+x)^(1/2)*(1+5^(1/2))^-1)+-1*(1+x)^-1*log(x+(1+x)^(1/2))^2+-1*log(1+5^(1/2)+2*(1+x)^(1/2))^2*(3//2+-1//2*5^(1/2))+-1*log(1+-1*5^(1/2)+2*(1+x)^(1/2))^2*(3//2+(1/2)*5^(1/2))+-1*(1+5^(1/2))*log(1+-1*5^(1/2)+2*(1+x)^(1/2))+-1*(1+-1*5^(1/2))*log(1+5^(1/2)+2*(1+x)^(1/2))+-1*(3+5^(1/2))*Polylogarithms.polylog(2,-1//10*5^(1/2)*(1+-1*5^(1/2)+2*(1+x)^(1/2)))+-1*(3+-1*5^(1/2))*Polylogarithms.polylog(2,1//10*5^(1/2)*(1+5^(1/2)+2*(1+x)^(1/2)))+-6*log((1+x)^(1/2))*log(x+(1+x)^(1/2))+2*(1+x)^-1//2*log(x+(1+x)^(1/2))+6*log((1+x)^(1/2))*log(1+2*(1+x)^(1/2)*(1+5^(1/2))^-1)+6*log(-1//2+(1/2)*5^(1/2))*log(1+-1*5^(1/2)+2*(1+x)^(1/2))+(3+5^(1/2))*log(x+(1+x)^(1/2))*log(1+-1*5^(1/2)+2*(1+x)^(1/2))+(3+-1*5^(1/2))*log(x+(1+x)^(1/2))*log(1+5^(1/2)+2*(1+x)^(1/2))+-1*(3+5^(1/2))*log(1//10*5^(1/2)*(1+5^(1/2)+2*(1+x)^(1/2)))*log(1+-1*5^(1/2)+2*(1+x)^(1/2))+-1*(3+-1*5^(1/2))*log(-1//10*5^(1/2)*(1+-1*5^(1/2)+2*(1+x)^(1/2)))*log(1+5^(1/2)+2*(1+x)^(1/2))+log(1+x))
@test integrate(x^-1*log(x+(1+x)^(1/2)), x) == :(-1*Polylogarithms.polylog(2,(1+5^(1/2))^-1*(2+2*(1+x)^(1/2)))+-1*Polylogarithms.polylog(2,(1+-1*5^(1/2))^-1*(2+2*(1+x)^(1/2)))+-1*Polylogarithms.polylog(2,(3+5^(1/2))^-1*(2+-2*(1+x)^(1/2)))+-1*Polylogarithms.polylog(2,(3+-1*5^(1/2))^-1*(2+-2*(1+x)^(1/2)))+log(1+(1+x)^(1/2))*log(x+(1+x)^(1/2))+log(-1+(1+x)^(1/2))*log(x+(1+x)^(1/2))+-1*log((3+5^(1/2))^-1*(1+5^(1/2)+2*(1+x)^(1/2)))*log(-1+(1+x)^(1/2))+-1*log((3+-1*5^(1/2))^-1*(1+-1*5^(1/2)+2*(1+x)^(1/2)))*log(-1+(1+x)^(1/2))+-1*log(-1*(1+5^(1/2))^-1*(1+-1*5^(1/2)+2*(1+x)^(1/2)))*log(1+(1+x)^(1/2))+-1*log(-1*(1+-1*5^(1/2))^-1*(1+5^(1/2)+2*(1+x)^(1/2)))*log(1+(1+x)^(1/2)))
@test integrate(arctan(2*tan(x)), x) == :(-1//4*Polylogarithms.polylog(2,1//3*exp(2*im*x))+1//4*Polylogarithms.polylog(2,3*exp(2*im*x))+x*arctan(2*tan(x))+(1/2)*im*x*log(1+-3*exp(2*im*x))+-1//2*im*x*log(1+-1//3*exp(2*im*x)))
@test integrate(x^-1*arctan(x)*log(x), x) == :((1/2)*im*Polylogarithms.polylog(3,im*x)+-1//2*im*Polylogarithms.polylog(3,-1*im*x)+(1/2)*im*Polylogarithms.polylog(2,-1*im*x)*log(x)+-1//2*im*Polylogarithms.polylog(2,im*x)*log(x))
@test integrate((1+x^2)^(1/2)*arctan(x)^2, x) == :(-1*Polylogarithms.polylog(3,-1*im*exp(im*arctan(x)))+-1*(1+x^2)^(1/2)*arctan(x)+im*Polylogarithms.polylog(2,-1*im*exp(im*arctan(x)))*arctan(x)+(1/2)*x*(1+x^2)^(1/2)*arctan(x)^2+-1*im*arctan(x)^2*arctan(exp(im*arctan(x)))+-1*im*Polylogarithms.polylog(2,im*exp(im*arctan(x)))*arctan(x)+Polylogarithms.polylog(3,im*exp(im*arctan(x)))+arcsinh(x))

max_stars_repo_path: test/0 Independent test suites/Bondarenko Problems.jl


text:
"""
    CategoricalLikelihood

Categorical likelihood is to be used if we assume that the
uncertainity associated with the data follows a Categorical distribution.
```math
    p(y|f_1, f_2, \\dots, f_{n-1}) = Categorical(y | softmax(f_1, f_2, \\dots, f_{n-1}, 0))
```
On calling, this would return a Categorical distribution with `f_i`
probability of `i` category.
"""
struct CategoricalLikelihood end

(l::CategoricalLikelihood)(f::AbstractVector{<:Real}) = Categorical(softmax(vcat(f, 0)))

(l::CategoricalLikelihood)(fs::AbstractVector) = Product(Categorical.(softmax.(vcat.(fs, 0))))