text:
using Base.Test
using NonlinearOptics
import Unitful: nm, ps, THz, m, c0, °, W

detected_tests = filter(name->startswith(name, "test_") && endswith(name, ".jl"),
                        readdir(@__DIR__))

@time @testset "NonlinearOptics" begin
    @time @testset "$(test[6:(end-3)])" for test in detected_tests
        include(test)
    end
end

ext: jl


text:
module LoadGraph

using DataStructures: DefaultDict
using DelimitedFiles: readdlm

export load_adjacency_matrix,
       load_edge_list,
       load_adjacency_list,
       AdjacencyList

include("adjacency_matrix.jl")
include("edge_list.jl")
include("adjacency_list.jl")
include("utils.jl")

end

max_stars_repo_name: quentin-burthier/GrandsGraphes

max_issues_repo_name: quentin-burthier/GrandsGraphes

max_forks_repo_name: quentin-burthier/GrandsGraphes


text:
#xortraintest-1.jl

#trains an xornet using noisy data.  Ten data inputs, top two xor'ed to get the
#correct values.

#encapsulate in a function so that the GC can take care of these variables.
function xortrain_1()

  println("working on noisy xor dataset")

  input_matrix = rand(Bool, 10, 500)
  expected_results = Array{Bool,2}(1,500)
  expected_results[:] = [input_matrix[1, col] $ input_matrix[2, col] for col in 1:size(input_matrix,2)]

  xornet = GenML.MLP.MultilayerPerceptron{Float64,(10,2,1)}(randn)

  GenML.Optimizers.gradientoptimize(xornet, input_matrix, expected_results, GenML.CF.crossentropy)

  #verify that the optimization has resulted in a good data set.
  wrongcount = 0
  for x = 1:50
    input_vector = rand(Bool, 10)
    wrongcount += (xornet(input_vector)[1] > 0.5) != (input_vector[1] $ input_vector[2])
  end
  println("incorrect responses for noisy xor: $wrongcount")

  return wrongcount
end

xortrain_1() == 0 || xortrain_1() == 0 || throw(ErrorException("two executions failed."))

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML


text:
abstract type ExtendedFloat{P} <: AbstractFloat end

abstract type Interval end

abstract type BallInterval  <: Interval end  # (midpoint, radius)
abstract type BoundInterval <: Interval end  # (lowerbound, upperbound)

struct Midpoint <: BallInterval end
struct Radius   <: BallInterval end

struct LowerBound <: BoundInterval end
struct UpperBound <: BoundInterval end

struct RealPart <: Number end
struct ImagPart <: Number end


text:
using Test
using TensorGames
using LinearAlgebra: norm
using Zygote: gradient, forwarddiff
using Random: MersenneTwister

include("test_equilibrium.jl")
include("test_derivatives.jl")


text:
module Tofu

export ◻

import LinearAlgebra
import REPL

togetfield(ex) = ex
togetfield(ex::Expr) =
    if ex.head == :. && ex.args[1] == :g
        @assert length(ex.args) == 2
        :($getfield(g, $(ex.args[2])))
    else
        Expr(ex.head, togetfield.(ex.args)...)
    end

"""
  @G ex

Convert `g.PROPERTY` in `ex` to `getfield(g, :PROPERTY)`.
"""
macro G(ex)
    esc(togetfield(ex))
end

# --- Call graph

abstract type Graph end
struct Hole <: Graph end
struct Call{F,A,K} <: Graph
    f::F
    args::A
    kwargs::K
end
call(f, args...; kwargs...) = Call(f, args, kwargs.data)

Base.getproperty(g::Graph, name::Symbol) = call(getproperty, g, name)
Base.getindex(g::Graph, idx...) = call(getindex, g, idx...)

const binop_symbols = [
    :(=>),
    :>, :<, :≥, :≤, :(==), :!=, :∈, :∉, :∋, :∌,
    :+, :-, :|,
    :*, :/, :÷, :%, :&, :⋅, :×, :\,
    ://,
    :<<, :>>, :>>>,
    :^,
]
const binop_map = Dict(
    if n in (:⋅, :×)
        getproperty(LinearAlgebra, n) => n
    else
        getproperty(Base, n) => n
    end
    for n in binop_symbols)
const binop_types = Union{typeof.(keys(binop_map))...}

for (f, name) in binop_map
    mod = parentmodule(f)
    @eval $mod.$name(x::Graph, y) = call($f, x, y)
    @eval $mod.$name(x, y::Graph) = call($f, x, y)
    @eval $mod.$name(x::Graph, y::Graph) = call($f, x, y)
end

# TODO: use https://github.com/JuliaDiff/DiffRules.jl
# (or maybe https://github.com/JuliaDiff/ChainRules.jl if it's ready)

const unaryops = [
    ~,
    adjoint,
]

for f in unaryops
    mod = parentmodule(f)
    name = nameof(f)
    @eval $mod.$name(x::Graph) = call($f, x)
end

Base.:(:)(x::Graph, y) = call(:, x, y)
Base.:(:)(x, y::Graph) = call(:, x, y)
Base.:(:)(x::Graph, y::Graph) = call(:, x, y)
Base.:(:)(x::Graph, y, z) = call(:, x, y, z)
Base.:(:)(x, y::Graph, z) = call(:, x, y, z)
Base.:(:)(x, y, z::Graph) = call(:, x, y, z)
Base.:(:)(x::Graph, y::Graph, z) = call(:, x, y, z)
Base.:(:)(x::Graph, y, z::Graph) = call(:, x, y, z)
Base.:(:)(x, y::Graph, z::Graph) = call(:, x, y, z)
Base.:(:)(x::Graph, y::Graph, z::Graph) = call(:, x, y, z)

# Disambiguation:
Base.:(:)(x::T, y::Graph, z::T) where {T <: Real} = call(:, x, y, z)
Base.:(:)(x::T, y::Graph, z::T) where {T} = call(:, x, y, z)
Base.:(:)(x::Real, y::Graph, z::Real) = call(:, x, y, z)

# --- Evaluation

(g::Hole)(x) = materialize(g, x)
(g::Call)(x) = materialize(g, x)

materialize(y, _) = y
materialize(g::Hole, x) = x
materialize(g::Call, x) = @G g.f(feed(x, g.args...)...; g.kwargs...)

feed(x) = ()
feed(x, a, args...) = (materialize(a, x), feed(x, args...)...)

# --- Show it like you build it

Base.show(io::IO, ::Hole) = print(io, "◻")
Base.show(io::IO, g::Call) = @G show_impl(io, g.f, g.args, g.kwargs)

function show_impl(io, f, args, kwargs)
    print(io, f, '(')
    if length(args) > 0
        show(io, args[1])
        for a in args[2:end]
            print(io, ", ")
            show(io, a)
        end
    end
    if length(kwargs) > 0
    end
    print(io, ')')
end

function show_impl(io, f::binop_types, args, kwargs)
    @assert length(kwargs) == 0
    if length(args) > 0
        show_term(io, f, args[1])
        for a in args[2:end]
            print(io, ' ', binop_map[f], ' ')
            show_term(io, f, a)
        end
    end
end

function show_term(io, f, a)
    if need_paren(f, a)
        print(io, '(')
        show(io, a)
        print(io, ')')
    else
        show(io, a)
    end
end

need_paren(f, g) = false
need_paren(f, g::Call) = @G g.f isa binop_types && g.f !== f

function show_impl(io, ::typeof(getproperty), args, kwargs)
    value, (name::Symbol) = args
    @assert length(args) == 2
    @assert length(kwargs) == 0
    show(io, value)
    print(io, '.')
    print(io, name)
    return
end

function show_impl(io, ::typeof(getindex), args, kwargs)
    @assert length(kwargs) == 0
    show(io, args[1])
    print(io, '[')
    if length(args) > 1
        show(io, args[2])
        for a in args[3:end]
            print(io, ", ")
            show(io, a)
        end
    end
    print(io, ']')
    return
end

function show_impl(io, ::Colon, args, kwargs)
    @assert length(kwargs) == 0
    if length(args) > 0
        show(io, args[1])
        for a in args[2:end]
            print(io, ":")
            show(io, a)
        end
    end
    return
end

# --- User interface

const ◻ = Hole()

function __init__()
    try
        REPL.REPLCompletions.latex_symbols["\\tofu"] = "◻"
    catch err
        @error "Incompatible REPL module; Disabling LaTeX command `\\tofu`" err
    end
end

end # module


text:
params = BasisMatrices.ChebParams(7, 1e-4, 10.0)

nod = BasisMatrices.nodes(params)
nod_end = BasisMatrices.nodes(params, 1) # extended nodes to endpoints
nod_l = BasisMatrices.nodes(params, 3) # Lobatto nodes

function manualeval(x, a, b, n)

    z = @. 2*(x-a)/(b-a)-1
    m = length(x)

    bas = Array{Float64}(undef, m, n)
    bas[:,1:2] = [ones(m, 1) z]

    for i in 3:n

        bas[:, i] .= 2.0.*z .* bas[:, i-1] .- bas[:, i-2]

    end

    return bas

end

mB = manualeval(nod, params.a[1], params.b[1], params.n[1])

#=
function manualderiv1(a,b,n)

    x = BasisMatrices.nodes(BasisMatrices.ChebParams(n-1,a,b))
    z = (2*(x-a)/(b-a)-1)''

    D = [zeros(length(z),1) ones(length(z),1) 4*z 12*z.^2-3 32*z.^3-16*z 80*z.^4-60*z.^2+5 192*z.^5-192*z.^3+36*z]

    return D

end

function manualint1(a,b,n)

    x = BasisMatrices.nodes(BasisMatrices.ChebParams(n+1,a,b))
    z = (2*(x-a)/(b-a)-1)''

    I = [z .5*z.^2 2/3*z.^3-z z.^4-3/2*z.^2 8/5*z.^5-8/3*z.^3+z 16/6*z.^6-5*z.^4+2.5*z.^2 32/7*z.^7-48/5*z.^5+6*z.^3-z]

    return I

end

D = manualderiv1(nod,params.a[1],params.b[1])
I = manualint1(nod,params.a[1],params.b[1])
=#

