text:
abstract type AbstractThrustCoefficientModel end

"""
    ThrustModelConstantCt(ct::Float)

Stores a constant ct value for wake calculations

# Arguments
- `ct::Float`: a constant ct value for computation
"""
struct ThrustModelConstantCt{TF} <: AbstractThrustCoefficientModel
    ct::TF
end

"""
    ThrustModelCtPoints(vel_points, ct_points)

Stores the thrust coefficient curve in terms of corresponding velocity and thrust
coefficient points. ct and velocity points should be in the same order and ordered from
lowest wind speed to highest wind speed.

# Arguments
- `inflow_velocity::Float`: inflow velocity of the wind turbine
- `thrust_model::ThrustModelCtPoints`: Struct containing ct and velocity points for ct curve
"""
struct ThrustModelCtPoints{ATF} <: AbstractThrustCoefficientModel
    vel_points::ATF
    ct_points::ATF
end

"""
    calculate_ct(model::ThrustModelConstantCt)

Calculate the thrust coefficient for a wind turbine based on a pre-determined constant ct

# Arguments
- `inflow_velocity::Float`: inflow velocity of the wind turbine (unused for const. ct)
- `thrust_model::ThrustModelConstantCt`: struct containing a constant ct value for computation
"""
function calculate_ct(inflow_velocity, thrust_model::ThrustModelConstantCt)
    return thrust_model.ct
end

"""
    calculate_ct(inflow_velocity, thrust_model::ThrustModelCtPoints)

    Calculate the thrust coefficient for a wind turbine based on a pre-determined ct curve
        with linear interpolation.

# Arguments
- `inflow_velocity::Float`: inflow velocity of the wind turbine
- `thrust_model::ThrustModelCtPoints`: Struct containing ct and velocity points for ct curve
"""
function calculate_ct(inflow_velocity, thrust_model::ThrustModelCtPoints)

    minv = thrust_model.vel_points[1]
    maxv = thrust_model.vel_points[end]
    minct = thrust_model.ct_points[1]
    maxct = thrust_model.ct_points[end]

    if inflow_velocity < minv
        ct = minct
    elseif inflow_velocity < maxv
        ct = linear(thrust_model.vel_points, thrust_model.ct_points, inflow_velocity)
    else
        ct = maxct
    end

    if ct > 1
        ct = 1.0
    end

    return ct
end

"""
    _ct_to_axial_ind_func(ct)

Calculate axial induction from the thrust coefficient

# Arguments
- `ct::Float`: thrust coefficient
"""
function _ct_to_axial_ind_func(ct)

    # initialize axial induction to zero
    axial_induction = 0.0

    # calculate axial induction
    if ct > 0.96  # Glauert condition
        axial_induction = 0.143 + sqrt(0.0203 - 0.6427*(0.889 - ct))
    else
        axial_induction = 0.5*(1.0 - sqrt(1.0 - ct))
    end

    return axial_induction

end


text:
using CRCBS
using GraphUtils
using TOML
using Test, Logging
using Parameters

"""
    FatPathsSolver{S}

A simple wrapper around a MAPF solver so that profile_solver! will dispatch
correctly (by transforming adding the fat paths cost to the problem before
solving it).
"""
@with_kw struct FatPathsSolver{S,M} <: SolverWrapper
    solver::S = CBSSolver(AStar{NTuple{2,Float64}}())
    fat_path_cost_model::M = FlatFPCost()
end
function CRCBS.profile_solver!(solver::FatPathsSolver,mapf)
    fp_mapf = CRCBS.init_fat_path_mapf(mapf,solver.fat_path_cost_model)
    CRCBS.profile_solver!(solver.solver,fp_mapf)
end
# CRCBS.get_logger(solver::FatPathsSolver) = get_logger(solver.solver)
# CRCBS.low_level(solver::FatPathsSolver) = low_level(solver.solver)

# Problem Instances
base_scen_path = joinpath(ENV["HOME"],"Repos/mapf_benchmarks/scenarios")
map_path = joinpath(ENV["HOME"],"Repos/mapf_benchmarks/maps/")
results_path = "/scratch/mapf_experiments"

EXPERIMENTS_DIR = "/scratch/mapf_experiments/fat_path_experiments"
PROBLEM_DIR = joinpath(EXPERIMENTS_DIR,"problems")
RESULTS_DIR = joinpath(EXPERIMENTS_DIR,"results")

config = (
    solver_configs = [
        (
            solver=CBSSolver(),
            results_path=joinpath(RESULTS_DIR,"CBSSolver")
        ),
        (
            solver=FatPathsSolver(),
            results_path=joinpath(RESULTS_DIR,"FatPathsSolver")
        ),
        (
            solver=FatPathsSolver(fat_path_cost_model=NormalizedFPCost()),
            results_path=joinpath(RESULTS_DIR,"NormalizedFatPathsSolver")
        ),
    ],
    problem_dir = PROBLEM_DIR,
    feats = [
        RunTime(),IterationCount(),SolutionCost(),NumConflicts(),RobotPaths(),
        RobotSeparation(),
        TimeOutStatus(),IterationMaxOutStatus(),
        # MemAllocs(),ByteCount()
    ]
)
for solver_config in config.solver_configs
    set_runtime_limit!(solver_config.solver,50)
    set_verbosity!(solver_config.solver,1)
    set_iteration_limit!(solver_config.solver,10000)
    set_iteration_limit!(low_level(solver_config.solver),1000)
end

# scen_paths = get_files_matching(base_scen_path,".scen",["Berlin_1_256","Paris_1_256"])
scen_paths = get_files_matching(base_scen_path,".scen",[
    "empty-8-8-even",
    # "empty-16-16-even",
    # "empty-32-32-even",
    # "empty-48-48-even",
    ])

BenchmarkInterface.generate_problem_files_from_moving_ai(
    scen_paths,
    # [joinpath(base_scen_path,"scen-even","empty-8-8-even-10.scen")],
    map_path,
    PROBLEM_DIR
)

loader = BenchmarkInterface.init_mapf_loader(PROBLEM_DIR)
BenchmarkInterface.profile_with_skipping!(config,loader)
# run_profiling(config,loader)

ext: jl


text:
using Test
using TestSetExtensions
using LinearAlgebra
using Qaintessent
using SparseArrays
using StatsBase

##==----------------------------------------------------------------------------------------------------------------------

isunitary(cg::CircuitGate) = (sparse_matrix(cg) * sparse_matrix(Base.adjoint(cg)) ≈ I)

@testset ExtendedTestSet "circuit gates" begin
    θ = 0.7 * π
    ϕ = 0.4 * π
    n = randn(3); n /= norm(n)

    # single qubit gates
    @testset "single qubit circuit gates" begin
        for g in [X, Y, Z, HadamardGate(), SGate(), TGate(), RxGate(θ), RyGate(θ), RzGate(θ), RotationGate(θ, n), PhaseShiftGate(ϕ)]
            cg = CircuitGate((2,), g)
            cgadj = adjoint(cg)
            Qaintessent.sparse_matrix(cgadj.gate) == adjoint(Qaintessent.sparse_matrix(cg.gate))
            @test LinearAlgebra.ishermitian(cg) == (Qaintessent.sparse_matrix(cg) == Qaintessent.sparse_matrix(adjoint(cg)))
        end

        cgs = circuit_gate.((2,), [X, Y, Z, HadamardGate(), SGate(), TGate(), RxGate(θ), RyGate(θ), RzGate(θ), RotationGate(θ, n), PhaseShiftGate(ϕ)])

        @test all(sparse_matrix(adjoint(cgs)) .≈ adjoint(sparse_matrix(cgs)))
    end

    # two qubit gates
    @testset "two qubit circuit gates" begin
        for g in [EntanglementXXGate(θ), EntanglementYYGate(θ), EntanglementZZGate(θ), controlled_not(), SwapGate()]
            cg = CircuitGate((2, 3), g)
            cgadj = adjoint(cg)
            Qaintessent.sparse_matrix(cgadj.gate) == adjoint(Qaintessent.sparse_matrix(cg.gate))
            @test LinearAlgebra.ishermitian(cg) == (Qaintessent.sparse_matrix(cg) == Qaintessent.sparse_matrix(adjoint(cg)))
        end
    end

    # Y acting on second wire
    @testset "apply circuit gate to second wire" begin
        cg = CircuitGate((2,), Y)
        @test Qaintessent.sparse_matrix(cg, 3) ≈ kron(Matrix(I, 2, 2), Qaintessent.matrix(Y), Matrix(I, 2, 2))
        @test isunitary(cg)
    end

    # flip control and target
    @testset "flip control and target circuit gate" begin
        cg = CircuitGate((2, 1), controlled_not())
        @test Qaintessent.sparse_matrix(cg) ≈ [1 0 0 0; 0 0 0 1; 0 0 1 0; 0 1 0 0]
        @test isunitary(cg)
    end

    # third qubit as control and first qubit as target
    @testset "shift control and target circuit gate" begin
        cg = circuit_gate(1, HadamardGate(), 3)
        @test Qaintessent.sparse_matrix(cg) ≈ [
            Matrix(I, 4, 4) fill(0, 4, 2) fill(0, 4, 2);
            fill(0, 2, 4) Qaintessent.matrix(HadamardGate()) fill(0, 2, 2);
            fill(0, 2, 6) Qaintessent.matrix(HadamardGate())]
        @test isunitary(cg)
    end

    @testset "circuit gate exceptions" begin
        H = HadamardGate()
        S = SwapGate()

        N = 2
        @test_throws ErrorException("SwapGate affects 2 wires but 0 wires, (), were passed.") CircuitGate{0,SwapGate}(NTuple{0,Int}(), S)
        @test_throws ErrorException("Wire indices must be unique.") CircuitGate{2,SwapGate}((1, 1), S)
        @test_throws ErrorException("Wire index cannot be smaller than 1.") CircuitGate{2,SwapGate}((1, -1), S)
    end
end

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "circuit gate isapprox" begin
    θ = 0.7 * π
    ϕ = 0.4 * π
    n = randn(3); n /= norm(n)
    ϵ = 3*sqrt(eps())
    sqg = [RxGate(θ), RyGate(θ), RzGate(θ), RotationGate(θ, n), PhaseShiftGate(ϕ)]
    sqḡ = [RxGate(θ + eps()), RyGate(θ + eps()), RzGate(θ + eps()), RotationGate(θ + eps(), n), PhaseShiftGate(ϕ + eps())]
    sqĝ = [RxGate(θ + ϵ), RyGate(θ + ϵ), RzGate(θ + ϵ), RotationGate(θ + ϵ, n), PhaseShiftGate(ϕ + ϵ)]

    for (i, g) in enumerate(sqg)
        cg1 = CircuitGate((2,), sqg[i])
        cg2 = CircuitGate((2,), sqḡ[i])
        cg3 = CircuitGate((2,), sqĝ[i])

        @test cg1 ≈ cg2
        @test !(cg1 ≈ cg3)
    end
end

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "circuit gate helper functions" begin
    N = 5
    @testset ExtendedTestSet "circuit gate single qubit helper function" begin
        iwire = rand(1:N)
        g = XGate()
        @test circuit_gate(iwire, g) ≈ CircuitGate((iwire,), g)
    end

    @testset "circuit gate two qubit helper function" begin
        iwire1 = rand(1:N)
        iwire2 = rand(vcat(1:iwire1 - 1..., iwire1 + 1:N...))
        g2 = SwapGate()
        @test circuit_gate(iwire1, iwire2, g2) ≈ CircuitGate((iwire1, iwire2), g2)
    end

    @testset "circuit gate controlled gate helper function" begin
        cntrl_iwire1, cntrl_iwire2, targ_iwire1, targ_iwire2 = sample(1:N, 4, replace=false)

        g = YGate()

        ref_cg = CircuitGate((targ_iwire1, cntrl_iwire1), ControlledGate(g, 1))

        @test circuit_gate(targ_iwire1,    g,  cntrl_iwire1)   ≈ ref_cg
        @test circuit_gate(targ_iwire1,    g, (cntrl_iwire1,)) ≈ ref_cg
        @test circuit_gate((targ_iwire1,), g,  cntrl_iwire1)   ≈ ref_cg
        @test circuit_gate((targ_iwire1,), g, (cntrl_iwire1,)) ≈ ref_cg

        g = SwapGate()

        ref_cg2 = CircuitGate((targ_iwire1, targ_iwire2, cntrl_iwire1), ControlledGate(g, 1))

        @test circuit_gate(targ_iwire1, targ_iwire2, g,  cntrl_iwire1)   ≈ ref_cg2
        @test circuit_gate(targ_iwire1, targ_iwire2, g, (cntrl_iwire1,)) ≈ ref_cg2
        @test circuit_gate((targ_iwire1, targ_iwire2), g,  cntrl_iwire1)   ≈ ref_cg2
        @test circuit_gate((targ_iwire1, targ_iwire2), g, (cntrl_iwire1,)) ≈ ref_cg2

        ref_cg3 = CircuitGate((targ_iwire1, targ_iwire2, cntrl_iwire1, cntrl_iwire2), ControlledGate(g, 2))

        @test circuit_gate(targ_iwire1, targ_iwire2, g,  cntrl_iwire1, cntrl_iwire2)   ≈ ref_cg3
        @test circuit_gate(targ_iwire1, targ_iwire2, g, (cntrl_iwire1, cntrl_iwire2)) ≈ ref_cg3
        @test circuit_gate((targ_iwire1, targ_iwire2), g,  cntrl_iwire1, cntrl_iwire2)   ≈ ref_cg3
        @test circuit_gate((targ_iwire1, targ_iwire2), g, (cntrl_iwire1, cntrl_iwire2)) ≈ ref_cg3
    end

    # test sparse_matrix
    @testset "circuit gates sparse matrix" begin
        cgs = [circuit_gate(1, X), circuit_gate(2, X), circuit_gate(3, X)]
        m = sparse_matrix(cgs)

        @test m ≈ sparse([8, 7, 6, 5, 4, 3, 2, 1], [1, 2, 3, 4, 5, 6, 7, 8], Float64[1, 1, 1, 1, 1, 1, 1, 1])
    end

    # test sparse_matrix
    @testset "circuit gates sparse matrix exceptions" begin
        cgs = [circuit_gate(1, X), circuit_gate(2, X), circuit_gate(3, X)]
        @test_throws ErrorException("Circuit size `2` too small; vector of CircuitGate requires 3 wires.") sparse_matrix(cgs, 2)
    end
