text:
using TrajectoryOptimization
using Random
using StaticArrays
using LinearAlgebra
using MeshCat
using TrajOptPlots
using Rotations
using BenchmarkTools
import TrajectoryOptimization: Rotation
const TO = TrajectoryOptimization

max_con_viol = 1.0e-3
T = Float64
verbose = true

opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,
    iterations=300)

opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,
    opts_uncon=opts_ilqr,
    iterations=40,
    # cost_tolerance=1.0e-5,
    cost_tolerance_intermediate=1.0e-4,
    # constraint_tolerance=max_con_viol,
    # penalty_scaling=10.,
    penalty_initial=10.)

opts_pn = ProjectedNewtonSolverOptions{T}(verbose=verbose,
    feasibility_tolerance=max_con_viol,
    solve_type=:feasible)

opts_altro = ALTROSolverOptions{T}(verbose=verbose,
    opts_al=opts_al,
    # R_inf=1.0e-4,
    resolve_feasible_problem=false,
    opts_pn=opts_pn,
    projected_newton=false,
    projected_newton_tolerance=1.0e-4)

model = Dynamics.Quadrotor2{UnitQuaternion{Float64,IdentityMap}}(use_rot=false)
n,m = 13,4

if !isdefined(Main,:vis)
    vis = Visualizer(); open(vis);
    set_mesh!(vis, model)
end

function gen_quad_maze(Rot; kwargs...)
    prob = Problems.QuadrotorMaze(Rot; kwargs...)
    solver = ALTROSolver(prob, opts_altro, infeasible=true)
    solver.solver_al.opts.verbose = false
    solver
end

# Initialize visualization
solver = gen_quad_maze(UnitQuaternion{Float64,IdentityMap}, costfun=:Quadratic)
solver = gen_quad_maze(RPY{Float64}, use_rot=false)
solver = gen_quad_maze(UnitQuaternion{Float64,CayleyMap}, costfun=:QuatLQR)
# solver = gen_quad_maze(UnitQuaternion{Float64,IdentityMap}, use_rot=:slack, normcon=true)
Z0 = deepcopy(get_trajectory(solver))
opts_ilqr.verbose = false
opts_altro.projected_newton = false
solver.opts.verbose = true
add_cylinders!(vis, solver, height=4, robot_radius=1.0)
visualize!(vis, solver)

# Solve
# opts_altro.R_inf = 1e-8
# opts_al.penalty_initial = 10
# opts_al.penalty_scaling = 10
opts_al.cost_tolerance_intermediate = 1e-2
solver = gen_quad_maze(UnitQuaternion{Float64,IdentityMap}, costfun=:Quadratic, use_rot=false)
solver = gen_quad_maze(UnitQuaternion{Float64,IdentityMap}, costfun=:Quadratic, use_rot=:slack, normcon=true)
solver = gen_quad_maze(UnitQuaternion{Float64,IdentityMap}, costfun=:Quadratic, use_rot=false, normcon=true)
# solver = gen_quad_maze(RPY{Float64}, use_rot=false)
solver = gen_quad_maze(UnitQuaternion{Float64,CayleyMap}, costfun=:QuatLQR)
Z0 = deepcopy(get_trajectory(solver))
initial_trajectory!(solver, Z0)
@time solve!(solver)
iterations(solver)
visualize!(vis, solver)
benchmark_solve!(solver)
conSet = get_constraints(solver)

# Plot waypoints
waypoints!(vis, get_model(solver), get_trajectory(solver), length=41)

# Step through solve
prob = Problems.QuadrotorMaze(Rot,use_rot=true, costfun=:Quadratic)
prob.constraints.constraints
opts_ilqr.verbose = false
opts_altro.projected_newton = false
opts_altro.R_inf = 1e-8
opts_al.penalty_initial = 1
solver = ALTROSolver(prob, opts_altro, infeasible=true)

al = solver.solver_al
initialize!(al)
TO.set_tolerances!(al, al.solver_uncon, 1)
visualize!(vis, al)

TO.step!(al)
max_violation(al)
visualize!(vis, al)

function compare_quats(gen_prob; kwargs...)
    # Initial data structure
    data = Dict{Symbol,Vector}(:name=>Symbol[], :costfun=>Symbol[], :iterations=>Int[],
        :rotation=>Symbol[], :time=>Float64[], :cost=>Float64[])

    # Treat quaternion as normal vector
    println("Quat")
    log_solve!(gen_prob(UnitQuaternion{Float64, IdentityMap}, use_rot=false),
        data; name=:Quat, kwargs...)

    # Re-normalize after disretization
    println("ReNorm")
    log_solve!(gen_prob(UnitQuaternion{Float64, ReNorm}, use_rot=false),
        data; name=:ReNorm, kwargs...)

    # Use Unit Norm Constraint
    println("NormCon")
    log_solve!(gen_prob(UnitQuaternion{Float64, IdentityMap}, use_rot=false,
        normcon=true), data; name=:NormCon, kwargs...)

    # Use Unit Norm Constrained with slack
    println("QuatSlack")
    log_solve!(gen_prob(UnitQuaternion{Float64, IdentityMap}, use_rot=:slack,
        normcon=true), data; name=:QuatSlack, kwargs...)

    for rmap in [ExponentialMap, CayleyMap, MRPMap, VectorPart]
        println(rmap)
        for costfun in [:QuatLQR, :Quadratic] #[:Quadratic, :QuatLQR, :ErrorQuad]
            log_solve!(gen_prob(UnitQuaternion{Float64, rmap}, costfun=costfun,
                normcon=false), data; kwargs...)
        end
    end
    return data
end

qdata_maze = compare_quats(gen_quad_maze, samples=1, evals=1)

# Run benchmarks
function run_all(gen_prob;kwargs...)
    # Initial data structure
    data = Dict{Symbol,Vector}(:name=>Symbol[], :costfun=>Symbol[], :iterations=>Int[], :rotation=>Symbol[],
        :time=>Float64[], :cost=>Float64[])

    # Baseline methods
    println("Baseline methods")
    log_solve!(gen_prob(UnitQuaternion{Float64, IdentityMap}, use_rot=false), data; kwargs...)
    # log_solve!(gen_prob(RodriguesParam{Float64}, use_rot=false), data; kwargs...)
    # log_solve!(gen_prob(MRP{Float64}, use_rot=false), data; kwargs...)
    log_solve!(gen_prob(RPY{Float64}, use_rot=false), data; kwargs...)

    # Quaternion Methods
    println("Quaternion methods")
    for rmap in [CayleyMap,] #[ExponentialMap, CayleyMap, MRPMap, VectorPart]
        for costfun in [:QuatLQR,] #[:Quadratic, :QuatLQR, :ErrorQuad]
            log_solve!(gen_prob(UnitQuaternion{Float64, rmap}, costfun=costfun), data; kwargs...)
        end
    end
    return data
end
data_maze = run_all(gen_quad_maze, samples=10, evals=1)

df = DataFrame(data_maze)
df.rots = string.(short_names.(df.rotation))
df.time_per_iter = df.time ./ df.iterations
quats = in([:ExponentialMap, :CayleyMap, :MRPMap, :VectorPart])
bases = in([:IdentityMap, :MRP, :RP, :RPY])
comp = in([:IdentityMap, :RPY, :CayleyMap])

############################################################################################
# NEW SIMPLIFIED PLOTS
############################################################################################
label = ["Quat","RPY","iMLQR"]
p = @pgf Axis(
    {
        ybar,
        ylabel="solve time (ms)",
        # "enlarge limits" = 0.20,
        legend_style =
        {
            at = Coordinate(0.5, -0.07),
            anchor = "north",
            legend_columns = -1
        },
        symbolic_x_coords=label,
        xtick = "data",
        nodes_near_coords,
        nodes_near_coords_align={vertical},
        "every node near coord/.append style={/pgf/number format/.cd, fixed,precision=0}"
    },
    Plot(Coordinates(label, df.time)),
    # Legend(["Quadratic", "Error Quadratic", "Geodesic"])
)
pgfsave("figs/maze_time_best.tikz", p, include_preamble=false)

############################################################################################
# PLOT 1: Iterations
############################################################################################

base = df[bases.(df.rotation),:]
quat = df[quats.(df.rotation),:]
qlqr  = quat[quat.costfun .== :QuatLQR,:]
equad = quat[quat.costfun .== :ErrorQuadratic,:]
quad  = quat[quat.costfun .== :Quadratic,:]
coord_base = Coordinates(base.rots, base.iterations)
coord1 = Coordinates(qlqr.rots, qlqr.iterations)
coord2 = Coordinates(equad.rots, equad.iterations)
coord3 = Coordinates(quad.rots, quad.iterations)

p = @pgf Axis(
    {
        ybar,
        ylabel="iterations",
        "enlarge x limits" = 0.20,
        legend_style =
        {
            at = Coordinate(0.5, -0.15),
            anchor = "north",
            legend_columns = -1
        },
        symbolic_x_coords=quad.rots,
        xtick = "data",
        nodes_near_coords,
        nodes_near_coords_align={vertical},
    },
    Plot(coord3),
    Plot(coord2),
    Plot(coord1),
    Legend(["Quadratic", "Error Quadratic", "Geodesic"])
)
pgfsave("figs/maze_quat_iters.tikz", p, include_preamble=false)

############################################################################################
# PLOT 2: Time
############################################################################################
df[:,[:rotation,:costfun,:time]]
base = df[bases.(df.rotation),:]
quat = df[quats.(df.rotation),:]
qlqr  = quat[quat.costfun .== :QuatLQR,:]
equad = quat[quat.costfun .== :ErrorQuadratic,:]
quad  = quat[quat.costfun .== :Quadratic,:]
quad_all = df[df.costfun .== :Quadratic,:]
best = [base; qlqr]
p = @pgf Axis(
    {
        ybar,
        ylabel="solve time (ms)",
        # "enlarge limits" = 0.20,
        legend_style =
        {
            at = Coordinate(0.5, -0.07),
            anchor = "north",
            legend_columns = -1
        },
        symbolic_x_coords=best.rots,
        xtick = "data",
        nodes_near_coords,
        nodes_near_coords_align={vertical},
        "every node near coord/.append style={/pgf/number format/.cd, fixed,precision=0}"
    },
    Plot(Coordinates(best.rots, best.time)),
    # Legend(["Quadratic", "Error Quadratic", "Geodesic"])
)
pgfsave("figs/maze_time_best.tikz", p, include_preamble=false)

coord_base = Coordinates(base.rots, base.time)
coord_all = Coordinates(quad_all.rots, quad_all.time)
coord1 = Coordinates(qlqr.rots, qlqr.time)
coord2 = Coordinates(equad.rots, equad.time)
coord3 = Coordinates(quad.rots, quad.time)

p = @pgf Axis(
    {
        ybar,
        width="20cm",
        ylabel="solve time (ms)",
        # "enlarge limits" = 0.20,
        legend_style =
        {
            at = Coordinate(0.5, -0.07),
            anchor = "north",
            legend_columns = -1
        },
        symbolic_x_coords=quad_all.rots,
        xtick = "data",
        nodes_near_coords,
        nodes_near_coords_align={vertical},
        "every node near coord/.append style={/pgf/number format/.cd, fixed,precision=0}"
    },
    Plot(coord_all),
    Plot(coord2),
    Plot(coord1),
    Legend(["Quadratic", "Error Quadratic", "Geodesic"])
)
pgfsave("figs/maze_time_all.tikz", p, include_preamble=false)

############################################################################################
# PLOT 4: Quat method comparison - time
############################################################################################

df = DataFrame(qdata_maze)
df.rots = string.(short_names.(df.name))
sub = df[in([:Quat,:ReNorm,:NormCon,:QuatSlack,:CayleyMap,:RPY]).(df.name),:]
quats = by(df, :rots, :time=>minimum)

p = @pgf Axis(
    {
        ybar,
        ylabel="solve time (ms)",
        # "enlarge limits" = 0.20,
        legend_style =
        {
            at = Coordinate(0.5, -0.07),
            anchor = "north",
            legend_columns = -1
        },
        symbolic_x_coords=quats.rots,
        xtick = "data",
        nodes_near_coords,
        nodes_near_coords_align={vertical},
        "every node near coord/.append style={/pgf/number format/.cd, fixed,precision=0}"
    },
    # HLine({"ultra thick", "white"}, 35),
    Plot(Coordinates(quats.rots, quats.time_minimum)),
    # Legend(["Quadratic", "Error Quadratic", "Geodesic"])
)
pgfsave("figs/maze_qcomp.tikz", p, include_preamble=false)


text:
# Transcode
# =========

"""
    transcode(::Type{C}, data::Vector{UInt8})::Vector{UInt8} where C<:Codec

Transcode `data` by applying a codec `C()`.

Note that this method does allocation and deallocation of `C()` in every call,
which is handy but less efficient when transcoding a number of objects.
`transcode(codec, data)` is a recommended method in terms of performance.

Examples
--------

```julia
julia> using CodecZlib

julia> data = b"abracadabra";

julia> compressed = transcode(ZlibCompressor, data);

julia> decompressed = transcode(ZlibDecompressor, compressed);

julia> String(decompressed)
"abracadabra"

```
"""
function Base.transcode(::Type{C}, data::ByteData) where C<:Codec
    codec = C()
    initialize(codec)
    try
        return transcode(codec, data)
    finally
        finalize(codec)
    end
end

"""
    transcode(codec::Codec, data::Vector{UInt8})::Vector{UInt8}

Transcode `data` by applying `codec`.

Note that this method does not initialize or finalize `codec`. This is
efficient when you transcode a number of pieces of data, but you need to call
[`TranscodingStreams.initialize`](@ref) and
[`TranscodingStreams.finalize`](@ref) explicitly.

Examples
--------

```julia
julia> using CodecZlib

julia> data = b"abracadabra";

julia> codec = ZlibCompressor();

julia> TranscodingStreams.initialize(codec)

julia> compressed = transcode(codec, data);

julia> TranscodingStreams.finalize(codec)

julia> codec = ZlibDecompressor();

julia> TranscodingStreams.initialize(codec)

julia> decompressed = transcode(codec, compressed);

julia> TranscodingStreams.finalize(codec)

julia> String(decompressed)
"abracadabra"

```
"""
function Base.transcode(codec::Codec, data::ByteData)
    input = Buffer(data)
    output = Buffer(initial_output_size(codec, buffermem(input)))
    error = Error()
    code = startproc(codec, :write, error)
    if code === :error
        @goto error
    end
    n = minoutsize(codec, buffermem(input))
    @label process
    makemargin!(output, n)
    Δin, Δout, code = process(codec, buffermem(input), marginmem(output), error)
    @debug(
        "called process()",
        code = code,
        input_size = buffersize(input),
        output_size = marginsize(output),
        input_delta = Δin,
        output_delta = Δout,
    )
    consumed!(input, Δin)
    supplied!(output, Δout)
    if code === :error
        @goto error
    elseif code === :end
        if buffersize(input) > 0
            if startproc(codec, :write, error) === :error
                @goto error
            end
            n = minoutsize(codec, buffermem(input))
            @goto process
        end
        resize!(output.data, output.marginpos - 1)
        return output.data
    else
        n = max(Δout, minoutsize(codec, buffermem(input)))
        @goto process
    end
    @label error
    if !haserror(error)
        set_default_error!(error)
    end
    throw(error[])