@testset "test Cheb" begin

    @testset "test nodes" begin
        @test (nod[1] == params.a[1],nod[end] == params.b[1]) == (false, false)
        @test ≈([nod_end[1],nod_end[end]], [params.a[1],params.b[1]], atol=1e-15)
        @test ≈([nod_l[1],nod_l[end]], [params.a[1],params.b[1]], atol=1e-15)
    end

    @testset "test derivative/integral" begin
        deriv, params_d = BasisMatrices.derivative_op(params, nod, 1)
        int, params_i = BasisMatrices.derivative_op(params, nod, -1)

        @test params_d.n  ==  params.n[1]-1
        @test params_i.n  ==  params.n[1]+1
        @test ≈([params_d.a, params_d.b], [params.a[1], params.b[1]], atol=1e-15)
        @test ≈([params_i.a, params_i.b], [params.a[1], params.b[1]], atol=1e-15)
    end

    @testset "test evalbase" begin
        B = @inferred BasisMatrices.evalbase(params, nod, 0, 0)
        @test ≈(B, mB, atol=1e-14)
    end

    @test BasisMatrices.family(ChebParams) == Cheb
    @test BasisMatrices.family_name(ChebParams) == "Cheb"
    @test !issparse(ChebParams)

    @test BasisMatrices.family(params) == Cheb
    @test BasisMatrices.family_name(params) == "Cheb"
    @test !issparse(params)

    @test_throws DimensionMismatch BasisMatrices.evalbasex!(rand(2, 7), rand(4), params, rand(4))
    @test_throws DimensionMismatch BasisMatrices.evalbasex!(rand(4, 6), rand(4), params, rand(4))

    @testset "non Float64 eltype" begin
        p = ChebParams(10, Float32(-3.0), Float32(1.0))
        @test Float32 == eltype(@inferred BasisMatrices.evalbase(p, Float32[0.25, 0.35]))
        @test Float32 == eltype(@inferred BasisMatrices.evalbase(p, Float32[0.25, 0.35], 1))
        @test Float64 == eltype(@inferred BasisMatrices.evalbase(p, [0.25, 0.35]))
        @test Float32 == eltype(@inferred nodes(p))
        @test Float32 == eltype(@inferred nodes(p, 1))
        @test Float32 == eltype(@inferred nodes(p, 2))
    end

end  # testset


text:
using StateMachines, Simulate
using Test

println("...... testing StateMachines.jl ......")
@testset "types.jl" begin
    include("test_types.jl")
end

@testset "box.jl" begin
    include("test_box.jl")
end

@testset "block.jl" begin
    include("test_block.jl")
end

@testset "statemachine.jl" begin
    include("test_statemachine.jl")
end

@testset "state.jl" begin
    include("test_state.jl")
end

@testset "event.jl" begin
    include("test_event.jl")
end

@testset "transition.jl" begin
    include("test_transition.jl")
end

@testset "action.jl" begin
    include("test_action.jl")
end


text:
export TR_Cub
function TR_Cub(h :: AbstractNLPModel,
                nlpstop :: AbstractStopping;
                t₀ = h.meta.x0, verbose = true, eps1 = 0.2, eps2 = 0.8,
                red = 0.5, aug = 2.0, Δ = 1.0)
    Δ  = [Δ]
    t  = t₀; iter = 0;
    fₖ = obj(h, t); gₖ = grad(h, t);

    s = t .- t₀; y = gₖ # just for first iteration

    dN = -gₖ # for the first iteration
    A = [0.5]; B = [0.0];

    q(d) = fₖ .+ gₖ .* d .+ A .* d .^ 2 .+ B .* d .^ 3

    OK  = update_and_start!(nlpstop, x = t, fx = fₖ, gx = gₖ, g0 = copy(gₖ))

    verbose && @printf(" iter  t         gₖ          Δ        pred         ared\n")
    verbose && @printf(" %4d %7.2e  %7.2e  %7.2e \n", iter, t[1], gₖ[1], Δ[1])

    while !OK

        if true in (q(Δ) .< q(-Δ))
            d = Δ
        else
            d = -Δ
        end

        cub(t) = fₖ .+ gₖ .* t .+ A .* t .^ 2 .+ B .* t .^ 3

        dR = roots(Poly([gₖ[1], 2.0 .* A[1], 3.0 .* B[1]]))

        if isreal(dR)
          dR = real(dR)
          dN2 = dR[1]
          if length(dR) > 1
            if true in ((cub(dR[2]) .< cub(dR[1])))
              dN2 = dR[2]
            end
          end
        else
          dN2 = -gₖ .* s ./ y
        end

        if true in ((abs.(dN2) .< Δ)) && true in ((q(d) .> q(dN2)))
            d = dN2
        end

        if true in (d .== 0.0)
            return OK, nlpstop  # still useful?
        end

        ftestTR = obj(h, t .+ d)
        gtestTR = grad(h, t .+ d)

        pred = gₖ .* d .+ A .* d .^ 2 .+ B .* d .^ 3

        if true in (pred .> - 1e-10)
            ared = (gₖ .+ gtestTR) .* d ./ 2
        else
            ared = ftestTR .- fₖ
        end

        ratio = ared ./ pred
        if true in ((ratio .< eps1)) && true in ((abs.(d) .== Δ))
            Δ = red * Δ
        else
          #two points memory
          tpred = t; gkm1 = gₖ; fkm1 = fₖ
          t = t .+ d; gₖ = gtestTR; fₖ = ftestTR

          #cubic two points interpolation
          s = t .- tpred; y = gₖ .- gkm1

          α = -s
          z = gₖ .+ gkm1 .+ 3. .* (fₖ .- fkm1) ./ α
          discr = z .^ 2 .- gₖ .* gkm1
          denom = gₖ .+ gkm1 .+ 2 .* z
          B = 1 ./ 3 .* (gₖ .+ gkm1 .+ 2. .* z) ./ (α .* α)
          A = -(gₖ .+ z) ./ α

          if true in (ratio .> eps2)
            Δ = aug * Δ
          end
        end

        OK  = update_and_stop!(nlpstop, x = t, fx = fₖ, gx = gₖ);

        iter += 1
        verbose && @printf(" %4d %7.2e  %7.2e  %7.2e %7.2e %7.2e\n", iter, t[1], gₖ[1], Δ[1], pred[1], ared[1])
    end

    optimal = nlpstop.meta.optimal
    return optimal, nlpstop
end


text:
module PythonOT

using PyCall: PyCall

export emd, emd2, sinkhorn, sinkhorn2, barycenter, sinkhorn_unbalanced, sinkhorn_unbalanced2

const pot = PyCall.PyNULL()

include("lib.jl")

function __init__()
    return copy!(pot, PyCall.pyimport_conda("ot", "pot", "conda-forge"))
end

end


text:
# Utilities

## repeat each element in a vector

function repeach{T}(x::AbstractVector{T}, n::Integer)
    k = length(x)
    r = Array(T, k * n)
    p = 0
    @inbounds for i = 1:k
        xi = x[i]
        for j = 1:n
            r[p += 1] = xi
        end
    end
    return r
end

function repeach{T}(x::AbstractVector{T}, ns::IntegerVector)
    k = length(x)
    length(ns) == k || throw(DimensionMismatch("length(ns) should be equal to k."))
    r = Array(T, sum(ns))
    p = 0
    @inbounds for i = 1:k
        xi = x[i]
        ni = ns[i]
        for j = 1:ni
            r[p += 1] = xi
        end
    end
    return r
end

## repeat each column in a matrix

function repeachcol{T}(x::DenseArray{T,2}, n::Integer)
    m = size(x, 1)
    k = size(x, 2)
    r = Array(T, m, k * n)
    p = 0
    @inbounds for i = 1:k
        xi = view(x, :, i)
        for j = 1:n
            r[:, p += 1] = xi
        end
    end
    return r
end

function repeachcol{T}(x::DenseArray{T,2}, ns::IntegerVector)
    m = size(x, 1)
    k = size(x, 2)
    r = zeros(T, m, sum(ns))
    p = 0
    @inbounds for i = 1:k
        xi = view(x, :, i)
        ni = ns[i]
        for j = 1:ni
            r[:, p += 1] = xi
        end
    end
    return r
end

## repeat each row in a matrix

function repeachrow{T}(x::DenseArray{T,2}, n::Integer)
    k = size(x, 1)
    m = size(x, 2)
    r = Array(T, k * n, m)
    p = 0
    @inbounds for icol = 1:m
        p = 0
        for i = 1:k
            xi = x[i, icol]
            for j = 1:n
                r[p += 1, icol] = xi
            end
        end
    end
    return r
end

function repeachrow{T}(x::DenseArray{T,2}, ns::IntegerVector)
    k = size(x, 1)
    m = size(x, 2)
    r = Array(T, sum(ns), m)
    @inbounds for icol = 1:m
        p = 0
        for i = 1:k
            xi = x[i, icol]
            ni = ns[i]
            for j = 1:ni
                r[p += 1, icol] = xi
            end
        end
    end
    return r
end

## return the unique values of A and a vector of vectors of indexes to reconstruct
## the original array

function unique_inverse(A)
    out = Array(eltype(A),0)
    out_idx = Array(Vector{Int}, 0)
    seen = Dict{eltype(A), Int}()
    for (idx, x) in enumerate(A)
        if !in(x, keys(seen))
            seen[x] = length(seen) + 1
            push!(out, x)
            push!(out_idx, Int[])
        end
        push!(out_idx[seen[x]], idx)
    end
    out, out_idx
end


text:
const project_root_dir = realpath(joinpath(@__DIR__, ".."))
unique!(push!(LOAD_PATH, realpath(joinpath(project_root_dir, "experiments/utils/MonteCarloStudy"))))
unique!(push!(LOAD_PATH, realpath(joinpath(project_root_dir, "experiments/utils"))))
unique!(push!(LOAD_PATH, realpath(joinpath(project_root_dir, "test/utils"))))

import Distributed
Distributed.@everywhere begin
    import Pkg
    Pkg.activate($project_root_dir)
    union!(LOAD_PATH, $LOAD_PATH)

    import MonteCarloStudy
    import CollisionAvoidanceGame
    import TestDynamics
    using PartiallyObservedInverseGames.ForwardGame: IBRGameSolver, KKTGameSolver
    using PartiallyObservedInverseGames.InverseGames:
        InverseKKTConstraintSolver, InverseKKTResidualSolver
end

import VegaLite
import Random
import Distributor
import CostHeatmapVisualizer
using PartiallyObservedInverseGames.TrajectoryVisualization: visualize_trajectory

# Utils
include("utils/misc.jl")
include("utils/simple_caching.jl")
load_cache_if_not_defined!("highway")

#==================================== Forward Game Formulation =====================================#

T = 20
Δt = 1.0
rng = Random.MersenneTwister(1)