end


text:
abstract type Estimator end

Base.iterate(estimator::Estimator, state = 0) = state > 0 ? nothing : (estimator, state + 1)
Base.length(estimator::Estimator) = 1

Base.show(io::IO, estimator::Estimator) = print(io, string(estimator))
Base.show(io::IO, ::MIME"application/prs.juno.inline", estimator::Estimator) = print(io, string(estimator))

(::Estimator)(x1::Integer, x2::Integer, design::AbstractDesign) = error("not implemented")
function estimate(estimator::TE, x1::TI, x2::TI, design::TD) where {TE<:Estimator,TI<:Integer,TD<:AbstractDesign}
    return estimator(x1, x2, design)
end

string(estimator::Estimator) = error("not implemented")

function bias(p::Real, estimator::Estimator, design::AbstractDesign)
    0 <= p <= 1 ? nothing : error("p must be between 0 and 1")
    XX = sample_space(design)
    x1, x2 = XX[:,1], XX[:,2]
    return (estimator.(x1, x2, design) .- p) .*
        pmf.(x2, n2.(design, x1), p) .*
        pmf.(x1, n1(design), p) |>
        sum
end

function mean_squared_error(p::Real, estimator::Estimator, design::AbstractDesign)
    0 <= p <= 1 ? nothing : error("p must be between 0 and 1")
    XX = sample_space(design)
    x1, x2 = XX[:,1], XX[:,2]
    return (estimator.(x1, x2, design) .- p).^2 .*
        pmf.(x2, n2.(design, x1), p) .*
        pmf.(x1, n1(design), p) |>
        sum
end

function mean_absolute_error(p::Real, estimator::Estimator, design::AbstractDesign)
    0 <= p <= 1 ? nothing : error("p must be between 0 and 1")
    XX = sample_space(design)
    x1, x2 = XX[:,1], XX[:,2]
    return abs.(estimator.(x1, x2, design) .- p) .*
        pmf.(x2, n2.(design, x1), p) .*
        pmf.(x1, n1(design), p) |>
        sum
end


text:
module Mod1DynamicDefinition
using ValidatedNumerics
using ..DynamicDefinition, ..Contractors

using ..DynamicDefinition: derivative

export Mod1Dynamic, preim, nbranches, plottable

"""
Defines a Dynamic on [0,1] as the Mod-1 quotient of a given map.

An alternative newer implementation relying on piecewise-defined functions is in `mod1_dynamic`
"""
struct Mod1Dynamic{FT} <: MarkovDynamic
	T::FT
	nbranches::Int
	orientation::Float64
	domain::Interval{Float64}
	is_full_branch::Bool
end

Mod1Dynamic(T::FT, nbranches = undef, domain = Interval{Float64}(0,1)) where {FT} = Mod1Dynamic{FT}(T, nbranches, domain)

function Mod1Dynamic{FT}(T, nbranches = undef, domain = Interval{Float64}(0,1)) where {FT}
	@assert domain == 0..1 # TODO: this only works for domain == 0..1, for now
	range_diff = T(@interval(1.))-T(@interval(0.))
	orientation = unique_sign(range_diff)
	nbranches = ceil(orientation * range_diff).hi
	is_full_branch = isinteger(T(0..0)) & isinteger(T(1..1))
	return Mod1Dynamic{FT}(T, nbranches, orientation, domain, is_full_branch)
end

DynamicDefinition.domain(S::Mod1Dynamic{FT}) where {FT} = S.domain

DynamicDefinition.nbranches(S::Mod1Dynamic{FT}) where {FT} =S.nbranches

DynamicDefinition.is_full_branch(S::Mod1Dynamic{FT}) where {FT} = S.is_full_branch

# TODO: serious doubts that this works if T(0) is not an integer...
function DynamicDefinition.preim(D::Mod1Dynamic{FT}, k, y, ϵ) where {FT}
	# we need to treat the case with the other orientation, 0 not fixed point...
	@assert 1 <= k <= D.nbranches
	f(x) = D.T(x)-D.T(0)-(y-D.T(0)+(k-1)*D.orientation)
	root(f, D.domain, ϵ)
end

DynamicDefinition.derivative(n, D::Mod1Dynamic{FT}, x) where {FT} = derivative(n, D.T, x)
DynamicDefinition.distorsion(D::Mod1Dynamic{FT}, x) where {FT} = distorsion(D.T, x)

DynamicDefinition.max_distorsion(D::Mod1Dynamic{FT}, tol=1e-3) where {FT} = maximise(x -> distorsion(D.T, x), domain(D), tol=tol)[1]
DynamicDefinition.expansivity(D::Mod1Dynamic{FT}, tol=1e-3) where {FT} = maximise(x -> abs(1/derivative(D, x)), domain(D), tol=tol)[1]

function DynamicDefinition.plottable(D::Mod1Dynamic{FT}, x) where {FT}
	@assert 0 <= x <= 1
	return mod(D.T(x), 1.)
end

using RecipesBase
@recipe f(::Type{Mod1Dynamic{FT}}, D::Mod1Dynamic{FT}) where {FT} = x -> plottable(D, x)

orientation(D::Mod1Dynamic, k) = D.orientation

end


text:
const berstandata = [
  Dict(
  "N" => 10,
  "y" => [0,1,0,0,0,0,0,0,0,1]
  )
]


text:
mutable struct MinimaxObjectiveCore
    points::Vector{Float64}
    targets::Vector{Float64}
    IIF::Matrix{Float64}

    MinimaxObjectiveCore() =
        new(Vector{Float64}(), Vector{Float64}(), Matrix{Float64}(undef, 0, 0))
end


text:
#!/usr/bin/env julia

# Reproduce the VTK_BEZIER_TETRA_quartic_solidSphereOctant.vtu file from the VTK
# test suite.

using WriteVTK

using Test

const VTK_BASENAME = "bezier_tetra_quartic_solidSphereOctant"

function main()
    cell_type = VTKCellTypes.VTK_BEZIER_TETRAHEDRON

    # Copied from VTK generated file.
    points_in = [
        1, 0, 0, 0, 1, 0,
        0, 0, 1, 0, 0, 0,
        1, 0.4226497411727905, 0, 0.7886751294136047, 0.7886751294136047, 0,
        0.4226497411727905, 1, 0, 0, 1, 0.4226497411727905,
        0, 0.7886751294136047, 0.7886751294136047, 0, 0.4226497411727905, 1,
        0.4226497411727905, 0, 1, 0.7886751294136047, 0, 0.7886751294136047,
        1, 0, 0.4226497411727905, 0.75, 0, 0,
        0.5, 0, 0, 0.25, 0, 0,
        0, 0.75, 0, 0, 0.5, 0,
        0, 0.25, 0, 0, 0, 0.75,
        0, 0, 0.5, 0, 0, 0.25,
        0.5, 0.25, 0, 0.25, 0.5, 0,
        0.25, 0.25, 0, 0, 0.25, 0.5,
        0, 0.25, 0.25, 0, 0.5, 0.25,
        0.5, 0, 0.25, 0.25, 0, 0.25,
        0.25, 0, 0.5, 1, 0.5893534421920776, 0.5893534421920776,
        0.5893534421920776, 0.5893534421920776, 1, 0.5893534421920776, 1, 0.5893534421920776,
        0.30000001192092896, 0.30000001192092896, 0.30000001192092896,
    ]
    points = reshape(points_in, 3, :)

    connectivity = 1:size(points, 2)

    rational_weights = [
        1, 1, 1, 1, 0.8365163037378083, 0.7886751345948131,
        0.8365163037378083, 0.8365163037378083, 0.7886751345948131, 0.8365163037378083, 0.8365163037378083, 0.7886751345948131,
        0.8365163037378083, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1,
        1, 0.6594659464934113, 0.6594659464934113, 0.6594659464934113, 1,
    ]

    cells = [MeshCell(cell_type, connectivity)]

    outfiles = vtk_grid(VTK_BASENAME, points, cells) do vtk
        vtk["RationalWeights"] = rational_weights
    end

    println("Saved:  ", join(outfiles, "  "))

    outfiles
end

main()


text:
using CairoMakie
using Colors
using Makie.GeometryBasics
using JLD2

training_data = [
    (-5e-4, 5e-8),
    (-3.5e-4, 5e-8),
    (-2e-4, 5e-8),
]

train_files = [
    "wind_-5e-4_diurnal_5e-8"
    "wind_-3.5e-4_diurnal_5e-8"
    "wind_-2e-4_diurnal_5e-8"
]

interpolating_data = [
    (-4.5e-4, 4e-8),
    (-4.5e-4, 2e-8),
    (-3e-4, 4e-8),
    (-3e-4, 2e-8),
]

interpolating_files = [
    "wind_-4.5e-4_diurnal_4e-8"
    "wind_-4.5e-4_diurnal_2e-8"
    "wind_-3e-4_diurnal_4e-8"
    "wind_-3e-4_diurnal_2e-8"
]

extrapolating_data_diurnal = [
    (-5.5e-4, 5.5e-8),
    (-1.5e-4, 5.5e-8),
]

extrapolating_diurnal_files = [
    "wind_-5.5e-4_diurnal_5.5e-8"
    "wind_-1.5e-4_diurnal_5.5e-8"
]

extrapolating_data = [
    (-1.5e-4, 3.5e-8),
    (-5.5e-4, 3.5e-8),
    (-5.5e-4, 0),
    (-5.5e-4, -3.5e-8),
    (-1.5e-4, -3.5e-8),
]

extrapolating_files = [
    "wind_-1.5e-4_cooling_3.5e-8"
    "wind_-5.5e-4_cooling_3.5e-8"

    "wind_-5.5e-4_new"

    "wind_-5.5e-4_heating_-3.5e-8"
    "wind_-1.5e-4_heating_-3.5e-8"
]

momentum_fluxes_training = [data[1] for data in training_data]
buoyancy_fluxes_training = [data[2] for data in training_data]

momentum_fluxes_interpolating = [data[1] for data in interpolating_data]
buoyancy_fluxes_interpolating = [data[2] for data in interpolating_data]

momentum_fluxes_extrapolating = [data[1] for data in extrapolating_data]
buoyancy_fluxes_extrapolating = [data[2] for data in extrapolating_data]

momentum_fluxes_extrapolating_diurnal = [data[1] for data in extrapolating_data_diurnal]
buoyancy_fluxes_extrapolating_diurnal = [data[2] for data in extrapolating_data_diurnal]

fig = CairoMakie.Figure(resolution=(1000, 650))

ax =fig[1,1] = CairoMakie.Axis(fig, xlabel="Buoyancy Flux / m² s⁻³", ylabel="Momentum Flux / m² s⁻²")
color_palette = distinguishable_colors(4, [RGB(1,1,1), RGB(0,0,0)], dropseed=true)

rectangle = CairoMakie.poly!(ax, Point2f0[(-5e-8, -2e-4), (-5e-8, -5e-4), (5e-8, -5e-4), (5e-8, -2e-4)], color=("paleturquoise3", 0.5))

training_lines = [CairoMakie.lines!(ax, [-buoyancy_fluxes_training[i], buoyancy_fluxes_training[i]],
                    [momentum_fluxes_training[i], momentum_fluxes_training[i]], color=color_palette[1]) for i in 1:length(training_data)]

interpolating_lines = [CairoMakie.lines!(ax, [-buoyancy_fluxes_interpolating[i], buoyancy_fluxes_interpolating[i]],
                    [momentum_fluxes_interpolating[i], momentum_fluxes_interpolating[i]], color=color_palette[2]) for i in 1:length(interpolating_data)]

extrapolating_lines = [CairoMakie.lines!(ax, [-buoyancy_fluxes_extrapolating_diurnal[i], buoyancy_fluxes_extrapolating_diurnal[i]],
                    [momentum_fluxes_extrapolating_diurnal[i], momentum_fluxes_extrapolating_diurnal[i]], color=color_palette[3]) for i in 1:length(extrapolating_data_diurnal)]

extrapolating_points = CairoMakie.scatter!(ax, buoyancy_fluxes_extrapolating, momentum_fluxes_extrapolating, color=color_palette[3])

# diurnal_line_1 = CairoMakie.lines!(ax, [-5.5e-8, 5.5e-8], [-1.5e-4, -1.5e-4], color=color_palette[4])
# diurnal_line_2 = CairoMakie.lines!(ax, [-5.5e-8, 5.5e-8], [-5.5e-4, -5.5e-4], color=color_palette[4])

# # diurnal_points = CairoMakie.scatter!(ax, buoyancy_fluxes_diurnal, momentum_fluxes_diurnal, color=color_palette[4])

legend = fig[2,1] = CairoMakie.Legend(fig, [training_lines[1], interpolating_lines[1], extrapolating_lines[1], extrapolating_points, rectangle],
                                            ["Training", "Interpolating", "Extrapolating", "Constant Fluxes", "Interpolation Region"], orientation=:horizontal)

rowsize!(fig.layout, 1, CairoMakie.Relative(0.95))
trim!(fig.layout)
fig
save("final_results/data_diurnal.png", fig, px_per_unit = 4)

T_loss(data) = sum(data) / 1153

losses_training = []