end

# Return the initial output buffer size.
function initial_output_size(codec::Codec, input::Memory)
    return max(
        minoutsize(codec, input),
        expectedsize(codec, input),
        8,  # just in case where both minoutsize and expectedsize are foolish
    )
end

ext: jl


text:
module NamedIndicesMeta

using NamedDims, ImageCore, ImageMetadata, ImageAxes, FieldProperties,
      AxisIndices, LightGraphs, SimpleWeightedGraphs, Reexport, MappedArrays
using Base: tail

import ImageAxes: timeaxis, timedim, colordim, checknames, isstreamedaxis

export
    NamedDimsArray,
    AxisIndicesArray,
    NIMArray,
    NIArray,
    IMArray,
    MetaArray,
    properties,
    dimnames,
    getdim,
    named_axes,
    @defdim

include("types.jl")
include("properties.jl")
include("data_format.jl")
include("defdim.jl")

include("time.jl")
using .TimeData

include("color.jl")
using .ColorData

include("observation.jl")
using .ObservationData

include("spatial.jl")
using .SpatialData

include("graphs.jl")
using .AxisIndicesGraphs

include("deprecations.jl")

@reexport using NamedIndicesMeta.TimeData

@reexport using NamedIndicesMeta.ColorData

@reexport using NamedIndicesMeta.ObservationData

@reexport using NamedIndicesMeta.SpatialData

@reexport using NamedIndicesMeta.AxisIndicesGraphs

end # module


text:
@testset "Testing Asap2" begin
    input_file = joinpath(@__DIR__, "aes.json")
    c = AsapMapper.PMConstructor(
        input_file,
        (
            architecture = "Array_Asap2",
            ruleset = AsapMapper.Asap2(),
        )
    )

    m = AsapMapper.build_map(c)
    m = AsapMapper.asap_pnr(m, move_attempts = 2000)
    AsapMapper.Mapper2.MapperCore.report_routing_stats(m)
end


text:
# define aliases
const SPREADSHEET_ALIASES = ["Spreadsheet", "spreadsheet", "ss", "Ss", "SS"]
const TRANSPOSE_ALIASES = ["Spreadsheet_transpose", "Spreadsheet_Transpose", "spreadsheet_transpose", "transpose", "Transpose", "sst", "Sst", "SST"]
const SECTIONWISE_ALIASES = ["sectionwise", "sw"]

"""
    curead(path::String; mode::String="auto")

**Keyword Arguments**
- `mode::String`: "sectionwise", "spreadsheet", "spreadsheet_transpose"

!!! note "COMSOL"
    Latest tested version: 5.3.0.223
"""
function curead(path::String; mode::String="auto")
    data = readdlm(path, ','; use_mmap=true) # read file @ path
    idx_first_value = _firstdataline(data) # get first line of values

    # auto-detect correct mode
    mode = mode == "auto" ? sum(data[idx_first_value-1,:] .== "") > 0 ? "sst" : "ss" : mode

    # switch to correct function depending on export mode
    return mode ∈ SPREADSHEET_ALIASES ? curead_spreadsheet(path, data, idx_first_value) :
           mode ∈ TRANSPOSE_ALIASES ? curead_spreadsheet_transpose(path, data, idx_first_value) :
           mode ∈ SECTIONWISE_ALIASES ? curead_sectionwise(path, data, idx_first_value) :
           error("Invalid argument for mode=\"$mode\".")

end

"""
    curead_spreadsheet(path::String, data::Matrix{Any}, idx_first_value::Int)
"""
function curead_spreadsheet(path::String, data::Matrix{Any}, idx_first_value::Int)
    meta = Dict{String, Any}() # meta data from file's header
    for i in 1:idx_first_value-2
        meta[replace(data[i,1], "% " => "")] = data[i,2]
    end
    fields = replace.(data[idx_first_value - 1, :], "% " => "") # field names
    values = Dict{String, Any}() # exported values of fields
    for i in 1:size(data, 2)
        values[fields[i]] = data[idx_first_value:end,i]
    end

    return CDataRaw(path, meta, fields, values)
end

"""
    _curead_spreadsheet_transpose(path::String)
"""
function curead_spreadsheet_transpose(path::String, data::Matrix{Any}, idx_first_value::Int)
    meta = Dict{String, Any}()  # meta data from file's header
    for i in 1:idx_first_value-2
        meta[replace(data[i,1], "% " => "")] = data[i,2]
    end
    fields = replace.(data[idx_first_value - 1:2:end, 1], "% " => "") # field names
    values = Dict{String, Any}() # values of the respective field
    for i in 1:length(fields)
        values[fields[i]] = data[idx_first_value + 2*(i-1),:]
    end

    return CDataRaw(path, meta, fields, values) # returns custom data-type
end

"""
    _curead_sectionwise(path::String)
"""
function curead_sectionwise(path::String)
    error("Sectionwise exported data not supported.")
end

## INTERNAL FUNCTIONS ======================================================= ##
function _firstdataline(data::Matrix{Any})
    for idx in 1:size(data, 1)
        if !(typeof(data[idx,1]) <: Union{String, SubString}) # could break if value is String
            return idx
        end
    end
end


text:
using SpinQubits, LinearAlgebra, Test