control_system = TestDynamics.ProductSystem([
    TestDynamics.Unicycle(Δt),
    TestDynamics.Unicycle(Δt),
    TestDynamics.Unicycle(Δt),
    TestDynamics.Unicycle(Δt),
    TestDynamics.Unicycle(Δt),
])

player_configurations = [
    # Vehicle on the right lane wishing to merge left to go faster
    (;
        initial_lane = 1.0,
        initial_progress = 0,
        initial_speed = 0.2,
        target_speed = 0.3,
        speed_cost = 1.0,
        target_lane = 0.0,
        prox_cost = 0.3,
    ),
    # Fast vehicle from the back that would like to maintain its speed.
    (;
        initial_lane = 0,
        initial_progress = -3.0,
        initial_speed = 0.4,
        target_speed = 0.4,
        target_lane = 0.0,
        speed_cost = 1.0,
        prox_cost = 0.3,
    ),
    # Slow truck on the right lane
    (;
        initial_lane = 1.0,
        initial_progress = 2,
        initial_speed = 0.10,
        target_speed = 0.10,
        speed_cost = 1.0,
        target_lane = 1.0,
        prox_cost = 0.05,
    ),
    # Slow truck on the right lane
    (;
        initial_lane = 1.0,
        initial_progress = 4,
        initial_speed = 0.10,
        target_speed = 0.10,
        speed_cost = 1.0,
        target_lane = 1.0,
        prox_cost = 0.05,
    ),
    # Fast vehicle on the left lane wishing to merge back on the right lane and slow down
    (;
        initial_lane = 0.0,
        initial_progress = 6,
        initial_speed = 0.3,
        target_speed = 0.2,
        speed_cost = 1.0,
        target_lane = 1.0,
        prox_cost = 0.3,
    ),
]

x0 = mapreduce(vcat, player_configurations) do player_config
    [
        player_config.initial_lane,
        player_config.initial_progress,
        player_config.initial_speed,
        deg2rad(90),
    ]
end

position_indices = mapreduce(vcat, eachindex(control_system.subsystems)) do subsystem_idx
    TestDynamics.state_indices(control_system, subsystem_idx)[1:2]
end

partial_state_indices = mapreduce(vcat, eachindex(control_system.subsystems)) do subsystem_idx
    TestDynamics.state_indices(control_system, subsystem_idx)[[1, 2, 4]]
end

player_cost_models_gt = map(Iterators.countfrom(1), player_configurations) do ii, player_config
    cost_model_p1 = CollisionAvoidanceGame.generate_player_cost_model(;
        player_idx = ii,
        control_system,
        T,
        goal_position = nothing,
        weights = (;
            state_proximity = player_config.prox_cost,
            state_velocity = player_config.speed_cost,
            control_Δv = 1,
            control_Δθ = 1,
        ),
        y_lane_center = player_config.target_lane,
        target_speed = player_config.target_speed,
    )
end

#===================================== Additional Visualization ====================================#

function visualize_highway(x; subsampling = 1, kwargs...)
    viz = let
        max_size = 500
        y_position_domain = [-4, 14]
        x_position_domain = [-1, 2]
        x_range = only(diff(extrema(x_position_domain) |> collect))
        y_range = only(diff(extrema(y_position_domain) |> collect))
        max_range = max(x_range, y_range)
        canvas = VegaLite.@vlplot(
            width = max_size * x_range / max_range,
            height = max_size * y_range / max_range
        )

        subsampled_taj = x[:, 1:subsampling:end]

        visualize_trajectory(
            control_system,
            subsampled_taj;
            x_position_domain,
            y_position_domain,
            canvas,
            kwargs...,
        )
    end
end

#======================================== Monte Carlo Study ========================================#

## Dataset Generation
n_observation_sequences_per_noise_level = 40

@run_cached forward_solution_gt, dataset = MonteCarloStudy.generate_dataset(;
    solve_args = (; solver = IBRGameSolver(), control_system, player_cost_models_gt, x0, T),
    noise_levels = unique([0:0.001:0.01; 0.01:0.005:0.03; 0.03:0.01:0.1]),
    n_observation_sequences_per_noise_level,
)

## Estimation
estimator_setup = (;
    dataset,
    control_system,
    player_cost_models = player_cost_models_gt,
    solver_attributes = (; print_level = 1),
)
estimator_setup_partial =
    merge(estimator_setup, (; expected_observation = x -> x[partial_state_indices, :]))
@run_cached estimates_conKKT =
    MonteCarloStudy.estimate(InverseKKTConstraintSolver(); estimator_setup...)
@run_cached estimates_conKKT_partial =
    MonteCarloStudy.estimate(InverseKKTConstraintSolver(); estimator_setup_partial...)
@run_cached estimates_resKKT =
    MonteCarloStudy.estimate(InverseKKTResidualSolver(); estimator_setup...)
@run_cached estimates_resKKT_partial =
    MonteCarloStudy.estimate(InverseKKTResidualSolver(); estimator_setup_partial...)

## Forward Solution Augmentation
augmentor_kwargs = (;
    solver = KKTGameSolver(),
    control_system,
    player_cost_models_gt,
    x0,
    T,
    match_equilibrium = (; forward_solution_gt.x),
    init = (; forward_solution_gt.x, forward_solution_gt.u),
    solver_attributes = (; print_level = 1),
)
@run_cached augmented_estimates_resKKT =
    MonteCarloStudy.augment_with_forward_solution(estimates_resKKT; augmentor_kwargs...)
@run_cached augmented_estimates_resKKT_partial =
    MonteCarloStudy.augment_with_forward_solution(estimates_resKKT_partial; augmentor_kwargs...)
estimates = [
    estimates_conKKT
    estimates_conKKT_partial
    augmented_estimates_resKKT
    augmented_estimates_resKKT_partial
]

## Error Ststistics Computation
demo_gt = merge((; player_cost_models_gt), forward_solution_gt)
@save_json errstats = map(estimates) do estimate
    MonteCarloStudy.estimator_statistics(estimate; dataset, demo_gt, position_indices)
end

## Visualization
frame = [-1.5n_observation_sequences_per_noise_level, 0]
@saveviz parameter_error_viz = errstats |> MonteCarloStudy.visualize_paramerr(; frame)
@saveviz position_error_viz = errstats |> MonteCarloStudy.visualize_poserr(; frame)
@saveviz highway_frontfig_cost1 = CostHeatmapVisualizer.cost_viz(
    1,
    player_configurations[1];
    x_sequence = forward_solution_gt.x[:, 1:1],
    control_system,
)
@saveviz highway_frontfig_cost5 = CostHeatmapVisualizer.cost_viz(
    5,
    player_configurations[5];
    x_sequence = forward_solution_gt.x[:, 1:1],
    control_system,
)
@saveviz highway_frontfig_gt_viz = visualize_highway(forward_solution_gt.x)
@saveviz highway_frontfig_observation_viz = visualize_highway(dataset[end].x; draw_line = false)


text:
using Primes: isprime
import Base.divgcd

const primetable = GrowingList([2, 3]; sizehint = 1024)
const factortable = GrowingList([UInt8[], UInt8[1], UInt8[0,1]]; sizehint = 4096)
const factorialtable = GrowingList([UInt32[], UInt32[1], UInt32[1,1]]; sizehint = 4096)
const bigprimetable = GrowingList([GrowingList([big(2)]; sizehint = 2048),
                                    GrowingList([big(3)]; sizehint = 1024)];
                                    sizehint = 1024)
const bigone = big(1)

# Make a prime iterator
struct PrimeIterator
end
primes() = PrimeIterator()

Base.IteratorSize(::Type{PrimeIterator}) = Base.IsInfinite()
Base.IteratorEltype(::Type{PrimeIterator}) = Base.HasEltype()
Base.eltype(::PrimeIterator) = Int

# Get the `n`th prime; store all primes up to the `n`th if not yet available
function prime(n::Int)
    k = min(length(primetable), length(bigprimetable))
    while k < n
        @inbounds p = primetable[k]
        p = p + 2
        while !isprime(p)
            p += 2
        end
        k += 1
        # these lines do not get but set new elements; provided no other task did so earlier
        get!(primetable, k, p)
        bp = big(p)
        bpf = GrowingList{BigInt}((big(p),); sizehint = 4)
        get!(bigprimetable, k, bpf)
        k = min(length(primetable), length(bigprimetable))
         # other threads might have inserted additional entries,
         # make sure they are finished with both primetable and bigprimetable
    end
    @inbounds return primetable[n]
end

Base.iterate(::PrimeIterator, n = 1) = prime(n), n+1

# get primes and their powers as `BigInt`, also cache all results
function bigprime(n::Integer, e::Integer=1)
    e == 0 && return bigone
    p = prime(n) # triggers computation of prime(n) if necessary
    @inbounds powerlist = bigprimetable[n]
    l = length(powerlist)
    @inbounds while l < e
        # compute next prime power as approximate square of existing results
        l += 1
        k = l>>1
        newpower = powerlist[k]*powerlist[l-k]
        get!(powerlist, l, newpower)
        l = length(powerlist) # other threads might have inserted more powers
    end
    @inbounds return powerlist[e]
end

# A custom `Integer` subtype to store an integer as its prime factorization
# mutable to allow in place update of sign
mutable struct PrimeFactorization{U<:Unsigned} <: Integer
    powers::Vector{U}
    sign::Int8
    PrimeFactorization{U}(powers::Vector, sign = one(Int8)) where {U<:Unsigned} =
        new{U}(convert(Vector{U}, powers), sign)
end
# convenience constructor: normalizes powers to have last entry nonzero
PrimeFactorization(powers::Vector{U}, sign = one(Int8)) where {U<:Unsigned} =
    PrimeFactorization{U}(_normalize_powers!(powers), sign)

function _normalize_powers!(v::Vector{<:Integer})
    i = findlast(!iszero, v)
    l = ifelse(i === nothing, 0, i)
    l < length(v) && resize!(v, l)
    return v
end

# define our own factor function, returning an instance of PrimeFactorization
function primefactor(n::Integer)
    iszero(n) && return PrimeFactorization{UInt8}(UInt8[], zero(Int8))
    sn = n < 0 ? -one(Int8) : one(Int8)
    n = abs(n)
    m = length(factortable)
    while m < abs(n)
        m += 1
        powers = UInt8[] # should be sufficient for all integers up to 2^255
        a = m
        for p in primes()
            f = 0
            anext, r = divrem(a, p)
            while r == 0
                f += 1
                a = anext
                anext, r = divrem(a, p)
            end
            push!(powers, f)
            a == 1 && break
        end
        get!(factortable, m, powers)
        m = length(factortable) # other threads may have inserted other entries
    end
    @inbounds return PrimeFactorization{UInt8}(factortable[n], sn)