for train_file in train_files
    file = jldopen("final_results/3sim_diurnal/train_$(train_file)/profiles_fluxes_oceananigans.jld2")
    data = file["NDE_profile"]
    close(file)

    loss = T_loss(data["T_losses"])
    loss_mpp = T_loss(data["T_losses_modified_pacanowski_philander"])
    loss_kpp = T_loss(data["T_losses_kpp"])

    loss_min = argmin([loss, loss_mpp, loss_kpp])

    loss_str = ["NDE", "mpp", "kpp"]
    push!(losses_training, loss_str[loss_min])
end

losses_training

losses_interpolating = []

for interpolating_file in interpolating_files
    file = jldopen("final_results/3sim_diurnal/test_$(interpolating_file)/profiles_fluxes_oceananigans.jld2")
    data = file["NDE_profile"]
    close(file)

    loss = T_loss(data["T_losses"])
    loss_mpp = T_loss(data["T_losses_modified_pacanowski_philander"])
    loss_kpp = T_loss(data["T_losses_kpp"])

    @show loss, loss_mpp, loss_kpp

    loss_min = argmin([loss, loss_mpp, loss_kpp])

    loss_str = ["NDE", "mpp", "kpp"]
    push!(losses_interpolating, loss_str[loss_min])
end

losses_interpolating

losses_extrapolating = []

for extrapolating_file in extrapolating_files
    file = jldopen("final_results/3sim_diurnal/test_$(extrapolating_file)/profiles_fluxes_oceananigans.jld2")
    data = file["NDE_profile"]
    close(file)

    loss = T_loss(data["T_losses"])
    loss_mpp = T_loss(data["T_losses_modified_pacanowski_philander"])
    loss_kpp = T_loss(data["T_losses_kpp"])

    @show loss, loss_mpp, loss_kpp

    loss_min = argmin([loss, loss_mpp, loss_kpp])

    loss_str = ["NDE", "mpp", "kpp"]
    push!(losses_extrapolating, loss_str[loss_min])
end

losses_extrapolating

losses_extrapolating_diurnal = []

for diurnal_file in extrapolating_diurnal_files
    file = jldopen("final_results/3sim_diurnal/test_$(diurnal_file)/profiles_fluxes_oceananigans.jld2")
    data = file["NDE_profile"]
    close(file)

    loss = T_loss(data["T_losses"])
    loss_mpp = T_loss(data["T_losses_modified_pacanowski_philander"])
    loss_kpp = T_loss(data["T_losses_kpp"])

    @show loss, loss_mpp, loss_kpp

    loss_min = argmin([loss, loss_mpp, loss_kpp])

    loss_str = ["NDE", "mpp", "kpp"]
    push!(losses_extrapolating_diurnal, loss_str[loss_min])
end

losses_extrapolating_diurnal

fig = CairoMakie.Figure(resolution=(1000, 650))

ax =fig[1,1] = CairoMakie.Axis(fig, xlabel="Buoyancy Flux / m² s⁻³", ylabel="Momentum Flux / m² s⁻²")
color_palette = distinguishable_colors(4, [RGB(1,1,1), RGB(0,0,0)], dropseed=true)

loss_colors = Dict(
    "NDE" => color_palette[1],
    "mpp" => color_palette[3],
    "kpp" => color_palette[4],
)

rectangle = CairoMakie.poly!(ax, Point2f0[(-5e-8, -2e-4), (-5e-8, -5e-4-3f-6), (5e-8, -5e-4-3f-6), (5e-8, -2e-4)], color=("paleturquoise3", 0.5))

NDE_line = CairoMakie.lines!(ax, [buoyancy_fluxes_training[1], buoyancy_fluxes_training[1]], [momentum_fluxes_training[1], momentum_fluxes_training[1]], color=loss_colors["NDE"])
mpp_line = CairoMakie.lines!(ax, [buoyancy_fluxes_training[1], buoyancy_fluxes_training[1]], [momentum_fluxes_training[1], momentum_fluxes_training[1]], color=loss_colors["mpp"])
kpp_line = CairoMakie.lines!(ax, [buoyancy_fluxes_training[1], buoyancy_fluxes_training[1]], [momentum_fluxes_training[1], momentum_fluxes_training[1]], color=loss_colors["kpp"])

training_lines = [CairoMakie.lines!(ax, [-buoyancy_fluxes_training[i], buoyancy_fluxes_training[i]],
                    [momentum_fluxes_training[i], momentum_fluxes_training[i]], color=loss_colors[losses_training[i]]) for i in 1:length(losses_training)]

interpolating_lines = [CairoMakie.lines!(ax, [-buoyancy_fluxes_interpolating[i], buoyancy_fluxes_interpolating[i]],
                    [momentum_fluxes_interpolating[i], momentum_fluxes_interpolating[i]], color=loss_colors[losses_interpolating[i]]) for i in 1:length(losses_interpolating)]

extrapolating_lines = [CairoMakie.lines!(ax, [-buoyancy_fluxes_extrapolating[i], buoyancy_fluxes_extrapolating[i]],
                    [momentum_fluxes_extrapolating[i], momentum_fluxes_extrapolating[i]], color=loss_colors[losses_extrapolating[i]]) for i in 1:length(losses_extrapolating_diurnal)]

extrapolating_points = [CairoMakie.scatter!(ax, [buoyancy_fluxes_extrapolating[i]], [momentum_fluxes_extrapolating[i]], color=loss_colors[losses_extrapolating[i]]) for i in 1:length(losses_extrapolating)]

# training_points = [CairoMakie.scatter!(ax, [buoyancy_fluxes_training[i]], [momentum_fluxes_training[i]], color=loss_colors[losses_training[i]]) for i in 1:length(losses_training)]
# interpolating_points = [CairoMakie.scatter!(ax, [buoyancy_fluxes_interpolating[i]], [momentum_fluxes_interpolating[i]], color=loss_colors[losses_interpolating[i]]) for i in 1:length(losses_interpolating)]
# extrapolating_points = [CairoMakie.scatter!(ax, [buoyancy_fluxes_extrapolating[i]], [momentum_fluxes_extrapolating[i]], color=loss_colors[losses_extrapolating[i]]) for i in 1:length(losses_extrapolating)]

# # extrapolating_points = CairoMakie.scatter!(ax, buoyancy_fluxes_extrapolating, momentum_fluxes_extrapolating, color=color_palette[3])

# diurnal_line_1 = CairoMakie.lines!(ax, [-5.5e-8, 5.5e-8], [-5.5e-4, -5.5e-4], color=loss_colors[losses_diurnal[1]])
# diurnal_line_2 = CairoMakie.lines!(ax, [-5.5e-8, 5.5e-8], [-1.5e-4, -1.5e-4], color=loss_colors[losses_diurnal[2]])

legend = fig[2,1] = CairoMakie.Legend(fig, [mpp_line, kpp_line, NDE_line],
                                ["Ri-based Diffusivity Only", "K-Profile Parameterisation", "NN Embedded in Oceananigans.jl"], orientation=:horizontal)

rowsize!(fig.layout, 1, CairoMakie.Relative(0.95))
trim!(fig.layout)
fig
save("final_results/data_diurnal_loss_results.png", fig, px_per_unit = 4)


text:
# using Pipe
using Test

_macroexpand(q) = macroexpand(Main, q)

rml! = Base.remove_linenums!
# performs linenum removal and temp variable replacing to avoid different names of temp variables in different julia versions
stringify_expr(e::Expr) = replace(string(rml!(e)), r"##\d{3}" => "##000")
pipe_equals(e1::Expr, e2::Expr) = stringify_expr(_macroexpand(e1)) == stringify_expr(e2)

@testset "`_` 1st location" begin
  @test _macroexpand(:(@pipe a |> b(x))) == :((b(a, x)))  #applying to function calls returning functions
  @test _macroexpand(:(@pipe a |> bb[2])) == :(bb[a, 2]) #Should work with RHS that is a array reference
  @test _macroexpand(:(@pipe a |> b(xb, _) |> c |> d(_, xd) |> e(xe) |> f(xf, _, yf) |> _[i])) == :((f(xf, e(d(c(b(xb, a)), xd), xe), yf))[i])
end

#Marked locations
@testset "Pipe main" begin
  #No change to nonpipes functionality
  @test _macroexpand(:(@pipe a)) == :a #doesn't change single inputs
  @test _macroexpand(:(@pipe b(a))) == :(b(a)) #doesn't change inputs that a function applications

  #Compatable with Julia 1.3 piping functionality
  @test _macroexpand(:(@pipe a |> b)) == :(b(a)) #basic
  @test _macroexpand(:(@pipe a |> b |> c)) == :(c(b(a)))  #Keeps chaining 3
  @test _macroexpand(:(@pipe a |> b |> c |> d)) == :(d(c(b(a)))) #Keeps chaining 4

  # @test _macroexpand( :(@pipe a|>b(x)) ) == :((b(x))(a))  #applying to function calls returning functions
  @test _macroexpand(:(@pipe a(x) |> b)) == :(b(a(x)))   #feeding functioncall results on wards

  @test _macroexpand(:(@pipe 1 |> a)) == :(a(1)) #Works with literals (int)
  @test _macroexpand(:(@pipe "foo" |> a)) == :(a("foo")) #Works with literal (string)
  # @test _macroexpand( :(@pipe a|>bb[2])) == :((bb[2])(a)) #Should work with RHS that is a array reference

  @test _macroexpand(:(@pipe a |> _)) == :(a) #Identity works
  @test _macroexpand(:(@pipe a |> _[b])) == :(a[b]) #Indexing works
  @test _macroexpand(:(@pipe a |> b(_))) == :(b(a)) #Marked location only
  @test _macroexpand(:(@pipe a |> b(x, _))) == :(b(x, a)) # marked 2nd (and last)
  @test _macroexpand(:(@pipe a |> b(_, x))) == :(b(a, x)) # marked first
  @test _macroexpand(:(@pipe a |> b(_, _))) == :(b(a, a)) # marked double (Not certain if this is a good idea)
  @test _macroexpand(:(@pipe a |> bb[2](x, _))) == :((bb[2])(x, a)) #Should work with RHS that is a array reference
end

#Macros and blocks
macro testmacro(arg, n)
  esc(:($arg + $n))
end
@test _macroexpand(:(@pipe a |> @testmacro _ 3)) == :(a + 3) # Can pipe into macros
@test _macroexpand(:(@pipe a |> begin
  b = _
  c + b + _
end)) == :(
  begin
    b = a
    c + b + a
  end)

#marked Unpacking
@test _macroexpand(:(@pipe a |> b(_...))) == :(b(a...)) # Unpacking
@test _macroexpand(:(@pipe a |> bb[2](_...))) == :((bb[2])(a...)) #Should work with RHS of arry ref and do unpacking

#Mixing modes
@test _macroexpand(:(@pipe a |> b |> c(_))) == :(c(b(a)))
@test _macroexpand(:(@pipe a |> b(x, _) |> c |> d(_, y))) == :(d(c(b(x, a)), y))
# @test _macroexpand( :(@pipe a|>b(xb,_)|>c|>d(_,xd)|>e(xe) |>f(xf,_,yf)|>_[i] ) ) == :(f(xf,(e(xe))(d(c(b(xb,a)),xd)),yf)[i]) #Very Complex

# broadcasting
@testset "Pipe broadcasting" begin
  vars = 1:10 .|> y -> gensym() # Julia < 1.3 changes how Symbols are stringified so we compute the representation here
  @test pipe_equals(:(@pipe 1:10 .|> _ * 2), :(1:10 .|> $(vars[1]) -> $(vars[1]) * 2))
  @test pipe_equals(:(@pipe 1:10 .|> fn), :(1:10 .|> $(vars[2]) -> fn($(vars[2]))))
  @test pipe_equals(:(@pipe a .|> fn .|> _ * 2), :(a .|> ($(vars[3]) -> fn($(vars[3]))) .|> ($(vars[4]) -> $(vars[4]) * 2)))
  @test pipe_equals(:(@pipe a .|> fn |> _ * 2), :((a .|> $(vars[5]) -> fn($(vars[5]))) * 2))
  @test pipe_equals(:(@pipe [1, 2, 2] |> atan.([10, 20, 30], _)), :(atan.([10, 20, 30], [1, 2, 2])))
  @test pipe_equals(:(@pipe [1, 2, 2] .|> atan.([10, 20, 30], _)), :([1, 2, 2] .|> $(vars[6]) -> atan.([10, 20, 30], $(vars[6]))))
  @test pipe_equals(:(@pipe fn |> _.(1:2)), :(fn.(1:2)))
  @test pipe_equals(:(@pipe fn .|> _.(1:2)), :(fn .|> $(vars[7]) -> $(vars[7]).(1:2)))

  @test pipe_equals(:(@pipe [true, false] .|> !), :([true, false] .|> $(vars[8]) -> !$(vars[8])))
  @test pipe_equals(:(@pipe [1, 2] |> .+(_, x)), :([1, 2] .+ x))
  @test pipe_equals(:(@pipe [1, 2] |> _ .+ x), :([1, 2] .+ x))
  @test pipe_equals(:(@pipe [1, 2] .|> .+(_, x)), :([1, 2] .|> $(vars[9]) -> $(vars[9]) .+ x))
  @test pipe_equals(:(@pipe [1, 2] .|> _ .+ x), :([1, 2] .|> $(vars[10]) -> $(vars[10]) .+ x))
end


text:
using ClusterTrees, BEAST, CompScienceMeshes

a, h = 1.0, 0.25
Γ1 = meshsphere(a, h)
Γ2 = CompScienceMeshes.translate(Γ1, point(2.1,0,0))

X1 = lagrangecxd0(Γ1)
X2 = lagrangecxd0(Γ2)

p = positions(X1)
q = positions(X2)

p, tp, permp = clustertree(p)
q, tq, permq = clustertree(q)

function adm(b) # p and q are passed in implicitly
    nmin = 20
    I = b[1][1].begin_idx : b[1][1].end_idx-1
    J = b[2][1].begin_idx : b[2][1].end_idx-1
    length(I) < nmin && return true
    length(J) < nmin && return true
    ll1, ur1 = ClusterTrees.boundingbox(p[I]); c1 = (ll1+ur1)/2;
    ll2, ur2 = ClusterTrees.boundingbox(q[J]); c2 = (ll2+ur2)/2;
    diam1 = norm(ur1-c1)
    diam2 = norm(ur2-c2)
    dist12 = norm(c2-c1)
    return dist12 >= η*max(diam1, diam2)