@testset verbose=true "All tests" begin
	@testset verbose=true "Comparison w/ Mathematica tests" begin
		@testset "Noiseless fidelities" begin
			@test calculateFidelities(4, 0.00, 0.0, zeros(3), 10, 3.0; verbose=false)[2][end] ≈ 0.0000148240466523175 rtol=1e-8
			@test calculateFidelities(4, 0.01, 0.0, zeros(3), 10, 3.0; verbose=false)[2][end] ≈ 0.0018048111455146731 rtol=1e-8
			@test calculateFidelities(4, 0.05, 0.0, zeros(3), 10, 3.0; verbose=false)[2][end] ≈ 0.0396766971366405 rtol=1e-8
		end
		@testset "Noisy fidelities" begin
			@test calculateFidelities(4, 0.01, 0.0, [0.01, 0.00, 0.00], 5000, 3.0; verbose=false)[2][end] ≈ 0.004784666103416213 rtol=5e-2
			@test calculateFidelities(4, 0.02, 0.0, [0.01, 0.00, 0.00], 5000, 3.0; verbose=false)[2][end] ≈ 0.00951159418717562 rtol=5e-2
		end
	end
	@testset "utils.jl" begin
		@testset "diagexp!" begin
			ringo = Diagonal(ones(3))
			SpinQubits.diagexp!(ringo)
			@test ringo == Diagonal(fill(exp(1.0),3))
			paul = Diagonal(zeros(3))
			SpinQubits.diagexp!(paul)
			@test paul == Diagonal(ones(3))
			for i in 1:10
				startDiag = randn(3)
				resultDiag = exp.(startDiag)
				testMat = Diagonal(startDiag)
				SpinQubits.diagexp!(testMat)
				@test testMat == Diagonal(resultDiag)
			end
		end
		@testset "zeros!" begin
			testMat = zeros(5,5)
			SpinQubits.zeros!(testMat)
			@test testMat == zeros(5,5)
			for i in 1:10
				testMat = randn(5,5)
				SpinQubits.zeros!(testMat)
				@test testMat == zeros(5,5)
			end
		end
	end

	@testset "spinors.jl" begin
		spinor1 = SpinQubits.Spinor([[1,0,0,1],[0,0,1,1]],[2.0,1.0])
		spinor2 = SpinQubits.Spinor([[1,0,0,1],[0,1,0,1]],[im,1.0])
		@testset "isequal" begin
			@test spinor1 != spinor2
			@test SpinQubits.isequal(spinor1,spinor1*1.0)
			@test spinor1 == 1.0*spinor1*1.0
			@test SpinQubits.Spinor([[1,0]],[1]) == SpinQubits.Spinor([[1,0]], [1])
		end
		@testset "+" begin
			@test spinor1 + spinor1 == SpinQubits.Spinor([[1,0,0,1],[0,0,1,1]],[4.0,2.0])
			@test spinor2 + spinor2 == SpinQubits.Spinor([[1,0,0,1],[0,1,0,1]],[2*im,2.0])
			@test spinor1 + spinor2 == SpinQubits.Spinor([[1,0,0,1],[0,0,1,1],[0,1,0,1]],[2.0 + 1.0im,1.0,1.0])
			@test spinor1 + spinor2 == spinor2 + spinor1
			@test_throws MethodError 3.0 + spinor1
		end
		@testset "*" begin
			@test spinor1 * spinor2 ≈ 2*im
			@test spinor2 * spinor1 ≈ -2*im
			@test spinor1 * spinor1 ≈ 5.0
			@test spinor2 * spinor2 ≈ 2.0
			@test 3.0 * spinor1 == SpinQubits.Spinor([[1,0,0,1],[0,0,1,1]],[6.0,3.0])
			@test im * spinor1 == SpinQubits.Spinor([[1,0,0,1],[0,0,1,1]],[2.0im,im])
			@test (spinor2 * -im) * spinor1 == spinor2 * (im * spinor1)
		end
		@testset "normalize" begin
			@test SpinQubits.normalize(spinor1) == SpinQubits.Spinor([[1,0,0,1],[0,0,1,1]],[2.0/sqrt(5),1.0/sqrt(5)])
			@test SpinQubits.normalize(spinor2) == SpinQubits.Spinor([[1,0,0,1],[0,1,0,1]],[im/sqrt(2),1.0/sqrt(2)])
			@test SpinQubits.normalize(spinor1 + spinor2) == SpinQubits.Spinor([[1,0,0,1],[0,0,1,1],[0,1,0,1]],[(2.0 + im)/sqrt(7), 1.0/sqrt(7), 1.0/sqrt(7)])
		end
		@testset "getInitKet" begin # Sanity check, the order of "basis" is ddd ddu dud duu udd udu uud uuu
			@test SpinQubits.getInitKet(2, true) ≈ [0.0, -1/sqrt(2), 1/sqrt(2), 0.0]
			@test SpinQubits.getInitKet(2, false) ≈ [0.0, 1.0, 0.0, 0.0]
			@test SpinQubits.getInitKet(3, true)  ≈ [0.0, 0.0, 0.0, -1/sqrt(2), 0.0, 1/sqrt(2), 0.0, 0.0]
			@test SpinQubits.getInitKet(3, false) ≈ [0.0, 0.0, 0.0,       1.0,  0.0,       0.0, 0.0, 0.0]
			@test SpinQubits.getInitKet(4, true) ≈ [0.0, 0.0, 0.0,        0.0, 0.0, 0.0, 0.0, -1/sqrt(2),
										  0.0, 0.0, 0.0, 1/sqrt(2), 0.0, 0.0, 0.0,       0.0]
			@test SpinQubits.getInitKet(4, false) ≈ [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
										  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
		end
	end

	@testset "tensors.jl" begin
		@testset "J tensors" begin
			@test sum(SpinQubits.getjtensor(2, 0.0),dims=3) ≈ [1.0  0.0  0.0 0.0;
																	0.0 -1.0  2.0 0.0;
																	0.0  2.0 -1.0 0.0;
																	0.0  0.0  0.0 1.0]
			@test sum(SpinQubits.getjtensor(2, 0.1),dims=3) ≈ [1.0  0.0  0.0 0.0;
																0.0 -1.0  2.0 0.0;
																0.0  2.0 -1.0 0.0;
																0.0  0.0  0.0 1.0]
			@test sum(SpinQubits.getjtensor(3, 0.01),dims=3) ≈ [2.01 0.0	0.0	0.0	0.0	0.0	0.0	0.0;
																0.0 -0.01 2.0	0.0	0.02 0.0 0.0 0.0;
																0.0 2 -1.99 0.0 2.0	0.0	0.0	0.0;
																0.0 0.0 0.0 -0.01	0.0 2 0.02	0.0;
																0.0 0.02	2.0 0.0 -0.01 0.0 0.0 0.0;
																0.0 0.0 0.0 2.0 0.0 -1.99 2.0 0.0;
																0.0 0.0 0.0 0.02 0.0 2.0 -0.01	0.0;
																0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.01]
		end
	end
end


text:
@testset "Test GenericStopping" begin
  include("rosenbrock.jl")
  #We build a first stopping: to test the change of tol_check function
  x0 = ones(6)
  state0 = GenericState(x0)
  stop0 = GenericStopping(
    rosenbrock,
    state0,
    tol_check = (atol, rtol, opt0) -> atol + rtol * opt0,
    list = ListofStates(state0),
  )

  io = IOBuffer()
  show(io, stop0)
  stop0.listofstates[1]

  meta = StoppingMeta(tol_check = (atol, rtol, opt0) -> atol + rtol * opt0)
  stop0_meta = GenericStopping(rosenbrock, meta, state0, list = ListofStates(state0))
  stop0_src = GenericStopping(rosenbrock, meta, cheap_stop_remote_control(), state0)
  stop0_src_without_meta = GenericStopping(
    rosenbrock,
    StopRemoteControl(),
    state0,
    tol_check = (atol, rtol, opt0) -> atol + rtol * opt0,
    list = ListofStates(state0),
  )
  @test isnan(elapsed_time(stop0))
  @test start!(stop0) == false
  @test status(stop0) == :Unknown
  #We now illustrate the impact of the choice of the norm for the unboundedness of the iterate
  stop0.meta.unbounded_x = sqrt(6)
  stop!(stop0)
  @test status(stop0, list = true) == [:Unknown] #ok as ||x||_\infty = 1 < sqrt(6)
  @test !isnan(elapsed_time(stop0))

  #We now test that stop! verifies that:
  #- there are no NaN in the score
  #- if the listofstates != nothing, stop! increases the list of states with the current_state.
  stop0.meta.optimality_check = (a, b) -> NaN
  #stop0.listofstates = ListofStates(state0)
  stop!(stop0)
  @test :DomainError in status(stop0, list = true)
  @test length(stop0.listofstates) == 4

  #Initialize a GenericStopping by default
  stop_def = GenericStopping(rosenbrock, x0, atol = 1.0)
  @test stop_def.current_state.x == x0
  @test stop_def.meta.atol == 1.0
  @test start!(stop_def) == false

  #We build a first stopping:
  x0 = ones(6)
  state = GenericState(x0)
  stop = GenericStopping(rosenbrock, state, max_time = 2.0, rtol = 0.0)
  #If rtol != 0, any point is a solution as optimality0 = Inf.

  io = IOBuffer()
  show(io, stop)

  @test start!(stop) == false
  @test stop.meta.start_time != NaN
  @test stop!(stop) == false

  #We build a substopping:
  x1 = zeros(6)
  state1 = GenericState(x1)
  # to avoid the stop by counting stop calls
  substop = GenericStopping(
    rosenbrock,
    state1,
    stop_remote = StopRemoteControl(),
    main_stp = stop,
    max_iter = typemax(Int),
    rtol = 0.0,
  )
  # substop.stop_remote = StopRemoteControl()
  #If rtol != 0, any point is a solution as optimality0 = Inf.

  io = IOBuffer()
  show(io, substop)

  @test start!(substop) == false
  @test stop!(substop) == false

  function infinite_algorithm(stp::AbstractStopping)
    x0 = stp.current_state.x
    smallest_f = stp.pb(x0) #stp.pb is a function here

    if !(typeof(stp.main_stp) <: VoidStopping)
      start!(stp.main_stp) #start the time counter of the upper loop
    end
    OK = start!(stp)

    while !OK
      x = 10 * rand(length(x0))
      smallest_f = min(smallest_f, stp.pb(x))

      OK = update_and_stop!(stp, x = x)
    end
    return smallest_f
  end

  #Solve the problem with substopping
  res = infinite_algorithm(substop)

  #The algorithm stopped because the time limit of stop is attained
  @test status(substop) == :ResourcesOfMainProblemExhausted
  @test status(stop) == :TimeLimit
  @test substop.current_state.x != x1

  #Reinitialize the stop
  #test the reinit!
  timesave, opt0save = stop.meta.start_time, stop.meta.optimality0
  reinit!(stop)
  @test timesave != stop.meta.start_time
  #@test opt0save != stop.meta.optimality0

  #Solve again the problem
  res2 = infinite_algorithm(stop)

  #The algorithm stopped as it attained the iteration limit (stop! call)
  @test status(stop) == :IterationLimit

  reinit!(stop)
  reinit!(substop)

  #
  # Test the triple sub-Stopping now:
  #
  subsubstop =
    GenericStopping(rosenbrock, state1, main_stp = substop, max_iter = typemax(Int), rtol = 0.0)
  #If rtol != 0, any point is a solution as optimality0 = Inf.

  #Solve again the problem
  start!(stop) #initialize here as infinite_algorithm has 2 "loops" only
  res3 = infinite_algorithm(subsubstop)

  @test status(stop) == :TimeLimit #stop because of the main main stopping.
  @test status(substop) == :ResourcesOfMainProblemExhausted
  @test status(subsubstop) == :ResourcesOfMainProblemExhausted

  stop.meta.infeasible = true
  @test status(stop, list = true) == [:TimeLimit, :Infeasible]

  try
    fill_in!(stop, zeros(5))
    @test false
  catch
    @test true
  end

  #Test with the cheap remote
  stop = GenericStopping(nothing, cheap_stop_remote_control(), GenericState(rand(2)), rtol = 0.0)
  @test stop.stop_remote.cheap_check
  OK = update_and_start!(stop, x = ones(2))
  @test !OK
  OK = update_and_stop!(stop, x = zeros(2))
  @test !OK

  #Test _inequality_check
  @test Stopping._inequality_check(0.0, 1.0, -1.0) == true
  @test Stopping._inequality_check(0.0, 1.0, 1.0) == false
  @test Stopping._inequality_check(Inf, Inf, -1.0) == true
  @test Stopping._inequality_check(NaN, NaN, -1.0) == false
  #works with tuple:
  @test Stopping._inequality_check((0.0, 1.0), 1.0, -1.0) == true
  @test Stopping._inequality_check((0.0, 1.5), 1.0, -1.0) == false
  @test Stopping._inequality_check(-ones(4), 1.0, -1.0) == true
  @test Stopping._inequality_check(vcat(0.0, 2 * ones(4)), 1.0, -1.0) == false
  #same types
  @test Stopping._inequality_check((0.0, 1.0), (1.0, 1.0), (-1.0, -1.0)) == true
  @test_throws ErrorException Stopping._inequality_check(zeros(2), ones(3), -ones(2))
end


text:
using AlphaZero
using ProgressMeter
using Statistics: mean

gspec = Examples.games["tictactoe"]
baseline = MinMax.Player(depth=5, amplify_rewards=true, τ=0.2)
mcts = MCTS.Env(gspec, MCTS.RolloutOracle(gspec))
mcts = MctsPlayer(mcts, niters=1000, τ=ConstSchedule(0.5))
player = TwoPlayers(mcts, baseline)

num_games = 200
bar = Progress(num_games)
rewards = map(1:num_games) do _
  trace = play_game(gspec, player)
  next!(bar)
  return total_reward(trace)
end

println("Average reward: $(mean(rewards))")


text:
using Plots
plotlyjs()

# function ploteye(ellipsoids::Dict{String,Ellipsoid}, l::Light, opticunits::Vector{OpticUnit}; n_rays::Int = 15)
function ploteye(;n_rays::Int = 15)
    include("main.jl")
    x, y = coordinates2d(ellipsoids)
    plot(x, y, aspect_ratio = :equal, leg=false, color = :black, reuse=false)
    r = Ray()
    for b in linspace(0, 1, n_rays), b2 in [0, 0.5]
        getRay!(r, l, b, b2)
        x = Float64[]
        y = Float64[]
        push!(x, r.orig[1])
        push!(y, r.orig[3])
        complete = true
        for i in opticunits
            if raytrace!(r, i)
                complete = false
                break
            end
            push!(x, r.orig[1])
            push!(y, r.orig[3])
        end
        complete && plot!(x, y, color = :grey)
    end
    plot!(ylims = (-500,100), xlims = (-300,300), xlabel = "μm", ylabel = "μm")
    gui()
end

# function plotpixels(opticunits::Vector{OpticUnit}, l::Light; n_rays::Int = 100000)
function plotpixels(;n_rays::Int = 100000)
    include("main.jl")
    psf = getpsf(opticunits, l, n_rays)
    ph = []
    for (k, v) in psf
        push!(ph, heatmap(collect(v), title = k, aspect_ratio = :equal, reuse = false))
    end
    plot(reuse = false, ph...)
    gui()
end

# function plotfwhm_aperture(opticunits::Vector{OpticUnit}, source_distance::Number, min_aperture::Number, max_aperture::Number; n_data::Int = 10, n_rays::Int = 10000, angle_step = 0.25°)
function plotfwhm_aperture(;n_data::Int = 10, n_rays::Int = 10000, angle_step = 0.25°)
    include("main.jl")
    apertures = linspace(min_aperture, max_aperture, n_data)
    fwhms = Dict(k => zeros(n_data) for k in retinas)
    for (i, aperture) in enumerate(apertures), (k, v) in getfwhm(opticunits, source_distance, aperture, n_rays, format(angle_step))
        fwhms[k][i] = get(v)
    end
    plot(reuse=false)
    [plot!(apertures, v, label = k) for (k, v) in fwhms]
    plot!(xlabel = "Aperture (μm)", ylabel = "FWHM (°)")
    gui()
end

# function plotfwhm_distance(opticunits::Vector{OpticUnit}, aperture::Number, min_distance::Number, max_distance::Number; n_data::Int = 10, n_rays::Int = 10000, angle_step = 0.25°, u = cm)
function plotfwhm_distance(;n_data::Int = 10, n_rays::Int = 10000, angle_step = 0.25°, u = cm)
    include("main.jl")
    viewdistances = logspace(log10(min_distance), log10(max_distance), n_data)
    fwhms = Dict(k => zeros(n_data) for k in retinas)
    for (i, viewdistance) in enumerate(viewdistances), (k, v) in getfwhm(opticunits, viewdistance, aperture, n_rays, format(angle_step))
        fwhms[k][i] = get(v)
    end
    plot(xscale = :log10, xlabel = "Source distance ($u)", ylabel = "FWHM (°)", reuse=false)
    [plot!(ustrip.(uconvert.(u, μm.*viewdistances)), v, label = k) for (k, v) in fwhms]
    gui()
end

# function plotfwhm_aperture_distance(opticunits::Vector{OpticUnit}, min_aperture::Number, max_aperture::Number, min_distance::Number, max_distance::Number; n_data::Int = 10, n_rays::Int = 10000, angle_step = 0.25°)
function plotfwhm_aperture_distance(;n_data::Int = 10, n_rays::Int = 10000, angle_step = 0.25°)
    include("main.jl")
    viewdistances = logspace(log10(min_distance), log10(max_distance), n_data)
    apertures = linspace(min_aperture, max_aperture, n_data)
    fwhms = Dict(k => zeros(n_data,n_data) for k in retinas)
    for (i, viewdistance) in enumerate(viewdistances), (j, aperture) in enumerate(apertures), (k, v) in getfwhm(opticunits, viewdistance, aperture, n_rays, format(angle_step))
        fwhms[k][i,j] = get(v)
    end
    ph = []
    for (k, v) in fwhms
        push!(ph, heatmap(apertures, ustrip.(uconvert.(cm, μm.*viewdistances)), v, title = k, xlabel = "Aperture (μm)", ylabel = "Source distance (cm)", reuse = false))
    end
    plot(ph..., yscale = :log10, zlims = (0,20), reuse=false)
    gui()
end

#=l = getLight(1e9, aperture, ellipsoids["cornea_distal"], .05)
psf = get2dpsf(opticunits, l, 100000)
plot()
for (k, v) in psf
    plot!(collect(v), label = k)
end
gui()=#

#=nt = 20
thetas = linspace(0, .5, nt)
ratios = Dict(k => zeros(nt) for k in retinas)
for (i, θ) in enumerate(thetas)
    l = getLight(1e9, aperture, ellipsoids["cornea_distal"], θ)
    psf = get2dpsf(opticunits, l, 100000)
    for k in keys(ratios)
        ratios[k][i] = maximum(psf[k])/psf[k][0]
    end
end
plot(ylims = (0,5))
for (k, v) in ratios
    plot!(thetas, v, label = k)
end
gui()=#


text:
using Dice
using Dice: num_flips, num_nodes, to_dice_ir

code = @dice begin
    x = flip(0.5)
    y = if x x else false
    if y y else false
end

# BDD analysis
bdd = compile(code)
num_flips(bdd), num_nodes(bdd)
infer(code, :bdd)

# IR analysis
println(to_dice_ir(code))
has_dice_binary() && rundice(code)
has_dice_binary() && infer(code, :ocaml)


text:
abstract type QueryStrategy end

function qs_score end

"""
get_query_objects(qs::QueryStrategy, data::Array{T, 2}, pools::Vector{Symbol}, global_indices::Vector{Int}, history::Vector{Int})

# Arguments
- `query_data`: Subset of the full data set
- `pools`: Labels for `query_data`
- `global_indices`: Indices of the observations in `query_data` relative to the full data set.
- `history`: Indices of previous queries
"""
function get_query_objects end

function initialize_qs(qs::DataType, model::SVDD.OCClassifier, data::Array{T, 2}, params::Dict{Symbol, <:Any})::QueryStrategy where T <: Real
    if qs <: BatchPQs
        if haskey(params, :SequentialStrategy)
            sequential_strategy_type = eval(params[:SequentialStrategy][:type])
            sequential_strategy_params = params[:SequentialStrategy][:param]
            sequential_strategy = initialize_qs(sequential_strategy_type, model, data, sequential_strategy_params)
            # use all params except SequentialStrategy for batch query initialization
            batch_params = filter(tuple -> first(tuple) != :SequentialStrategy, params)
            if qs <: ExtendingBatchQs
                return qs(sequential_strategy; batch_params...)
            else
                return qs(model, sequential_strategy; batch_params...)
            end
        else
            return qs(; params...)
        end
    elseif qs <: SequentialPQs
        if qs <: HybridPQs
            return qs(model, data; params...)
        elseif qs <: ModelBasedPQs
            return qs(model; params...)
        elseif qs <: DataBasedPQs
            kernel = SVDD.get_kernel(model)
            if typeof(kernel) == MLKernels.GaussianKernel
                return qs(data, bw_method=kernel.alpha.value.x; params...)
            else
                return qs(data; params...)
            end
        else
            return qs(; params...)
        end
    elseif qs <: QuerySynthesisStrategy
        if qs <: HybridQss
            return qs(model, data; params...)
        elseif qs <: ModelBasedQss
            return qs(model; params...)
        elseif qs <: DataBasedQss
            kernel = get_kernel(model)
            if typeof(kernel) == MLKernels.GaussianKernel
                return qs(data, bw_method=kernel.alpha.value.x; params...)
            else
                return qs(data; params...)
            end
        else
            return qs(; params...)
        end
    elseif qs <: SubspaceQueryStrategy
        return qs(model, data; params...)
    elseif qs <: QueryStrategy
        return qs(; params...)
    end
    throw(ErrorException("Unknown query strategy of type $(qs)."))
end


text:
using Empirikos
using Test

_target = MarginalDensity(StandardNormalSample(3.0))
@test Base.extrema(_target)[2] == pdf(Normal(),0)

@test Base.extrema(MarginalDensity(NormalSample(3.0, 0.5)))[2] == pdf(Normal(0,0.5),0.0)


text:
mutable struct Context <: AContext
  cache::Union{IdDict{Any,Any},Nothing}
  globals::Union{Dict{GlobalRef,Any},Nothing}
end

Context() = Context(nothing, nothing)

cache(cx::Context) = cx.cache === nothing ? (cx.cache = IdDict()) : cx.cache
globals(cx::Context) = cx.globals === nothing ? (cx.globals = Dict{GlobalRef,Any}()) : cx.globals

struct Pullback{S,T}
  t::T
end

Pullback{S}(x) where S = Pullback{S,typeof(x)}(x)

struct CompileError
  T
  e
end

function Base.showerror(io::IO, e::CompileError)
  print(io, "Compiling $(e.T): ")
  showerror(io, e.e)
end

# interface2.jl

# Wrappers

_forward(f, args...) = _forward(Context(), f, args...)

tailmemaybe(::Nothing) = nothing
tailmemaybe(x::Tuple) = Base.tail(x)

function forward(f, args...)
  y, back = _forward(f, args...)
  y, Δ -> tailmemaybe(back(Δ))
end

sensitivity(y::Number) = one(y)
sensitivity(y::Vector) = _vec(y)
sensitivity(y::Complex) = error("Output is complex, so the gradient is not defined.")
sensitivity(y) = error("Output should be scalar; gradients are not defined for output $(repr(y))")

function gradient(f, args...)
  y, back = forward(f, args...)
  return back(sensitivity(y))
end

Base.adjoint(f::Function) = x -> gradient(f, x)[1]

# Param-style wrappers

# TODO store ids only
struct Params
  order::Vector{Any}
  params::IdSet{Any}
  Params() = new([], IdSet())
end

@forward Params.order Base.iterate, Base.length

function Base.push!(ps::Params, x)
  if !(x in ps.params)
    push!(ps.order, x)
    push!(ps.params, x)
  end
  return ps
end

Base.push!(ps::Params, x...) = (foreach(x -> push!(ps, x), x); ps)

Params(xs) = push!(Params(), xs...)

function Base.show(io::IO, ps::Params)
  print(io, "Params([")
  join(io, ps.order, ", ")
  print(io, "])")
end

struct Grads
  grads::IdDict{Any,Any}
end

Base.show(io::IO, ps::Grads) = print(io, "Grads(...)")

@forward Grads.grads Base.getindex, Base.haskey

function Base.getindex(gs::Grads, x)
  isbits(x) && error("Only reference types can be differentiated with `Params`.")
  return gs.grads[x]
end

function forward(f, ps::Params)
  cx = Context()
  y, back = _forward(cx, f)
  y, function (Δ)
    for p in ps
      cache(cx)[p] = nothing
    end
    back(Δ)
    Grads(cx.cache) # TODO make a copy
  end
end

# Code Reflection

using InteractiveUtils
using InteractiveUtils: typesof
using Core: Typeof

function code_ir(f, T)
  m = meta(Tuple{Typeof(f),T.parameters...})
  return IR(m)
end

function code_irm(ex)
  isexpr(ex, :call) || error("@code_ir f(args...)")
  f, args = ex.args[1], ex.args[2:end]
  :(code_ir($(esc(f)), typesof($(esc.(args)...))))
end

macro code_ir(ex)
  code_irm(ex)
end

macro code_adjoint(ex)
  :(Adjoint($(code_irm(ex)), varargs = varargs($(esc(:($InteractiveUtils.@which $ex))), length(($(esc.(ex.args)...),)))))
end

function _vec(f)
  println("Privet")
  println("---------------------------------")
  println("F: ", f)
  N = length(f)
  res = zeros(N)
  res[1] = f[2] - f[1]
  res[end] = f[end] - f[end-1]
  for i=2:N-1
    #Central differences
    res[i] = (f[i+1] - f[i-1]) / 2
  end
  println("Result: ", res)
  println("---------------------------------")
  return res
end


text:
using JLD, Serialization
using BenchmarkTools
using TypedTables

res = if endswith(ARGS[1], ".jld")
    JLD.load(ARGS[1])
elseif endswith(ARGS[1], ".jls")
    deserialize(ARGS[1])
else
    error("Unknown file type")
end

serial_results = filter(x->!occursin("dagger", x[1]), res["results"])
@assert length(keys(serial_results)) > 0 "No serial results found"
dagger_results = filter(x->occursin("dagger", x[1]), res["results"])
@assert length(keys(dagger_results)) > 0 "No Dagger results found"

scale_set = sort([key=>parse(Int, lstrip(last(split(key, ':')), ' ')) for key in keys(first(serial_results)[2])]; by=x->x[2])
nw_set = sort([key=>parse(Int, lstrip(last(split(key, ':')), ' ')) for key in keys(first(dagger_results)[2][first(first(scale_set))])]; by=x->x[2])
raw_table = NamedTuple[]
for bset_key in keys(res["results"])
    bset = res["results"][bset_key]
    if typeof(bset[first(first(scale_set))]) <: BenchmarkGroup
        procs = parse(Int, lstrip(last(split(first(first(bset[first(first(scale_set))])), ':')), ' '))
        for nw in nw_set
            for i in 1:length(scale_set)
                set_times = [minimum(bset[scale][nw[1]]).time/(10^9) for scale in first.(scale_set)]
                push!(raw_table, (name=bset_key, time=set_times[i], scale=last.(scale_set)[i], procs=nw[2]))
            end
        end
    else
        set_times = [minimum(bset[scale]).time/(10^9) for scale in first.(scale_set)]
        procs = 8  # default for OpenBLAS
        for i in 1:length(set_times)
            push!(raw_table, (name=bset_key, time=set_times[i], scale=last.(scale_set)[i], procs=procs))
        end
    end
end
table = Table(raw_table)

btable = copy(table[map(x->!x, occursin.(Ref("dagger"), table.name))])
dtable = copy(table[occursin.(Ref("dagger"), table.name)])

#table = Table(name=[:Base for _ in 1:3], time=serial_times, scale=last.(scale_set), procs=[8 for _ in 1:3])

#btable = copy(table)

#=
for (nw,nw_val) in nw_set
    dagger_times = [minimum(dagger_results[scale][nw]).time/(10^9) for scale in first.(scale_set)]
    t = Table(name=[:Dagger for _ in 1:3], time=dagger_times, scale=last.(scale_set), procs=[parse(Int,split(nw, ":")[2]) for _ in 1:3])
    append!(table, t)
end
=#

#dtable = table[table.name .== :Dagger]

# Plotting

using CairoMakie
using Colors, ColorSchemes

# 1. Strong scaling analysis
#   - Keep the problem size constant
#   - Increase amount of compute

fig = Figure(resolution = (1200, 800))
strong_scaling(t1, tn, N) = t1/(N*tn)

ssp = fig[1, 1] = Axis(fig, title = "Strong Scaling")
ssp.xlabel = "Number of processes"
ssp.ylabel = "Scaling efficiency"

line_plots = Any[]
legend_names = String[]

scales = unique(dtable.scale)

colors = distinguishable_colors(length(scales), ColorSchemes.seaborn_deep.colors)

for (i, scale) in enumerate(scales)
    stable = dtable[dtable.scale .== scale]
    t1 = first(stable[stable.procs .== minimum(dtable.procs)].time)
    ss_efficiency = strong_scaling.(t1, stable.time, stable.procs)
    push!(line_plots, lines!(ssp, stable.procs, ss_efficiency, linewidth=3.0, color = colors[i]))
    push!(legend_names, "scale = $scale")
end

legend = fig[1, 2] = Legend(fig, line_plots, legend_names)

save("strong_scaling.png", fig)

# 2. Weak scaling
#   - Increase work and available compute
# Q: Are we actually scaling the work correctly?
#    too little data

fig = Figure(resolution = (1200, 800))
weak_scaling(t1, tn, p_prime, p) = t1/((p_prime/p)*tn)

t1 = first(dtable[map(row->(row.scale == 10) && (row.procs == 1), dtable)]).time

fig = Figure(resolution = (1200, 800))
perf = fig[1, 1] = Axis(fig, title = "Weak Scaling")
perf.xlabel = "Number of processes"
perf.ylabel = "Scaling efficiency"

line_plots = Any[]
legend_names = String[]

wstable = similar(dtable, 0)
for pair in [(10,1),(35,4),(85,8)]
    append!(wstable, dtable[map(row->(row.scale == pair[1]) && (row.procs == pair[2]), rows(dtable))])
end
push!(line_plots, lines!(perf, wstable.procs, weak_scaling.(t1, wstable.time, wstable.procs .* 10, wstable.scale), linewidth=3.0))
push!(legend_names, "cpu+dagger")

legend = fig[1, 2] = Legend(fig, line_plots, legend_names)
save("weak_scaling.png", fig)

# 3. Comparision against Base

fig = Figure(resolution = (1200, 800))
perf = fig[1, 1] = Axis(fig, title = "Dagger vs Base")
perf.xlabel = "Scaling factor"
perf.ylabel = "time (s)"

line_plots = Any[]
legend_names = String[]

procs = unique(dtable.procs)

colors = distinguishable_colors(length(procs) + 1, ColorSchemes.seaborn_deep.colors)

for (i, nproc) in enumerate(procs)
    stable = dtable[dtable.procs .== nproc]
    push!(line_plots, lines!(perf, stable.scale, stable.time, linewidth=3.0, color = colors[i]))
    push!(legend_names, "Dagger (nprocs = $nproc)")
end

push!(line_plots, lines!(perf, btable.scale, btable.time, linewidth=3.0, color = colors[end]))
push!(legend_names, "Base (threads = 8)")

legend = fig[1, 2] = Legend(fig, line_plots, legend_names)
save("raw_timings.png", fig)

# 4. Speedup
fig = Figure(resolution = (1200, 800))
speedup = fig[1, 1] = Axis(fig, title = "Speedup vs. 1 processor")
speedup.xlabel = "Number of processors"
speedup.ylabel = "Speedup"

line_plots = Any[]
legend_names = String[]

colors = distinguishable_colors(length(procs), ColorSchemes.seaborn_deep.colors)

t1 = sort(dtable[dtable.scale .== 10]; by=r->r.procs)

for (i, scale) in enumerate(unique(dtable.scale))
    stable = dtable[dtable.scale .== scale]
    sort!(stable, by=r->r.procs)
    push!(line_plots, lines!(speedup, stable.procs, stable.time ./ t1.time, linewidth=3.0, color = colors[i]))
    push!(legend_names, "Dagger (scale = $scale)")
end

legend = fig[1, 2] = Legend(fig, line_plots, legend_names)
save("speedup.png", fig)


text:
# included from Potentials.jl
# part of the module JuLIP.Potentials

using JuLIP: JVec, JMat, neighbourlist
using LinearAlgebra: I
using JuLIP.Chemistry: atomic_number
using NeighbourLists

export ZeroPairPotential, PairSitePotential, ZBLPotential,
         LennardJones, lennardjones,
         Morse, morse

## TODO: kill this one?
grad(V::PairPotential, r::Real, R::JVec) = ((@D V(r)) / r) * R

evaluate!(tmp, V::PairPotential, r::Union{Number, JVec}) = V(r)
evaluate_d!(tmp, V::PairPotential, r::Union{Number, JVec}) = @D V(r)
evaluate_dd!(tmp, V::PairPotential, r::Union{Number, JVec}) = @DD V(r)

function evaluate!(tmp, V::PairPotential, R::AbstractVector{JVec{T}}) where {T}
   Es = zero(T)
   for i = 1:length(R)
      Es += T(0.5) * evaluate!(tmp, V, norm(R[i]))
   end
   return Es
end

function evaluate_d!(dEs, tmp, V::PairPotential, R::AbstractVector{JVec{T}}) where {T}
   for i = 1:length(R)
      r = norm(R[i])
      dEs[i] = (T(0.5) * evaluate_d!(tmp, V, r) / r) * R[i]
   end
   return dEs
end

function evaluate_dd!(hEs, tmp, V::PairPotential, R::AbstractVector{<: JVec})
   n = length(R)
   for i = 1:n
      hEs[i,i] = 0.5 * _hess!(tmp, V, norm(R[i]), R[i])
   end
   return hEs
end

function _hess!(tmp, V::PairPotential, r::Number, R::JVec)
   R̂ = R/r
   P = R̂ * R̂'
   dV = evaluate_d!(tmp, V, r) / r
   ddV = evaluate_dd!(tmp, V, r)
   return (ddV - dV) * P + dV * I
end

function precon!(hEs, tmp, V::PairPotential, R::AbstractVector{<: JVec}, innerstab=T(0.0))
   n = length(R)
   for i = 1:n
      hEs[i,i] = precon!(tmp, V, norm(R[i]), R[i], innerstab)
   end
   return hEs
end

# an FF preconditioner for pair potentials
function precon!(tmp, V::PairPotential, r::T, R::JVec{T}, innerstab=T(0.1)
                 ) where {T <: Number}
   r = norm(R)
   dV = evaluate_d!(tmp, V, r)
   ddV = evaluate_dd!(tmp, V, r)
   R̂ = R/r
   return (1-innerstab) * (abs(ddV) * R̂ * R̂' + abs(dV / r) * (I - R̂ * R̂')) +
             innerstab  * (abs(ddV) + abs(dV / r)) * I