end

function primefactorial(n::Integer)
    n < 0 && throw(DomainError(n,"primefactorial only works for non-negative numbers"))
    n <= 1 && return PrimeFactorization{UInt32}(UInt32[], one(Int8))
    m = length(factorialtable)
    @inbounds while m < n
        prevfactorial = factorialtable[m]
        m += 1
        f = primefactor(m).powers
        if length(f) > length(prevfactorial) # can at most be 1 larger
            powers = similar(prevfactorial, length(f))
            powers[1:end-1] = prevfactorial
            powers[end] = 0
        else
            powers = copy(prevfactorial)
        end
        for k = 1:length(f)
            powers[k] += f[k]
        end
        get!(factorialtable, m, powers)
        m = length(factorialtable)  # other threads may have inserted other entries
    end
    @inbounds return PrimeFactorization{UInt32}(factorialtable[n])
end

# Methods for PrimeFactorization:
Base.copy(a::PrimeFactorization) = PrimeFactorization(copy(a.powers), a.sign)
function Base.copy!(c::PrimeFactorization, a::PrimeFactorization)
    c.sign = a.sign
    copy!(c.powers, a.powers)
    return c
end

Base.one(::Type{PrimeFactorization{U}}) where {U<:Unsigned} =
    PrimeFactorization{U}(Vector{U}(), one(Int8))
Base.zero(::Type{PrimeFactorization{U}}) where {U<:Unsigned} =
    PrimeFactorization{U}(Vector{U}(), zero(Int8))

one!(c::PrimeFactorization) = (c.sign = one(Int8); empty!(c.powers); return c)
zero!(c::PrimeFactorization) = (c.sign = zero(Int8); empty!(c.powers); return c)

Base.promote_rule(P::Type{<:PrimeFactorization},::Type{<:Integer}) = P
Base.promote_rule(::Type{<:PrimeFactorization},::Type{BigInt}) = BigInt
Base.promote_rule(::Type{PrimeFactorization{U1}},
                    ::Type{PrimeFactorization{U2}}) where {U1<:Unsigned, U2<:Unsigned} = PrimeFactorization{promote_type(U1, U2)}

Base.convert(P::Type{<:PrimeFactorization}, n::Integer) = convert(P, primefactor(n))
function _convert!(x::BigInt, a::PrimeFactorization)
    MPZ.set!(x, bigone)
    for (n, e) in enumerate(a.powers)
        if !iszero(e)
            MPZ.mul!(x, bigprime(n, e))
        end
    end
    return a.sign < 0 ? MPZ.neg!(x) : x
end
Base.convert(::Type{BigInt}, a::PrimeFactorization) = _convert!(one(BigInt), a)
Base.convert(::Type{PrimeFactorization{U}}, a::PrimeFactorization{U}) where {U<:Unsigned} =
    a
Base.convert(::Type{PrimeFactorization{U}}, a::PrimeFactorization) where {U<:Unsigned} =
    PrimeFactorization{U}(convert(Vector{U}, a.powers), a.sign)

Base.:(==)(a::PrimeFactorization, b::PrimeFactorization) =
    a.powers == b.powers && a.sign == b.sign
function Base.:<(a::PrimeFactorization, b::PrimeFactorization)
    if a.sign != b.sign
        return a.sign < b.sign
    elseif a.sign < 0
        return <(-b, -a)
    else
        ag, bg = divgcd(a, b)
        ag == bg && return false
        if length(ag.powers) <= length(bg.powers)
                all(k->ag.powers[k]<bg.powers[k], 1:length(ag.powers))
            return true
        else
            return convert(BigInt, ag) < convert(BigInt, bg)
        end
    end
end

# Methods for PrimeFactorization: only fast multiplication, and lcm and gcd.
# Addition and subtraction will require conversion to BigInt
Base.sign(a::PrimeFactorization) = a.sign

neg!(a::PrimeFactorization) = (a.sign = -a.sign; return a)

function mul!(c::PrimeFactorization, a::PrimeFactorization, b::PrimeFactorization)
    if a.sign == 0 || b.sign == 0
        zero!(c)
    else
        c.sign = a.sign * b.sign
        la = length(a.powers)
        lb = length(b.powers)
        lc = max(la, lb)
        lc === length(c.powers) || resize!(c.powers, lc)
        @inbounds for k = 1:min(la,lb)
            c.powers[k] = +(a.powers[k], b.powers[k])
        end
        if c !== a
            @inbounds for k = lb+1:la
                c.powers[k] = a.powers[k]
            end
        end
        @inbounds for k = la+1:lb
            c.powers[k] = b.powers[k]
        end
    end
    return c
end
# unlike div, this one errors if the a is not divisible by b
function divexact!(c::PrimeFactorization, a::PrimeFactorization, b::PrimeFactorization)
    if iszero(a.sign)
        zero!(c)
    elseif iszero(b.sign)
        throw(DivideError())
    else
        c.sign = a.sign * b.sign
        la = length(a.powers)
        lb = length(b.powers)
        if lb > la
            throw(DivideError())
        end
        lc = la
        if lb == lc
            while lc > 0 && a.powers[lc] == b.powers[lc]
                lc -= 1
            end
        end
        lc == length(c.powers) || resize!(c.powers, lc)
        @inbounds for k = 1:min(lb, lc)
            if b.powers[k] > a.powers[k]
                throw(DivideError())
            end
            c.powers[k] = a.powers[k] - b.powers[k]
        end
        if c !== a
            @inbounds for k = lb+1:lc
                c.powers[k] = a.powers[k]
            end
        end
    end
    return c
end
function gcd!(c::PrimeFactorization, a::PrimeFactorization, b::PrimeFactorization)
    if a.sign == 0
        copy!(c.powers, b.powers)
    elseif b.sign ==0
        copy!(c.powers, a.powers)
    else
        c.sign = one(Int8)
        la = length(a.powers)
        lb = length(b.powers)
        lc = min(la, lb)
        lc === length(c.powers) || resize!(c.powers, lc)
        @inbounds for k = 1:lc
            c.powers[k] = min(a.powers[k], b.powers[k])
        end
    end
    c.sign = one(Int8)
    return c
end
function lcm!(c::PrimeFactorization, a::PrimeFactorization, b::PrimeFactorization)
    if a.sign == 0 || b.sign == 0
        return zero!(c)
    else
        c.sign = one(Int8)
        la = length(a.powers)
        lb = length(b.powers)
        lc = max(la, lb)
        lc === length(c.powers) || resize!(c.powers, lc)
        @inbounds for k = 1:min(la,lb)
            c.powers[k] = max(a.powers[k], b.powers[k])
        end
        if c !== a
            @inbounds for k = lb+1:la
                c.powers[k] = a.powers[k]
            end
        end
        @inbounds for k = la+1:lb
            c.powers[k] = b.powers[k]
        end
    end
    c.sign = one(Int8)
    return c
end
function divgcd!(a::PrimeFactorization, b::PrimeFactorization)
    af, bf = a.powers, b.powers
    for k = 1:min(length(af), length(bf))
        gk = min(af[k], bf[k])
        af[k] -= gk
        bf[k] -= gk
    end
    _normalize_powers!(a.powers)
    _normalize_powers!(b.powers)
    return a, b
end

mul!(a::PrimeFactorization, b::PrimeFactorization) = mul!(a, a, b)
divexact!(a::PrimeFactorization, b::PrimeFactorization) = divexact!(a, a, b)
gcd!(a::PrimeFactorization, b::PrimeFactorization) = gcd!(a, a, b)
lcm!(a::PrimeFactorization, b::PrimeFactorization) = lcm!(a, a, b)

Base.:-(a::PrimeFactorization) = neg!(copy(a))
function Base.:*(a::PrimeFactorization, b::PrimeFactorization)
    P = promote_type(typeof(a), typeof(b))
    if length(a.powers) >= length(b.powers)
        return typeof(a) == P ? mul!(copy(a), b) : mul!(convert(P, a), b)
    else
        return typeof(b) == P ? mul!(copy(b), a) : mul!(convert(P, b), a)
    end
end
function Base.lcm(a::PrimeFactorization, b::PrimeFactorization)
    P = promote_type(typeof(a), typeof(b))
    if length(a.powers) >= length(b.powers)
        return typeof(a) == P ? lcm!(copy(a), b) : lcm!(convert(P, a), b)
    else
        return typeof(b) == P ? lcm!(copy(b), a) : lcm!(convert(P, b), a)
    end
end
function Base.gcd(a::PrimeFactorization, b::PrimeFactorization)
    P = promote_type(typeof(a), typeof(b))
    if length(a.powers) <= length(b.powers)
        return typeof(a) == P ? lcm!(copy(a), b) : lcm!(convert(P, a), b)
    else
        return typeof(b) == P ? lcm!(copy(b), a) : lcm!(convert(P, b), a)
    end
end
Base.divgcd(a::PrimeFactorization, b::PrimeFactorization) = divgcd!(copy(a), copy(b))
# no promotion necessary, should be smaller than a
divexact(a::PrimeFactorization, b::PrimeFactorization) = divexact!(copy(a), b)

# split `a::PrimeFactorization` into a square `s` and a remainder `r`, such that
# `a = s^2 * r` and the powers in the prime factorization of `r` are zero or one
function splitsquare(a::PrimeFactorization)
    r = PrimeFactorization(map(p->convert(UInt8, isodd(p)), a.powers), a.sign)
    s = PrimeFactorization(map(p->(p>>1), a.powers))
    return s, r
end

# given a list of numerators and denominators, compute the common denominator and
# the rescaled numerator after putting all fractions at the same common denominator
function commondenominator!(nums::Vector{P}, dens::Vector{P}) where {P<:PrimeFactorization}
    isempty(nums) && return one(P)
    # accumulate lcm of denominator
    den = copy(dens[1])
    for i = 2:length(dens)
        lcm!(den, dens[i])
    end
    # rescale numerators
    for i = 1:length(nums)
        divexact!(mul!(nums[i], den), dens[i])
    end
    return den
end