end

blocktree, η = (tp,tq), 2.0
P = admissable_partition(blocktree, adm)

# sanity check on the tree
for (τ,σ) in P
    tree_size = 0
    depthfirst(τ) do c,l
        tree_size += 1
    end
    @assert τ[1].num_children+1 == tree_size
end

depthfirst(tp) do τ,l
    I = τ[1].begin_idx : τ[1].end_idx-1
    for c in children(τ)
        J = c[1].begin_idx : c[1].end_idx-1
        @assert J ⊆ I
    end
end

# gather all the observer clusters that participate in an interaction pair
Q = Vector{typeof(P[1][1])}()
for (i,p) in enumerate(P)
    τ1, _ = p[1], p[2]
    I = τ1[1].begin_idx : τ1[1].end_idx-1
    descendant = false
    for q in P
        τ2, _ = q[1], q[2]
        J = τ2[1].begin_idx : τ2[1].end_idx-1
        # Can τ2 be reached from τ1?
        descendant = ((J != I) && (J ⊆ I))
        !isempty(intersect(I,J)) && !(J ⊆ I) && !(I ⊆ J) && (@show I J)
        descendant && break
    end
    !descendant && push!(Q,τ1)
end
Q = unique(Q)
length(Q)
length(P)

balance = zeros(numfunctions(X1))
for q in Q
    I = q[1].begin_idx : q[1].end_idx-1
    balance[permp[I]] .+= 1
end
@show extrema(balance)

using MATLAB
mat"figure()"
mat"hold("on")"
V = [v[i] for v in Γ1.vertices, i in 1:3]
F = [f[i] for f in Γ1.faces, i in 1:3]
@mput V F
for (r,q) in enumerate(Q)

    I = q[1].begin_idx : q[1].end_idx-1
    I = permp[I]
    @mput I r
    mat"patch("Vertices",V,"Faces",F[I,:],"FaceColor",rand(3,1))"
end

## Another sanity check: everything can be reached from the root
for (i,τ) in enumerate(tp)
    reached = false
    depthfirst(tp) do st,l
        st[1] == τ && (reached = true)
    end
    if !reached
        @show τ
        @show i
        error()
    end
end


text:
c1 = Patience(1)
c2 = InvalidValue()
c3 = TimeLimit(t=100)

@test Disjunction(c1) == c1

d = c1 + c2 + Never() + c3 + c1
show(d)

# codecov:
@test zero(typeof(c1)) == Never()

@test sum(Disjunction[]) == Never()
@test Never() in c1

@testset "_criteria" begin
    criteria = EarlyStopping._criteria(d)
    @test length(criteria) == 3
    @test issubset([c1, c2, c3], criteria)
end

@testset "stoppping times" begin
    d2 = Patience(3) + InvalidValue()
    @test stopping_time(d2, [12.0, 10.0, 11.0, 12.0, 13.0, NaN]) == 5
    @test stopping_time(d2, [NaN, 12.0, 10.0, 11.0, 12.0, 13.0]) == 1
end

@testset "message" begin
    state = EarlyStopping.update(d, NaN)
    @test EarlyStopping.message(d, state) ==
        "Stopping early as `NaN`, "*
        "`Inf` or `-Inf` encountered. "
end

state = EarlyStopping.update(d, 1.0)
state = EarlyStopping.update(d, 2.0, state)
@test EarlyStopping.message(d, state) ==
    "Stop triggered by Patience(1) stopping criterion. "


text:
export FT_Init_FreeType
export FT_Done_FreeType
export FT_New_Face
export FT_New_Memory_Face
export FT_Open_Face
export FT_Attach_File
export FT_Attach_Stream
export FT_Reference_Face
export FT_Done_Face
export FT_Select_Size
export FT_Request_Size
export FT_Set_Char_Size
export FT_Set_Pixel_Sizes
export FT_Load_Glyph
export FT_Load_Char
export FT_Set_Transform
export FT_Render_Glyph
export FT_Get_Kerning
export FT_Get_Track_Kerning
export FT_Get_Glyph_Name
export FT_Get_Postscript_Name
export FT_Select_Charmap
export FT_Set_Charmap
export FT_Get_Charmap_Index
export FT_Get_Char_Index
export FT_Get_First_Char
export FT_Get_Next_Char
export FT_Get_Name_Index
export FT_Get_SubGlyph_Info
export FT_Get_FSType_Flags
export FT_Face_GetCharVariantIndex
export FT_Face_GetCharVariantIsDefault
export FT_Face_GetVariantSelectors
export FT_Face_GetVariantsOfChar
export FT_Face_GetCharsOfVariant
export FT_MulDiv
export FT_DivFix
export FT_RoundFix
export FT_CeilFix
export FT_FloorFix
export FT_Vector_Transform
export FT_Library_Version
export FT_Face_CheckTrueTypePatents
export FT_Face_SetUnpatentedHinting
export FT_Outline_Decompose

export FT_Int16
export FT_UInt16
export FT_Int32
export FT_UInt32
export FT_Fast
export FT_UFast
export FT_Memory
export FT_Alloc_Func
export FT_Free_Func
export FT_Realloc_Func
export FT_Stream
export FT_Stream_IoFunc
export FT_Stream_CloseFunc
export FT_Pos
export FT_Pixel_Mode
export FT_Outline_MoveToFunc
export FT_Outline_LineToFunc
export FT_Outline_ConicToFunc
export FT_Outline_CubicToFunc
export FT_Glyph_Format
export FT_Raster
export FT_SpanFunc
export FT_Raster_BitTest_Func
export FT_Raster_BitSet_Func
export FT_Raster_NewFunc
export FT_Raster_DoneFunc
export FT_Raster_ResetFunc
export FT_Raster_SetModeFunc
export FT_Raster_RenderFunc
export FT_Bool
export FT_FWord
export FT_UFWord
export FT_Char
export FT_Byte
export FT_Bytes
export FT_Tag
export FT_String
export FT_Short
export FT_UShort
export FT_Int
export FT_UInt
export FT_Long
export FT_ULong
export FT_F2Dot14
export FT_F26Dot6
export FT_Fixed
export FT_Error
export FT_Pointer
export FT_Offset
export FT_PtrDist
export FT_Generic_Finalizer
export FT_ListNode
export FT_List
export FT_Library
export FT_Module
export FT_Driver
export FT_Renderer
export FT_Face
export FT_Size
export FT_GlyphSlot
export FT_CharMap
export FT_Encoding
export FT_Face_Internal
export FT_Size_Internal
export FT_SubGlyph
export FT_Slot_Internal
export FT_Size_Request_Type
export FT_Size_Request
export FT_Render_Mode
export FT_Kerning_Mode

export FT_LOAD_DEFAULT
export FT_LOAD_TARGET_NORMAL
export FT_LOAD_NO_SCALE
export FT_LOAD_NO_HINTING
export FT_LOAD_RENDER
export FT_LOAD_NO_BITMAP
export FT_LOAD_VERTICAL_LAYOUT
export FT_LOAD_FORCE_AUTOHINT
export FT_LOAD_CROP_BITMAP
export FT_LOAD_PEDANTIC
export FT_LOAD_ADVANCE_ONLY
export FT_LOAD_IGNORE_GLOBAL_ADVANCE_WIDTH
export FT_LOAD_NO_RECURSE
export FT_LOAD_IGNORE_TRANSFORM
export FT_LOAD_MONOCHROME
export FT_LOAD_LINEAR_DESIGN
export FT_LOAD_SBITS_ONLY
export FT_LOAD_NO_AUTOHINT
export FT_LOAD_TARGET_LIGHT
export FT_LOAD_TARGET_MONO
export FT_LOAD_TARGET_LCD
export FT_LOAD_TARGET_LCD_V
export FT_LOAD_COLOR

export FT_FACE_FLAG_SCALABLE
export FT_FACE_FLAG_FIXED_SIZES
export FT_FACE_FLAG_FIXED_WIDTH
export FT_FACE_FLAG_SFNT
export FT_FACE_FLAG_HORIZONTAL
export FT_FACE_FLAG_VERTICAL
export FT_FACE_FLAG_KERNING
export FT_FACE_FLAG_FAST_GLYPHS
export FT_FACE_FLAG_MULTIPLE_MASTERS
export FT_FACE_FLAG_GLYPH_NAMES
export FT_FACE_FLAG_EXTERNAL_STREAM
export FT_FACE_FLAG_HINTER
export FT_FACE_FLAG_CID_KEYED
export FT_FACE_FLAG_TRICKY
export FT_FACE_FLAG_COLOR

export FT_RENDER_MODE_NORMAL
export FT_RENDER_MODE_LIGHT
export FT_RENDER_MODE_MONO
export FT_RENDER_MODE_LCD
export FT_RENDER_MODE_LCD_V
export FT_RENDER_MODE_MAX

export FT_PIXEL_MODE_NONE
export FT_PIXEL_MODE_MONO
export FT_PIXEL_MODE_GRAY
export FT_PIXEL_MODE_GRAY2
export FT_PIXEL_MODE_GRAY4
export FT_PIXEL_MODE_LCD
export FT_PIXEL_MODE_LCD_V
export FT_PIXEL_MODE_BGRA
export FT_PIXEL_MODE_MAX


text:
module GatedLinearNetworks

export GaussianGGLN, train!, predict

using Flux
using Statistics

"""
Normalize input using a online estimate of mean and variance using Welford's
algorithm.
"""
struct NormalizationLayer{VF <: AbstractVector{<:Real}, VI <: AbstractVector{Int}}
    count::VI # [1]
    mean::VF # [data_dim]
    m2::VF # [data_dim]
end

function NormalizationLayer(n::Int)
    return NormalizationLayer(
        Int[0],
        zeros(Float32, n),
        zeros(Float32, n))
end

"""
Apply standardization to the input, and update parameters if `training` is true.
This is basically batch normalization.

x is [data_dim, batch_dim]
"""
function forward(layer::NormalizationLayer, x::AbstractMatrix, training::Bool)
    if training
        for j in 1:size(x, 2)
            layer.count .+= 1
            xj = x[:,j]
            delta = xj - layer.mean
            layer.mean .+= delta ./ layer.count
            delta2 = xj - layer.mean
            layer.m2 .+= delta .* delta2
        end
    end

    @assert layer.count[1] > 0

    sd = sqrt.(layer.m2 ./ layer.count)
    return (x .- layer.mean) ./ sd
end

function inverse(layer::NormalizationLayer, z::AbstractArray)
    sd = sqrt.(layer.m2 ./ layer.count)
    return (z .* sd) .+ layer.mean
end

function inverse(layer::NormalizationLayer, μ::AbstractArray, σ2::AbstractArray)
    sd = sqrt.(layer.m2 ./ layer.count)

    μ_inv = (μ .* sd) .+ layer.mean
    σ2_inv = σ2 .* sd.^2

    return μ_inv, σ2_inv
end

abstract type GaussianWeakLearner end

"""
Bayesian gaussian linear regression with known precision. Used as a weak
learner input to the GGLN.
"""
struct BasicGaussianLinearRegression{
        VF <: AbstractVector{<:Real},
        MF <: AbstractMatrix{<:Real},
        VI <: AbstractVector{Int}} <: GaussianWeakLearner
    # accumulated sufficient statistics
    xy::MF # [predictor_dim, prediction_dim]
    xx::VF # [predictor_dim]
    x::VF # [prediction_dim]
    y::VF # [prediction_dim]
    count::VI # [1]

    # known prior precision
    τ0::VF # [prediction_dim]
    τ::VF # [prediction_dim]
end

"""
Construct a bayesian linear regression.

  * `predictor_dim`: dimensionality of predictors
  * `prediction_dim`: dimensionality of predictions
"""
function BasicGaussianLinearRegression(predictor_dim::Int, prediction_dim::Int)
    return BasicGaussianLinearRegression(
        zeros(Float32, (predictor_dim, prediction_dim)),
        zeros(Float32, predictor_dim),
        zeros(Float32, predictor_dim),
        zeros(Float32, prediction_dim),
        zeros(Int, 1),
        Float32[1f0],
        Float32[1f0])
end

"""
Forward pass for the linear regression. Train when `y` is given, otherwise
only output posterior predictive parameters.
"""
function forward(
        layer::BasicGaussianLinearRegression,
        x::AbstractMatrix, y::Union{Nothing, AbstractMatrix})

    predictor_dim = size(layer.xy, 1)
    prediction_dim = size(layer.xy, 2)
    batch_dim = size(x, 2)

    @assert size(x, 1) == predictor_dim
    @assert y === nothing || size(y, 1) == prediction_dim
    @assert y === nothing || size(y, 2) == batch_dim

    # [predictor_dim]
    c1 = (layer.τ0 .+ layer.τ .* layer.count) .*
        (layer.τ0 .+ layer.τ .* layer.xx) .- (layer.τ .* layer.x).^2

    # [predictor_dim]
    c2 = layer.τ ./ c1

    # [predictor_dim, prediction_dim]
    uw = (layer.τ0 .+ layer.τ .* layer.count) .* layer.xy

    vw = layer.τ .* reshape(layer.x, (predictor_dim, 1)) .*
        reshape(layer.y, (1, prediction_dim))

    # [predictor_dim, prediction_dim]
    μ_w = c2 .* (uw .- vw)

    # [predictor_dim, prediction_dim]
    ub = (layer.τ0 .+ layer.τ .* reshape(layer.xx, (predictor_dim, 1))) .*
        reshape(layer.y, (1, prediction_dim))

    vb = layer.τ .* layer.x .* layer.xy

    # [predictor_dim, predction_dim]
    μ_b = c2 .* (ub .- vb)

    # [predictors_dim, batch_dim]
    σ2 = inv.(layer.τ) .+ inv.(c1) .* (
        (layer.τ0 .+ layer.τ .* layer.count) .* x.^2 .-
        2f0 .* layer.τ .* x .* layer.x .+
        layer.τ0 .+
        layer.τ .* layer.xx)

    # repeat this across the prediction_dim, since this is an isotropic gaussian

    # [predictior_dim, prediction_dim, batch_dim]
    σ2 = repeat(reshape(σ2, (predictor_dim, 1, batch_dim)), 1, prediction_dim, 1)

    # [predictors_dim, predictions_dim, batch_dim]
    μ = reshape(x, (predictor_dim, 1, batch_dim)) .*
        reshape(μ_w, (predictor_dim, prediction_dim, 1)) .+
        reshape(μ_b, (predictor_dim, prediction_dim, 1))

    if y !== nothing
        layer.xx .+= dropdims(sum(x.*x, dims=2), dims=2)
        layer.xy .+=
            dropdims(sum(reshape(x, (predictor_dim, 1, batch_dim)) .*
            reshape(y, (1, prediction_dim, batch_dim)), dims=3), dims=3)

        layer.x .+= dropdims(sum(x, dims=2), dims=2)
        layer.y .+= dropdims(sum(y, dims=2), dims=2)
        layer.count .+= size(x, 2)
    end

    # [predictor_dim, predicton_dim, batch_dim]
    return μ, σ2