end

"""
`LennardJones(σ, e0):` constructs the 6-12 Lennard-Jones potential [wiki](https://en.wikipedia.org/wiki/Lennard-Jones_potential)

   e0 * 4 * ( (σ/r)¹² - (σ/r)⁶ )

Constructor with kw arguments: `LennardJones(; kwargs...)`

* `e0, σ` : standard LJ parameters, default `e0 = 1.0, σ = 1.0`
* `r0` : equilibrium distance - if `r0` is specified then `σ` is ignored
* `a0` : FCC lattice parameter, if `a0` is specified then `σ` is ignored

(`r0, a0` cannot both be specified at the same time)
"""
LennardJones(σ, e0) = (@analytic r -> e0 * 4.0 * ((σ/r)^(12) - (σ/r)^(6)))

function ljparams(; σ=1.0, e0=1.0, r0 = nothing, a0 = nothing)
   if r0 != nothing && a0 != nothing
      error("`LenndardJones`: cannot specify both `r0` and `a0`")
   end
   if a0 != nothing    # r0 = nn-dist = a0/sqrt(2) in FCC
      r0 = a0 / sqrt(2)
   end
   if r0 != nothing    # standard LJ is minimised at r = 2^(1/6)
      σ =  r0 / 2^(1/6)
   end
   return σ, e0
end

LennardJones(; kwargs...) = LennardJones(ljparams(;kwargs...)...)

"""
`lennardjones(; kwargs...)`

simplified constructor for `LennardJones` (note this is type unstable!)

In addition to the `kwargs` of `LennardJones`, this accepts also

* `rcut` : default `:auto` which gives `rcut = (1.9*σ, 2.7*σ)`. Use
`nothing` or `Inf` to specify no cutoff, or specify a tuple or
array with two elements specifying the lower and upper cut-off radii to be
used with `SplineCutoff`.
"""
function lennardjones(; rcut = :auto, kwargs...)
   σ, e0 = ljparams(; kwargs...)
   if (rcut == nothing || rcut == Inf)
      return LennardJones(σ, e0)
   elseif rcut == :auto
      rcut = (1.9*σ, 2.7*σ)
   end
   return SplineCutoff(rcut[1], rcut[2]) * LennardJones(σ, e0)
end

"""
`Morse(A, e0, r0)` or `Morse(;A=4.0, e0=1.0, r0=1.0)`: constructs a
`PairPotential` for
```
   e0 ( exp( -2 A (r/r0 - 1) ) - 2 exp( - A (r/r0 - 1) ) )
```
"""
Morse(A, e0, r0) = @analytic(
   r -> e0 * ( exp(-(2.0*A) * (r/r0 - 1.0)) - 2.0 * exp(-A * (r/r0 - 1.0)) ) )
Morse(;A=4.0, e0=1.0, r0=1.0) = Morse(A, e0, r0)

"""
`morse(A=4.0, e0=1.0, r0=1.0, rcut=(1.9*r0, 2.7*r0))`

simplified constructor for `Morse` (type unstable)
"""
morse(;A=4.0, e0=1.0, r0=1.0, rcut=(1.9*r0, 2.7*r0)) = (
   (rcut == nothing || rcut == Inf)
         ?  Morse(A, e0, r0)
         :  SplineCutoff(rcut[1], rcut[2]) * Morse(A, e0, r0) )

"""
`ZeroPairPotential()`: creates a potential that just returns zero
"""
struct ZeroPairPotential <: PairPotential end

@pot ZeroPairPotential