# auxiliary function to compute sums of a list of PrimeFactorizations as quickly as possible
function sumlist!(list::Vector{<:PrimeFactorization}, ind = 1:length(list))
    # first compute gcd to take out common factors
    g = copy(list[ind[1]])
    for k in ind
        gcd!(g, list[k])
    end
    for k in ind
        divexact!(list[k], g)
    end
    L = length(ind)
    i = big(1)
    if L > 32
        l = L >> 1
        s = sumlist!(list, first(ind).+(0:l-1))
        s = MPZ.add!(s, sumlist!(list, first(ind).+(l:L-1)))
    else # do sum, add to s
        s = big(0)
        for k in ind
            MPZ.add!(s, _convert!(i, list[k]))
        end
    end
    return MPZ.mul!(s, _convert!(i, g))
end


text:
module Changepoints

using Distributions
using Distributions: Sampleable
using Base.Meta
import Base.rand

export PELT, @PELT, BS, @BS, CROPS, @segment_cost, NormalVarSegment, NormalMeanSegment, NormalMeanVarSegment, PoissonSegment, BetaSegment, ExponentialSegment, ChangepointSampler, GammaShapeSegment, GammaRateSegment, NonparametricSegment, OLSSegment, @changepoint_sampler

include("segment_costs.jl")
include("PELT.jl")
include("CROPS.jl")
include("BS.jl")
include("sim_type.jl")
include("macros.jl")

# This approach to loading supported plotting packages is taken directly from the "KernelDensity" package
macro glue(pkg)
    path = joinpath(dirname(@__FILE__),"glue",string(pkg,".jl"))
    init = Symbol(string(pkg,"_init"))
    quote
        $(esc(init))() = Base.include($path)
        isdefined(Main,$(QuoteNode(pkg))) && $(esc(init))()
    end
end

@glue Winston
@glue Gadfly

end # module


text:
using LinearAlgebra, DiffEqOperators, Random, Test, BandedMatrices, SparseArrays

# 4th derivative, 4th order
function fourth_deriv_approx_stencil(N)
    A = zeros(N,N+2)
    A[1,1:8] = [3.5 -56/3 42.5 -54.0 251/6 -20.0 5.5 -2/3]
    A[2,1:8] = [2/3 -11/6 0.0 31/6 -22/3 4.5 -4/3 1/6]

    A[N-1,N-5:end] = reverse([2/3 -11/6 0.0 31/6 -22/3 4.5 -4/3 1/6], dims=2)
    A[N,N-5:end] = reverse([3.5 -56/3 42.5 -54.0 251/6 -20.0 5.5 -2/3], dims=2)

    for i in 3:N-2
        A[i,i-2:i+4] = [-1/6 2.0 -13/2 28/3 -13/2 2.0 -1/6]
    end
    return A
end

function second_deriv_fourth_approx_stencil(N)
    A = zeros(N,N+2)
    A[1,1:6] = [5/6 -15/12 -1/3 7/6 -6/12 5/60]
    A[N,N-3:end] = [1/12 -6/12 14/12 -4/12 -15/12 10/12]
    for i in 2:N-1
        A[i,i-1:i+3] = [-1/12 4/3 -5/2 4/3 -1/12]
    end
    return A
end

function second_derivative_stencil(N)
  A = zeros(N,N+2)
  for i in 1:N, j in 1:N+2
      (j-i==0 || j-i==2) && (A[i,j]=1)
      j-i==1 && (A[i,j]=-2)
  end
  A
end

@testset "Test Constructor, Multiplication, and Concretization" begin
    # Generate random parameters
    al = rand()
    bl = rand()
    cl = rand()

    ar = rand()
    br = rand()
    cr = rand()
    dx = rand()

    Q = RobinBC((al, bl, cl), (ar, br, cr), dx)
    N = 20
    L = CenteredDifference(4,4, 1.0, N)
    L2 = CenteredDifference(2,4, 1.0, N)

    function coeff_func(du,u,p,t)
        du .= u
    end

    cL = coeff_func*L
    coeffs = rand(N)
    DiffEqOperators.update_coefficients!(cL,coeffs,nothing,0.0)
    cA = cL*Q
    DiffEqOperators.update_coefficients!(cA,coeffs,nothing,0.0)
    @test cL.coefficients == coeffs == cA.L.coefficients

    A = L*Q
    c2A = coeff_func*A
    DiffEqOperators.update_coefficients!(c2A,coeffs,nothing,0.0)
    @test coeffs == c2A.L.coefficients

    # Test GhostDerivativeOperator constructor by *
    u = rand(N)
    A = L*Q
    # Test for consistency of GhostDerivativeOperator*u with L*(Q*u)
    @test A*u ≈ L*(Q*u)

    # Test for consistency of GhostDerivativeOperator*M with L*(Q*M)
    M = rand(N,10)
    LQM = zeros(N,10)
    for i in 1:10
        mul!(view(LQM,:,i), L, Q*M[:,i])
    end
    ghost_LQM = A*M
    @test ghost_LQM ≈ LQM

    u = rand(22)
    @test (L + L2) * u ≈ convert(AbstractMatrix,L + L2) * u ≈ (BandedMatrix(L) + BandedMatrix(L2)) * u

    # Test concretization
    @test Array(A)[1] ≈ (Array(L)*Array(Q,N)[1], Array(L)*Array(Q,N)[2])[1]
    @test Array(A)[2] ≈ (Array(L)*Array(Q,N)[1], Array(L)*Array(Q,N)[2])[2]
    @test SparseMatrixCSC(A)[1] ≈ (SparseMatrixCSC(L)*SparseMatrixCSC(Q,N)[1], SparseMatrixCSC(L)*SparseMatrixCSC(Q,N)[2])[1]
    @test SparseMatrixCSC(A)[2] ≈ (SparseMatrixCSC(L)*SparseMatrixCSC(Q,N)[1], SparseMatrixCSC(L)*SparseMatrixCSC(Q,N)[2])[2]
    @test sparse(A)[1] ≈ (sparse(L)*sparse(Q,N)[1], sparse(L)*sparse(Q,N)[2])[1]
    @test sparse(A)[2] ≈ (sparse(L)*sparse(Q,N)[1], sparse(L)*sparse(Q,N)[2])[2]
    # BandedMatrix not implemeted for boundary operator
    @test_broken BandedMatrix(A)[1] ≈ (BandedMatrix(L)*BandedMatrix(Q,N)[1], BandedMatrix(L)*BandedMatrix(Q,N)[2])[1]
    @test_broken BandedMatrix(A)[2] ≈ (BandedMatrix(L)*BandedMatrix(Q,N)[1], BandedMatrix(L)*BandedMatrix(Q,N)[2])[2]

    # Test that concretization works with multiplication
    u = rand(20)
    @test Array(A)[1]*u + Array(A)[2] ≈ L*(Q*u) ≈ A*u
    @test sparse(A)[1]*u + sparse(A)[2] ≈ L*(Q*u) ≈ A*u
end

@testset "Test Left Division L2 (second order)" begin
    # Test \ Homogenous BC
    # f(x) = -x^2 + x
    # f''(x) = -2.0
    # f'(0) = f'(1) = 0
    f(x) = -x^2 + x
    f2(x) = -2.0

    dx = 0.01
    x = dx:dx:1.0-dx
    N = length(x)

    L = CenteredDifference(2, 2, dx, N)
    Q = RobinBC((1.0, 0.0, 0.0), (1.0, 0.0, 0.0), dx)
    A = L*Q

    analytic_L = second_derivative_stencil(N) ./ dx^2
    analytic_QL = [transpose(zeros(N)); Diagonal(ones(N)); transpose(zeros(N))]
    analytic_AL = analytic_L*analytic_QL

    # No affine component to the this system
    analytic_f = analytic_AL \ f2.(x)
    ghost_f = A \ f2.(x)

    # Check that A\f2.(x) is consistent with analytic_AL \ f2.(x)
    @test analytic_f ≈ ghost_f

    # Additionally test that A\f2.(x) ≈ f.(x)
    @test f.(x) ≈ ghost_f ≈ analytic_f

    # Check ldiv!
    f_temp = zeros(N)
    ldiv!(f_temp, A, f2.(x))
    @test f_temp ≈ ghost_f ≈ analytic_f

    # Check that left division with matrices works
    ghost_fM = A \ [f2.(x) f2.(x)]
    analytic_fM = analytic_AL \ [f2.(x) f2.(x)]
    @test ghost_fM ≈ analytic_fM

    fM_temp = zeros(N,2)
    ldiv!(fM_temp, A, [f2.(x) f2.(x)])
    @test fM_temp ≈ analytic_fM

    # Test \ Inhomogenous BC
    # f(x) = -x^2 + x + 4.0
    # f''(x) = -2.0
    # f'(0) = f'(1) = 0
    f(x) = -x^2 + x + 4.0
    f2(x) = -2.0

    dx = 0.01
    x = dx:dx:1.0-dx
    N = length(x)

    L = CenteredDifference(2, 2, dx, N)
    Q = RobinBC((1.0, 0.0, 4.0), (1.0, 0.0, 4.0), dx)
    A = L*Q

    analytic_L = second_derivative_stencil(N) ./ dx^2
    analytic_QL = [transpose(zeros(N)); Diagonal(ones(N)); transpose(zeros(N))]
    analytic_Qb = [4.0; zeros(N); 4.0]
    analytic_AL = analytic_L*analytic_QL
    analytic_Ab = analytic_L*analytic_Qb

    analytic_f = analytic_AL \ (f2.(x) - analytic_Ab)
    ghost_f = A \ f2.(x)

    # Check that A\f2.(x) is consistent with analytic_AL \ f2.(x)
    @test analytic_f ≈ ghost_f

    # Additionally test that A\f2.(x) ≈ f.(x)
    @test f.(x) ≈ ghost_f ≈ analytic_f

    # Check ldiv!
    f_temp = zeros(N)
    ldiv!(f_temp, A, f2.(x))
    @test f_temp ≈ ghost_f ≈ analytic_f

    # Check \ for Matrix
    M2 = [f2.(x) 2.0*f2.(x) 10.0*f2.(x)]
    analytic_M = analytic_AL \ (M2 .- analytic_Ab)
    ghost_M = A \ M2
    @test analytic_M ≈ ghost_M

    # Check ldiv! for Matrix
    M_temp = zeros(N,3)
    ldiv!(M_temp, A, M2)
    @test M_temp ≈ analytic_M ≈ ghost_M

    # Additionally test that A\M2 ≈ [f, 2.0(f-4.0)+4.0, 10.0(f-4.0)+4.0]
    M = [f.(x) 2.0*(f.(x) .- 4.0).+4.0 10.0*(f.(x) .- 4.0).+4.0]
    @test M ≈ M_temp ≈ analytic_M ≈ ghost_M
end