end

"""
Basically MinHash.
"""
struct LocalitySensitiveHash{
        MF <: AbstractMatrix{<:Real},
        VF <: AbstractVector{<:Real},
        VI <: AbstractVector{Int32}}

    output_dim::Int
    context_dim::Int

    # context function parameters
    hyperplanes::MF
    hyperplanes_bias::VF

    # used to compute weight indexes
    bit_offsets::VI
    k_offsets::VI
end

"""
Construct a locality sensitive hash function.

  * `output_dim`: number of outputs (i.e. number of units in this layer)
  * `context_dim`: number of context functions (inducing 2^context_dim weight vectors)
  * `predictor_dim`: dimensionality of predictors
"""
function LocalitySensitiveHash(output_dim::Int, context_dim::Int, predictor_dim::Int)
    hyperplanes = randn(Float32, (output_dim*context_dim, predictor_dim))
    hyperplanes_bias = randn(Float32, (output_dim*context_dim))

    bit_offsets = collect(Int32, 0:context_dim-1)
    k_offsets = collect(Int32, 1:output_dim)

    return LocalitySensitiveHash(
        output_dim, context_dim, hyperplanes, hyperplanes_bias, bit_offsets, k_offsets)
end

"""
Apply the context functions for this layer, mapping (standardized) input vectors
to indexes in [1, 2^context_dim].
"""
function (lsh::LocalitySensitiveHash)(X::AbstractMatrix)
    bits = (lsh.hyperplanes * X) .> lsh.hyperplanes_bias
    batch_dim = size(X, 2)
    bits_reshape = reshape(bits, (lsh.context_dim, lsh.output_dim, batch_dim))
    cs = dropdims(sum(bits_reshape .<< lsh.bit_offsets, dims=1), dims=1) .+ lsh.k_offsets

    return cs
end

"""
A single GGLN layer with arbitrary number of units.
"""
struct GGLNLayer{
        MF <: AbstractMatrix{<:Real},
        VF <: AbstractVector{<:Real},
        VI <: AbstractVector{Int32}}

    input_dim::Int
    output_dim::Int
    context_dim::Int
    predictor_dim::Int
    prediction_dim::Int

    learning_rate::Float64

    lsh::LocalitySensitiveHash{MF, VF, VI}

    weights::MF
end

"""
Construct a single GGLN layer.

  * `input_dim`: number of inputs (i.e. number of units in the prev layer)
  * `output_dim`: number of outputs (i.e. number of units in this layer)
  * `context_dim`: number of context functions (inducing 2^context_dim weight vectors)
  * `predictor_dim`: dimensionality of predictors
  * `prediction_dim`: dimensionality of predictions
  * `learning_rate`: controls step size
"""
function GGLNLayer(
        input_dim::Int, output_dim::Int,
        context_dim::Int, predictor_dim::Int, prediction_dim::Int,
        learning_rate::Float64=1e-2)

    lsh = LocalitySensitiveHash(output_dim, context_dim, predictor_dim)
    weights = fill(log(1f0/input_dim), (input_dim, output_dim*(2^context_dim)))

    return GGLNLayer(
        input_dim,
        output_dim,
        context_dim,
        predictor_dim,
        prediction_dim,
        learning_rate,
        lsh,
        weights)
end

"""
Forward pass for a single GGLN layer, training if `y` is given, predicting if not.

input_μ: [input_dim, prediction_dim, batch_dim]
input_σ2: [input_dim, prediction_dim, batch_dim]
z: [predictor_dim, batch_dim]
y: [prediction_dim, batch_dim]
"""
function forward(
        layer::GGLNLayer, input_μ::AbstractArray, input_σ2::AbstractArray,
        z::AbstractMatrix, y::Union{Nothing, AbstractMatrix})

    # input_μ, input_σ should be [input_dim, batch_dim]
    batch_dim = size(input_μ, 3)

    # [output_dim, batch_dim]
    cs = layer.lsh(z)

    function predict(weights)
        # [input_dim, output_dim, 1, batch_dim]
        weights_cs = weights[:,cs]
        penalty = 1f-5 * sum(weights_cs.^2)

        weights_ = exp.(reshape(
            weights_cs, (layer.input_dim, layer.output_dim, 1, batch_dim)))

        # [input_dim, 1, prediction_dim, batch_dim]
        input_μ_ = reshape(input_μ, (layer.input_dim, 1, layer.prediction_dim, batch_dim))
        input_σ2_ = reshape(input_σ2, (layer.input_dim, 1, layer.prediction_dim, batch_dim))

        # [output_dim, prediction_dim, batch_dim]
        σ2 = inv.(dropdims(sum(weights_ ./ input_σ2_, dims=1), dims=1))

        # [output_dim, prediction_dim, batch_dim]
        μ = σ2 .* dropdims(sum(weights_ .* input_μ_ ./ input_σ2_, dims=1), dims=1)

        return μ, σ2, penalty
    end

    if y !== nothing
        function loss(weights)
            # [output_dim, prediction_dim, batch_dim]
            μ, σ2, penalty = predict(weights)
            σ = sqrt.(σ2)

            y_ = reshape(y, (1, layer.prediction_dim, batch_dim))

            # negative normal log-pdf
            neg_ll = -sum(.- log.(σ) .- 0.5 * log.(2*π) .- 0.5 .* ((y_ .- μ) ./ σ).^2)

            return neg_ll + penalty

            # return neg_ll
        end

        dloss_dweights = gradient(loss, layer.weights)[1]
        layer.weights .-= layer.learning_rate .* dloss_dweights
        clamp!(layer.weights, -10f0, 10f0)
    end

    return predict(layer.weights)
end

"""
Complete GGLN model.
"""
struct GaussianGGLN{AF <: AbstractArray{<:Real,3}}
    x_norm::NormalizationLayer
    y_norm::NormalizationLayer
    weak_learner::BasicGaussianLinearRegression
    layers::Vector{GGLNLayer}

    # [bias_dim, prediction_dim, 1]
    μ_bias::AF
    σ2_bias::AF
end

"""
Construct an untrained gaussian gated linear network.

  * `predictor_dim`: dimensionality of predictors (x)
  * `prediction_dim`: dimensionality of predictions (y)
  * `num_layers`: number of layers in the model
  * `layer_width`: width of each layer in the model
  * `context_dim`: number of context functions (inducing 2^context_dim weight vectors)
  * `bias`: values for `bias` units.
"""
function GaussianGGLN(
        predictor_dim::Int, prediction_dim::Int,
        num_layers::Int, layer_width::Int, context_dim::Int,
        bias::Float32=5.0f0)

    bias_dim = 2*prediction_dim

    layers = GGLNLayer[]

    if num_layers > 0
        push!(
            layers, GGLNLayer(
                predictor_dim + bias_dim, layer_width, context_dim,
                predictor_dim, prediction_dim))
    end

    for i in 1:num_layers-1
        push!(
            layers, GGLNLayer(
                layer_width + bias_dim, layer_width, context_dim,
                predictor_dim, prediction_dim))
    end

    # compute the bias arrays
    μ_bias = Array{Float32}(undef, (bias_dim, prediction_dim, 1))
    μ_bias = zeros(Float32, (bias_dim, prediction_dim, 1))

    # We want a bias unit for +/- bias for every prediction dimension.
    for i in 1:prediction_dim
        μ_bias[i, i, 1] = -bias
        μ_bias[prediction_dim+i, i, 1] = bias
    end
    σ2_bias = ones(Float32,  (bias_dim, prediction_dim, 1))

    return GaussianGGLN(
        NormalizationLayer(predictor_dim),
        NormalizationLayer(prediction_dim),
        BasicGaussianLinearRegression(predictor_dim, prediction_dim),
        layers, μ_bias, σ2_bias)
end

"""
Forward pass for the GGLN, training if y is given, predicting otherwise.
"""
function forward(ggln::GaussianGGLN, x::AbstractMatrix, y::Union{Nothing, AbstractMatrix})

    batch_dim = size(x, 2)

    @assert y === nothing || size(x, 2) == size(y, 2)
    @assert y === nothing || size(y, 1) == size(ggln.μ_bias, 2)

    xz = forward(ggln.x_norm, x, y !== nothing)
    yz = y === nothing ? nothing : forward(ggln.y_norm, y, true)

    μ, σ2 = forward(ggln.weak_learner, xz, yz)

    for layer in ggln.layers
        # concatenate bias inputs
        μ = cat(repeat(ggln.μ_bias, 1, 1, batch_dim), μ, dims=1)
        σ2 = cat(repeat(ggln.σ2_bias, 1, 1, batch_dim), σ2, dims=1)

        # [output_dim, prediction_dim, batch_dim]
        μ, σ2 = forward(layer, μ, σ2, xz, yz)
    end

    return inverse(
        ggln.y_norm,
        dropdims(mean(μ, dims=1), dims=1),
        dropdims(mean(σ2, dims=1), dims=1))
end

"""
Train the model one step using a the batch of predictors `x` and predictions
`y`.

`x` should have shape [predictors dimension, batch size]
`y` should have shape [predictions dimension, batch size]
"""
function train!(ggln::GaussianGGLN, x::AbstractMatrix, y::AbstractMatrix)
    return forward(ggln, x, y)
end

"""
Make predictions using a trained model for predictors `x`.

`x` should have shape [predictors dimension, batch size]
"""
function predict(ggln::GaussianGGLN, x::AbstractMatrix)
    return forward(ggln, x, nothing)
end

end # module


text:
using Plots
# plotly()

@testset "plot with no errors" begin
    res0 = [get_interval(
        [3., 2., 2.1],
        i,
        f_3p_1im_dep,
        :LIN_EXTRAPOL;
        local_alg = :LN_NELDERMEAD,
        theta_bounds = [(-12., 12.), (-12., 12.), (-12., 12.)],
        loss_tol = 1e-3,
        loss_crit = 9.,
        silent = true
    ) for i in 1:3]
    update_profile_points!(res0[1]; max_recursions=1)
    a_grid_1 = LikelihoodProfiler.get_grid(res0[1])
    @test length(a_grid_1[2]) > 0
    p = plot(res0[1])
    @test p isa Plots.Plot
end


text:
module Contractions

# export
export fix_contractions