evaluate(p::ZeroPairPotential, r::T) where {T <: Number} = T(0.0)
evaluate_d(p::ZeroPairPotential, r::T) where {T <: Number} = T(0.0)
evaluate_dd(p::ZeroPairPotential, r::T) where {T <: Number} = T(0.0)
cutoff(p::ZeroPairPotential) = Bool(0) # the weakest number type

# ------------------------------------------------------------------------

# TODO: write more docs + tests for ZBL

"""
Implementation of the ZBL potential to model close approach.
"""
struct ZBLPotential{TV} <: PairPotential
   Z1::Int
   Z2::Int
   V::TV   # analytic
end

@pot ZBLPotential

evaluate(V::ZBLPotential, r::Number) = evaluate(V.V, r::Number)
evaluate_d(V::ZBLPotential, r::Number) = evaluate_d(V.V, r::Number)
cutoff(::ZBLPotential) = Inf

ZBLPotential(Z1::Integer, Z2::Integer) =
   let Z1=Z1, Z2=Z2
      au = 0.8854 * 0.529 / (Z1^0.23 + Z2^0.23)
      ϵ0 = 0.00552634940621
      C = Z1*Z2/(4*π*ϵ0)
      E1, E2, E3, E4 = 0.1818, 0.5099, 0.2802, 0.02817
      A1, A2, A3, A4 = 3.2/au, 0.9423/au, 0.4028/au, 0.2016/au
      V = @analytic(r -> C * (E1*exp(-A1*r) + E2*exp(-A2*r) +
                              E3*exp(-A4*r) + E4*exp(-A4*r) ) / r)
      ZBLPotential(Z1, Z2, V)
   end

ZBLPotential(Z::Integer) = ZBLPotential(Z, Z)
ZBLPotential(s1::Symbol, s2::Symbol) = ZBLPotential(atomic_number(s1), atomic_number(s2))
ZBLPotential(s::Symbol) = ZBLPotential(s, s)

Dict(V::ZBLPotential) = Dict("__id__" => "JuLIP_ZBLPotential",
                             "Z1" => V.Z1,
                             "Z2" => Z.Z2)
ZBLPotential(D::Dict) = ZBLPotential(D["Z1"], D["Z2"])
Base.convert(::Val{:JuLIP_ZBLPotential}, D::Dict) = ZBLPotential(D)

# ====================================================================
#   A product of two pair potentials: primarily used for cutoff mechanisms

"product of two pair potentials"
mutable struct ProdPot{P1, P2} <: PairPotential
   p1::P1
   p2::P2
end

@pot ProdPot

import Base.*
*(p1::PairPotential, p2::PairPotential) = ProdPot(p1, p2)
@inline evaluate(p::ProdPot, r::Number) = p.p1(r) * p.p2(r)
evaluate_d(p::ProdPot, r::Number) = (p.p1(r) * (@D p.p2(r)) + (@D p.p1(r)) * p.p2(r))
evaluate_dd(p::ProdPot, r::Number) = (p.p1(r) * (@DD p.p2(r)) +
              2 * (@D p.p1(r)) * (@D p.p2(r)) + (@DD p.p1(r)) * p.p2(r))
cutoff(p::ProdPot) = min(cutoff(p.p1), cutoff(p.p2))

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

"""
`struct WrappedPairPotential`

wraps a pairpotential using `FunctionWrappers` in order to allow
type-stable storage of multiple potentials. This is the main technique
required at the moment to work with multi-component systems.
Otherwise, this is not advisable since it disables a range of
possible compiler optimisations.
"""
struct WrappedPairPotential <: SimplePairPotential
   f::F64fun
   f_d::F64fun
   f_dd::F64fun
   rcut::Float64
end

@pot WrappedPairPotential

cutoff(V::WrappedPairPotential) = V.rcut
# evaluate, etc are all derived from SimplePairPotential

function WrappedPairPotential(V::AnalyticFunction, rcut)
   @assert (0 < rcut < Inf)
   f, f_d, f_dd = let V=V, rc = rcut
      (F64fun(r -> evaluate(V, r) * (r<rc)),
              F64fun(r -> evaluate_d(V, r) * (r<rc)),
              F64fun(r -> evaluate_dd(V, r) * (r<rc)))
   end
   return WrappedPairPotential(f, f_d, f_dd, cutoff(V))
end

function WrappedPairPotential(V::PairPotential)
   @assert (0 < cutoff(V) < Inf)
   f, f_d, f_dd = let V=V, rc = cutoff(V)
      (F64fun(r -> evaluate(V, r) * (r<rc)),
              F64fun(r -> evaluate_d(V, r) * (r<rc)),
              F64fun(r -> evaluate_dd(V, r) * (r<rc)))
   end
   return WrappedPairPotential(f, f_d, f_dd, cutoff(V))
end


text:
using Luna
import Test: @test, @testset
import HCubature: hquadrature
import SpecialFunctions: besselj
import GSL: sf_bessel_zero_Jnu

@testset "On-axis intensity statistics" begin
    # Manually normalise the field distribution of HE11 to find scaling factor between
    # power and intensity
    a = 13e-6
    unm = sf_bessel_zero_Jnu(0, 1)
    E(r) = besselj(0, unm*r/a)
    norm, err = hquadrature(0, a) do r
        2π*r * abs2(E(r))
    end

    energy = 1e-6
    gas = :Ar
    pres = 5
    τ = 30e-15
    λ0 = 800e-9
    grid = Grid.RealGrid(15e-2, 800e-9, (160e-9, 3000e-9), 1e-12)
    m = Capillary.MarcatiliMode(a, gas, pres, loss=false)
    aeff(z) = Modes.Aeff(m, z=z)
    energyfun, energyfunω = Fields.energyfuncs(grid)

    dens0 = PhysData.density(gas, pres)
    densityfun(z) = dens0
    responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
    linop, βfun!, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0)
    inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6)
    Eω, transform, FT = Luna.setup(
        grid, densityfun, responses, inputs, βfun!, aeff)

    statsfun = Stats.default(grid, Eω, m, linop, transform; gas=gas, onaxis=true)
    output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
    Luna.run(Eω, grid, linop, transform, FT, output, status_period=5)

    @test all(output["stats"]["peakintensity"] .≈ output["stats"]["peakpower"]/norm)
end


text:
using Revise
using Flipping
##
Directory_path = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/";
Exp_type = "Photometry";
Exp_name = "GCaMP_close_pokes"
training_days = Date(2017,5,15):Day(1):Date(2017,5,26)
issues_days = Date(2016,12,31):Day(1):Date(2017,1,1)
bad_days = vcat(training_days, issues_days)
Camera_path = "/Users/dariosarra/Google Drive/Flipping/run_task_photo/DatcreGCaMPAutomaticBarr/Cam/"
Behavior_path = "/Users/dariosarra/Google Drive/Flipping/run_task_photo/raw_data"
Mice_suffix = "DG";
DataIndex = ProcessPhotometry.create_photometry_DataIndex(Directory_path, Exp_type,Exp_name, Mice_suffix;bad_days = bad_days)
names(DataIndex)
#
dformat = Dates.DateFormat("yyyymmdd")
mask = DataIndex[:Day].<=Date("20190717",dformat)
# mask = DataIndex[:MouseID].== "DG2"
DataIndex = DataIndex[mask,:]
# names(DataIndex)
##
pokes, cam_dict = save_bhv_photo(DataIndex)

DataIndex[.!ispath.(DataIndex[:Log_Path]),:Log_Path]
size(DataIndex,1)
union(column(pokes,:MouseID))
union(column(pokes,:Session))
union(DataIndex[:Session])
process_pokes(DataIndex[56,:Bhv_Path])[:MouseID]

Flipping.session_info(DataIndex[56,:Bhv_Path])


text:
using NeuralVerification, LazySets, LinearAlgebra
using Test

import NeuralVerification: ReLU, Id

macro no_error(ex)
    quote
        try $(esc(ex))
            true
        catch e
            @error(e)
            false
        end
    end
end

include("identity_network.jl")
include("relu_network.jl")
include("inactive_relus.jl")
if Base.find_package("Flux") != nothing
    include("flux.jl")
end
include("complements.jl")
include("splitting.jl")
include("mnist_1000.jl")
include("fully_split.jl")
include("write_nnet_test.jl")


text:
using Base.Test
using Lora

println("    Testing BasicUnvVariableState constructors...")

v = Float64(1.21)
s = BasicUnvVariableState(v)

@test eltype(s) == Float64
@test s.value == v

println("    Testing BasicMuvVariableState constructors...")

v = Float32[1.5, 4.1]
s = BasicMuvVariableState(v)

@test eltype(s) == Float32
@test s.value == v
@test s.size == length(v)

ssize = 3
s = BasicMuvVariableState(ssize, Float16)

@test eltype(s) == Float16
@test s.size == ssize

println("    Testing BasicMavVariableState constructors...")

v = BigFloat[3.11 7.34; 9.7 6.72; 1.18 8.1]
s = BasicMavVariableState(v)

@test eltype(s) == BigFloat
@test s.value == v
@test s.size == size(v)

ssize = (3, 5)
s = BasicMavVariableState(ssize, Float16)

@test eltype(s) == Float16
@test s.size == ssize


text:
# Matern 5/2 ARD covariance Function

"""
    Mat52Ard <: MaternARD

ARD Matern 5/2 kernel (covariance)
```math
k(x,x') = σ²(1 + √5|x-x'|/L + 5|x-x'|²/(3L²))\\exp(- √5|x-x'|/L)
```
with length scale ``ℓ = (ℓ₁, ℓ₂, …)`` and signal standard deviation ``σ`` where
``L = diag(ℓ₁, ℓ₂, …)``.
"""
mutable struct Mat52Ard <: MaternARD
    "Inverse squared length scale"
    iℓ2::Vector{Float64}
    "Signal variance"
    σ2::Float64
    "Priors for kernel parameters"
    priors::Array

    """
        Mat52Ard(ll::Vector{Float64}, lσ::Float64)

    Create `Mat52Ard` with length scale `exp.(ll)` and signal standard deviation `exp(lσ)`.
    """
    Mat52Ard(ll::Vector{Float64}, lσ::Float64) = new(exp.(-2 .* ll), exp(2 * lσ), [])
end

function set_params!(mat::Mat52Ard, hyp::VecF64)
    length(hyp) == num_params(mat) || throw(ArgumentError("Mat52 kernel only has $(num_params(mat)) parameters"))
    @views @. mat.iℓ2 = exp(-2 * hyp[1:(end-1)])
    mat.σ2 = exp(2 * hyp[end])
end

get_params(mat::Mat52Ard) = [-log.(mat.iℓ2) / 2; log(mat.σ2) / 2]
get_param_names(mat::Mat52Ard) = [get_param_names(mat.iℓ2, :ll); :lσ]
num_params(mat::Mat52Ard) = length(mat.iℓ2) + 1

Statistics.cov(mat::Mat52Ard, r::Float64) =
    (s = √5 * r; mat.σ2 * (1 + s + s^2 / 3) * exp(-s))

dk_dll(mat::Mat52Ard, r::Float64, wdiffp::Float64) =
    (s = √5 * r; 5 / 3 * mat.σ2 * wdiffp *(1 + s) * exp(-s))


text:
using DataDeps