@testset "Test Left Division L4 (fourth order)" begin

    # Test \ homogenous and inhomogenous BC
    dx = 0.01
    x = 0.01:dx:0.2
    N = length(x)
    u = sin.(x)

    L = CenteredDifference(4, 4, dx, N)
    Q = RobinBC((1.0, 0.0, sin(0.0)), (1.0, 0.0, sin(0.2+dx)), dx)
    A = L*Q

    analytic_L = fourth_deriv_approx_stencil(N) ./ dx^4
    analytic_QL = [transpose(zeros(N)); Diagonal(ones(N)); transpose(zeros(N))]
    analytic_AL = analytic_L*analytic_QL
    analytic_Qb = [zeros(N+1); sin(0.2+dx)]
    analytic_Ab = analytic_L*analytic_Qb

    analytic_u = analytic_AL \ (u - analytic_Ab)
    ghost_u = A \ u

    # Check that A\u.(x) is consistent with analytic_AL \ u.(x)
    @test analytic_u ≈ ghost_u

    # Check ldiv!
    u_temp = zeros(N)
    ldiv!(u_temp, A, u)
    @test u_temp ≈ ghost_u ≈ analytic_u

    # Check \ for Matrix
    M2 = [u 2.0*u 10.0*u]
    analytic_M = analytic_AL \ (M2 .- analytic_Ab)
    ghost_M = A \ M2
    @test analytic_M ≈ ghost_M

    # Check ldiv! for Matrix
    M_temp = zeros(N,3)
    ldiv!(M_temp, A, M2)
    @test M_temp ≈ ghost_M ≈ analytic_M
end


text:
@testset "Attitude" begin
    @testset "Static Attitude" begin
        next_attitude = @inferred GNSSSimulator.propagate(
            RotXYZ(1.0, 2.0, 3.0),
            1ms,
            Random.GLOBAL_RNG
        )
        @test next_attitude == RotXYZ(1.0, 2.0, 3.0)

        noisy_static_attitude = @inferred NoisyStaticAttitude(
            RotXYZ(1.0, 2.0, 3.0),
            0.1,
            0.2,
            0.3
        )
        rng = MersenneTwister(1234)
        next_noisy_static_attitude = @inferred GNSSSimulator.propagate(
            noisy_static_attitude,
            1ms,
            rng
        )
        @test next_noisy_static_attitude.base_attitude == RotXYZ(1.0, 2.0, 3.0)
        rng = MersenneTwister(1234)
        @test get_attitude(next_noisy_static_attitude) == RotXYZ(
            1.0 + randn(rng) * 0.1,
            2.0 + randn(rng) * 0.2,
            3.0 + randn(rng) * 0.3
        )

    end

    @testset "Dynamic Attitude" begin
        attitudes = [RotXYZ(0.1*i, 0.1*i+0.1, 0.1*i+0.2) for i = 0:10]
        dynamic_attitude = @inferred DynamicAttitude(attitudes, 1Hz, 0.0s)
        next_dynamic_attitude = @inferred GNSSSimulator.propagate(
            dynamic_attitude,
            1s,
            Random.GLOBAL_RNG
        )
        forward_dynamic_attitude = @inferred GNSSSimulator.propagate(
            dynamic_attitude,
            20s,
            Random.GLOBAL_RNG
        )
        linear_dynamic_attitude = @inferred LinearDynamicAttitude(
            RotXYZ(0.0, 0.1, 0.2),
            0.1rad/s,
            0.1rad/s,
            0.1rad/s
        )
        next_linear_dynamic_attitude = @inferred GNSSSimulator.propagate(
            linear_dynamic_attitude,
            1s,
            Random.GLOBAL_RNG
        )

        @test next_dynamic_attitude == DynamicAttitude(attitudes, 1Hz, 1.0s)
        @test get_attitude(next_dynamic_attitude) ≈ RotXYZ(0.1, 0.2, 0.3)
        @test get_attitude(forward_dynamic_attitude) == RotXYZ(1.0, 1.1, 1.2)
        @test next_linear_dynamic_attitude.attitude ≈ RotXYZ(0.1, 0.2, 0.3)
        @test get_attitude(next_linear_dynamic_attitude) ≈ RotXYZ(0.1, 0.2, 0.3)

        dynamic_attitude = @inferred DynamicAttitude(attitudes, 1Hz, time = 1.0s)
        @test dynamic_attitude.attitudes == attitudes
        @test dynamic_attitude.sample_freq == 1Hz
        @test dynamic_attitude.time == 1.0s

        attitudes2 = [RotXYZ(0.0, 0.1, 0.2) for i = 0:10]
        noisy_dynamic_attitude = @inferred NoisyDynamicAttitude(
            attitudes2,
            1Hz,
            0.1,
            0.2,
            0.3,
            time = 0.0s
        )

        rng = MersenneTwister(1234)
        next_noisy_dynamic_attitude = @inferred GNSSSimulator.propagate(
            noisy_dynamic_attitude,
            1s,
            rng
        )
        @test next_noisy_dynamic_attitude.attitudes == attitudes2
        rng = MersenneTwister(1234)
        @test get_attitude(next_noisy_dynamic_attitude) == RotXYZ(
            0.0 + randn(rng) * 0.1,
            0.1 + randn(rng) * 0.2,
            0.2 + randn(rng) * 0.3
        )
        noisy_linear_dynamic_attitude = @inferred NoisyLinearDynamicAttitude(
            RotXYZ(0.0, 0.1, 0.2),
            0.1rad/s,
            0.2rad/s,
            0.3rad/s,
            0.1,
            0.2,
            0.3
        )
        rng = MersenneTwister(1234)
        next_noisy_linear_dynamic_attitude = @inferred GNSSSimulator.propagate(
            noisy_linear_dynamic_attitude,
            1s,
            rng
        )
        @test next_noisy_linear_dynamic_attitude.base_attitude ≈ RotXYZ(0.1, 0.3, 0.5)
        rng = MersenneTwister(1234)
        @test get_attitude(next_noisy_linear_dynamic_attitude) ≈ RotXYZ(
            0.1 + randn(rng) * 0.1,
            0.3 + randn(rng) * 0.2,
            0.5 + randn(rng) * 0.3
        )
    end
end


text:
# Inference results
[
(
    method = "solve",
    mu = 0.2,

    ds = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.4, 0.4, 0.4, 0.4,
          0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7,
          0.8, 0.8, 0.8, 0.8, 0.8, 0.8],

    dfits = [0.109, 0.229, 0.028, 0.075, 0.184, 0.115, 0.105, 0.15, 0.035, 0.157, 0.275, 0.331, 0.228, 0.39, 0.209, 0.268, 0.279, 0.416, 0.129, 0.275, 0.447, 0.397, 0.123, 0.436,
             0.361, 0.6, 0.209, 0.594, 0.426, 0.328, 0.296, 0.613, 0.507, 0.559, 0.466, 0.341, 0.672, 0.62, 0.564, 0.427, 0.489, 0.557, 0.609, 0.611, 0.798, 0.85, 0.504, 0.63,
             0.719, 0.415, 0.887, 0.83, 0.698, 0.886, 0.675, 0.49, 0.756, 0.575, 0.63, 0.47, 0.918, 0.822, 0.869, 0.796, 0.586, 0.729, 0.912, 0.78, 0.796, 0.602, 0.907, 0.825,
             0.897, 0.717],

    mufits = [0.163, 0.259, 0.347, 0.373, 0.241, 0.258, 0.134, 0.213, 0.386, 0.203, 0.159, 0.154, 0.114, 0.092, 0.272, 0.203, 0.219, 0.079, 0.303, 0.243, 0.132, 0.176, 0.374,
              0.142, 0.146, 0.031, 0.319, 0.073, 0.185, 0.219, 0.194, 0.043, 0.176, 0.103, 0.138, 0.279, 0.111, 0.118, 0.129, 0.154, 0.162, 0.178, 0.132, 0.139, 0.033, 0.018,
              0.205, 0.119, 0.088, 0.282, 0.016, 0.053, 0.124, 0.033, 0.186, 0.298, 0.09, 0.173, 0.171, 0.297, 0.03, 0.097, 0.078, 0.107, 0.28, 0.172, 0.038, 0.144, 0.166, 0.335,
              0.083, 0.13, 0.085, 0.252]
),

(
    method = "solve",
    d = 0.4,

    mus = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.4, 0.4, 0.4,
           0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7,
           0.7, 0.7, 0.7, 0.7, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8],

    dfits = [0.206, 0.38, 0.158, 0.363, 0.368, 0.488, 0.32, 0.347, 0.602, 0.492, 0.35, 0.506, 0.227, 0.509, 0.287, 0.48, 0.505, 0.628, 0.52, 0.323, 0.436, 0.427, 0.433, 0.801,
             0.506, 0.663, 0.481, 0.558, 0.829, 0.551, 0.524, 0.441, 0.497, 0.406, 0.693, 0.365, 0.419, 0.468, 0.547, 0.716, 0.459, 0.432, 0.552, 0.416, 0.501, 0.367, 0.452,
             0.563, 0.388, 0.482, 0.506, 0.438, 0.488, 0.363, 0.585, 0.358, 0.529, 0.373, 0.573, 0.636, 0.507, 0.349, 0.46, 0.439, 0.558, 0.566, 0.645, 0.476, 0.387, 0.334, 0.439,
             0.56, 0.571, 0.588, 0.451, 0.556, 0.47, 0.623, 0.518, 0.471],

    mufits = [0.328, 0.14, 0.12, 0.091, 0.08, 0.085, 0.025, 0.127, 0.023, 0.068, 0.153, 0.101, 0.37, 0.105, 0.256, 0.183, 0.176, 0.102, 0.169, 0.288, 0.284, 0.27, 0.29, 0.076,
              0.23, 0.144, 0.25, 0.202, 0.036, 0.184, 0.271, 0.366, 0.335, 0.366, 0.186, 0.403, 0.393, 0.366, 0.263, 0.149, 0.41, 0.412, 0.32, 0.436, 0.375, 0.504, 0.418, 0.301,
              0.481, 0.405, 0.414, 0.486, 0.424, 0.553, 0.354, 0.55, 0.411, 0.547, 0.365, 0.297, 0.436, 0.589, 0.494, 0.516, 0.403, 0.397, 0.322, 0.484, 0.57, 0.63, 0.539, 0.415,
              0.411, 0.398, 0.533, 0.427, 0.514, 0.365, 0.455, 0.512]
),