"""
fix_contractions: Fix Contractions in a given text

params:
------
	- text:: AbstractString , given text

usage:
------
    >>> fix_contractions(mytext)

"""
function fix_contractions(text::AbstractString)
    text = replace.(text,"ain't"=>"am not")
    text = replace.(text,"aren't"=>"are not")
    text = replace.(text,"can't"=>"cannot")
    text = replace.(text,"can't've"=>"cannot have")
    text = replace.(text,"'cause"=>"because")
    text = replace.(text,"could've"=>"could have")
    text = replace.(text,"couldn't"=>"could not")
    text = replace.(text,"couldn't've"=>"could not have")
    text = replace.(text,"didn't"=>"did not")
    text = replace.(text,"doesn't"=>"does not")
    text = replace.(text,"don't"=>"do not")
    text = replace.(text,"hadn't"=>"had not")
    text = replace.(text,"hadn't've"=>"had not have")
    text = replace.(text,"hasn't"=>"has not")
    text = replace.(text,"haven't"=>"have not")
    text = replace.(text,"he'd"=>"he would")
    text = replace.(text,"he'd've"=>"he would have")
    text = replace.(text,"he'll"=>"he will")
    text = replace.(text,"he'll've"=>"he will have")
    text = replace.(text,"he's"=>"he is")
    text = replace.(text,"how'd"=>"how did")
    text = replace.(text,"how'd'y"=>"how do you")
    text = replace.(text,"how'll"=>"how will")
    text = replace.(text,"how's"=>"how is")
    text = replace.(text,"I'd"=>"I would")
    text = replace.(text,"I'd've"=>"I would have")
    text = replace.(text,"I'll"=>"I will")
    text = replace.(text,"I'll've"=>"I will have")
    text = replace.(text,"I'm"=>"I am")
    text = replace.(text,"I've"=>"I have")
    text = replace.(text,"isn't"=>"is not")
    text = replace.(text,"it'd"=>"it had")
    text = replace.(text,"it'd've"=>"it would have")
    text = replace.(text,"it'll"=>"it will")
    text = replace.(text,"it'll've"=>"it will have")
    text = replace.(text,"it's"=>"it is")
    text = replace.(text,"let's"=>"let us")
    text = replace.(text,"ma'am"=>"madam")
    text = replace.(text,"mayn't"=>"may not")
    text = replace.(text,"might've"=>"might have")
    text = replace.(text,"mightn't"=>"might not")
    text = replace.(text,"mightn't've"=>"might not have")
    text = replace.(text,"must've"=>"must have")
    text = replace.(text,"mustn't"=>"must not")
    text = replace.(text,"mustn't've"=>"must not have")
    text = replace.(text,"needn't"=>"need not")
    text = replace.(text,"needn't've"=>"need not have")
    text = replace.(text,"o'clock"=>"of the clock")
    text = replace.(text,"oughtn't"=>"ought not")
    text = replace.(text,"oughtn't've"=>"ought not have")
    text = replace.(text,"shan't"=>"shall not")
    text = replace.(text,"sha'n't"=>"shall not")
    text = replace.(text,"shan't've"=>"shall not have")
    text = replace.(text,"she'd"=>"she would")
    text = replace.(text,"she'd've"=>"she would have")
    text = replace.(text,"she'll"=>"she will")
    text = replace.(text,"she'll've"=>"she will have")
    text = replace.(text,"she's"=>"she is")
    text = replace.(text,"should've"=>"should have")
    text = replace.(text,"shouldn't"=>"should not")
    text = replace.(text,"shouldn't've"=>"should not have")
    text = replace.(text,"so've"=>"so have")
    text = replace.(text,"so's"=>"so is")
    text = replace.(text,"that'd"=>"that would")
    text = replace.(text,"that'd've"=>"that would have")
    text = replace.(text,"that's"=>"that is")
    text = replace.(text,"there'd"=>"there had")
    text = replace.(text,"there'd've"=>"there would have")
    text = replace.(text,"there's"=>"there is")
    text = replace.(text,"they'd"=>"they would")
    text = replace.(text,"they'd've"=>"they would have")
    text = replace.(text,"they'll"=>"they will")
    text = replace.(text,"they'll've"=>"they will have")
    text = replace.(text,"they're"=>"they are")
    text = replace.(text,"they've"=>"they have")
    text = replace.(text,"to've"=>"to have")
    text = replace.(text,"wasn't"=>"was not")
    text = replace.(text,"we'd"=>"we had")
    text = replace.(text,"we'd've"=>"we would have")
    text = replace.(text,"we'll"=>"we will")
    text = replace.(text,"we'll've"=>"we will have")
    text = replace.(text,"we're"=>"we are")
    text = replace.(text,"we've"=>"we have")
    text = replace.(text,"weren't"=>"were not")
    text = replace.(text,"what'll"=>"what will")
    text = replace.(text,"what'll've"=>"what will have")
    text = replace.(text,"what're"=>"what are")
    text = replace.(text,"what's"=>"what is")
    text = replace.(text,"what've"=>"what have")
    text = replace.(text,"when's"=>"when is")
    text = replace.(text,"when've"=>"when have")
    text = replace.(text,"where'd"=>"where did")
    text = replace.(text,"where's"=>"where is")
    text = replace.(text,"where've"=>"where have")
    text = replace.(text,"who'll"=>"who will")
    text = replace.(text,"who'll've"=>"who will have")
    text = replace.(text,"who's"=>"who is")
    text = replace.(text,"who've"=>"who have")
    text = replace.(text,"why's"=>"why is")
    text = replace.(text,"why've"=>"why have")
    text = replace.(text,"will've"=>"will have")
    text = replace.(text,"won't"=>"will not")
    text = replace.(text,"won't've"=>"will not have")
    text = replace.(text,"would've"=>"would have")
    text = replace.(text,"wouldn't"=>"would not")
    text = replace.(text,"wouldn't've"=>"would not have")
    text = replace.(text,"y'all"=>"you all")
    text = replace.(text,"y'alls"=>"you alls")
    text = replace.(text,"y'all'd"=>"you all would")
    text = replace.(text,"y'all'd've"=>"you all would have")
    text = replace.(text,"y'all're"=>"you all are")
    text = replace.(text,"y'all've"=>"you all have")
    text = replace.(text,"you'd"=>"you had")
    text = replace.(text,"you'd've"=>"you would have")
    text = replace.(text,"you'll"=>"you you will")
    text = replace.(text,"you'll've"=>"you you will have")
    text = replace.(text,"you're"=>"you are")
    text = replace.(text,"you've"=>"you have")
    return text
end

end


text:
using TensorToolbox
include("LTR.jl")
include("methods/NTD.jl")
include("methods/lraSNTD.jl")
include("methods/LD/LD.jl")

function tucker_to_cprank(tucker_rank)
    if all(isequal(first(tucker_rank)),tucker_rank)
        cp = tucker_rank[1]
    else
        error("rank $tucker_rank cannot convert to CP rank")
    end
    return cp
end

function calc(X, method, reqrank)
    if method == "LT1R"
        @assert tucker_to_cprank(reqrank) == 1
        return LT1R(X)
    elseif method == "LTR"
        return LTR(X, reqrank)
    elseif method == "NTD_KL"
        _, _, Xr = NTD(X, reqrank, cost="KL", init_method=NTD_init, verbose=false)
        return Xr
    elseif method == "NTD_LS"
        _, _, Xr = NTD(X, reqrank, cost="LS", init_method=NTD_init, verbose=false)
        return Xr
    elseif method == "lraSNTD"
        return lraSNTD(X, reqrank)
    else
        error("calc method error")
    end
end

max_stars_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_issues_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_forks_repo_name: gkazunii/Legendre-tucker-rank-reduction


text:
#!/usr/local/bin/julia

using MultivariatePolynomials
using JuMP
using PolyJuMP
using SumOfSquares

using ArgParse

function try_import(name::Symbol)
    try
        @eval import $name
        return true
    catch e
        return false
    end
end

scs = try_import(:SCS)
csdp = try_import(:CSDP)

isscs(solver) = contains(string(typeof(solver)),"SCSSolver")

# Semidefinite solvers
sdp_solvers = Any[]

# Need 54000 iterations for sosdemo3 to pass on Linux 64 bits
# With 55000, sosdemo3 passes for every platform except Windows 64 bits on AppVeyor
scs && push!(sdp_solvers, SCS.SCSSolver(eps=1e-6, max_iters=60000, verbose=0))

csdp && push!(sdp_solvers, CSDP.CSDPSolver(printlevel=1, maxiter=60000, axtol=1.0e-8, atytol=1.0e-8))

function main(args)
  global sdp_solvers

  s = ArgParseSettings("Example 4 for argparse.jl: " *
                         "more tweaking of the arg fields: " *
                         "dest_name, metvar, range_tested, " *
                         "alternative actions.")

  @add_arg_table s begin
    "--opt1"
      action = :append_const   # appends 'constant' to 'dest_name'
      arg_type = ByteString    # the only utility of this is restricting the dest array type
      constant = "O1"
      dest_name = "O_stack"    # this changes the destination
      help = "append O1"
  end

  s.epilog = """
    freeSir.jl
    """

  parsed_args = parse_args(args, s)
  println("Parsed args:")
  for (key,val) in parsed_args
    println("  $key  =>  $(repr(val))")
  end

  println("calling the solvers")
  for solver in sdp_solvers
    @polyvar x[1:2]
    μ = 0.2
    β = 0.5
    γ = 0.6
    s = 1.
    i = 0.
    R0 = β / ( μ + γ )

    # Constructing the vector field dx/dt = f

    f = [μ - β * (x[1] + s) * (x[2] + i) - μ * (x[1] + s),
        β * (x[1] + s) * (x[2] + i) - (μ + γ) * (x[2] + i)]

    #f = [0.06 - 0.01 * (x[1] + 1.) * x[2] - 0.06 * (x[1] + 1.),
    #    0.01 * (x[1] + 1.) * x[2] - (0.06 + 0.5) * x[2]]

    #println(f)

    m = Model(solver = solver)

    # The Lyapunov function V(x):
    Z = x.^2
    @polyvariable m V Z

    @polyconstraint m V >= 1e-16 * sum(x.^2)

    # dV/dx*f <= 0
    P = dot(differentiate(V, x), f)
    @polyconstraint m P <= 0
    @polyconstraint m x[1] * ( x[1] + 1) <= 0
    @polyconstraint m x[2] * ( x[2] - 1) <= 0

    println(x)

    println("este es el modelo final")
    print(m)
    println("valor de la funcion objetivo", getobjectivevalue(m))

    status = solve(m);

    flyapunov = getvalue(V)
    println("Lyapunov function using $(typeof(solver))")

    println(flyapunov)
    println("Derivada Orbital")
    println(dot(differentiate(flyapunov, x), f))
    println(status)

    println("Basic Reproductive Number")
    println(R0)

      #status --> :Optimal

    #removemonomials(getvalue(V), Z) --> zero(Polynomial{true, Float64})
  end
end

main(ARGS)

max_stars_repo_name: polislizarralde/jsos

max_issues_repo_name: polislizarralde/jsos

max_forks_repo_name: polislizarralde/jsos


text:
using DrWatson
@quickactivate("HyperspectraWithNeXL")
# Load the necessary libraries
using HyperspectraWithNeXL
using NeXLSpectrum
using Unitful
using FileIO, ImageIO

# Load the HyperSpectrum from disk
lt = 0.72*4.0*18.0*3600.0/(1024*1024) # 18.0 hours on 4 detectors
hs = NeXLSpectrum.compress(HyperSpectrum(
    LinearEnergyScale(0.0,10.0),
    Dict{Symbol,Any}(
      :TakeOffAngle => deg2rad(35.0),
      :ProbeCurrent => 1.0,
      :LiveTime => lt,
      :BeamEnergy => 20.0e3,
      :Name => "Mn Nodule"
    ),
    readrplraw(joinpath(datadep"MnNodule","map[15]")),
    fov = [ 4.096u"mm", 4.096u"mm"], offset= [ 0.0u"mm", 0.0u"mm" ]
))
outpath = joinpath(papersdir(),"Figures","Figure 5")
mkpath(outpath)
cp(joinpath(datadep"MnNodule","Image[0][[1]].png"), joinpath(outpath,"Mn_nodule_BSED.png"), force=true)

# Independent normalization (to brightest pixel)
outpath = joinpath(papersdir(),"Figures","Figure 6")
mkpath(outpath)
FileIO.save(File{format"PNG"}(joinpath(outpath,"I[C K-L2].png")), hs[n"Mn K-L3"])
FileIO.save(File{format"PNG"}(joinpath(outpath,"I[Mn K-L3].png")), hs[n"Mn K-L3"])
FileIO.save(File{format"PNG"}(joinpath(outpath,"I[Fe K-L3].png")), hs[n"Fe K-L3"])
FileIO.save(File{format"PNG"}(joinpath(outpath,"I[O K-L3].png")), hs[n"O K-L3"])
# Normalized as a set (normalize to sum of intensities)

outpath=plotsdir()
mkpath(outpath)
cxrs = [n"Mn K-L3", n"O K-L3", n"Fe K-L3" ]
imgs = hs[ cxrs ]
for (cxr, img) in zip(cxrs, imgs)
  FileIO.save(File{format"PNG"}(joinpath(outpath,"I_rel[$cxr].png")), img)
end
# RGB
outpath = joinpath(papersdir(),"Figures","Figure 7")
mkpath(outpath)
img = colorize(hs, cxrs, :All)
FileIO.save(File{format"PNG"}(joinpath(outpath,"colorized[Mn,O,Fe,All].png")), img)
img = colorize(hs, cxrs, :Each)
FileIO.save(File{format"PNG"}(joinpath(outpath,"colorized[Mn,O,Fe,Each].png")), img)


text:
@testset "Levenberg-Marquardt: NIST Dataset (BoxBOD)" begin
    print_head("Levenberg-Marquardt: NIST Dataset (BoxBOD)")

    obj = optimizer(nd=2, ny=6, method="lm")
    fun = BoxBOD()
    sol = Vector{Float64}(undef, 2)

    parμ = [10.0, 5.5]
    βvec = zeros(Float64, 6)
    βmat = zeros(Float64, 6, 6)

    print_body("\033[1m\033[34mLikelihood Precision Matrix: Identity\033[0m")

    @simd for i in eachindex(βvec)
        @inbounds βvec[i] = 1.0
    end

    @simd for i in eachindex(βvec)
        @inbounds βmat[i,i] = 1.0
    end

    ans = [213.809409, 0.54723748]

    print_body(string("Precision Const.: ", optimize!(sol, fun, parμ, 1.0, obj)))
    @inbounds for i in eachindex(ans)
        @test sol[i] ≈ ans[i] rtol=1e-7
    end

    print_body(string("Precision Vector: ", optimize!(sol, fun, parμ, βvec, obj)))
    @inbounds for i in eachindex(ans)
        @test sol[i] ≈ ans[i] rtol=1e-7
    end

    print_body(string("Precision Matrix: ", optimize!(sol, fun, parμ, βmat, obj)))
    @inbounds for i in eachindex(ans)
        @test sol[i] ≈ ans[i] rtol=1e-7
    end

    print_body("\033[1m\033[34mLikelihood Precision Matrix: Custom\033[0m")

    @simd for i in eachindex(βvec)
        @inbounds βvec[i] = 1.0 + 0.1 * (i-1)
    end

    @simd for i in eachindex(βvec)
        @inbounds βmat[i,i] = 1.0 + 0.1 * (i-1)
    end

    @inbounds ans[1], ans[2] = 216.132114, 0.52157598

    print_body(string("Precision Vector: ", optimize!(sol, fun, parμ, βvec, obj)))
    @inbounds for i in eachindex(ans)
        @test sol[i] ≈ ans[i] rtol=1e-7
    end

    print_body(string("Precision Matrix: ", optimize!(sol, fun, parμ, βmat, obj)))
    @inbounds for i in eachindex(ans)
        @test sol[i] ≈ ans[i] rtol=1e-7
    end
end