register(DataDep(
  "pyrim",
  "https://www.openml.org/d/217",
  "https://www.openml.org/data/get_csv/3654/dataset_2203_pyrim.arff",
  "a2625a60c6a144ddcf2334df9e2f8cd087f7e95908be3b167e17d3c9587282b3",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=28,
        feature_indices=1:27,
        header=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
# from https://github.com/JuliaStats/MLBase.jl/blob/master/src/crossval.jl

abstract type CrossValGenerator end

struct Kfold <: CrossValGenerator
    permseq::Vector{Int}
    k::Int
    coeff::Float64

    function Kfold(n::Int, k::Int)
        2 <= k <= n || error("DTM: (Kfold CV Sampler): The value of k must be in [2, length(a)]. You provided k=$(k)")
        new(Random.randperm(n), k, n / k)
    end
end

length(c::Kfold) = c.k

struct KfoldState
    i::Int      # the i-th of the subset
    s::Int      # starting index
    e::Int      # ending index
end

start(c::Kfold) = KfoldState(1, 1, round.(Integer, c.coeff))
next(c::Kfold, s::KfoldState) = (i = s.i + 1; (setdiff(1:length(c.permseq), c.permseq[s.s:s.e]), KfoldState(i, s.e + 1, round.(Integer, c.coeff * i))))
done(c::Kfold, s::KfoldState) = (s.i > c.k)
@inline iterate(x::Kfold) = next(x, start(x))
@inline iterate(x::Kfold, i) = done(x, i) ? nothing : next(x, i)

# Repeated random sub-sampling

struct RandomSub <: CrossValGenerator
    n::Int    # total length
    sn::Int   # length of each subset
    k::Int    # number of subsets
end

length(c::RandomSub) = c.k

start(c::RandomSub) = 1
next(c::RandomSub, s::Int) = (sort!(StatsBase.sample(1:c.n, c.sn; replace=false)), s + 1)
done(c::RandomSub, s::Int) = (s > c.k)
@inline iterate(x::RandomSub) = next(x, start(x))
@inline iterate(x::RandomSub, i) = done(x, i) ? nothing : next(x, i)

# Kfold, but disjoing sets ('negation' of Kfold sets)
struct KfoldDisjoint <: CrossValGenerator
    permseq::Vector{Int}
    k::Int
    coeff::Float64

    function KfoldDisjoint(n::Int, k::Int)
        2 <= k <= n || error("DTM: (KfoldDisjoint CV Sampler): The value of k must be in [2, length(a)]. You provided k=$(k)")
        # 2 <= k <= n || error("The value of k must be in [2, length(a)].")
        new(Random.randperm(n), k, n / k)
    end
end

length(c::KfoldDisjoint) = c.k

start(c::KfoldDisjoint) = KfoldState(1, 1, round.(Integer, c.coeff))
next(c::KfoldDisjoint, s::KfoldState) =
    (i = s.i + 1; (sort(c.permseq[s.s:s.e]), KfoldState(i, s.e + 1, round.(Integer, c.coeff * i))))
done(c::KfoldDisjoint, s::KfoldState) = (s.i > c.k)
@inline iterate(x::KfoldDisjoint) = next(x, start(x))
@inline iterate(x::KfoldDisjoint, i) = done(x, i) ? nothing : next(x, i)

#=

#testing

collect(KfoldDisjoint(10,3))
collect(RandomSub(10,4,20))
collect(Kfold(10,3))

    nn=10
    kk=3
    sampled_kfold=Kfold(nn,kk)
    for i in sampled_kfold
        @show i
        @show setdiff(1:nn,i)
    end =#


text:
include("../roblib.jl")

function horizion_borders(r::Robot)

    while isborder(r, Left)==false || isborder(r, Down)==false
        moves!(r, Down)
        moves!(r, Left)
    end

    num_steps = get_num_movements!(r, Right) # определяем размер поля
    num_borders = 0

    while isborder(r, Up) == false || isborder(r, Right) == false # пока не дойдёт до верхней границы поля
        for side in(Left, Right)
            if isborder(r, Up) == true && isborder(r, Right) == true # дошёл до верхней границы  поля
                break
            end
            num_borders += search_border(r, side, num_steps)
            if isborder(r, Up)==false
                move!(r, Up)
            end
        end
    end

    print(num_borders)
end

function search_border(r::Robot, side::HorizonSide, numsteps::Int)
    num_borders = 0
    num_steps = 0
    state = false

    for _ in 1:numsteps
        if isborder(r,side)==false
            move!(r,side)
            if isborder(r,Up)==true && state==false
                num_borders += 1
                state = true
            else
                state=false
            end
        else # обходит вертикальную перегородку
            while isborder(r,side)==true
                move!(r, Down)
                num_steps += 1
            end
            move!(r, side)
            for _ in 1:num_steps
                move!(r, Up)
            end
            num_steps=0
        end
    end

    return(num_borders)
end


text:
@testset "Basic reading" begin
    anonymous_auth = FHIRClient.AnonymousAuth()
    oauth2_auth = FHIRClient.OAuth2()
    FHIRClient.set_token!(oauth2_auth, "helloworld")
    jwt_auth = FHIRClient.JWTAuth()
    FHIRClient.set_token!(jwt_auth, "helloworld")
    username_password_auth_1 = FHIRClient.UsernamePassAuth("helloworld")
    FHIRClient.set_password!(username_password_auth_1, "helloworld")
    username_password_auth_2 = FHIRClient.UsernamePassAuth(; username = "helloworld")
    FHIRClient.set_password!(username_password_auth_2, "helloworld")
    all_auths = [
        anonymous_auth,
        oauth2_auth,
        jwt_auth,
        username_password_auth_1,
        username_password_auth_2,
    ]
    for auth in all_auths
        fhir_version = FHIRClient.R4()
        base_url = FHIRClient.BaseURL("https://hapi.fhir.org/baseR4")
        client = FHIRClient.Client(fhir_version, base_url, auth)
        @test FHIRClient.get_fhir_version(client) == fhir_version
        @test FHIRClient.get_base_url(client) == base_url
        search_request_path = "/Patient?given=Jason&family=Argonaut"
        json_response_search_results_bundle = FHIRClient.request_json(client, "GET", search_request_path)
        patient_id = json_response_search_results_bundle.entry[1].resource.id
        patient_request = "/Patient/$(patient_id)"
        patients = [
            FHIRClient.request(FHIRClient.R4Types.Patient, client, "GET", patient_request),
            FHIRClient.request(FHIRClient.R4Types.Patient, client, "GET", patient_request; body = JSON3.read("{}")),
            FHIRClient.request(FHIRClient.R4Types.Patient, client, "GET", patient_request; query = Dict{String, String}()),
            FHIRClient.request(FHIRClient.R4Types.Patient, client, "GET", patient_request; body = JSON3.read("{}"), query = Dict{String, String}()),
            # FHIRClient.request(FHIRClient.R4Types.AbstractResource, client, "GET", patient_request),
            # FHIRClient.request(FHIRClient.R4Types.AbstractResource, client, "GET", patient_request; body = JSON3.read("{}")),
            # FHIRClient.request(FHIRClient.R4Types.AbstractResource, client, "GET", patient_request; query = Dict{String, String}()),
            # FHIRClient.request(FHIRClient.R4Types.AbstractResource, client, "GET", patient_request; body = JSON3.read("{}"), query = Dict{String, String}()),
        ]
        for patient in patients
            @test patient isa FHIRClient.R4Types.AbstractResource
            @test patient isa FHIRClient.R4Types.Patient
            @test only(patient.name).use == "usual"
            @test only(patient.name).text == "Jason Argonaut"
            @test only(patient.name).family == "Argonaut"
            @test only(only(patient.name).given) == "Jason"
            @test patient.birthDate == Dates.Date("1985-08-01")
        end
        Base.shred!(auth)
        Base.shred!(client)
    end
    for i in 1:length(all_auths)
        Base.shred!(all_auths[i])
    end
end


text:
const SMALL_K = 0.075

abstract type AbstractStructureFactor <: Function end

struct StructureFactor{L, S} <: AbstractStructureFactor
    liquid::L
    scheme::S
end

function StructureFactor(liquid::L, T::Type; kw...) where {L<:AbstractLiquid}
    scheme = T(liquid; kw...)
    return StructureFactor{L, typeof(scheme)}(liquid, scheme)
end

(S::StructureFactor)(k) = structure_factor(S.liquid, S.scheme, k)

"""
    structure_factor(liquid, scheme, k)

Returns the static structure factor of a `liquid` using the approximation
defined by `scheme` at the wavenumber `k`.
"""
function structure_factor(liquid, scheme, k)
    Ck = Ĉ(liquid, scheme, k)
    return 1 / (1 - Ck)
end

function structure_factor(liquid::DipolarHardSpheres, scheme::MSA, k)
    C₀₀k, C₁₀k, C₁₁k = Ĉ(liquid, scheme, k)

    S₀₀k = 1 / (1 - C₀₀k)
    S₁₀k = 1 / (1 - C₁₀k)
    S₁₁k = 1 / (1 - C₁₁k)

    return (S₀₀k, S₁₀k, S₁₁k)
end

function Ĉ(liquid::HardDisks, scheme::RosenfeldFMT, k′)
    η = volume_fraction(liquid)
    σ = length_scale(potential(liquid))

    k  = σ * k′
    k² = k * k
    J₀ = besselj0(k / 2)
    J₁ = besselj1(k / 2)
    J₁′ = besselj1(k)

    A, B, G = scheme.A, scheme.B, scheme.G

    smallk = k < SMALL_K # for small wavenumber use a Taylor series approximation

    C = smallk ? (A * (1 -  k² / 16) / 4 + B * (1 - 3k² / 32) / 2 +
                  G * (1 -  k² / 8 )) :
                 (A * (2J₁ / k)^2 + B * 2J₀ * J₁ / k + G * 2J₁′ / k)

    return -4η * C
end

"""
    Ĉ(liquid::HardSpheres, scheme::PercusYevick, k)

Returns the product of the bulk density and the Fourier transform of the direct
correlation function for a hard-spheres liquid using the Percus-Yevick
approximation.

`k` is the wavenumber.
"""
function Ĉ(liquid::HardSpheres, scheme::PercusYevick, k′)
    η = volume_fraction(liquid)
    σ = length_scale(potential(liquid))

    k  = σ * k′
    k² = k * k
    k³ = k² * k
    k⁴ = k² * k²
    k⁶ = k³ * k³
    sink, cosk = sin(k), cos(k)

    α, β, δ = scheme.α, scheme.β, scheme.δ

    smallk = k < SMALL_K # for small wavenumber use a Taylor series approximation

    C₀ = smallk ? (1 - k² / 10) / 3 :
                  (sink - k * cosk) / k³
    C₁ = smallk ? (1 - k² / 9 ) / 4 :
                  (2k * sink - (k² - 2) * cosk - 2) / k⁴
    C₃ = smallk ? (1 - k² / 8 ) / 6 :
                  ((4k³ - 24k) * sink - (k⁴ - 12k² + 24) * cosk + 24) / k⁶

    return 24η * (α * C₀ + β * C₁ + δ * C₃)
end

"""
    Ĉ(::HardSpheres, ::VerletWeis, k)

Returns the product of the bulk density and the Fourier transform of the direct
correlation function for a hard-spheres liquid using the Percus-Yevick
approximation with the Verlet-Weis correction.

`k` is the wavenumber.
"""
Ĉ(liquid::HardSpheres, scheme::VerletWeis, k) =
    Ĉ(scheme.coreliquid, scheme.subscheme, scheme.α * k)

"""
    Ĉ(::AttractiveHardSpheres{<:YHSPotential}, ::SharmaSharma, k)

Returns the product of the bulk density and the Fourier transform (FT) of the
direct correlation function for an attractive Yukawa hard-spheres liquid using
the Sharma and Sharma approximation, that is, the sum of the FT of the direct
correlation function for the hard spheres core and ``\\beta\\cdot\\hat{u}(k)``,
where ``\\hat{u}(k)`` is the FT of the Yukawa potential:

```math
    u(r) = -\\epsilon\\exp(-zr/\\sigma) / (r/\\sigma)
```

`k` is the wavenumber.
"""
function Ĉ(liquid::AttractiveHardSpheres{P, U}, scheme::SharmaSharma, k′) where
    {P, U<:YHSPotential}

    η = volume_fraction(liquid)
    T = temperature(liquid)
    σ, z = length_scale(potential(liquid))

    k  = σ * k′
    k² = k * k
    k⁴ = k² * k²
    sink, cosk = sin(k), cos(k)

    smallk = k < SMALL_K # for small wavenumber use a Taylor series approximation

    C = smallk ? ((1 + z) - (3 + z) / 6 * k² + (5 + z) / 120 * k⁴) :
                  (k * cosk + z * sink) / k
    C = C / (k² + z * z)

    C₀ = Ĉ(scheme.coreliquid, scheme.subscheme, k′)

    return C₀ + 24η / T * C
end

"""
    Ĉ(::AttractiveHardSpheres{<:SWHSPotential}, ::SharmaSharma, k)

Returns the product of the bulk density and the Fourier transform (FT) of the
direct correlation function for an attractive square-well hard-spheres liquid
using the Sharma and Sharma approximation, that is, the sum of the FT of the
direct correlation function for the hard spheres core and
``\\beta\\cdot\\hat{u}(k)``, where ``\\hat{u}(k)`` is the FT of the square well
potential:

```math
    u(r) = -\\epsilon\\sigma, \\qquad \\sigma < r \\le \\lambda\\sigma
```

`k` is the wavenumber.
"""
function Ĉ(liquid::AttractiveHardSpheres{P, U}, scheme::SharmaSharma, k′) where
    {P, U<:SWHSPotential}

    η = volume_fraction(liquid)
    T = temperature(liquid)
    σ, λ = length_scale(potential(liquid))

    λ³ = λ * λ * λ
    λ⁵ = λ³ * λ * λ

    k  = σ * k′
    k² = k * k
    k³ = k² * k
    sink, cosk = sin(k), cos(k)
    sinλk, cosλk = sin(λ * k), cos(λ * k)

    smallk = k < SMALL_K # for small wavenumber use a Taylor series approximation

    C = smallk ? (λ³ - 1) / 3 - (λ⁵ - 1) / 30 * k² :
                 (cosk - λ * cosλk) / k² + (sinλk - sink) / k³

    C₀ = Ĉ(scheme.coreliquid, scheme.subscheme, k′)

    return C₀ + 24η / T * C
end

"""
    Ĉ(::AttractiveHardSpheres{<:SWHSPotential}, ::NonLinearSharmaSharma, k)

Returns the product of the bulk density and the Fourier transform (FT) of the
direct correlation function for an attractive square-well hard-spheres liquid
using the Non-Linear Sharma and Sharma approximation, that is, the sum of the FT of the
direct correlation function for the hard spheres core and
``exp[\\beta\\cdot\\u(r)]-1``, where ``u(r)`` is the potential:

```math
    u(r) = -\\epsilon\\sigma, \\qquad \\sigma < r \\le \\lambda\\sigma
```

`k` is the wavenumber.
"""
function Ĉ(liquid::AttractiveHardSpheres{P, U}, scheme::NonLinearSharmaSharma, k′) where
    {P, U<:SWHSPotential}

    η = volume_fraction(liquid)
    T = temperature(liquid)
    σ, λ = length_scale(potential(liquid))
    x = exp(1 / T) - 1

    λ³ = λ * λ * λ
    λ⁵ = λ³ * λ * λ

    k = σ * k′
    k² = k * k
    k³ = k² * k
    sink, cosk = sin(k), cos(k)
    sinλk, cosλk = sin(λ * k), cos(λ * k)

    smallk = k < 0.075

    C = smallk ? (λ³ - 1) / 3 - (λ⁵ - 1) / 30 * k² :
                 (cosk - λ * cosλk) / k² + (sinλk - sink) / k³

    C₀ = Ĉ(scheme.coreliquid, scheme.subscheme, k′)

    return C₀ + 24η * x * C
end

"""
    Ĉ(liquid::DipolarHardSpheres, scheme::MSA, k)

Returns the product of the bulk density and the Fourier transform of
projections of the direct correlation function `(C₀₀, C₁₀, C₁₁)` for a
dipolar hard-spheres liquid using the MSA approximation.

`k` is the wavenumber.
"""
function Ĉ(liquid::DipolarHardSpheres, scheme::MSA, k′)
    η = volume_fraction(liquid)
    σ = length_scale(potential(liquid))

    k  = σ * k′
    k² = k * k
    k⁶ = k² * k² * k²
    sink, cosk = sin(k), cos(k)

    α₀, α₁, α₂, α₃ = scheme.α₀, scheme.α₁, scheme.α₂, scheme.α₃
    β₀, β₁, β₂, β₃ = scheme.β₀, scheme.β₁, scheme.β₂, scheme.β₃
    a₁, a₃, b₁, b₃ = scheme.a₁, scheme.a₃, scheme.b₁, scheme.b₃

    smallk = k < SMALL_K # for small wavenumber use a Taylor series approximation

    C₁₀ = smallk ? (α₀ + 2β₀ - (α₁ + 2β₁) * k²) :
                   (  a₃ - a₁ * k² + (-a₃ + (α₃ + 2β₃) * k² / 3) * k * sink +
                    (-a₃ + ((a₁ + a₃ / 2) - (α₂ + 2β₂) * k² / 3) * k²) * cosk
                   ) / k⁶
    C₁₁ = smallk ? (α₀ -  β₀ - (α₁ -  β₁) * k²) :
                   (  b₃ - b₁ * k² + (-b₃ + (α₃ -  β₃) * k² / 3) * k * sink +
                    (-b₃ + ((b₁ + b₃ / 2) - (α₂ -  β₂) * k² / 3) * k²) * cosk
                   ) / k⁶

    C₀₀ = Ĉ(scheme.coreliquid, scheme.subscheme, k)
    C₁₀ = 24η * C₁₀
    C₁₁ = 24η * C₁₁

    return (C₀₀, C₁₀, C₁₁)
end


text:
using JSON
using Reactive
using Compat
import Compat.String

import Base: writemime
import Interact: update_view, Slider, Widget, InputWidget, Latex, HTML, recv_msg, statedict,
                 Progress, Checkbox, Button, ToggleButton, Textarea, Textbox, Options

export mimewritable, writemime

const ijulia_js = readstring(joinpath(dirname(@__FILE__), "ijulia.js"))

if displayable("text/html")
    display("text/html", """
     <div id="interact-js-shim">
         <script charset="utf-8">$(ijulia_js)</script>
         <script>
             window.interactLoadedFlag = true
            \$("#interact-js-shim").bind("destroyed", function () {
                if (window.interactLoadedFlag) {
                    console.warn("JavaScript required by Interact will be removed if you remove this cell or run using Interact more than once.")
                }
            })
            \$([IPython.events]).on("kernel_starting.Kernel kernel_restarting.Kernel", function () { window.interactLoadedFlag = false })
        </script>
     </div>""")
end

import IJulia
import IJulia: metadata, display_dict
using  IJulia.CommManager
import IJulia.CommManager: register_comm
import Base: show, mimewritable

const comms = Dict{Signal, Comm}()

function get_data_dict(value, mimetypes)
    dict = Dict{Compat.ASCIIString, Compat.String}()
    for m in mimetypes
        if mimewritable(m, value)
            dict[m] = stringmime(m, value)
        elseif m == "text/latex" && mimewritable("application/x-latex", value)
            dict[string("text/latex")] =
                stringmime("application/x-latex", value)
        else
            warn("IPython seems to be requesting an unavailable mime type")
        end
    end
    return dict
end

function init_comm(x::Signal)
    if !haskey(comms, x)
        subscriptions = Dict{Compat.ASCIIString, Int}()
        function handle_subscriptions(msg)
            if haskey(msg.content, "data")
                action = get(msg.content["data"], "action", "")
                if action == "subscribe_mime"
                    mime = msg.content["data"]["mime"]
                    subscriptions[mime] = get(subscriptions, mime, 0) + 1
                elseif action == "unsubscribe_mime"
                    mime = msg.content["data"]["mime"]
                    subscriptions[mime] = get(subscriptions, mime, 1) - 1
                end
            end
        end
        # One Comm channel per signal object
        comm = Comm(:Signal)
        comms[x] = comm   # Backend -> Comm
        # Listen for mime type registrations
        comm.on_msg = handle_subscriptions
        # prevent resending the first time?
        function notify(value)
            mimes = keys(filter((k,v) -> v > 0, subscriptions))
            if length(mimes) > 0
                send_comm(comm, @compat Dict(:value =>
                                 get_data_dict(value, mimes)))
            end
        end
        preserve(map(notify, x))
    else
        comm = comms[x]
    end

    return comm
end

function metadata(x::Signal)
    comm = init_comm(x)
    return @compat Dict("reactive"=>true,
                        "comm_id"=>comm.id)
end

function IJulia.display_dict(x::Signal)
    IJulia.display_dict(value(x))
end

# Render the value of a signal.
mimewritable(m::MIME, s::Signal) =
    mimewritable(m, s.value)

# fixes ambiguity warnings
@compat function Base.show(io::IO, m::MIME"text/plain", s::Signal)
    Base.show(io, m, s.value)
end

@compat function Base.show(io::IO, m::MIME"text/csv", s::Signal)
    Base.show(io, m, s.value)
end

@compat function Base.show(io::IO, m::MIME"text/tab-separated-values", s::Signal)
    Base.show(io, m, s.value)
end

@compat function Base.show(io::IO, m::MIME, s::Signal)
    Base.show(io, m, s.value)
end

@compat function Base.show(io::IO, ::MIME"text/html", w::InputWidget)
    create_view(w)
end

@compat function Base.show(io::IO, ::MIME"text/html", w::Widget)
    create_view(w)
end

@compat function Base.show{T<:Widget}(io::IO, ::MIME"text/html", x::Signal{T})
    create_widget_signal(x)
end

## This is for our own widgets.
function register_comm(comm::Comm{:InputWidget}, msg)
    w_id = msg.content["data"]["widget_id"]
    comm.on_msg = (msg) -> recv_msg(w, msg.content["data"]["value"])
end

JSON.print(io::IO, s::Signal) = JSON.print(io, s.value)

##################### IPython IPEP 23: Backbone.js Widgets #################

## ButtonView ✓
## CheckboxView ✓
## DropdownView ✓
## FloatSliderView ✓
## FloatTextView ✓
## IntSliderView ✓
## IntTextView ✓
## ProgressView
## RadioButtonsView ✓
## SelectView ✓
## TextareaView ✓
## TextView ✓
## ToggleButtonsView ✓
## ToggleButtonView ✓
## AccordionView W
## ContainerView W
## PopupView W
## TabView W

# Interact -> IJulia view names
widget_class(::HTML) = "HTML"
widget_class(::Latex) = "LaTeX"
widget_class(::Progress) = "Progress"
widget_class{T<:Integer}(::Slider{T}) = "IntSlider"
widget_class(::Button) = "Button"
widget_class(::Textarea) = "Textarea"
widget_class{T<:AbstractFloat}(::Slider{T}) = "FloatSlider"
widget_class{T<:Integer}(::Textbox{T}) = "IntText"
widget_class(::Checkbox) = "Checkbox"
widget_class(::ToggleButton) = "ToggleButton"
widget_class{T<:AbstractFloat}(::Textbox{T}) = "FloatText"
widget_class(::Textbox) = "Text"
widget_class{view}(::Options{view}) = string(view)
widget_class(w, suffix) = widget_class(w) * suffix
view_name(w) = widget_class(w, "View")
model_name(w) = widget_class(w, "Model")

function metadata{T <: Widget}(x::Signal{T})
    Dict()
end

function add_ipy3_state!(state)
    for attr in ["color" "background" "width" "height" "border_color" "border_width" "border_style" "font_style" "font_weight" "font_size" "font_family" "padding" "margin" "border_radius"]
        state[attr] = ""
    end
end

function add_ipy4_state!(state)
    state["_view_module"] = "jupyter-js-widgets"
    state["_model_module"] = "jupyter-js-widgets"
end

const widget_comms = Dict{Widget, Comm}()
function update_view(w; src=w)
    send_comm(widget_comms[w], view_state(w, src=src))
end

function view_state(w::InputWidget; src::InputWidget=w)
    msg = Dict()
    msg["method"] = "update"
    state = Dict()
    state["msg_throttle"] = 3
    state["_view_name"] = view_name(src)
    state["_model_name"] = model_name(src)
    state["model_name"] =  model_name(src)
    state["description"] = w.label
    state["visible"] = true
    state["disabled"] = false
    state["readout"] = true
    add_ipy3_state!(state)
    add_ipy4_state!(state)
    msg["state"] = merge(state, statedict(src))
    msg
end

function view_state(w::Widget; src::Widget=w)
    msg = Dict()
    msg["method"] = "update"
    state = Dict()
    state["msg_throttle"] = 3
    state["_view_name"] = view_name(src)
    state["_model_name"] = model_name(src)
    state["model_name"] = model_name(src)
    state["description"] = w.label
    state["visible"] = true
    state["disabled"] = false
    add_ipy3_state!(state)
    add_ipy4_state!(state)

    msg["state"] = merge(state, statedict(src))
    msg
end

function create_view(w::Widget)
    if haskey(widget_comms, w)
        comm = widget_comms[w]
    else
        comm = Comm("jupyter.widget", data=merge(Dict{AbstractString, Any}([
            ("model_name", model_name(w)),
            ("_model_name", model_name(w)), # Jupyter 4.0 missing (https://github.com/ipython/ipywidgets/pull/84)
            ("_view_module", "jupyter-js-widgets"),
            ("_model_module", "jupyter-js-widgets"),
        ]), view_state(w)))
        widget_comms[w] = comm
        # Send a full state update message.
        update_view(w) # This is redundant on 4.0 but keeps it working on Jupyter 3.0

        # dispatch messages to widget's handler
        comm.on_msg = msg -> handle_msg(w, msg)
        nothing # display() nothing
    end

    send_comm(comm, @compat Dict("method"=>"display"))
end

function create_widget_signal(s)
    create_view(s.value)
    local target = s.value
    preserve(map(x->update_view(target, src=x), s, init=nothing))
end

include("statedict.jl")
include("handle_msg.jl")


text:
module ReportReasons

"""
Abstract type to represent the reason why equality saturation halted
"""
abstract type ReportReason end

struct Saturated <: ReportReason end
Base.show(io::IO, x::Saturated) = print(io, "EGraph Saturated")
export Saturated

struct Timeout <: ReportReason end
Base.show(io::IO, x::Timeout) = print(io, "Iteration Timeout")
export Timeout

struct ConditionSatisfied <: ReportReason end
Base.show(io::IO, x::ConditionSatisfied) = print(io, "Halting condition satisfied")
export ConditionSatisfied

struct Contradiction <: ReportReason end
Base.show(io::IO, x::Contradiction) = print(io, "Equality Contradiction detected")
export Contradiction

struct EClassLimit <: ReportReason
    limit::Int
end
Base.show(io::IO, x::EClassLimit) = print(io, "Limit of $(x.limit) EClasses Exceeded")
export EClassLimit

struct NodeLimit <: ReportReason end
Base.show(io::IO, x::NodeLimit) = print(io, "ENode Limit Exceeded")
export NodeLimit

end


text:
import Base.convert
import Base.size

struct ForwardSim
    post_mean::Array{Float64,3}
    post_sd::Array{Float64,3}
    post_pd::Array{Float64,3}
    n::Matrix{Int}
    param::Matrix{Float64}
end

function get_interim_time(d::Int, data::AbstractData)
    Y = data.entry .+ data.time
    ev = Y[data.status]

    success = true

    if d > length(ev)
        @debug length(data.time), length(ev), d
        @debug "Requested number of events exceeds total number of events!"
        success = false
    end

    if length(ev) == 0
        @debug "No events"
        maximum(Y), false
    else
        events = sort(ev) ## ordered event times
        events[min(d, length(events))], success
    end
end

function append_data(x::T, y::T; add_ids=false) where T <: Union{BiPHWeibullData, BiPHSplineData}
    if add_ids
        id = vcat(x.id, maximum(x.id) .+ y.id)
    else
        id = vcat(x.id, y.id)
    end
    T(id, vcat(x.entry, y.entry), vcat(x.time, y.time), vcat(x.status, y.status), vcat(x.b, y.b), vcat(x.z, y.z))
end

function append_data(x::T, y::T; add_ids=false) where T <: WeibullPHData
    if add_ids
        id = vcat(x.id, maximum(x.id) .+ y.id)
    else
        id = vcat(x.id, y.id)
    end
    T(id, vcat(x.entry, y.entry), vcat(x.time, y.time), vcat(x.status, y.status), vcat(x.z, y.z))
end

function slice_data(x::T, sel::BitArray{1}) where T <: Union{BiPHWeibullData, BiPHSplineData}
    T(x.id[sel], x.entry[sel], x.time[sel], x.status[sel], x.b[sel], x.z[sel])
end

function select_ia_data(nevents, max_data::WeibullPHData; n_max::Int=size(max_data))
    at, success = get_interim_time(nevents, max_data)
    at_max = sort(max_data.entry)[min(length(max_data.entry),n_max)]
    at = min(at, at_max)
    at, select_data(at, max_data)
end

function select_ia_data(nevents, max_data::DoubleWeibullData; n_max::Int=size(max_data))
    at1, success = get_interim_time(sum(nevents), max_data.x)
    at2, success = get_interim_time(sum(nevents), max_data.y)

    et1 = sort(max_data.x.entry)[min(length(max_data.x.entry),n_max)]
    et2 = sort(max_data.y.entry)[min(length(max_data.y.entry),n_max)]

    max(at1, at2), DoubleWeibullData(select_data(at1, max_data.x; entry_max=min(et1, at1)),
                                      select_data(at2, max_data.y; entry_max=min(et2, at2)))
end

function select_ia_data(nevents, max_data::Union{BiPHWeibullData, BiPHSplineData}; n_max::Int=size(max_data))
    max_data_Sc = slice_data(max_data, max_data.b .== 1)
    max_data_S = slice_data(max_data, max_data.b .== 2)
    select_ia_data(nevents, max_data_Sc, max_data_S; n_max=n_max)
end

function select_ia_data(nevents, max_data_Sc, max_data_S; n_max=size(max_data_Sc)+size(max_data_S))
    at_Sc, s1 = get_interim_time(nevents[1], max_data_Sc)
    at_S, s2 = get_interim_time(nevents[2], max_data_S)

    if !s1 || !s2
        @debug length(max_data_Sc.time), length(max_data_S.time)
        @debug nevents, sum(max_data_Sc.status), sum(max_data_S.status)
        @debug at_Sc, at_S
    end

    at = max(at_Sc, at_S)

    new_data_Sc = select_data(at, max_data_Sc; entry_stop=at_Sc)
    new_data_S = select_data(at, max_data_S; entry_stop=at_S)
    new_data = append_data(new_data_Sc, new_data_S)

    if size(new_data) > n_max
        e_max = sort(new_data.entry)[min(length(new_data.entry),n_max)]
        new_data = select_data(new_data, e_max)
    end

    at, new_data
end

function extract_subpop(bm_pos::Bool, x::BiPHWeibullData)
    b = bm_pos ? 2 : 1
    sel = (x.b .== b)
    WeibullPH.Data(x.id[sel], x.entry[sel], x.time[sel], x.status[sel], x.z[sel,:])
end

get_par(param::Union{BiPHWeibullParam, BiPHSplineParam}) = vcat(collect(param.beta), param.prev)
get_par(param::DoubleWeibullParam) = vcat(param.x.beta, param.y.beta, 0.0)

convert(::Type{DoubleWeibullParam}, x::BiPHWeibullParam) = DoubleWeibullParam(WeibullPHParam(x.shape[1], x.scale[1], [x.beta[1]]),
                                                                              WeibullPHParam(x.shape[2], x.scale[2], [x.beta[2]]))

convert(::Type{DoubleWeibullData}, x::BiPHWeibullData) = DoubleWeibullData(extract_subpop(false, x), extract_subpop(true, x))

size(data::Union{BiPHWeibullData, BiPHSplineData}) = length(data.time)
size(data::WeibullPHData) = length(data.time)
size(data::DoubleWeibullData) = size(data.x) + size(data.y)

nparams(pr::Union{BiPHWeibullPrior, BiPHSplinePrior}) = length(pr.beta)
nparams(pr::DoubleWeibullPrior) = length(pr.x.beta) + length(pr.y.beta)

get_model_data(::Type{BiPHWeibullParam}, data::BiPHWeibullData) = data
get_model_data(::Type{BiPHSplineParam}, data::BiPHSplineData) = data
get_model_data(::Type{DoubleWeibullParam}, data::BiPHWeibullData) = convert(DoubleWeibullData, data)

## forward simulation of future IAs
function forward_simulate(J::Int, ia_time::Float64, delta, max_ev, data::D,
    m::Model{T,S}, x::F, prior0::P; B=10000) where {D <: AbstractData, F <: AbstractFit, P <: AbstractPrior, T <: AbstractParam, S <: AbstractFixedParam}

    B0 = size(x.beta, 1)
    K = nparams(prior0)

    n = SharedArray{Int,2}((B, J))
    post_mean = SharedArray{Float64,3}((B, J, K))
    post_sd = SharedArray{Float64,3}((B, J, K))
    post_pd = SharedArray{Float64,3}((B, J, K))
    par = SharedArray{Float64,2}((B,K+1))

    nmiss = SharedArray{Int,1}((J))
    nmiss .= 0

    @sync @distributed for b in 1:B
        param = slice(sample(1:B0), x)

        ## convert BiPHWeibullParam to DoubleWeibullParam
        model_param = convert(T, param)
        model_data = get_model_data(T, data)
        max_data = predict_data(model_param, m.fparam, model_data, ia_time, max_ev[end])

        for js in 1:J
            at, new_data = select_ia_data(max_ev[js], max_data)
            par[b,:] = get_par(param)

            if nevents(new_data) < sum(max_ev[js])
                @debug nevents(new_data), max_ev[js]
                nmiss[js] += 1
            end

            beta_sample = get_approx(model_param, new_data, prior0)
            post_mean[b, js, :] = mean.(beta_sample)
            post_sd[b, js, :] = sqrt.(var.(beta_sample))
            if K == 2
                post_pd[b, js, :] = [cdf(beta_sample[1], delta), cdf(beta_sample[2], delta)]
            else
                post_pd[b, js, :] = [cdf(beta_sample[1], delta)]
            end

            ##end
            n[b,js] = size(new_data)
        end
    end

    @debug nmiss

    ForwardSim(post_mean, post_sd, post_pd, n, par)
end

function fsim_model(K, B, M, warmup, m::Model, data::AbstractData, p::NTuple{3,AbstractPrior}, delta, nevents)
    x = fit(M, data; prior=p[1], warmup=warmup)
    x0 = fit(M, data; prior=p[2], warmup=warmup)
    x1 = fit(M, data; prior=p[3], warmup=warmup)

    fsim = forward_simulate(K, 0.0, delta, nevents, data, m, x, x.prior; B=B)
    fsim0 = forward_simulate(K, 0.0, delta, nevents, data, m, x0, x.prior; B=B)
    fsim1 = forward_simulate(K, 0.0, delta, nevents, data, m, x1, x.prior; B=B)

    fsim, fsim0, fsim1
end


text:
using FactCheck
using Distributions
using PDMats

using VinDsl
srand(12357)

facts("Flattening matrices") do
    context("Basic flattens") do
        @fact flatten(5) --> 5
        @fact flatten([[1, 2, 3] [4, 5, 6]]) --> [1, 2, 3, 4, 5, 6]
    end
    context("Triangulars") do
        A = rand(5, 5)
        U = UpperTriangular(A)
        L = LowerTriangular(A)
        @fact UpperTriangular(flatten(U)) --> U
        @fact LowerTriangular(flatten(L)) --> L
    end
    context("PDMats") do
        A = Float64[[10, 2, 3] [2, 40, 5] [3, 5, 60]]
        P = PDMat(A)
        @fact flatten(P) --> Float64[10, 2, 40, 3, 5, 60]
    end
end

facts("Checking normal_from_unconstrained") do
    context("Univariate") do
        x = randn(2)
        d = VinDsl.normal_from_unconstrained(x)
        @fact isa(d, Normal) --> true
        @fact d.μ --> x[1]
        @fact d.σ --> roughly(exp(x[2]))
    end
    context("Multivariate: mean field") do
        p = 7
        l = 2p
        x = randn(l)
        d = VinDsl.mvnormal_from_unconstrained(x)
        @fact isa(d, MvNormal) --> true
        @fact d.μ --> x[1:p]
        @fact d.Σ.diag --> roughly(exp(2 .* x[p+1:end]))
    end
    context("Multivariate: full covariance") do
        p = 7
        l = p * (p + 3) ÷ 2
        x = randn(l)
        d = VinDsl.mvnormal_from_unconstrained(x, true)
        @fact isa(d, MvNormal) --> true
        @fact d.μ --> x[1:p]
        @fact unconstrain(RCovMat(p), d.Σ) --> roughly(x[p+1:end])
    end
end


text:
using SequentialPointProcesses
using Distributions: logpdf
using Test

# Test that things run
@testset "SequentialPointProcesses.jl" begin
    M1(theta) = Softcore(d -> (2.3/d)^(theta))
    pp = rand(M1(4.5), (x=(0,10), y=(0,10)), 10)
    @test length(pp) == 10
    l1 = logpdf(M1(4.5), pp, 100)
    l2 = logpdf(M1(4.5), pp, SequentialPointProcesses.Options(nx=100, parallel=:threads))
    @test isfinite(l1)
    @test l1 ≈ l2 rtol=10eps(l1)

    M2(theta, prob2) = Mixture(M1(theta), SequentialPointProcesses.Uniform(), prob2)
    pp = rand(M2(4.5, 0.1), (x=(0,10), y=(0,10)), 10)
    @test length(pp) == 10
    l1 = logpdf(M2(4.5, 0.1), pp, 100)
    l2 = logpdf(M2(4.5, 0.1), pp, SequentialPointProcesses.Options(nx=100, parallel=:threads))
    @test isfinite(l1)
    @test l1 ≈ l2 rtol=10eps(l1)
end

max_stars_repo_name: mikkoku/SequentialPointProcesses

max_issues_repo_name: mikkoku/SequentialPointProcesses

max_forks_repo_name: mikkoku/SequentialPointProcesses


text:
using LinearAlgebra

using SparseArrays

using IterativeSolvers

using DelimitedFiles

using Revise

includet("Operators.jl")
using ..Operators

function avgq(velovec,p,u,v,nx,ny,dx,dy)
    cuti = (nx-1)*ny

    AVG = spzeros(length(velovec),length(velovec))

        for i = 1:nx
            for j = 1:ny

                if i != 1 && i != nx-1
                    AVG[u[j,i],u[j,i]] =      1/2
                    AVG[u[j,i],u[j,i-1]] =   -1/2
                end

                if i == 1

                end
            end
        end
end

using Plots
pyplot()

nx = 100
ny = 100

Lx = 1
Ly = 1

dx = Lx/nx
dy = Ly/ny

cuti  = (nx-1)*ny

Re = 1000
nu = 1/Re
CFL = 0.8
Δt = CFL*dx
tend = 100.0
tspan = 0:Δt:tend

px = range(0+0.5Lx/nx,Lx-0.5Lx/nx,length=nx)'  # note ': this is a row vector
py = range(0+0.5Ly/ny,Ly-0.5Ly/ny,length=ny)
# pvalue =  @. sin(2*pi*px) * cos(2*pi*py);
pvalue =  @. 0*px + 0*py;

ux = range(0+Lx/nx,Lx-Lx/nx,length=nx-1)'
uy = range(0+0.5Ly/ny,Ly-0.5Ly/ny,length=ny)
# uvalue =  @. 0*cos(2*pi*ux) + 0*sin(2*pi*uy);
uvalue =  @. 0*ux + 0*uy;

vx = range(0+0.5Lx/nx,Lx-0.5Lx/nx,length=nx)'  # note ': this is a row vector
vy = range(0+Ly/ny,Ly-Ly/ny,length=ny-1)
# vvalue=  @.  0*sin(2*pi*vx) + sin(2*pi*vy);
vvalue=  @.  0*vx + 0*vy;

uvec = vcat(uvalue...)
vvec = vcat(vvalue...)

velovec = [uvec;vvec]

pvec = vcat(pvalue...)

vec = [velovec;pvec];

p = reshape([i for i = 1:length(pvalue[:])],ny,nx)
u = reshape([i for i = 1:length(uvalue[:])],ny,nx-1)
v = reshape([i for i = length(uvalue[:])+1:length(uvalue[:])+length(vvalue[:])],ny-1,nx);

u_up = @. 0*ux + 1
u_bottom = @. 0*ux
u_left = @. 0*uy'
u_right = @. 0*uy'

v_up = @. 0*vx
v_bottom = @. 0*vx
v_left = @. 0*vy'
v_right = @. 0*vy'

vx

u_up

G = grad(velovec,p,u,v,pvec,nx,ny,dx,dy)
D = -transpose(G)
bcD = bcdiv(v_up,v_bottom,u_left,u_right,pvec,p,uy,vx,dx,dy)
L = lapq(velovec,p,u,v,nx,ny,dx,dy);
bcL = bclap(velovec,u_up,u_bottom,u_left,u_right,v_up,v_bottom,v_left,v_right,p,u,v,nx,ny,dx,dy)
A = adv(velovec,u_up,u_bottom,u_left,u_right,v_up,v_bottom,v_left,v_right,p,u,v,nx,ny,dx,dy);

R = sparse(I,size(L)) - Δt*nu/2 * L
invR = sparse(I,size(L)) + Δt*nu/2 * L + (Δt*nu/2 * L)^2;
S = sparse(I,size(L)) + Δt*nu/2 * L;

cgA = -D*invR*G

# initializing n-1, n, and n+1 (n=2)
Ahist = []
push!(Ahist,A)

# bcDhist = []
# push!(bcDhist,bcD)
# push!(bcDhist,bcD)
# push!(bcDhist,bcD);

# bcLhist = []
# push!(bcLhist,bcL)
# push!(bcLhist,bcL);

velovechist = []
push!(velovechist,velovec)
push!(velovechist,velovec);

pvechist = []
push!(pvechist,pvec)
push!(pvechist,pvec);

bcLhist = []
push!(bcLhist,bcL)
push!(bcLhist,bcL);
for i = 3:length(tspan)+1
    push!(bcLhist,bcL);
end

bcDhist = []
push!(bcDhist,bcD)
push!(bcDhist,bcD);
for i = 3:length(tspan)
    push!(bcDhist,bcD);
end

Pnew

Pnew = zero.(pvec)
un = velovechist[2]
# constant B.C.
@time for n = 2:length(tspan)
# @time for n = 2:20
#      un = velovechist[2]

      A = adv(un,u_up,u_bottom,u_left,u_right,v_up,v_bottom,v_left,v_right,p,u,v,nx,ny,dx,dy);

#     uf = invR*(S*un + Δt/2*(3*Ahist[n] - Ahist[n-1]) + Δt*nu/2*(bcLhist[n+1] + bcLhist[n]))
      i1,uf = cgsolv(R,(S*un + Δt/2*(3*A - Ahist[n-1]) + Δt*nu/2*(bcLhist[n+1] + bcLhist[n])),un,length(un),1e-6)

      cgb = -Array((bcDhist[n]+D*Array(uf)))/Δt

      i2,Pnew = cgsolv(cgA,cgb,Pnew,length(cgb),1e-6)

      unew = uf - Δt * invR * G * Pnew

      un = deepcopy(unew)

      push!(Ahist,A);

     if (mod(n,100) == 2) || n == length(tspan)

         push!(velovechist,unew);
         push!(pvechist,Pnew)
        end

#            print("n=",n," ")
end

uresult = []
vresult = []
for i = 1:length(velovechist)
    utemp = reshape(velovechist[i][1:cuti],ny,nx-1)
    push!(uresult,utemp)

    vtemp = reshape(velovechist[i][cuti+1:end],ny-1,nx)
    push!(vresult,vtemp)
end

presult = []
for i = 1:length(pvechist)
    ptemp = reshape(pvechist[i],ny,nx)
    push!(presult,ptemp)
end

@gif for i = 3:length(uresult)
    contour(uresult[i],aspect_ratio=:equal,levels=50,line_z=true)
    end every 1

@gif for i = 3:length(vresult)
    contour(vresult[i],aspect_ratio=:equal,levels=50,line_z=true)
    end every 1

function vq_to_vp(v_bottom,vresult,v_up)

    v_q = vcat(v_bottom,vresult,v_up)
    v_p = zeros(size(v_q,1)-1,size(v_q,2))

    for i = 1:size(v_q,2)
        v_p[i,:] = (v_q[i,:] + v_q[i+1,:])/2
    end

    return v_p
end

function uq_to_up(u_left,uresult,u_right)

    u_q = hcat(u_left',uresult,u_right')
    u_p = zeros(size(u_q,1),size(u_q,2)-1)

    for i = 1:size(u_q,1)
        u_p[i,:] = (u_q[:,i] + u_q[:,i+1])/2
    end

    return u_p'
end

v_p = vq_to_vp(v_bottom,vresult[end],v_up);
u_p = uq_to_up(u_left,uresult[end],u_right);

open

open(string("matlabx",Re,".txt"), "w") do io
           writedlm(io, px, '\n')
       end

open(string("matlaby",Re,".txt"), "w") do io
           writedlm(io, py, '\n')
       end

open(string("matlabup",Re,".txt"), "w") do io
           writedlm(io, u_p, '\n')
       end

open(string("matlabvp",Re,".txt"), "w") do io
           writedlm(io, v_p, '\n')
       end

open(string("matlabp",Re,".txt"), "w") do io
        writedlm(io, presult[end], '\n')
       end

open(string("matlabu",Re,".txt"), "w") do io
        writedlm(io, uresult[end], '\n')
       end

open(string("matlabv",Re,".txt"), "w") do io
        writedlm(io, vresult[end], '\n')
       end

ualog = [0.0;-0.03717;-0.04192;-0.04775;-0.06434;-0.10150;-0.15662;-0.21090;-0.20581;-0.13641;0.00332;0.23151;0.68717;0.73722;0.78871;0.84123;1.0]
ualogy = 1/128*([1;8;9;10;14;23;37;59;65;80;95;110;123;124;125;126;129] .- 1)


text:
# Commonly used constants

const q  = 1.6021765e-19    # electron charge, [C]
const me = 9.10938356e-31   # electron mass, [kg]
const mp = 1.6726219e-27    # proton mass, [kg]
const qₑ = -q
const mₑ = me
const qᵢ = q
const mᵢ = mp
const amu= 1.66054e-27      # Atomic Mass Unit, 1/12 of C12, [kg]
const ϵ0 = 8.8542e-12       # vacuum permittivity, [F/m]
const ϵ₀ = ϵ0
const c  = 3e8              # speed of light, [m/s]
const μ₀ = 4π*1e-7          # Vacuum permeability, [H/m]

const Rₑ = 6.38e6           # Earth radius, [m]
const Rm = 2444000.         # Mercury's radius, [m]

# Solar
const Tsolarrotate = 25.05*24*3600 # [s]
const ωsun = 2π / Tsolarrotate # [rad/s]
const Rsun = 696340e3 # solar radius, [m]
const AU = 1.496e11 # sun-Earth distance, [m]

# Distance from sun
# https://www.enchantedlearning.com/subjects/astronomy/planets/
const dMercury = 0.39  # [AU], 57.91e9 [m]
const dVenus   = 0.723 # [AU], 108.2e9 [m]
const dMars    = 1.524 # [AU], 227.9e9 [m]
const dJupiter = 5.203 # [AU], 778.5e9 [m]
const dSaturn  = 9.539 # [AU], 1434e9  [m]
const dUranus  = 19.18 # [AU], 2871e9  [m]
const dNeptune = 30.06 # [AU], 4495e9  [m]