(
    method = "fit",
    mu = 0.2,

    ds = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.4, 0.4, 0.4, 0.4,
          0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7,
          0.7, 0.7, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8],

    dfits = [0.109, 0.229, 0.028, 0.075, 0.184, 0.115, 0.105, 0.15, 0.035, 0.157, 0.275, 0.331, 0.228, 0.39, 0.209, 0.268, 0.279, 0.416, 0.129, 0.275, 0.447, 0.397, 0.123, 0.436,
             0.361, 0.6, 0.209, 0.594, 0.426, 0.328, 0.296, 0.613, 0.507, 0.559, 0.466, 0.341, 0.672, 0.62, 0.564, 0.427, 0.489, 0.557, 0.609, 0.611, 0.798, 0.85, 0.504, 0.63,
             0.719, 0.415, 0.887, 0.83, 0.698, 0.886, 0.675, 0.49, 0.756, 0.575, 0.63, 0.47, 0.918, missing, 0.822, 0.869, 0.796, 0.586, 0.729, 0.912, 0.78, missing, 0.796, 0.602,
             0.907, missing, 0.825, 0.897, missing, missing, 0.717, missing],

    mufits = [0.218, 0.456, 0.183, 0.349, 0.391, 0.304, 0.157, 0.301, 0.25, 0.298, 0.248, 0.304, 0.17, 0.303, 0.296, 0.293, 0.312, 0.289, 0.206, 0.35, 0.254, 0.269, 0.169, 0.263,
              0.213, 0.266, 0.238, 0.315, 0.302, 0.239, 0.117, 0.19, 0.276, 0.227, 0.188, 0.233, 0.324, 0.285, 0.253, 0.193, 0.164, 0.229, 0.217, 0.223, 0.244, 0.268, 0.22, 0.221,
              0.239, 0.223, 0.215, 0.215, 0.195, 0.249, 0.251, 0.22, 0.201, 0.169, 0.202, 0.218, 0.194, 0.21, 0.196, 0.221, 0.179, 0.21, 0.204, 0.206, 0.217, 0.24, missing,
              missing, missing, missing, missing, missing, missing, missing, missing, missing]
),

(
    method = "fit",
    d = 0.4,

    mus = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.4, 0.4, 0.4,
           0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7,
           0.7, 0.7, 0.7, 0.7, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8],

    dfits = [0.206, 0.38, 0.158, 0.363, 0.368, 0.488, 0.32, 0.347, 0.602, 0.492, 0.35, 0.506, 0.227, 0.509, 0.287, 0.48, 0.505, 0.628, 0.52, 0.323, 0.436, 0.427, 0.433, 0.801,
             0.506, 0.663, 0.481, 0.558, 0.829, 0.551, 0.524, 0.441, 0.497, 0.406, 0.693, 0.365, 0.419, 0.468, 0.547, 0.716, 0.459, 0.432, 0.552, 0.416, 0.501, 0.367, 0.452,
             0.563, 0.388, 0.482, 0.506, 0.438, 0.488, 0.363, 0.585, 0.358, 0.529, 0.373, 0.573, 0.636, 0.507, 0.349, 0.46, 0.439, 0.558, 0.566, 0.645, 0.476, 0.387, 0.334, 0.439,
             0.56, 0.571, 0.588, 0.451, 0.556, 0.47, 0.623, 0.518, 0.471],

    mufits = [0.147, 0.147, 0.0, 0.094, 0.08, 0.143, 0.0, 0.142, 0.161, 0.14, 0.111, 0.187, 0.222, 0.192, 0.176, 0.255, 0.266, 0.277, 0.265, 0.224, 0.33, 0.308, 0.336, 0.414,
              0.328, 0.369, 0.329, 0.343, 0.382, 0.314, 0.405, 0.425, 0.438, 0.39, 0.462, 0.391, 0.432, 0.445, 0.404, 0.434, 0.496, 0.469, 0.488, 0.476, 0.491, 0.502, 0.494,
              0.474, 0.5, 0.509, 0.551, 0.556, 0.543, 0.553, 0.566, 0.546, 0.572, 0.56, 0.573, 0.554, 0.584, 0.582, 0.649, 0.649, 0.652, 0.655, 0.716, 0.654, 0.655, 0.723, 0.834,
              0.947, 0.835, 0.95, 0.946, 0.828, 0.842, 0.947, 0.832, 0.84]
)

]

max_stars_repo_name: aangaji/Turnover

max_issues_repo_name: aangaji/Turnover

max_forks_repo_name: aangaji/Turnover


text:
using DistributionsAD

using Combinatorics
using Distributions
using FiniteDifferences
using ForwardDiff
using PDMats
using ReverseDiff
using Tracker
using Zygote

using Random, LinearAlgebra, Test

using Distributions: meanlogdet
using DistributionsAD: TuringUniform, TuringMvNormal, TuringMvLogNormal,
    TuringPoissonBinomial
using ForwardDiff: Dual
using StatsBase: entropy
using StatsFuns: binomlogpdf, logsumexp

const FDM = FiniteDifferences
const GROUP = get(ENV, "GROUP", "All")

if GROUP == "All" || GROUP == "Others"
    include("others.jl")
end

if GROUP == "All" || GROUP == "AD"
    include("ad/utils.jl")
    include("ad/distributions.jl")
end


text:
# PoC on jacobians for a factor

using IncrementalInference
using TensorCast
using Manifolds
using Test

# overloading with new dispatch
import IncrementalInference: getSample, getManifold

##

@testset "test manual call to gradient lambda utilities" begin
##

pp = LinearRelative(MvNormal([10;0],[1 0; 0 1]))
measurement = [10.0;0.0]
varTypes = (ContinuousEuclid{2}, ContinuousEuclid{2})
pts = ([0;0.0], [9.5;0])

##

J__, λ_fncs, λ_sizes = IIF._prepFactorGradientLambdas(pp, measurement, varTypes, pts; h=1e-4);

##

λ_fncs[1][1]()
λ_fncs[1][2]()
λ_fncs[2][1]()
λ_fncs[2][2]()

##

J__

@test norm( J__ - [0 0 1 0; 0 0 0 1; 1 0 0 0; 0 1 0 0] ) < 1e-4

## build new functor container

gradFct = FactorGradientsCached!(pp, varTypes, measurement, pts);

## test grad calc for current values

J_c = gradFct()

@test norm( J_c - [0 0 1 0; 0 0 0 1; 1 0 0 0; 0 1 0 0] ) < 1e-4

##

J = gradFct(measurement, pts...)

##

@test norm( J - [0 0 1 0; 0 0 0 1; 1 0 0 0; 0 1 0 0] ) < 1e-4

## check on transmitted info per coords

ret = calcPerturbationFromVariable(gradFct, [1=>[1;1]])

# the fromVar itself should be zero
@test length(ret[1]) == 2
@test isapprox( ret[1], [0;0], atol=1e-6 )

# the other variable
@test length(ret[2]) == 2
@test isapprox( ret[2], [1;1], atol=1e-6 )

##
end

##

struct TestPartialRelative2D{B <: SamplableBelief} <: IIF.AbstractRelativeMinimize
  Z::B
  partial::Tuple{Int}
end
# standard helper with partial set
TestPartialRelative2D(z::SamplableBelief) = TestPartialRelative2D(z, (2,))

# imported earlier for overload
getManifold(fnc::TestPartialRelative2D) = TranslationGroup(2)
getSample(cf::CalcFactor{<:TestPartialRelative2D}) = rand(cf.factor.Z, 1)

# currently requires residual to be returned as a tangent vector element
(cf::CalcFactor{<:TestPartialRelative2D})(z, x1, x2) = x2[2:2] - (x1[2:2] + z[1:1])

##

@testset "test a partial, binary relative factor perturbation (a new user factor)" begin
##

tpr = TestPartialRelative2D(Normal(10,1))

measurement = [10.0;]
pts = ([0;0.0], [0;10.0])
varTypes = (ContinuousEuclid{2}, ContinuousEuclid{2})

## construct the lambdas
gradients = FactorGradientsCached!(tpr, varTypes, measurement, pts);

## calculate new gradients
J = gradients(measurement, pts...)

## check on transmitted info per coords

ret = calcPerturbationFromVariable(gradients, [1=>[1;1]])

# the fromVar itself should be zero
@test length(ret[1]) == 2
@test isapprox( ret[1], [0;0], atol=1e-6 )

# the other variable only affects the second coordinate dimension
@test length(ret[2]) == 2
@test isapprox( ret[2], [0;1], atol=1e-6 )

## check the reverse perturbation

ret = calcPerturbationFromVariable(gradients, [2=>[1;1]])

# only the first coordinate dimension is affected
@test length(ret[1]) == 2
@test isapprox( ret[1], [0;1], atol=1e-6 )

## the fromVar itself should be zero
@test length(ret[2]) == 2
@test isapprox( ret[2], [0;0], atol=1e-6 )

##
end

#


text:
export dropout

type DropoutDesc
  ptr::Ptr{Void}
end

function DropoutDesc(dropout::Float64)
  p = Ptr{Void}[0]
  cudnnCreateDropoutDescriptor(p)
  cudnnSetDropoutDescriptor(p[1], h)

  desc = new(p[1])
  finalizer(desc, cudnnDestroyDropoutDescriptor)
  desc
end

Base.unsafe_convert(::Type{Ptr{Void}}, desc::DropoutDesc) = desc.ptr

function dropout()
  p = Ptr{Void}[0]
  cudnnCreateDropoutDescriptor(p)
  dropoutdesc = p[1]

  xdesc = TensorDesc(x)
  ydesc = TensorDesc(y)

  cudnnDropoutForward(h, dropoutdesc, xdesc, x, ydesc, y, reservespace, reservesize)

  cudnnDestroyDropoutDescriptor(dropoutdesc)
end


text:
@doc """
  Splits a floating-point number `a` into two parts `a = h + l` such that the
  halves `h` and `l` are non-overlapping and `|x| ≥ |y|`. That is, both halves
  have at most `s-1` nonzero bits, where `s = ceil(p/2)` and `p` is the
  precision (`p = 24` for `Float32` and `p = 53` for `Float64`).

  References:

    * T.J. Dekker, A floating-point technique for extending the available
      precision, Numerische Mathematik 18, pp. 224–242, 1971.
  """ ->
# TODO: Extend the definition to floating-point arrays
for T in (Float32, Float64)
  @eval begin
    function fsplit(a::$T)
      # Float32: factor = 2^12 + 1 = 4097
      # Float64: factor = 2^27 + 1 = 134217729
      f = $(convert(T, 2^ceil(precision(T)/2)+1));
      c = a * f
      h = c - (c-a)
      l = a - h
      h, l
    end
  end