text:
using DataFrames  # readtable function
# using Requests
type Bus
   nodeID::Int
   root::Int
   Pd::Float64
   Qd::Float64
   Vmax::Float64
   Vmin::Float64
   B::Float64
   R::Float64
   X::Float64
   children::Vector{Int}
   ancestor::Vector{Int}
   genids::Int
   function Bus(nodeID, root, Pd, Qd, B, R, X, Vmax, Vmin)
      b = new(nodeID, root, Pd, Qd)
      b.Vmax = Vmax
      b.Vmin = Vmin
      b.R = R
      b.X = X
      b.B = B
      b.children = Int[]
      b.ancestor = Int[]
      # b.genids = 0
      return b
   end
end
#################################################################
type Generator
   genID::Int
   busidx::Int
   Pgmax::Float64
   Pgmin::Float64
   Qgmax::Float64
   Qgmin::Float64
   cost::Float64
   function Generator( busidx, Pgmax, Pgmin, Qgmax, Qgmin, cost)
      g = new(busidx)
      g.busidx = busidx
      g.cost = cost
      g.Pgmax = Pgmax
      g.Pgmin = Pgmin
      g.Qgmax = Qgmax
      g.Qgmin = Qgmin
      return g
   end
end

##################################################################
type Line
   arcID::Int
   tail::Int # the "to" node
   head::Int # the "from" node
   r::Float64 # the resistance value
   x::Float64 # the reactance value
   u::Float64 # the capacity of the line
   function Line(arcID, tail, head, r, x, u)
      line = new(arcID, tail, head, r, x)
      line.u = u
      return line
   end
end
#########################################
## storage data
type Storage
   ID::Int
   busidx::Int
   Pkmax::Float64
   Ekmax::Float64
   Aleph::Float64
   Co::Float64
   Cp::Float64
   Ce::Float64
   function Storage(ID, busidx, Pkmax, Ekmax, Aleph, Co, Cp, Ce)
      s = new(ID)
      s.ID = ID
      s.busidx = busidx
      s.Pkmax = Pkmax
      s.Ekmax = Ekmax
      s.Aleph = Aleph
      s.Co = Co
      s.Cp = Cp
      s.Ce = Ce
      return s
   end
end
#########################################
## station data
type Station
   ID::Int
   busidx::Int
   Cf::Float64
   function Station(ID, busidx, Cf)
      st = new(ID)
      st.ID = ID
      st.busidx = busidx
      st.Cf = Cf
      return st
   end
end
######################################################################################################

function DataImport(filename_Node, filename_Generator, filename_Line, filename_Storage, filename_Station, filename_SMP)

   ######################################################################################################
   # Bus/Node Data
   # busmat = readtable(filename_Node)
   busmat = readcsv(filename_Node, header=true)[1]

   buses = Bus[]
   busIDmap = Dict()
   for i in 1:size(busmat,1)

      nodeID = i
      busIDmap[busmat[i,1]] = i

      if i==1
         root = busIDmap[busmat[i,1]]
      else
         root = 0
      end

      Pd = busmat[i,2]
      Qd = busmat[i,3]
      R = busmat[i,6]
      X = busmat[i,7]

      Vmax = busmat[i,4]
      Vmin = busmat[i,5]
      B = busmat[i,8]

      b = Bus(nodeID, root, Pd, Qd, B, R, X, Vmax, Vmin)
      push!(buses, b)
   end

   #######################################################################################################
   ## generator data
   generatorlist = Int[]
   generators = Generator[]
   # genmat = readtable(filename_Generator)
   genmat = readcsv(filename_Generator, header=true)[1]

   for i in 1:size(genmat,1)
      busidx = genmat[i,1]
      Pgmax = genmat[i,2]
      Pgmin = genmat[i,3]
      Qgmax = genmat[i,4]
      Qgmin = genmat[i,5]
      cost = genmat[i,6]

      g = Generator(busidx, Pgmax, Pgmin, Qgmax, Qgmin, cost)
      push!(generators, g)
      # setg(buses[busidx], i)
   end
   #for g in 1:length(generators)
   #   generators[g].cost = genmat[g,4]
   #end
   #################################################################
   ## branch data
   # branchmat = readtable(filename_Line)
   branchmat = readcsv(filename_Line, header=true)[1]

   lines = Line[]
   for i in 1:size(branchmat,1)
      fbus = busIDmap[branchmat[i,2]]
      tbus = busIDmap[branchmat[i,1]]
      abus = busIDmap[branchmat[i,1]]
      x = branchmat[i,4]
      r = branchmat[i,3]
      u = branchmat[i,6]  # flow limit
      push!(buses[tbus].children, fbus)#children
      push!(buses[fbus].ancestor, abus)#ancestor
      l = Line(i, tbus, fbus, r, x, u)
      push!(lines,l)
   end
   #########################################
   ## storage data
   # storagemat = readtable(filename_Storage)
   storagemat = readcsv(filename_Storage, header=true)[1]

   storages = Storage[]
   for i in 1:size(storagemat,1)
      storageID = i
      busidx = 1
      Pkmax = storagemat[i,2]
      Ekmax = storagemat[i,3]
      Aleph = storagemat[i,4]
      Co = storagemat[i,5]
      Cp = storagemat[i,6]
      Ce = storagemat[i,7]

      temp = Storage(storageID, busidx, Pkmax, Ekmax, Aleph, Co, Cp, Ce)
      push!(storages,temp)
   end

   #########################################
   ## charging/discharing station data
   # stationmat = readtable(filename_Station)
   stationmat = readcsv(filename_Station, header=true)[1]
   stations = Station[]
   for i in 1:size(stationmat,1)
      stationID = i
      busidx = stationmat[i,2]
      Cf = stationmat[i,3]

      temp = Station(stationID, busidx, Cf)
      push!(stations,temp)
   end

   #########################################
   SMP_raw = readcsv(filename_SMP, header=true)[1]
   SMP = SMP_raw[3:end]

   return buses, generators,lines, storages, stations, SMP
end


text:
# SPDX-License-Identifier: X11
# 2020-12-07
# Day 2, Part 1

const reg = r"^(\d*)-(\d*) ([a-z]): ([a-z]*)"

function main()
  input = readlines()

  cnt = 0
  for x ∈ input
    validpass(x) && (cnt += 1)
  end

  println(cnt)
end

function validpass(str::AbstractString)::Bool
  m = match(reg, str)

  min  = parse.(Int, m.captures[1])
  max  = parse.(Int, m.captures[2])
  chr  = m.captures[3][1]
  pass = m.captures[end]

  cnt = 0
  for c ∈ pass
    c == chr && (cnt += 1)
  end

  return cnt >= min && cnt <= max
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
using Docile, Lexicon, JuliaFEM

const api_directory = "api"

"""
Searches recursively all the modules from packages. As documentation grows, it's a bit
troublesome to add all the new modules manually, so this function searches all the modules
automatically.

Parameters
----------
    module_: Module
        Module where we want to search modules inside
    append_list: Array{Module, 1}
        Array, where we append Modules as we find them

Returns
-------
    None. Void function, which manipulates the append_list
"""
function search_modules!(module_::Module, append_list::Array{Module, 1})
    all_names = names(module_, true)
    for each in all_names
        inner_module = module_.(each)
        if (typeof(inner_module) == Module) && !(inner_module in append_list)
            push!(append_list, inner_module)
            search_modules!(inner_module, append_list)
        end
    end
end

append_list = Array(Module, 0)
search_modules!(JuliaFEM, append_list)

const modules = append_list

# main_folder = dirname(dirname(@__FILE__))
# this_folder = dirname(@__FILE__)

# file_ = "README.md"
# run(`cp $main_folder/$file_ $this_folder`)

cd(dirname(@__FILE__)) do
    # Run the doctests *before* we start to generate *any* documentation.
#    for m in modules
#        failures = failed(doctest(m))
#        if !isempty(failures.results)
#            println("\nDoctests failed, aborting commit.\n")
#            display(failures)
#            exit(1) # Bail when doctests fail.
#        end
#    end
    # Generate and save the contents of docstrings as markdown files.
    index  = Index()
    for mod in modules
        Lexicon.update!(index, save(joinpath(api_directory, "$(mod).rst"), mod))
    end
    save(joinpath(api_directory, "index.rst"), index)

    # Add a reminder not to edit the generated files.
  #  open(joinpath(api_directory, "README.md"), "w") do f
  #      print(f, """
  #      Files in this directory are generated using the `build.jl` script. Make
  #      all changes to the originating docstrings/files rather than these ones.
  #      """)
  #  end
    # save(joinpath(api_directory, "index.rst"), index; md_subheader = :category)
    # info("Adding all documentation changes in $(api_directory) to this commit.")
    # success(`git add $(api_directory)`) || exit(1)

end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