end

@doc """
  Computes the error-free multiplication `a * b = x + y`, where `x = fl(a * b)`,
  for floating-point numbers `a` and `b`. This algorithm uses all cross products
  between the halves generated by `fsplit` from `a` and `b`. It requires 17
  flops.

  References:

    * T.J. Dekker, A floating-point technique for extending the available
      precision, Numerische Mathematik 18, pp. 224–242, 1971.
  """ ->
# TODO: Extend the definition to floating-point arrays
function err_mul{T<:IEEE754}(a::T, b::T)
  x = a * b
  ah, al = fsplit(a)
  bh, bl = fsplit(b)
  y = al*bl - (((x - ah*bh) - al*bh) - ah*bl)
  x, y
end

@doc """
  Computes the error-free transformation `a * b = x + y`, where `x = fl(a * b)`,
  for floating-point numbers `a` and `b`. By using the fused multiply-add (FMA)
  operation, it only requires 2 flops instead of the 17 flops required by
  `err_mul` but the improvement in execution time typically a less than that.

  References:

    * T. Ogita, S.M. Rump, and S. Oishi, Accurate sum and dot product, SIAM
      Journal on Scientific Computing 26, pp. 1955–1988, 2005.
  """ ->
# TODO: Extend the definition to floating-point arrays
function err_fast_mul{T<:IEEE754}(a::T, b::T)
  x = a * b
  y = fma(a, b, -x)
  x, y
end


text:
function Zksτ(k, two_s, two_τ, two_λs, two_λs′, σs, tbs; refζs=(1,2,3,1))
    i,j = ij_from_k(k)
    #
    ms² = tbs.ms^2
    two_js = tbs.two_js
    two_j0 = two_js[4]
    two_λ0 = two_λs[4]
    #
    w0 = wr(k,refζs[4],0);
    wi = wr(k,refζs[i],i);
    wj = wr(k,refζs[j],j);
    wk = wr(k,refζs[k],k);
    #
    cosθ = cosθij(k,σs,ms²)
    #
    val =   phase(tbs.two_js[k]-two_λs′[k]) * # particle-2 convention
        sqrt(two_s+1) * wignerd_doublearg(two_s, two_τ, two_λs′[i]-two_λs′[j], cosθ) *
            phase(tbs.two_js[j]-two_λs′[j]) # particle-2 convention
    #
    val *= two_j0==0 ? (two_τ == two_λs′[k]) : wignerd_doublearg_sign(two_j0, two_λ0, two_τ-two_λs′[k], cosζ(w0, σs, ms²), ispositive(w0))
    val *= two_js[i]==0 ? 1 : wignerd_doublearg_sign(two_js[i], two_λs′[i], two_λs[i], cosζ(wi, σs, ms²), ispositive(wi))
    val *= two_js[j]==0 ? 1 : wignerd_doublearg_sign(two_js[j], two_λs′[j], two_λs[j], cosζ(wj, σs, ms²), ispositive(wj))
    val *= two_js[k]==0 ? 1 : wignerd_doublearg_sign(two_js[k], two_λs′[k], two_λs[k], cosζ(wk, σs, ms²), ispositive(wk))
    #
    return val
end

wignerd_doublearg_sign(two_j,two_λ1, two_λ2, cosθ, ispositive) =
    (ispositive ? 1 : x"-1"^div(two_λ1-two_λ2, 2)) *
        wignerd_doublearg(two_j, two_λ1, two_λ2, cosθ)

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
# Copyright 2020 Google LLC

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

#     https://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
raw_input = open(f->read(f, String), "day04input.txt")

passport_strings = split(raw_input, "\n\n")

function parsePassport(passport)
    fields = Dict{String, String}()
    words = split(passport)
    for word ∈ words
        field, value = split(word, ':')
        fields[field] = value
    end
    return fields
end

passports = map(parsePassport, passport_strings)

function hasRequiredFields(passport)
    haskey(passport, "byr") && haskey(passport, "iyr") && haskey(passport, "eyr") && haskey(passport, "hgt") && haskey(passport, "hcl") && haskey(passport, "ecl") && haskey(passport, "pid")
end

part1(passports) = count(hasRequiredFields, passports)

function byr_valid(byr)
    byr = parse(Int, byr)
    if byr < 1920 || byr > 2002
        return false
    end
    return true
end

function iyr_valid(iyr)
    iyr = parse(Int, iyr)
    if iyr < 2010 || iyr > 2020
        return false
    end
    return true
end

function eyr_valid(eyr)
    eyr = parse(Int, eyr)
    if eyr < 2020 || eyr > 2030
        return false
    end
    return true
end

function hgt_valid(hgt)
    hgtmatch = match(r"^([0-9]+)(cm|in)$", hgt)
    if hgtmatch === nothing
        return false
    end
    hgt = parse(Int, hgtmatch.captures[1])
    if hgtmatch.captures[2] == "cm"
        if hgt < 150 || hgt > 193
            return false
        end
    else
        if hgt < 59 || hgt > 76
            return false
        end
    end
    return true
end

function fieldsValidate(passport)
    if !byr_valid(passport["byr"])
        return false
    end
    if !iyr_valid(passport["iyr"])
        return false
    end
    if !eyr_valid(passport["eyr"])
        return false
    end
    if !hgt_valid(passport["hgt"])
        return false
    end
    hcl = match(r"^#[0-9a-f]{6}$", passport["hcl"])
    if hcl === nothing
        return false
    end
    ecl = match(r"^(amb|blu|brn|gry|grn|hzl|oth)$", passport["ecl"])
    if ecl === nothing
        return false
    end
    pid = match(r"^[0-9]{9}$", passport["pid"])
    if pid === nothing
        return false
    end
    return true
end

count(p->hasRequiredFields(p) && fieldsValidate(p), passports)


text:
__precompile__(true)

"""
# HavenOnDemand

Julia package to access HPE Haven OnDemand API.

* https://www.havenondemand.com
"""
module HavenOnDemand

export set_api_key,
       HODException,

       # Missing
       # list_stuff
       # job_status
       # job_result

       # Speech Recognition
       recognize_speech,

       # Connectors
       cancel_connector_schedule,
       connector_history,
       connector_status,
       create_connector,
       delete_connector,
       retrieve_config,
       start_connector,
       stop_connector,
       update_connector,

       # Format Conversion
       expand_container,
       ocr_document,
       store_object,
       extract_text,
       view_document,

       # Image Analysis
       recognize_barcodes,
       detect_faces,
       recognize_images,

       # Graph Analysis
       get_common_neighbors,
       get_neighbors,
       get_nodes,
       get_shortest_path,
       get_sub_graph,
       suggest_links,
       summarize_graph,

       # Policy
       create_classification_objects,
       create_policy_objects,
       delete_classification_objects,
       delete_policy_objects,
       retrieve_classification_objects,
       retrieve_policy_objects,
       update_classification_objects,
       update_policy_objects,

       # Prediction
       predict,
       recommend,
       train_predictor,

       # Query Profile and Manipulation
       create_query_profile,
       delete_query_profile,
       retrieve_query_profile,
       update_query_profile,

       # Search
       find_related_concepts,
       find_similar,
       get_content,
       get_parametric_values,
       query_text_index,
       retrieve_index_fields,

       # Text Analysis
       auto_complete,
       classify_document,
       extract_concepts,
       categorize_document,
       entity_extraction,
       expand_terms,
       highlight_text,
       identify_language,
       analyze_sentiment,
       tokenize_text,

       # Unstructured Text Indexing
       add_to_text_index,
       create_text_index,
       delete_text_index,
       delete_from_text_index,
       index_status,
       list_resources,
       restore_text_index

using DataFrames: DataFrame, eachrow, readtable
using HttpCommon: Response
using   Requests: json, post

import Base: showerror

typealias String AbstractString

"`HavenOnDemand` Exception."
immutable HODException <: Exception
    response::Response
end

function showerror(io::IO, ex::HODException)
    endpoint = match(r"sync/(.+)/", get(ex.response.request).resource)[1]
    println(io, typeof(ex), ": $endpoint")
    for (k, v) in json(ex.response)
        print(io, ucfirst(k), ": ")
        println(io, v)
    end
    println(io, "For more information, visit:")
    println(io, "* https://dev.havenondemand.com/apis/$(endpoint)")
end

"Main `HavenOnDemand` function, used to call **Haven OnDemand** API."
function call_HOD(
        endpoint        :: String,
        async           :: Bool;
        api_url         :: String = "https://api.havenondemand.com",
        version         :: Int    = 1,
        default_version :: Int    = 1,
        options         :: Dict   = Dict()
    )
    try
        options["apikey"] = _HOD_API_KEY
    catch
        error("Use `HavenOnDemand.set_api_key(api_key::AbstractString)` first.")
    end
    sync_str = async ? "async" : "sync"
    r = post(
        "$(api_url)/$(version)/api/$(sync_str)/$(endpoint)/v$(default_version)",
        data = options
    )
    return r.status == 200 ? json(r) : throw(HODException(r))
end

"""
Sets the **HavenOn Demand** API key, this function must be called after
`using HavenOnDemand` with a valid key, in order to be able to use the Haven
OnDemand *API*.
"""
function set_api_key(api_key::String)
    global _HOD_API_KEY
    const _HOD_API_KEY = api_key
    return nothing
end

"""
`DataFrame` that holds the `HavenOnDemand` data used to metaprogram a wrapper
function for each endpoint in the API.
"""
const HOD_API = readtable(
    joinpath(Pkg.dir("HavenOnDemand"), "src", "api.data"),
    separator = ' ',
)

# Meta wrap the API:
for row in eachrow(HOD_API::DataFrame)
    func_name   = row[:Func_Name]
    endpoint    = row[:Endpoint]
    async_only  = row[:Async_Only]
    description = row[:Description]
    title = join([ucfirst(s) for s in split(func_name, '_')], ' ')
    docstring = (
        """
        **HPE Haven OnDemand: $(title)**

        `$(func_name)([kwargs...])`

        $description

        All the arguments are optional and they must be supplied
        as keyword arguments, non valid keyword names are ignored.

        For information about valid arguments, visit:

        * https://dev.havenondemand.com/apis/$(endpoint)
        """
    )
    @eval begin
        @doc $docstring ->
        function $(symbol(func_name))(; kwargs...)
            return call_HOD($endpoint, $async_only, options = Dict(kwargs))
        end
    end
end

end