# ValueType describes the semantics and measurement units of a value.
mutable struct ValueType <: ProtoType
    _type::Int64  # Index into string table.
    unit::Int64  # Index into string table.
    ValueType(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ValueType

mutable struct Label <: ProtoType
    key::Int64  # Index into string table
    # At most one of the following must be present
    str::Int64  # Index into string table
    num::Int64
    # Should only be present when num is present.
    # Specifies the units of num.
    # Use arbitrary string (for example, "requests") as a custom count unit.
    # If no unit is specified, consumer may apply heuristic to deduce the unit.
    # Consumers may also  interpret units like "bytes" and "kilobytes" as memory
    # units and units like "seconds" and "nanoseconds" as time units,
    # and apply appropriate unit conversions to these.
    num_unit::Int64  # Index into string table
    Label(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Label

# Each Sample records values encountered in some program
# context. The program context is typically a stack trace, perhaps
# augmented with auxiliary information like the thread-id, some
# indicator of a higher level request being handled etc.
mutable struct Sample <: ProtoType
    # The ids recorded here correspond to a Profile.location.id.
    # The leaf is at location_id[0].
    location_id::Base.Vector{UInt64}
    # The type and unit of each value is defined by the corresponding
    # entry in Profile.sample_type. All samples must have the same
    # number of values, the same as the length of Profile.sample_type.
    # When aggregating multiple samples into a single sample, the
    # result has a list of values that is the elemntwise sum of the
    # lists of the originals.
    value::Base.Vector{Int64}
    # label includes additional context for this sample. It can include
    # things like a thread id, allocation size, etc
    label::Base.Vector{Label}
    Sample(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Sample
const __pack_Sample = Symbol[:location_id,:value]
meta(t::Type{Sample}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_Sample, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct Mapping <: ProtoType
    # Unique nonzero id for the mapping.
    id::UInt64
    # Address at which the binary (or DLL) is loaded into memory.
    memory_start::UInt64
    # The limit of the address range occupied by this mapping.
    memory_limit::UInt64
    # Offset in the binary that corresponds to the first mapped address.
    file_offset::UInt64
    # The object this entry is loaded from.  This can be a filename on
    # disk for the main binary and shared libraries, or virtual
    # abstractions like "[vdso]".
    filename::Int64  # Index into string table
    # A string that uniquely identifies a particular program version
    # with high probability. E.g., for binaries generated by GNU tools,
    # it could be the contents of the .note.gnu.build-id field.
    build_id::Int64  # Index into string table
    # The following fields indicate the resolution of symbolic info.
    has_functions::Bool
    has_filenames::Bool
    has_line_numbers::Bool
    has_inline_frames::Bool
    Mapping(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Mapping

mutable struct Line <: ProtoType
    # The id of the corresponding profile.Function for this line.
    function_id::UInt64
    # Line number in source code.
    line::Int64
    Line(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Line

# Describes function and line table debug information.
mutable struct Location <: ProtoType
    # Unique nonzero id for the location.  A profile could use
    # instruction addresses or any integer sequence as ids.
    id::UInt64
    # The id of the corresponding profile.Mapping for this location.
    # It can be unset if the mapping is unknown or not applicable for
    # this profile type.
    mapping_id::UInt64
    # The instruction address for this location, if available.  It
    # should be within [Mapping.memory_start...Mapping.memory_limit]
    # for the corresponding mapping. A non-leaf address may be in the
    # middle of a call instruction. It is up to display tools to find
    # the beginning of the instruction if necessary.
    address::UInt64
    # Multiple line indicates this location has inlined functions,
    # where the last entry represents the caller into which the
    # preceding entries were inlined.
    #
    # E.g., if memcpy() is inlined into printf:
    #    line[0].function_name == "memcpy"
    #    line[1].function_name == "printf"
    line::Base.Vector{Line}
    # Provides an indication that multiple symbols map to this location's
    # address, for example due to identical code folding by the linker. In that
    # case the line information above represents one of the multiple
    # symbols. This field must be recomputed when the symbolization state of the
    # profile changes.
    is_folded::Bool
    Location(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Location

mutable struct Function <: ProtoType
    # Unique nonzero id for the function.
    id::UInt64
    # Name of the function, in human-readable form if available.
    name::Int64  # Index into string table
    # Name of the function, as identified by the system.
    # For instance, it can be a C++ mangled name.
    system_name::Int64  # Index into string table
    # Source file containing the function.
    filename::Int64  # Index into string table
    # Line number in source file.
    start_line::Int64
    Function(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Function

mutable struct Profile <: ProtoType
    # A description of the samples associated with each Sample.value.
    # For a cpu profile this might be:
    #   [["cpu","nanoseconds"]] or [["wall","seconds"]] or [["syscall","count"]]
    # For a heap profile, this might be:
    #   [["allocations","count"], ["space","bytes"]],
    # If one of the values represents the number of events represented
    # by the sample, by convention it should be at index 0 and use
    # sample_type.unit == "count".
    sample_type::Base.Vector{ValueType}
    # The set of samples recorded in this profile.
    sample::Base.Vector{Sample}
    # Mapping from address ranges to the image/binary/library mapped
    # into that address range.  mapping[0] will be the main binary.
    mapping::Base.Vector{Mapping}
    # Useful program location
    location::Base.Vector{Location}
    # Functions referenced by locations
    _function::Base.Vector{Function}
    # A common table for strings referenced by various messages.
    # string_table[0] must always be "".
    string_table::Base.Vector{AbstractString}
    # frames with Function.function_name fully matching the following
    # regexp will be dropped from the samples, along with their successors.
    drop_frames::Int64  # Index into string table.
    # frames with Function.function_name fully matching the following
    # regexp will be kept, even if it matches drop_functions.
    keep_frames::Int64  # Index into string table.

    # The following fields are informational, do not affect
    # interpretation of results.

    # Time of collection (UTC) represented as nanoseconds past the epoch.
    time_nanos::Int64
    # Duration of the profile, if a duration makes sense.
    duration_nanos::Int64
    # The kind of events between sampled ocurrences.
    # e.g [ "cpu","cycles" ] or [ "heap","bytes" ]
    period_type::ValueType
    # The number of events between sampled occurrences.
    period::Int64
    # Freeform text associated to the profile.
    comment::Base.Vector{Int64}  # Indices into string table.
    # Index into the string table of the type of the preferred sample
    # value. If unset, clients should default to the last sample value.
    default_sample_type::Int64
    Profile(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Profile
const __pack_Profile = Symbol[:comment]
meta(t::Type{Profile}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_Profile, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

export Profile, ValueType, Sample, Label, Mapping, Location, Line, Function


text:
using JLD2, DataFrames, Serialization, LinearAlgebra, Distributions

import GEAS:breeding, sim_QTL, summarize, create_storage

function gpar(nqtl)
    Dict(
        :nSire=> 100,
        :nDam => 200,   # male:female = 1:2
        :nSib => 40,
        :nC7e => 30,    # number of sib to be challenged
        :nG8n => 11,    # number of generations
        :nQTL => [nqtl, nqtl],
        :h²   => [.5, .5],
        :p8e  => .5,    # percentage of death after challenge
        :e19e => 1.,    # edit_successsful_rate
        :w4t  => 1,     # weight on the binary trait EBV
        :nK3n => 10,    # number of known QTL on binary trait
        :dd   => 0.01,  # A small value to be added to diagonals
        :edit => false,
        :fix  => false
    )
end

function t_2021_09_12(nrpt = 10)
    BLAS.set_num_threads(12)
    nqtl = [100, 500]
    isdir("dat/tmp") || mkpath("dat/tmp")
    @load "dat/run/base.jld2" base
    qdist = [Normal(), Laplace()]
    rst = DataFrame()
    ng8n = 11

    for nq in nqtl
        @info "With n_QTL = $nq"
        dpr = gpar(nq)  # parameter dictionary
        for qd in qdist
            @info "QTL distribution $(string(qd))"
            for rpt in 1:nrpt
                @info "Repeat $rpt of $nrpt"
                @info "Preparing common base + generation 1"
                dpr[:fix], dpr[:edit] = false, false
                par = (; dpr...)
                qtl = sim_QTL(base, par.nQTL..., d=qd)
                snp, ped = create_storage(base, par, qtl)
                serialize("dat/tmp/snp.ser", snp)
                serialize("dat/tmp/ped.ser", ped)

                @info "Selection with SNP-BLUP"
                par = (; dpr...)
                breeding(ped, snp, base, qtl, par)
                df = summarize(ped.prd, snp.prd, qtl[2])
                df.repeat = repeat([rpt], ng8n)
                df.dist = repeat([string(qd)[1:3]], ng8n)
                df.nq = repeat([nq], ng8n)
                df.method = repeat(["SBLP"], ng8n)
                append!(rst, df)

                @info "Selection with emphasis on top $(par.nK3n) QTL"
                dpr[:fix], dpr[:edit] = true, false
                par = (; dpr...)
                snp = deserialize("dat/tmp/snp.ser")
                ped = deserialize("dat/tmp/ped.ser")
                breeding(ped, snp, base, qtl, par)
                df = summarize(ped.prd, snp.prd, qtl[2])
                df.repeat = repeat([rpt], ng8n)
                df.dist = repeat([string(qd)[1:3]], ng8n)
                df.nq = repeat([nq], ng8n)
                df.method = repeat(["Fix"], ng8n)
                append!(rst, df)

                dpr[:fix], dpr[:edit] = false, true
                par = (; dpr...)
                snp = deserialize("dat/tmp/snp.ser")
                ped = deserialize("dat/tmp/ped.ser")
                breeding(ped, snp, base, qtl, par)
                df = summarize(ped.prd, snp.prd, qtl[2])
                df.repeat = repeat([rpt], ng8n)
                df.dist = repeat([string(qd)[1:3]], ng8n)
                df.nq = repeat([nq], ng8n)
                df.method = repeat(["Edit"], ng8n)
                append!(rst, df)
                serialize("dat/rst/3w-cmp.ser", rst)
            end
        end
    end
end

function shortcut_2021_09_12()
    BLAS.set_num_threads(12)
    dpr = gpar(100)
    par = (; dpr...)
    @load "dat/run/base.jld2" base
    qtl = sim_QTL(base, par.nQTL...)
    snp, ped = create_storage(base, par, qtl)
    breeding(ped, snp, base, qtl, par)
    serialize("dat/tmp/sample.ser", Dict(:snp=>snp, :ped=>ped, :qtl=>qtl))
end


text:
"""
T Macrina
160310

Change to directory with the following files:
* segmentation.h5
* classification.csv
* class_description.csv

Run
`julia semantic.jl`
"""

using HDF5

"""
`CREATE_SEMANTIC_MASK` - create indexed image for semantic class of every voxel

Args:

* indexed_matrix: matrix indexed with segment IDs
* segment_class: dictionary, keys are segment IDs, values are the class
* uncertain_id: class index for a segment ID that was not classified

Returns:

* semantic_mask: matrix indexed with class IDs
"""
function create_semantic_mask(indexed_matrix, segment_class, uncertain_id=0)
	semantic_mask = ones(UInt8, size(indexed_matrix)...)
	for i in 1:length(indexed_matrix)
		k = indexed_matrix[i]
		if haskey(segment_class, k)
			semantic_mask[i] = UInt8(segment_class[k])
		else
			semantic_mask[i] = uncertain_id
		end
	end
	return semantic_mask
end

"""
`CHECK_SEGMENT_IDS` - Check that all segments have been classified
"""
function check_segment_ids(indexed_matrix, segment_class)
	segments_marked = Set(unique(indexed_matrix))
	segments_classified = Set(keys(segment_class))
	println("Marked but not classified: ", length(setdiff(segments_marked,
														segments_classified)))
	println("Classified but not marked: ", length(setdiff(segments_classified,
															segments_marked)))
	return assert(Set(segment_ids) == Set(keys(segment_class)))
end

"""
`WRITE_SEMANTIC_MASK` - create indexed image for semantic class of every voxel

Args:

* dir: path to the folder containing the following files
** segmentation.h5: the segment indexed volume with attribute "main"
** classification.csv: 2 column table - segment id, class id
** class_description.csv: 2 column table - class id, class description

Returns:

Writes out H5 file, semantic_mask - class id indexed image
* "main": matrix indexed with class ids &
* "class_id": table of class ids
* "class_description": table of class descriptions (match class ids)
"""
function write_semantic_mask(dir)
	segmentation_fn = joinpath(dir, "segmentation.h5")
	classification_fn = joinpath(dir, "classification.csv")
	class_description_fn = joinpath(dir, "class_description.csv")
	semantic_mask_fn = joinpath(dir, "semantic_mask.h5")

	indexed_matrix = h5read(segmentation_fn, "main")
	segment_class = convert_table_to_dict(readdlm(classification_fn, Int))
	semantic_mask = create_semantic_mask(indexed_matrix, segment_class)
	classes = readdlm(class_description_fn)
	class_id = map(UInt8, classes[:,1])
	class_description = map(String, classes[:,2])

    f = h5open(semantic_mask_fn, "w")
    f["main"] = semantic_mask
    f["class_id"] = class_id
    f["class_description"] = class_description
    close(f)
end

function convert_table_to_dict(table)
	d = Dict()
	for i in 1:size(table,1)
		d[table[i,1]] = table[i,2]
	end
	return d
end

if !isinteractive()
	write_semantic_mask(pwd())
end


text:
using DataFrames
using CSV
using FileIO
using Statistics

include("experiments/experimentalutils.jl")

resultsFolder = mainfolder * "experiments/twostagesvae/"
files = readdir(resultsFolder)

results = []
for f in files
    if isfile(resultsFolder * f)
        push!(results, CSV.read(resultsFolder * f))
    end
end
results = vcat(results...)

aggres = []
for d in unique(results[:dataset])
    ddf = results[results[:dataset] .== d, :]
    mean_auc_pxv = mean(ddf[:auc_pxv])
    mean_auc_pz = mean(ddf[:auc_pz])
    push!(aggres, DataFrame(dataset = d, auc_pxv = mean_auc_pxv, auc_pz = mean_auc_pz))
end
aggres = vcat(aggres...)


text:
module DataFittingBasics

import DataFitting: AbstractDomain, Domain_1D, Domain_2D,
    Parameter, AbstractComponent, AbstractComponentData,
    cdata, evaluate!

include("OffsetSlope.jl")
include("Polynomial.jl")
include("Gaussian.jl")
include("Lorentzian.jl")
end


text:
#=
Visual illustration of the law of large numbers.

@author : Spencer Lyon <spencer.lyon@nyu.edu>
          Victoria Gregory <victoria.gregory@nyu.edu>

References
----------

Based off the original python file illustrates_lln.py
=#
using Plots
pyplot()
using Distributions
using LaTeXStrings

n = 100
srand(42)  # reproducible results

# == Arbitrary collection of distributions == #
distributions = Dict("student's t with 10 degrees of freedom" => TDist(10),
                 "beta(2, 2)" => Beta(2.0, 2.0),
                 "lognormal LN(0, 1/2)" => LogNormal(0.5),
                 "gamma(5, 1/2)" => Gamma(5.0, 2.0),
                 "poisson(4)" => Poisson(4),
                 "exponential with lambda = 1" => Exponential(1))

num_plots = 3
dist_data = zeros(num_plots, n)
sample_means = []
dist_means = []
titles = []
for i = 1:num_plots
    dist_names = collect(keys(distributions))
    # == Choose a randomly selected distribution == #
    name = dist_names[rand(1:length(dist_names))]
    dist = pop!(distributions, name)

    # == Generate n draws from the distribution == #
    data = rand(dist, n)

    # == Compute sample mean at each n == #
    sample_mean = Array(Float64, n)
    for j=1:n
        sample_mean[j] = mean(data[1:j])
    end

    m = mean(dist)

    dist_data[i, :] = data'
    push!(sample_means, sample_mean)
    push!(dist_means, m*ones(n))
    push!(titles, name)

end

# == Plot == #
N = repmat(reshape(repmat(linspace(1, n, n), 1, num_plots)', 1, n*num_plots), 2, 1)
heights = [zeros(1,n*num_plots); reshape(dist_data, 1, n*num_plots)]
plot(N, heights, layout=(3, 1), label="", color=:grey, alpha=0.5)
plot!(1:n, dist_data', layout=(3, 1), color=:grey, markershape=:circle,
        alpha=0.5, label="", linewidth=0)
plot!(1:n, sample_means, linewidth=3, alpha=0.6, color=:green, legend=:topleft,
      layout=(3, 1), label=[LaTeXString("\$\\bar{X}_n\$") "" ""])
plot!(1:n, dist_means, color=:black, linewidth=1.5, layout=(3, 1),
      linestyle=:dash, grid=false, label=[LaTeXString("\$\\mu\$") "" ""])
plot!(title=titles')

max_stars_repo_name: mwhchen/quantecon

max_issues_repo_name: mwhchen/quantecon

max_forks_repo_name: mwhchen/quantecon


text:
mutable struct DynamicBuffer{O,A} <: AbstractBuffer
    observations::Vector{O}
    actions::Vector{A}
    rewards::Vector{Float64}
    terminals::Vector{Bool}
    episodes::Int

    DynamicBuffer{O,A}() where {O,A} = new(Vector{O}(), A[], Float64[], Bool[], 0)
end

DynamicBinaryBuffer{A} = DynamicBuffer{Vector{Bool},A}

function add!(buffer::DynamicBuffer{O,A}, observation::O) where {O,A}
    (length(buffer.observations) != length(buffer.actions)) && !buffer.terminals[end] && error(
      "Incomplete transition"
    )

    push!(buffer.observations, observation)
    push!(buffer.terminals, false)
    nothing
end

function add!(buffer::DynamicBuffer{O,A}, action::A, reward::Float64, next_obs::O, terminal::Bool) where {O,A}
    length(buffer.observations) == length(buffer.actions) && error(
      "Incomplete transition"
    )

    push!(buffer.actions, action)
    push!(buffer.rewards, reward)
    push!(buffer.observations, next_obs)
    push!(buffer.terminals, terminal)

    terminal && (buffer.episodes += 1)
    nothing
end

function reset!(buffer::DynamicBuffer)
    empty!(buffer.observations)
    empty!(buffer.actions)
    empty!(buffer.rewards)
    empty!(buffer.terminals)

    buffer.episodes = 0

    nothing
end

function Base.iterate(buffer::DynamicBuffer, state=1)
    state > length(buffer.actions) && return nothing

    return ((buffer.observations[state],
        buffer.actions[state],
        buffer.rewards[state],
        buffer.observations[state + 1],
        buffer.terminals[state + 1]),
      state + 1)
end

function Base.iterate(rbuffer::Base.Iterators.Reverse{DynamicBuffer{O,A}}, state=length(rbuffer.itr.actions)) where {O,A}
    state < 1 && return nothing

    return ((rbuffer.itr.observations[state],
        rbuffer.itr.actions[state],
        rbuffer.itr.rewards[state],
        rbuffer.itr.observations[state + 1],
        rbuffer.itr.terminals[state + 1]),
      state - 1)
end

Base.length(buffer::DynamicBuffer) = length(buffer.actions)
Base.eltype(buffer::DynamicBuffer{O,A}) where {O,A} = Tuple{O,A,Float64,O,Bool}