text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestLP

import MathOptInterface
using Test

const MOI = MathOptInterface
const LP = MOI.FileFormats.LP
const LP_TEST_FILE = "test.lp"

function test_show()
    @test sprint(show, LP.Model()) == "A .LP-file model"
    return
end

function test_comprehensive_write()
    model = LP.Model()
    MOI.Utilities.loadfromstring!(
        model,
        """
variables: a, x, y, z
minobjective: x
x >= -1.0
x <= 2.0
y == 3.0
z in Interval(4.0, 5.0)
c5: 1.1x + 0.0 <= 5.1
c6: 1.3x + -1.4 >= -0.1
c7: 1.5a + 1.6 == 0.2
c8: 1.7a + 1.8 in Interval(0.3, 0.4)
x in ZeroOne()
y in Integer()
c11: [x, y, z] in SOS1{Float64}([1.0, 2.0, 3.0])
c12: [x, y, z] in SOS2{Float64}([3.3, 1.1, 2.2])
""",
    )
    MOI.write_to_file(model, LP_TEST_FILE)
    @test read(LP_TEST_FILE, String) ==
          "minimize\n" *
          "obj: x\n" *
          "subject to\n" *
          "c5: 1.1 x <= 5.1\n" *
          "c6: -1.4 + 1.3 x >= -0.1\n" *
          "c7: 1.6 + 1.5 a = 0.2\n" *
          "c8: 0.3 <= 1.8 + 1.7 a <= 0.4\n" *
          "Bounds\n" *
          "x <= 2\n" *
          "x >= -1\n" *
          "y = 3\n" *
          "4 <= z <= 5\n" *
          "a free\n" *
          "General\n" *
          "y\n" *
          "Binary\n" *
          "x\n" *
          "SOS\n" *
          "c11: S1:: x:1.0 y:2.0 z:3.0\n" *
          "c12: S2:: x:3.3 y:1.1 z:2.2\n" *
          "End\n"

    @test !MOI.is_empty(model)
    MOI.empty!(model)
    @test MOI.is_empty(model)
end

function test_name_sanitization_start()
    for starting_letter in [".", "0", "E", "e"]
        model = LP.Model()
        MOI.Utilities.loadfromstring!(
            model,
            """
variables: x
minobjective: x
c1: 1.0 * x >= -1.0
""",
        )
        x = MOI.get(model, MOI.VariableIndex, "x")
        MOI.set(model, MOI.VariableName(), x, starting_letter * "x")
        c = MOI.get(model, MOI.ConstraintIndex, "c1")
        MOI.set(model, MOI.ConstraintName(), c, starting_letter * "c1")
        MOI.write_to_file(model, LP_TEST_FILE)
        @test read(LP_TEST_FILE, String) ==
              "minimize\n" *
              "obj: _$(starting_letter)x\n" *
              "subject to\n" *
              "_$(starting_letter)c1: 1 _$(starting_letter)x >= -1\n" *
              "Bounds\n" *
              "_$(starting_letter)x free\n" *
              "End\n"
    end
end

function test_name_sanitization_illegal()
    for illegal_letter in ["[", "]", "*", "^"]
        model = LP.Model()
        MOI.Utilities.loadfromstring!(
            model,
            """
variables: x
minobjective: x
c1: 1.0 * x >= -1.0
""",
        )
        x = MOI.get(model, MOI.VariableIndex, "x")
        MOI.set(model, MOI.VariableName(), x, "x$(illegal_letter)y")
        c = MOI.get(model, MOI.ConstraintIndex, "c1")
        MOI.set(model, MOI.ConstraintName(), c, "c$(illegal_letter)d")
        MOI.write_to_file(model, LP_TEST_FILE)
        @test read(LP_TEST_FILE, String) ==
              "minimize\n" *
              "obj: x_y\n" *
              "subject to\n" *
              "c_d: 1 x_y >= -1\n" *
              "Bounds\n" *
              "x_y free\n" *
              "End\n"
    end
end

function test_name_sanitization_duplicate()
    model = LP.Model()
    MOI.Utilities.loadfromstring!(
        model,
        """
variables: a, b, c, d
minobjective: a + b + c + d
c1: a + b + c + d >= -1.0
""",
    )
    variables = MOI.get(model, MOI.ListOfVariableIndices())
    MOI.set.(model, MOI.VariableName(), variables, ["a[", "a]", "a*", "a^"])
    MOI.write_to_file(model, LP_TEST_FILE)
    @test read(LP_TEST_FILE, String) ==
          "minimize\n" *
          "obj: 1 a_ + 1 a__1 + 1 a__2 + 1 a__3\n" *
          "subject to\n" *
          "c1: 1 a_ + 1 a__1 + 1 a__2 + 1 a__3 >= -1\n" *
          "Bounds\n" *
          "a_ free\n" *
          "a__1 free\n" *
          "a__2 free\n" *
          "a__3 free\n" *
          "End\n"
end

function test_name_sanitization_length()
    model = LP.Model(maximum_length = 3)
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "abcdefg")
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    MOI.write_to_file(model, LP_TEST_FILE)
    @test read(LP_TEST_FILE, String) ==
          "minimize\n" *
          "obj: abc\n" *
          "subject to\n" *
          "Bounds\n" *
          "abc free\n" *
          "End\n"
end

function test_name_sanitization_too_long()
    model = LP.Model(maximum_length = 3)
    MOI.Utilities.loadfromstring!(
        model,
        """
variables: abcd, abce
minobjective: abcd + abce
c1: abcd + abce >= -1.0
""",
    )
    MOI.write_to_file(model, LP_TEST_FILE)
    @test read(LP_TEST_FILE, String) ==
          "minimize\n" *
          "obj: 1 abc + 1 abc_1\n" *
          "subject to\n" *
          "c1: 1 abc + 1 abc_1 >= -1\n" *
          "Bounds\n" *
          "abc free\n" *
          "abc_1 free\n" *
          "End\n"
end

function test_name_sanitization_other()
    model = LP.Model()
    MOI.Utilities.loadfromstring!(
        model,
        """
variables: x
maxobjective: 2.0 * x + -1.0
""",
    )
    MOI.write_to_file(model, LP_TEST_FILE)
    @test read(LP_TEST_FILE, String) ==
          "maximize\n" *
          "obj: -1 + 2 x\n" *
          "subject to\n" *
          "Bounds\n" *
          "x free\n" *
          "End\n"
end

function test_free_variables()
    model = LP.Model()
    MOI.Utilities.loadfromstring!(
        model,
        """
variables: x, y, z
maxobjective: x
x in ZeroOne()
y in Integer()
""",
    )
    MOI.write_to_file(model, LP_TEST_FILE)
    @test read(LP_TEST_FILE, String) ==
          "maximize\n" *
          "obj: x\n" *
          "subject to\n" *
          "Bounds\n" *
          "y free\n" *
          "z free\n" *
          "General\n" *
          "y\n" *
          "Binary\n" *
          "x\n" *
          "End\n"
end

function test_quadratic_objective()
    model = LP.Model()
    @test_throws(
        MOI.UnsupportedAttribute,
        MOI.Utilities.loadfromstring!(
            model,
            """
variables: x
minobjective: 1.0*x*x
""",
        )
    )
end

###
### Read tests
###

function test_read_invalid()
    models = joinpath(@__DIR__, "models")
    for filename in filter(f -> startswith(f, "invalid_"), readdir(models))
        model = LP.Model()
        @test_throws(
            ErrorException,
            MOI.read_from_file(model, joinpath(models, filename)),
        )
    end
    return
end

function test_read_example_lo1()
    model = LP.Model()
    MOI.read_from_file(model, joinpath(@__DIR__, "models", "example_lo1.lp"))
    @test MOI.get(model, MOI.NumberOfVariables()) == 4
    constraints = MOI.get(model, MOI.ListOfConstraintTypesPresent())
    @test (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) in
          constraints
    @test (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) in
          constraints
    @test (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) in
          constraints
    @test (MOI.VariableIndex, MOI.GreaterThan{Float64}) in constraints
    @test (MOI.VariableIndex, MOI.Interval{Float64}) in constraints
    io = IOBuffer()
    write(io, model)
    seekstart(io)
    file = read(io, String)
    @test occursin("maximize", file)
    @test occursin("obj: 3 x1 + 1 x2 + 5 x3 + 1 x4", file)
    @test occursin("c1: 3 x1 + 1 x2 + 2 x3 = 30", file)
    @test occursin("c2: 2 x1 + 1 x2 + 3 x3 + 1 x4 >= 15", file)
    @test occursin("c3: 2 x2 + 3 x4 <= 25", file)
    @test occursin("x1 >= 0", file)
    @test occursin("0 <= x2 <= 10", file)
    @test occursin("x3 >= 0", file)
    @test occursin("x4 >= 0", file)
    return
end

function test_read_model1_tricky()
    model = LP.Model()
    MOI.read_from_file(model, joinpath(@__DIR__, "models", "model1_tricky.lp"))
    @test MOI.get(model, MOI.NumberOfVariables()) == 8
    var_names = MOI.get.(model, MOI.VariableName(), MOI.VariableIndex.(1:8))
    @test Set(var_names) ==
          Set(["Var4", "V5", "V1", "V2", "V3", "V6", "V7", "V8"])
    io = IOBuffer()
    write(io, model)
    seekstart(io)
    file = read(io, String)
    @test occursin("maximize", file)
    @test occursin("obj: -1 Var4 + 1 V5", file)
    @test occursin("CON3: 1 V3 <= 2.5", file)
    @test occursin("CON4: 1 V5 + 1 V6 + 1 V7 <= 1", file)
    @test occursin("CON1: 1 V1 >= 0", file)
    @test occursin("R1: 1 V2 >= 2", file)
    @test occursin("V1 <= 3", file)
    @test occursin("Var4 >= 5.5", file)
    @test occursin("V3 >= -3", file)
    @test occursin("V5 = 1", file)
    @test occursin("0 <= V2 <= 3", file)
    @test occursin("0 <= V7 <= 1", file)
    @test occursin("0 <= V8 <= 1", file)
    @test occursin("V6 free", file)
    @test occursin("\nVar4\n", file)
    @test occursin("\nV5\n", file)
    @test occursin("\nV6\n", file)
    @test occursin("Binary\nV8\n", file)
    @test occursin("sos1: S1:: V1:1.0 V2:2.0 V3:3.0", file)
    @test occursin("sos2: S2:: V1:8.5 V2:10.2 V3:18.3", file)
    return
end

function test_read_model1()
    model = LP.Model()
    MOI.read_from_file(model, joinpath(@__DIR__, "models", "model1.lp"))
    constraints = MOI.get(model, MOI.ListOfConstraintTypesPresent())
    @test (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) in
          constraints
    @test (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) in
          constraints
    @test (MOI.VariableIndex, MOI.GreaterThan{Float64}) in constraints
    @test (MOI.VariableIndex, MOI.Interval{Float64}) in constraints
    @test (MathOptInterface.VariableIndex, MathOptInterface.Integer) in
          constraints
    @test (MathOptInterface.VariableIndex, MathOptInterface.ZeroOne) in
          constraints
    @test (
        MathOptInterface.VectorOfVariables,
        MathOptInterface.SOS1{Float64},
    ) in constraints
    @test (
        MathOptInterface.VectorOfVariables,
        MathOptInterface.SOS2{Float64},
    ) in constraints
    return
end

function test_read_model2()
    model = LP.Model()
    MOI.read_from_file(model, joinpath(@__DIR__, "models", "model2.lp"))
    @test MOI.get(model, MOI.NumberOfVariables()) == 8
    constraints = MOI.get(model, MOI.ListOfConstraintTypesPresent())
    @test (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) in
          constraints
    @test (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) in
          constraints
    @test (MOI.VariableIndex, MOI.GreaterThan{Float64}) in constraints
    @test (MOI.VariableIndex, MOI.Interval{Float64}) in constraints
    @test (MathOptInterface.VariableIndex, MathOptInterface.Integer) in
          constraints
    @test (MathOptInterface.VariableIndex, MathOptInterface.ZeroOne) in
          constraints
    # Adicionar testes dos bounds de V8
    @test MOI.get(model, MOI.VariableName(), MOI.VariableIndex(8)) == "V8"
    @test model.variables.lower[8] == -Inf
    @test model.variables.upper[8] == -3
    obj_type = MOI.get(model, MOI.ObjectiveFunctionType())
    obj_func = MOI.get(model, MOI.ObjectiveFunction{obj_type}())
    @test obj_func.constant == 2.5
    return
end

function test_read_objective_sense()
    model = LP.Model()
    cases = Dict(
        "max" => MOI.MAX_SENSE,
        "maximize" => MOI.MAX_SENSE,
        "maximise" => MOI.MAX_SENSE,
        "maximum" => MOI.MAX_SENSE,
        "min" => MOI.MIN_SENSE,
        "minimize" => MOI.MIN_SENSE,
        "minimise" => MOI.MIN_SENSE,
        "minimum" => MOI.MIN_SENSE,
    )
    for (sense, result) in cases
        LP._set_objective_sense(LP._KW_OBJECTIVE, model, sense)
        @test MOI.get(model, MOI.ObjectiveSense()) == result
    end
    return
end

function test_read_nonempty_model()
    filename = joinpath(@__DIR__, "models", "model2.lp")
    model = LP.Model()
    MOI.read_from_file(model, filename)
    @test_throws(
        ErrorException("Cannot read in file because model is not empty."),
        MOI.read_from_file(model, filename),
    )
    return
end

function test_read_maximum_length_error()
    filename = joinpath(@__DIR__, "models", "model2.lp")
    model = LP.Model(; maximum_length = 1)
    @test_throws(
        ErrorException("Name exceeds maximum length: V4"),
        MOI.read_from_file(model, filename),
    )
    return
end

function test_default_bound()
    io = IOBuffer()
    write(io, "minimize\nobj: x + y")
    seekstart(io)
    model = LP.Model()
    MOI.read!(io, model)
    x = MOI.get(model, MOI.ListOfVariableIndices())
    F, S = MOI.VariableIndex, MOI.GreaterThan{Float64}
    c = [MOI.ConstraintIndex{F,S}(xi.value) for xi in x]
    sets = MOI.get.(model, MOI.ConstraintSet(), c)
    @test all(s -> s == MOI.GreaterThan(0.0), sets)
    @test length(sets) == 2
    return
end

function test_default_bound_double_bound()
    io = IOBuffer()
    write(io, "minimize\nobj: x\nsubject to\nbounds\n x <= -1\n x >= -2")
    seekstart(io)
    model = LP.Model()
    MOI.read!(io, model)
    x = first(MOI.get(model, MOI.ListOfVariableIndices()))
    F = MOI.VariableIndex
    @test MOI.get(
        model,
        MOI.ConstraintSet(),
        MOI.ConstraintIndex{F,MOI.GreaterThan{Float64}}(x.value),
    ) == MOI.GreaterThan(-2.0)
    @test MOI.get(
        model,
        MOI.ConstraintSet(),
        MOI.ConstraintIndex{F,MOI.LessThan{Float64}}(x.value),
    ) == MOI.LessThan(-1.0)
    return
end

function runtests()
    for name in names(@__MODULE__, all = true)
        if startswith("$(name)", "test_")
            @testset "$name" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    sleep(1.0)  # Allow time for unlink to happen.
    rm(LP_TEST_FILE, force = true)
    return
end

end

TestLP.runtests()

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
# Baseline version:
module MatlabPort

# matlab compat functions
_land(a, b) = a != 0 && b != 0
function land(a, b)
    map(_land, a, b)
end

ndgrid(v::AbstractVector) = copy(v)
function ndgrid(v1::AbstractVector, v2::AbstractVector)
    m, n = length(v1), length(v2)
    v1 = reshape(v1, m, 1)
    v2 = reshape(v2, 1, n)
    (repmat(v1, 1, n), repmat(v2, m, 1))
end

function ndgrid_fill(a, v, s, snext)
    for j = 1:length(a)
        a[j] = v[div(rem(j-1, snext), s)+1]
    end
end
function ndgrid{T}(vs::AbstractVector{T}...)
    n = length(vs)
    sz = map(length, vs)
    out = ntuple(i->Array{T}(sz), n)
    s = 1
    for i=1:n
        a = out[i]::Array
        v = vs[i]
        snext = s*size(a,i)
        ndgrid_fill(a, v, s, snext)
        s = snext
    end
    out
end

type ISF
    px::Array{Float64,3}; py::Array{Float64,3}; pz::Array{Float64,3}          # coordinates of grid points
    ix::UnitRange{Int}; iy::UnitRange{Int}; iz::UnitRange{Int}          # 1D index array
    iix::Array{Int, 3}; iiy::Array{Int, 3}; iiz::Array{Int, 3}         # 3D index array
    dx::Float64; dy::Float64; dz::Float64          # edge length
    sizex::Int; sizey::Int; sizez::Int # size of grid
    resx::Int; resy::Int; resz::Int    # number of grid points in each dimension

    hbar::Float64             # reduced Planck constant
    dt::Float64               # time step
    mask::Array{Complex{Float64},3} # Fourier coefficient for solving Schroedinger eq

    function ISF(vol_size::NTuple{3}, vol_res::NTuple{3}, hbar, dt)
        obj = new()
        obj.sizex, obj.sizey, obj.sizez = vol_size
        obj.resx = round(Int, vol_res[1])
        obj.resy = round(Int, vol_res[2])
        obj.resz = round(Int, vol_res[3])
        obj.dx = obj.sizex / obj.resx
        obj.dy = obj.sizey / obj.resy
        obj.dz = obj.sizez / obj.resz
        obj.ix = 1:obj.resx
        obj.iy = 1:obj.resy
        obj.iz = 1:obj.resz
        obj.iix, obj.iiy, obj.iiz = ndgrid(obj.ix,obj.iy,obj.iz)
        obj.px = (obj.iix-1) * obj.dx
        obj.py = (obj.iiy-1) * obj.dy
        obj.pz = (obj.iiz-1) * obj.dz
        obj.hbar = hbar
        obj.dt = dt
        BuildSchroedinger(obj)
        obj
    end

end

function BuildSchroedinger(obj::ISF)
    nx=obj.resx; ny=obj.resy; nz=obj.resz;
    fac = -4*pi^2*obj.hbar;
    kx = (obj.iix .- 1 .- nx ./ 2) ./ (obj.sizex);
    ky = (obj.iiy .- 1 .- ny ./ 2) ./ (obj.sizey);
    kz = (obj.iiz .- 1 .- nz ./ 2) ./ (obj.sizez);
    lambda = fac .* (kx .^ 2 .+ ky .^ 2 .+ kz .^ 2);
    obj.mask = exp.(1.0im*lambda*obj.dt/2.);
end
"""
solves Schroedinger equation for dt time.
"""
function SchroedingerFlow(obj, psi1, psi2)
    psi1 = fftshift(fft(psi1)); psi2 = fftshift(fft(psi2));
    psi1 = psi1.*obj.mask;
    psi2 = psi2.*obj.mask;
    psi1 = ifft(fftshift(psi1)); psi2 = ifft(fftshift(psi2));
    psi1, psi2
end

"""
 For a 1-form v compute the corresponding vector field `v^sharp` as
 a staggered vector field living on edges
"""
function StaggeredSharp(obj,vx,vy,vz)
    ux = vx/obj.dx
    uy = vy/obj.dy
    uz = vz/obj.dz
    ux,uy,uz
end

function VelocityOneForm(obj, psi1, psi2, hbar=1.0)
    ixp = mod.(obj.ix, obj.resx) + 1;
    iyp = mod.(obj.iy, obj.resy) + 1;
    izp = mod.(obj.iz, obj.resz) + 1;
    vx = angle.(
        conj.(psi1).*view(psi1, ixp,:,:) .+
        conj.(psi2).*view(psi2, ixp,:,:)
    );
    vy = angle.(
        conj.(psi1).*view(psi1,:,iyp,:) .+
        conj.(psi2).*view(psi2,:,iyp,:)
    )
    vz = angle.(
        conj.(psi1).*view(psi1,:,:,izp) .+
        conj.(psi2).*view(psi2,:,:,izp)
    )
    vx = vx*hbar;
    vy = vy*hbar;
    vz = vz*hbar;
    vx,vy,vz
end

function PoissonSolve(obj, f)
    f = fft(f)
    sx = sin.(pi .* (obj.iix .- 1) ./ obj.resx) ./ obj.dx
    sy = sin.(pi .* (obj.iiy .- 1) ./ obj.resy) ./ obj.dy
    sz = sin.(pi .* (obj.iiz .- 1) ./ obj.resz) ./ obj.dz
    denom = sx.^2 + sy.^2 + sz.^2
    fac = -0.25./denom
    fac[1,1,1] = 0.0
    f = f .* fac
    ifft(f)
end

function GaugeTransform(psi1, psi2, q)
    eiq = exp.(1.0im .* q);
    psi1 = psi1 .* eiq
    psi2 = psi2 .* eiq
    psi1, psi2
end

function PressureProject(obj, psi1, psi2)
    vx,vy,vz = VelocityOneForm(obj, psi1, psi2)
    div = Div(obj, vx,vy,vz)
    q = -PoissonSolve(obj, div)
    x = GaugeTransform(psi1, psi2, q)
    x
end

function Normalize(psi1, psi2)
    psi_norm = sqrt.(abs.(psi1) .^ 2 .+ abs.(psi2) .^ 2)
    psi1 = psi1 ./ psi_norm
    psi2 = psi2 ./ psi_norm
    psi1, psi2
end

function Div(obj, vx, vy, vz)
    ixm = mod.(obj.ix-2, obj.resx) .+ 1
    iym = mod.(obj.iy-2, obj.resy) .+ 1
    izm = mod.(obj.iz-2, obj.resz) .+ 1
    f =     (vx .- view(vx, ixm,:,:)) ./ (obj.dx .^ 2)
    f = f .+ (vy .- view(vy, :,iym,:)) ./ (obj.dy .^ 2)
    f = f .+ (vz .- view(vz, :,:,izm)) ./ (obj.dz .^ 2)
    f
end

"""
 Particles - class of particle which can be advected by staggered velocity
 field on a TorusDEC grid, using RK4 method.
 Velocities are trilinearly interpolated.
"""
type Particles
    x::Vector{Float32}
    y::Vector{Float32}
    z::Vector{Float32} # array of positions
    length::Int
    active::Vector{Int}
end
Base.length(p::Particles) = p.length

function Base.append!(p::Particles, x, y, z)
    @assert length(x) + p.length < length(p.x)
    range = (p.length+1):(p.length+length(x))
    p.x[range] = x
    p.y[range] = y
    p.z[range] = z
    p.length   = (p.length+length(x))
    append!(p.active, range)
    nothing
end
"""
 advect particle positions using RK4 in a grid torus with
 staggered velocity vx,vy,vz, for dt period of time
"""
function StaggeredAdvect(p, torus, vx, vy, vz, dt)
    px, py, pz = p.x[p.active], p.y[p.active], p.z[p.active]
    k1x,k1y,k1z = StaggeredVelocity(
        px,py,pz,
        torus,vx,vy,vz
    )
    k2x,k2y,k2z = StaggeredVelocity(
        px .+ k1x .* dt ./ 2,py .+ k1y .* dt/2,pz .+ k1z .* dt ./ 2,
        torus,vx,vy,vz
    )
    k3x,k3y,k3z = StaggeredVelocity(
        px .+ k2x .* dt ./ 2, py .+ k2y .* dt ./ 2, pz .+ k2z .* dt ./2,
        torus, vx, vy, vz
    )
    k4x,k4y,k4z = StaggeredVelocity(
        px .+ k3x .* dt,py .+ k3y .* dt,pz .+ k3z .* dt,
        torus,vx,vy,vz
    )
    p.x[p.active] .= px .+ dt ./ 6 .* (k1x .+ 2 .* k2x .+ 2 .* k3x .+ k4x)
    p.y[p.active] .= py .+ dt ./ 6 .* (k1y .+ 2 .* k2y .+ 2 .* k3y .+ k4y)
    p.z[p.active] .= pz .+ dt ./ 6 .* (k1z .+ 2 .* k2z .+ 2 .* k3z .+ k4z)
end
"""
for removing particles
"""
function Keep(particle, ind)
    particle.x = particle.x[ind]
    particle.y = particle.y[ind]
    particle.z = particle.z[ind]
end

"""
 evaluates velocity at (px,py,pz) in the grid torus with staggered
 velocity vector field vx,vy,vz
"""
function StaggeredVelocity(pxf,pyf,pzf,torus,vx,vy,vz)
    px = mod.(pxf, torus.sizex)
    py = mod.(pyf, torus.sizey)
    pz = mod.(pzf, torus.sizez)

    ix = floor.(Int, px/torus.dx) .+ 1
    iy = floor.(Int, py/torus.dy) .+ 1
    iz = floor.(Int, pz/torus.dz) .+ 1
    ixp = mod.(ix,torus.resx) .+ 1
    iyp = mod.(iy,torus.resy) .+ 1
    izp = mod.(iz,torus.resz) .+ 1

    res = Base.RefValue((torus.resx,torus.resy,torus.resz))

    ind0 = sub2ind.(res, ix, iy, iz)
    indxp = sub2ind.(res, ixp, iy, iz)
    indyp = sub2ind.(res, ix, iyp, iz)
    indzp = sub2ind.(res, ix, iy, izp)
    indxpyp = sub2ind.(res, ixp, iyp, iz)
    indypzp = sub2ind.(res, ix, iyp, izp)
    indxpzp = sub2ind.(res, ixp, iy, izp)

    wx = px .- (ix .- 1) .* torus.dx
    wy = py .- (iy .- 1) .* torus.dy
    wz = pz .- (iz .- 1) .* torus.dz

    ux = ((1-wz).*((1-wy).*vx[ind0 ]+wy.*vx[indyp  ]) .+
            wz .*((1-wy).*vx[indzp]+wy.*vx[indypzp]))

    uy = ((1-wz).*((1-wx).*vy[ind0 ]+wx.*vy[indxp  ]) .+
            wz .*((1-wx).*vy[indzp]+wx.*vy[indxpzp]))

    uz = ((1-wy).*((1-wx).*vz[ind0 ]+wx.*vz[indxp  ]) .+
            wy .*((1-wx).*vz[indyp]+wx.*vz[indxpyp]))

    ux,uy,uz
end

end

ext: jl


text:
import Base.keys
import Base.values

## These are the containers that can be looped over
## The prefix SDM is for SortedDict and SortedMultiDict
## The prefix SS is for SortedSet.  The prefix SA
## is for all sorted containers.
## The following two definitions now appear in tokens2.jl

# const SDMContainer = Union{SortedDict, SortedMultiDict}
# const SAContainer = Union{SDMContainer, SortedSet}

@inline extractcontainer(s::SAContainer) = s

## This holds an object describing an exclude-last
## iteration.
abstract type AbstractExcludeLast{ContainerType <: SAContainer} end

struct SDMExcludeLast{ContainerType <: SDMContainer} <:
                              AbstractExcludeLast{ContainerType}
    m::ContainerType
    first::Int
    pastlast::Int
end

struct SSExcludeLast{ContainerType <: SortedSet} <:
                              AbstractExcludeLast{ContainerType}
    m::ContainerType
    first::Int
    pastlast::Int
end

@inline extractcontainer(s::AbstractExcludeLast) = s.m
eltype(s::AbstractExcludeLast) = eltype(s.m)

## This holds an object describing an include-last
## iteration.

abstract type AbstractIncludeLast{ContainerType <: SAContainer} end

struct SDMIncludeLast{ContainerType <: SDMContainer} <:
                               AbstractIncludeLast{ContainerType}
    m::ContainerType
    first::Int
    last::Int
end

struct SSIncludeLast{ContainerType <: SortedSet} <:
                               AbstractIncludeLast{ContainerType}
    m::ContainerType
    first::Int
    last::Int
end

@inline extractcontainer(s::AbstractIncludeLast) = s.m
eltype(s::AbstractIncludeLast) = eltype(s.m)

## The basic iterations are either over the whole sorted container, an
## exclude-last object or include-last object.

const SDMIterableTypesBase = Union{SDMContainer,
                                   SDMExcludeLast,
                                   SDMIncludeLast}

const SSIterableTypesBase = Union{SortedSet,
                                  SSExcludeLast,
                                  SSIncludeLast}

const SAIterableTypesBase = Union{SAContainer,
                                  AbstractExcludeLast,
                                  AbstractIncludeLast}

## The compound iterations are obtained by applying keys(..) or values(..)
## to the basic iterations of the SDM.. type.
## Furthermore, semitokens(..) can be applied
## to either a basic iteration or a keys/values iteration.

struct SDMKeyIteration{T <: SDMIterableTypesBase}
    base::T
end

eltype(s::SDMKeyIteration) = keytype(extractcontainer(s.base))
length(s::SDMKeyIteration) = length(extractcontainer(s.base))

struct SDMValIteration{T <: SDMIterableTypesBase}
    base::T
end

eltype(s::SDMValIteration) = valtype(extractcontainer(s.base))
length(s::SDMValIteration) = length(extractcontainer(s.base))

struct SDMSemiTokenIteration{T <: SDMIterableTypesBase}
    base::T
end

eltype(s::SDMSemiTokenIteration) = Tuple{IntSemiToken,
                                         keytype(extractcontainer(s.base)),
                                         valtype(extractcontainer(s.base))}

struct SSSemiTokenIteration{T <: SSIterableTypesBase}
    base::T
end

eltype(s::SSSemiTokenIteration) = Tuple{IntSemiToken,
                                        eltype(extractcontainer(s.base))}

struct SDMSemiTokenKeyIteration{T <: SDMIterableTypesBase}
    base::T
end

eltype(s::SDMSemiTokenKeyIteration) = Tuple{IntSemiToken,
                                            keytype(extractcontainer(s.base))}

struct SAOnlySemiTokensIteration{T <: SAIterableTypesBase}
    base::T
end

eltype(::SAOnlySemiTokensIteration) = IntSemiToken

struct SDMSemiTokenValIteration{T <: SDMIterableTypesBase}
    base::T
end

eltype(s::SDMSemiTokenValIteration) = Tuple{IntSemiToken,
                                            valtype(extractcontainer(s.base))}

const SACompoundIterable = Union{SDMKeyIteration,
                                 SDMValIteration,
                                 SDMSemiTokenIteration,
                                 SSSemiTokenIteration,
                                 SDMSemiTokenKeyIteration,
                                 SDMSemiTokenValIteration,
                                 SAOnlySemiTokensIteration}

@inline extractcontainer(s::SACompoundIterable) = extractcontainer(s.base)

const SAIterable = Union{SAIterableTypesBase, SACompoundIterable}

## All the loops maintain a state which is an object of the
## following type.

struct SAIterationState
    next::Int
    final::Int
end

@inline start(m::SAContainer) = SAIterationState(nextloc0(m.bt,1), 2)
@inline start(e::SACompoundIterable) = start(e.base)

function start(e::AbstractExcludeLast)
    (!(e.first in e.m.bt.useddatacells) || e.first == 1 ||
        !(e.pastlast in e.m.bt.useddatacells)) &&
        throw(BoundsError())
    if compareInd(e.m.bt, e.first, e.pastlast) < 0
        return SAIterationState(e.first, e.pastlast)
    else
        return SAIterationState(2, 2)
    end
end

function start(e::AbstractIncludeLast)
    (!(e.first in e.m.bt.useddatacells) || e.first == 1 ||
        !(e.last in e.m.bt.useddatacells) || e.last == 2) &&
        throw(BoundsError())
    if compareInd(e.m.bt, e.first, e.last) <= 0
        return SAIterationState(e.first, nextloc0(e.m.bt, e.last))
    else
        return SAIterationState(2, 2)
    end
end

## All the loops share the same `iterate` method
@inline function iterate(s::SAIterable, state::SAIterationState=start(s))
    state.next == state.final && return nothing
    return next(s, state)
end

@inline exclusive(m::T, ii::(Tuple{IntSemiToken,IntSemiToken})) where {T <: SDMContainer} =
    SDMExcludeLast(m, ii[1].address, ii[2].address)

@inline exclusive(m::T, ii::(Tuple{IntSemiToken,IntSemiToken})) where {T <: SortedSet} =
    SSExcludeLast(m, ii[1].address, ii[2].address)

@inline exclusive(m::T, i1::IntSemiToken, i2::IntSemiToken) where {T <: SAContainer} =
    exclusive(m, (i1,i2))

@inline inclusive(m::T, ii::(Tuple{IntSemiToken,IntSemiToken})) where {T <: SDMContainer} =
    SDMIncludeLast(m, ii[1].address, ii[2].address)

@inline inclusive(m::T, ii::(Tuple{IntSemiToken,IntSemiToken})) where {T <: SortedSet} =
    SSIncludeLast(m, ii[1].address, ii[2].address)

@inline inclusive(m::T, i1::IntSemiToken, i2::IntSemiToken) where {T <: SAContainer} =
    inclusive(m, (i1,i2))

# Next definition needed to break ambiguity with keys(AbstractDict) from Dict.jl

@inline keys(ba::SortedDict{K,D,Ord}) where {K, D, Ord <: Ordering} = SDMKeyIteration(ba)
@inline keys(ba::T) where {T <: SDMIterableTypesBase} = SDMKeyIteration(ba)

in(k, keyit::SDMKeyIteration{SortedDict{K,D,Ord}}) where {K,D,Ord <: Ordering} =
    haskey(extractcontainer(keyit.base), k)

in(k, keyit::SDMKeyIteration{SortedMultiDict{K,D,Ord}}) where {K,D,Ord <: Ordering} =
    haskey(extractcontainer(keyit.base), k)

# Next definition needed to break ambiguity with values(AbstractDict) from Dict.jl
@inline values(ba::SortedDict{K,D,Ord}) where {K, D, Ord <: Ordering} = SDMValIteration(ba)
@inline values(ba::T) where {T <: SDMIterableTypesBase} = SDMValIteration(ba)
@inline semitokens(ba::T) where {T <: SDMIterableTypesBase} = SDMSemiTokenIteration(ba)
@inline semitokens(ba::T) where {T <: SSIterableTypesBase} = SSSemiTokenIteration(ba)
@inline semitokens(ki::SDMKeyIteration{T}) where {T <: SDMIterableTypesBase} =
                   SDMSemiTokenKeyIteration(ki.base)
@inline semitokens(vi::SDMValIteration{T}) where {T <: SDMIterableTypesBase} =
                   SDMSemiTokenValIteration(vi.base)
@inline onlysemitokens(ba::T) where {T <: SAIterableTypesBase} = SAOnlySemiTokensIteration(ba)

## The 'next' function returns different objects depending on whether
## it is a basic iteration, a key iteration, a values iterations,
## a semitokens/basic iteration, a semitokens/key iteration, or semitokens/values
## iteration.

@inline function next(u::SAOnlySemiTokensIteration, state::SAIterationState)
    sn = state.next
    (sn < 3 || !(sn in extractcontainer(u).bt.useddatacells)) && throw(BoundsError())
    IntSemiToken(sn),
    SAIterationState(nextloc0(extractcontainer(u).bt, sn), state.final)
end

@inline function nexthelper(u, state::SAIterationState)
    sn = state.next
    (sn < 3 || !(sn in extractcontainer(u).bt.useddatacells)) && throw(BoundsError())
    extractcontainer(u).bt.data[sn], sn,
    SAIterationState(nextloc0(extractcontainer(u).bt, sn), state.final)
end

@inline function next(u::SDMIterableTypesBase, state::SAIterationState)
    dt, t, ni = nexthelper(u, state)
    (dt.k => dt.d), ni
end

@inline function next(u::SSIterableTypesBase, state::SAIterationState)
    dt, t, ni = nexthelper(u, state)
    dt.k, ni
end

@inline function next(u::SDMKeyIteration, state::SAIterationState)
    dt, t, ni = nexthelper(u, state)
    dt.k, ni
end

@inline function next(u::SDMValIteration, state::SAIterationState)
    dt, t, ni = nexthelper(u, state)
    dt.d, ni
end

@inline function next(u::SDMSemiTokenIteration, state::SAIterationState)
    dt, t, ni = nexthelper(u, state)
    (IntSemiToken(t), dt.k, dt.d), ni
end

@inline function next(u::SSSemiTokenIteration, state::SAIterationState)
    dt, t, ni = nexthelper(u, state)
    (IntSemiToken(t), dt.k), ni
end

@inline function next(u::SDMSemiTokenKeyIteration, state::SAIterationState)
    dt, t, ni = nexthelper(u, state)
    (IntSemiToken(t), dt.k), ni
end

@inline function next(u::SDMSemiTokenValIteration, state::SAIterationState)
    dt, t, ni = nexthelper(u, state)
    (IntSemiToken(t), dt.d), ni
end

eachindex(sd::SortedDict) = keys(sd)
eachindex(sdm::SortedMultiDict) = onlysemitokens(sdm)
eachindex(ss::SortedSet) = onlysemitokens(ss)
eachindex(sd::SDMExcludeLast{SortedDict{K,D,Ord}}) where {K,D,Ord <: Ordering} = keys(sd)
eachindex(smd::SDMExcludeLast{SortedMultiDict{K,D,Ord}}) where {K,D,Ord <: Ordering} =
     onlysemitokens(smd)
eachindex(ss::SSExcludeLast) = onlysemitokens(ss)
eachindex(sd::SDMIncludeLast{SortedDict{K,D,Ord}}) where {K,D,Ord <: Ordering} = keys(sd)
eachindex(smd::SDMIncludeLast{SortedMultiDict{K,D,Ord}}) where {K,D,Ord <: Ordering} =
     onlysemitokens(smd)
eachindex(ss::SSIncludeLast) = onlysemitokens(ss)

empty!(m::SAContainer) =  empty!(m.bt)
@inline length(m::SAContainer) = length(m.bt.data) - length(m.bt.freedatainds) - 2
@inline isempty(m::SAContainer) = length(m) == 0


text:
@with_kw struct PoissonParameter
    rate::SNNFloat = 1Hz
end

@with_kw mutable struct Poisson
    param::PoissonParameter = PoissonParameter()
    N::SNNInt = 100
    randcache::Vector{SNNFloat} = rand(N)
    fire::Vector{Bool} = zeros(Bool, N)
    records::Dict = Dict()
end

function integrate!(p::Poisson, param::PoissonParameter, dt::SNNFloat)
    prob = rate * dt
    rand!(randcache)
    @inbounds for i = 1:N
        fire[i] = randcache[i] < prob
    end
end


text:
mutable struct BinaryNode{T}
    data::T
    left::BinaryNode{T}
    right::BinaryNode{T}

    BinaryNode{T}(data) where T = new{T}(data)
    BinaryNode{T}(data, left, right) where T = new{T}(data, left, right)
end

BinaryNode(data::T) where T = BinaryNode{typeof(data)}(data)

BinaryNode(data::T, left::BinaryNode{T}, right::BinaryNode{T}) where T =
BinaryNode{T}(data, left, right)

function inordertraversal(root::BinaryNode)
    acc = Array{BinaryNode,1}()
    function f(acc, root)
        push!(acc, root)
        if isdefined(root, :left)
            f(acc, root.left)
        end

        if isdefined(root, :right)
            f(acc, root.right)
        end
        return nothing
    end
    f(acc, root)
    return acc
end


text:
#=

Calculate recursions from paper http://eprints.lancs.ac.uk/745/1/online_chpt4.pdf

We are going to calculate an approximation to the posterior - approximate because particle filtering
is involved to increase the speed of the algorithm but reduce the computtional cost, this is an online
algorithm but it is suited at looking at retrospective

=#

function recursions( data::Array{Real}, LOS::DiscreteUnivariateDistribution,  segment_ML::Function )
	n=length(data)
	# weights and last cpt locations
	weights = Array{Array{Float64} }(undef, n)
	cpt_locations =  Array{Array{Int64}}(undef, n)
	weights[1] =
	cpt_locations[1] = [0,1]
	for t in 2:n

		for j in cpt_locations[t-1]
			logccdf(LOS , t-j) - logccdf(LOS, t-j-1)
		end

		if t > 20

		else
		# resample

		end

	end

end


text:
# Model: Anand, McFarland, Kihm, Wong (Optimization of Aerosol Penetration through Transport Lines, 1992)
# Based on: straightTubePenetrationAnand1992.m

if "./src" in LOAD_PATH
    # print("in src")
else
    dir = "./src"
    push!(LOAD_PATH,dir)
    # print("not in src")
end

module bronchiole

import physicalConstants: k, g, p, MM

struct inclinedTube
    # CarrierFluidProfile
    T :: Float64          # Carrier fluid temperature (T) SI units: K
    ro :: Float64         # Carrier fluid density (ro) SI units: kg m-3
    Q :: Float64          # Carrier fluid flow rate (Q) SI units: m3 s-1
    mu :: Float64         # Carrier fluid dynamic viscosity (mu) SI units: N s m-2

    # Droplet
    ro_d :: Float64         # Droplet density (ro_d) SI units: kg m-3
    # mu_d :: Float64       # Droplet dynamic viscosity (mu_d) SI units: N s m-2
    d_d :: Float64          # Droplet diameter (dd) SI units: m
    md :: Float64           # Droplet mass (md) SI units: kg

    # Tube
    d :: Float64            # Tube inner diameter (d) SI units: m
    phi :: Float64          # Inclination angle in radians (phi) SI units: rad
    L :: Float64            # Distance from the tube inlet (L) SI units: m

    # Velocities
    Vb :: Float64           # Brownian diffusion velocity (Vb) SI units: m s-1
    U :: Float64            # Fluid flow mean velocity (U) SI units: m s-1
    Vf :: Float64           # Friction velocity (Vf) SI units: m s-1
    Vs :: Float64           # Dimensionless depositional velocity (Vs)
    Vr :: Float64           # Dimensionless depositional velocity (Vr)
    Vt :: Float64           # Turbulent diffusional velocity (Vt)
    Vg :: Float64           # Gravitational settling velocity (Vg)
    Vd :: Float64           # Sum of Vt and Vb is defined as Vd (Vd)
    Ve :: Float64           # Ve is the effective depositional velocity, given by the vector sum of velocities due to turbulent diffusion (Vt)

    # FluidProperties
    Re :: Float64           # Reynolds number (Re)
    f :: Float64            # Fanning friction factor, Blasius equation for turbulent pipe flow (f) dimensionless
    lambda :: Float64       # Model for the mean free path of the carrier fluid (lambda)
    C :: Float64            # Cunningham slip correction (C)
    tau :: Float64          # Droplet relaxation time (tau)

    # Misc.
    extrapl :: Integer      # Flag to test crossed boundaries
    Rplus :: Float64        # Parameter Rplus (Rplus)
    ParamS :: Float64       # Parameter S (ParamS)
    Splus :: Float64        # Parameter Splus (Splus)
    TETAc :: Float64        # Critical Angle (TETAc)

    # Penetration
    P :: Float64            # P = penetration = C/C0 = ratio between concentration at distance x and initial concentration

    function inclinedTube(T, ro, Q, mu, ro_d, d_d, d, phi, L)
        extrapl = 0
        K = 5e-5;

        md = ro_d*(4/3)*(pi*(d_d/2)^3)          # Droplet mass (md) SI units: kg
        Vb = sqrt((k*T)/(2*pi*md))              # Brownian diffusion velocity (Vb) SI units: m s-1
        U=4*Q/(pi*d^2)                          # Fluid flow mean velocity (U) SI units: m s-1
        Re=(U*d*ro)/mu                          # Reynolds number (Re)
        f=0.316/(4*Re^0.25)                     # Fanning friction factor, Blasius equation for turbulent pipe flow (f) dimensionless

        if (Re .< 3000 || Re .> 100000)
            extrapl=extrapl+1
        end

        Vf = U*sqrt(f/2)                        # Friction velocity (Vf) SI units: m s-1
        lambda = (mu/p)*sqrt((pi*k*T)/(2*MM))
        C = 1 + (lambda/d_d)*(2.34+1.05*exp(-.39*d_d/lambda))    # Cunningham slip correction (C)
        tau = (C*ro_d*d_d^2)/(18*mu)            # Droplet relaxation time (tau)
        Rplus = (d_d*Vf*ro)/(2*mu)              # Parameter Rplus (Rplus)
        ParamS = 0.9*tau*Vf + d_d/2             # Parameter S (ParamS)
        Splus = (ParamS*Vf*ro)/(mu)             # Parameter Splus (Splus)

        if (Splus >= 0 && Splus <= 10)
            Vs = 0.05*Splus
        else
            Vs = 0.5 + 0.0125*(Splus - 10)
            extrapl=extrapl+1
        end

        if (Rplus >= 0 && Rplus <= 10)
            Vr = 0.05*Rplus
        else
            Vr = 0.5 + 0.0125*(Rplus - 10)
            extrapl=extrapl+1
        end

        Vt = (Vf*(Vs+Vr))/(4)                  # Turbulent diffusional velocity (Vt)
        Vg = tau*g*cos(phi)                    # Gravitational settling velocity (Vg)
        Vd = Vt + Vb*K                         # Sum of Vt and Vb is defined as Vd (Vd)

        if (Vd>Vg)                           # If Vd>Vg, there is no need to use TETAc, since the integral that originated the Ve equation can be evaluated from 0 to 2pi. Also, arcsine(x) for x>1 is not defined in the real domain, therefore, TETAc would be imaginary.
            Ve=Vd
            print("\n Vd > Vg")
        else
            TETAc = asin(Vd/Vg)                    # Critical angle (TETAc)
            Ve=Vd*TETAc/pi + Vd/2 + Vg*cos(TETAc)/pi        # Ve is the effective depositional velocity, given by the vector sum of velocities due to turbulent diffusion (Vt),
            print("\n Vd < Vg")
        end

        # print(Re)

        print(Ve)

        P = exp(-(pi*d*Ve*L)/(Q))             # P = penetration = C/C0 = ratio between concentration at distance x and initial concentration

        # results = [md, Vb, U, Re, f, extrapl, Vf, lambda, tau, Rplus, ParamS, Splus, Vs, Vr, Vt, Vg, Vd, TETAc, Ve, P]
        # labels = ["md", "Vb", "U", "Re", "f", "extrapl", "Vf", "lambda", "tau", "Rplus", "ParamS", "Splus", "Vs", "Vr", "Vt", "Vg", "Vd", "TETAc", "Ve", "P"]

        # for i in 1:20
        #     print(labels[i], ": ", results[i], '\n')
        # end

        # results = [md, Vb, U, Re, f, extrapl, Vf, lambda, tau, Rplus, ParamS, Splus, Vs, Vr, Vt, Vg, Vd, Ve, P]
        # labels = ["md", "Vb", "U", "Re", "f", "extrapl", "Vf", "lambda", "tau", "Rplus", "ParamS", "Splus", "Vs", "Vr", "Vt", "Vg", "Vd", "Ve", "P"]

        # for i in 1:19
        #     print(labels[i], ": ", results[i], '\n')
        # end

    end

end

dataset1 = [5.4444E-07, 5.8472E-07, 6.2656E-07, 6.6727E-07, 7.2656E-07, 7.6727E-07, 8.3492E-07, 9.0000E-07, 9.6627E-07, 1.0557E-06, 1.1621E-06, 1.2616E-06, 1.3686E-06, 1.4746E-06, 1.5741E-06, 1.6736E-06, 1.7828E-06, 1.8863E-06, 1.9858E-06, 2.1567E-06, 2.3432E-06, 2.5136E-06, 2.6939E-06, 2.8643E-06, 3.0672E-06, 3.3113E-06, 3.5555E-06, 3.8230E-06]
# d_d = 5.4444E-07

# straightTube(T, ro, Q, mu, ro_d, d_d, d, phi, L)
T = 309.5 # Temperature of Air
ro = 1.12027 # Density of Air (http://www.anaesthesia.med.usyd.edu.au/resources/lectures/humidity_clt/humidity.html#:~:text=While%20nose%20breathing%20at%20rest,%25%20to%2070%25%20relative%20humidity.)
# Q = 0.001 # Based on 60L/min flowrate of ventilators
# Q = 1e-4 # Based on 6L/min of average adult male (Guyton and Hall Textbook of Medical Physiology 12th Edition)
Q = 0.00050038/2 # Chapter 3
mu = 1.895e-5 # Dynamic viscosity of air (https://www.engineersedge.com/physics/viscosity_of_air_dynamic_and_kinematic_14483.htm#:~:text=The%20viscosity%20of%20air%20depends,10%2D5%20Pa%C2%B7s%20.&text=At%2025%20%C2%B0C%2C%20the,the%20kinematic%20viscosity%2015.7%20cSt.)
ro_d = 1000 # Droplet density of water
d = 0.018 # Diameter of first inclined tube
# phi = (pi/2) # Inclination angle in radians (phi) SI units: rad
IBA = 1.040216 # Interbronchial angle, average for all ages, 59.6 degrees SI Units: radians Haskin Goodman et. al. "Normal Tracheal Bifurcation Angle"
phi = pi-(pi/2)-(IBA/2) # Inclination angle in radians (phi) SI units: rad
L = 0.118 # Avg length of trachea (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5900092/#:~:text=On%20average%2C%20the%20length%20of,to%20be%20shorter%20in%20females.&text=Structure%20of%20the%20trachea.)

post_inclinedTube0 = []

for d_d = dataset1
    penetration = inclinedTube(T, ro, Q, mu, ro_d, d_d, d, phi, L)
    push!(post_inclinedTube0, penetration)
end

print(post_inclinedTube0)

# print(inclinedTube(T, ro, Q, mu, ro_d, d_d, d, phi, L))

# Dummy values
# straightTube(309.5, 1, 2, 3, 4, 5, 6, 7, 8)

end


text:
n = 100
d = 10
x = rand(Float64, (n, d))
obs_names = ["obs_$i" for i in 1:n]
ad = AnnData(X=x, obs_names=obs_names)

@testset "create anndata" begin
    @test size(ad) == (n, d)
end

function test_row_slice(ad, i1, n, d, x)
    ad1, ad2, ad3, ad4 = ad[i1, :], ad[collect(i1), :], ad[ad.obs_names[i1], :], ad[1:n .∈ (i1,), :]

    @test size(ad1) == size(ad2) == size(ad3) == size(ad4) == (length(i1), d)
    @test ad1.X == ad2.X == ad3.X == ad4.X == x[i1, :]
    @test ad1.obs_names == ad2.obs_names == ad3.obs_names == ad4.obs_names == ad.obs_names[i1]
end

function test_ad_slicing(ad, n, d, x)
    @testset "row slice $(1:i)" for i in 2:n
        test_row_slice(ad, 1:i, n, d, x)
    end
    @testset "row slice $(i:(i + 1))" for i in 1:(n - 1)
        test_row_slice(ad, i:(i + 1), n, d, x)
    end
end

@testset "slicing anndata" begin
    test_ad_slicing(ad, n, d, x)
end

@testset "views" begin
    i, j = (26:75, [string(i) for i in 3:7])
    adview = @view ad[i, j]
    subad = ad[i, j]
    @test adview.X == x[26:75, 3:7]
    @test adview.obs_names == subad.obs_names == ad.obs_names[i]
    @test adview.var_names == subad.var_names == ad.var_names[3:7]
    @test copy(adview).X == subad.X
    test_ad_slicing(subad, 50, 5, x[i, 3:7])
end


text:
@testset "1021.remove-outermost-parentheses.jl" begin
    @test remove_outer_parentheses("(()())(())(()(()))") == "()()()()(())"
    @test remove_outer_parentheses("(()())(())") == "()()()"
    @test remove_outer_parentheses("()()") == ""
end


text:
using TomoForward
using Images
using Plots
using XfromProjections
using StaticArrays
using PyCall
using Logging
using Printf

include("../../phantoms/simple_phantoms.jl")
include("../tv_primaldual_flow.jl")
include("../optical_flow.jl")
include("../../phantoms/sperm_phantom.jl")

H,W = 609, 609
function radon_operator(height, width, detcount, angles)
    proj_geom = ProjGeom(0.5, detcount, angles)
    A = fp_op_parallel2d_line(proj_geom, height, width, -38.0,38.0, -38.0,38.0)
    return A
end

angles = [0.0,π/3,2*π/3]
#Create phantom at smaller gridsize than reconstruction
grid = collect(detmin:0.1:detmax)
r(s) = 1.0
frames, tracks = get_sperm_phantom(2,r,grid)
detcount = Int(floor(H*1.4))
#match size of input image (generating data)
As = map(t -> radon_operator(size(frames[:,:,1],1),size(frames[:,:,1],2),detcount, angles),1:size(frames)[3])
bsex = zeros(size(As[1])[1],size(frames)[3])

map(t -> bsex[:,t] = As[t]*vec(frames[:,:,t]), 1:size(frames)[3])
@info "adding gaussion noise at level 0.01"
rho = 0.01
e = randn(size(bsex));
e = rho*norm(bsex)*e/norm(e);
bs = bsex + e;

u0s = zeros(H,W,size(frames)[3])

As  = map(t -> radon_operator(H,W, detcount, angles),1:size(u0s)[3])
cwd = @__DIR__
path = normpath(joinpath(@__DIR__, "result"))
!isdir(path) && mkdir(path)
cd(path)

c=10.0
w_tvs = [0.0001,0.001,0.01,0.05,0.1,0.2,0.3,0.5]
w_flows = [0.0001,0.001,0.01,0.05,0.1,0.2,0.3,0.5]
p = Progress(size(w_tvs,1)*size(w_flows,1),60, "Testing parameters")

for w_tv = w_tvs
    niter=450
    us_tv = zeros(H,W,size(frames)[3])
    for t = 1:size(frames)[3]
        A = As[t]
        p = bs[:,t]
        u0 = u0s[:,:,t]
        @suppress begin
            us_tv[:,:,t] .= recon2d_tv_primaldual!(us_tv[:,:,t], A, p, niter, w_tv, c)
        end
    end
    for w_flow = w_flows
        niter=50
        @suppress begin
            us_flow = recon2d_tv_primaldual_flow(As, bs, deepcopy(us_tv), niter, 20, w_tv, w_flow, c)
            filenamepng = @sprintf "flow_2ms_%f_%f_.png" w_tv w_flow;
            #PLOT
            plot(Gray.(us_flow[:,:,1]), aspect_ratio=:equal, framestyle=:none, legend=false, yflip=false)
            savefig(filenamepng)
            next!(p)
        end
    end
end

angles = rand(1:0.001:π,3)
#Create phantom at smaller gridsize than reconstruction
frames = get_sperm_phantom(2,r, grid)
detcount = Int(floor(H*1.4))
#match size of input image (generating data)
As = map(t -> radon_operator(size(frames[:,:,1],1),size(frames[:,:,1],2),detcount, angles[t]),1:size(frames)[3])
bsex = zeros(size(As[1])[1],size(frames)[3])

map(t -> bsex[:,t] = As[t]*vec(frames[:,:,t]), 1:size(frames)[3])
@info "adding gaussion noise at level 0.01"
rho = 0.01
e = randn(size(bsex));
e = rho*norm(bsex)*e/norm(e);
bs = bsex + e;

As  = map(t -> radon_operator(H,W, detcount, angles[t]),1:size(u0s)[3])

for w_tv = w_tvs
    niter=450
    us_tv = zeros(H,W,size(frames)[3])
    for t = 1:size(frames)[3]
        A = As[t]
        p = bs[:,t]
        u0 = u0s[:,:,t]
        @suppress begin
            us_tv[:,:,t] .= recon2d_tv_primaldual!(us_tv[:,:,t], A, p, niter, w_tv, c)
        end
    end
    for w_flow = w_flows
        niter=50
        @suppress begin
            us_flow = recon2d_tv_primaldual_flow(As, bs, deepcopy(us_tv), niter, 20, w_tv, w_flow, c)
            filenamepng = @sprintf "flow_radn_ang_2ms_%f_%f_.png" w_tv w_flow;
            #PLOT
            plot(Gray.(us_flow[:,:,1]), aspect_ratio=:equal, framestyle=:none, legend=false, yflip=false)
            savefig(filenamepng)
            next!(p)
        end
    end
end
cd(cwd)


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using Distributed
using Random

function _run_benchmark_sample(;
    case::String,
    method::SolutionMethod,
    formulation::Formulation,
    solution_filename::String,
    optimizer,
)::Nothing
    total_time = @elapsed begin
        @info "Reading: $case"
        time_read = @elapsed begin
            instance = read_benchmark(case)
        end
        @info @sprintf("Read problem in %.2f seconds", time_read)
        BLAS.set_num_threads(Threads.nthreads())
        model = build_model(
            instance = instance,
            formulation = formulation,
            optimizer = optimizer,
            variable_names = true,
        )
        @info "Optimizing..."
        BLAS.set_num_threads(1)
        optimize!(model, method)
    end
    @info @sprintf("Total time was %.2f seconds", total_time)

    @info "Writing solution: $solution_filename"
    solution = UnitCommitment.solution(model)
    write("$solution_filename", solution)

    @info "Verifying solution..."
    validate(instance, solution)
    return
end

function _run_benchmark_combination(
    case::String,
    optimizer_name::String,
    optimizer,
    method_name::String,
    method::SolutionMethod,
    formulation_name::String,
    formulation::Formulation,
    trial,
)
    dirname = "results/$optimizer_name/$method_name/$formulation_name/$case"
    function info(msg)
        @info @sprintf(
            "%-8s %-16s %-16s %-16s %-8s %s",
            msg,
            optimizer_name,
            method_name,
            formulation_name,
            trial,
            case
        )
    end
    mkpath(dirname)
    trial_filename = @sprintf("%s/%03d.json", dirname, trial)
    if isfile(trial_filename)
        info("skip")
        return
    end
    info("run")
    open("$trial_filename.log", "w") do file
        redirect_stdout(file) do
            redirect_stderr(file) do
                return _run_benchmark_sample(
                    case = case,
                    method = method,
                    formulation = formulation,
                    solution_filename = trial_filename,
                    optimizer = optimizer,
                )
            end
        end
    end
    return info("done")
end

function _run_benchmarks(;
    cases::Vector{String},
    optimizers::Dict,
    formulations::Dict,
    methods::Dict,
    trials,
)
    combinations = [
        (c, s.first, s.second, m.first, m.second, f.first, f.second, t) for
        c in cases for s in optimizers for f in formulations for
        m in methods for t in trials
    ]
    shuffle!(combinations)
    if nworkers() > 1
        @printf("%24s", "")
    end
    @info @sprintf(
        "%-8s %-16s %-16s %-16s %-8s %s",
        "STATUS",
        "SOLVER",
        "METHOD",
        "FORMULATION",
        "TRIAL",
        "CASE"
    )
    @sync @distributed for c in combinations
        _run_benchmark_combination(c...)
    end
end


text:
#= Implements Monte-Carlo tree search algorithm.

A key assumption is that both the action space and the state space are finite.
Otherwise, we will keep selecting unexplored actions, and no node of depth
higher than one would be added. Thus, after the state and action spaces are
discretized, the algorithm only works with these discrete states and actions.
For transition, we will map the given current state and action pair to the
closest discrete state and action. Additionally, the next state transitioned
into (via sampling) will also be mapped to the closest discrete state.
=#

using
  GridInterpolations,
  StatsBase

export
  SerialMCTS

const NIter = 50
const MaxDepth = 20
const ExploreExploit = 3.0
const Discount = 0.99

typealias ActualState Vector

type StateNode

  n::Vector{Int64}  # number of visits at the node for each action
  qval::Vector{Float64}  # estimated value for each action

  StateNode(nactions::Real) =
    new(zeros(Int64, Int64(nactions)), zeros(Float64, Int64(nactions)))

end

type LazyDiscrete

  varname::AbstractString
  step::Float64

  LazyDiscrete(varname::AbstractString, step::Float64) = new(varname, step)

end

abstract MCTS <: Solver

type SerialMCTS <: MCTS

  niter::Int64  # number of iterations during each action selection
  maxdepth::Int64  # max depth of search tree
  exex::Float64  # exploration-exploitation constant
  discount::Float64  # simulation/rollout discount factor

  tree::Dict{ActualState, StateNode}
  statemap::Dict{AbstractString, LazyDiscrete}
  actionmap::Dict{AbstractString, LazyDiscrete}
  stategrid::RectangleGrid
  actiongrid::RectangleGrid

  SerialMCTS(;
      niter::Real=NIter,
      maxdepth::Real=MaxDepth,
      exex::Real=ExploreExploit,
      discount::Real=Discount) =
    new(
      Int64(niter),
      Int64(maxdepth),
      Float64(exex),
      Float64(discount),
      Dict{ActualState, StateNode}(),
      Dict{AbstractString, LazyDiscrete}(),
      Dict{AbstractString, LazyDiscrete}(),
      RectangleGrid(),
      RectangleGrid())

end

function discretize_statevariable!(mcts::MCTS, varname::AbstractString, step::Real)
  if haskey(mcts.statemap, varname)
    warn(string(
      "state variable ", varname, " already discretized in ValueIteration object, ",
      "replacing existing discretization scheme"))
  end
  mcts.statemap[varname] = LazyDiscrete(varname, Float64(step))
end

function discretize_actionvariable!(mcts::MCTS, varname::AbstractString, step::Real)
  if haskey(mcts.actionmap, varname)
    warn(string(
      "action variable ", varname, " already discretized in ValueIteration object, ",
      "replacing existing discretization scheme"))
  end
  mcts.actionmap[varname] = LazyDiscrete(varname, Float64(step))
end

# returned policy is "live" and grows the search tree
function getpolicy(mdp::MDP, smcts::SerialMCTS)
  function policy(actualstate...)
    return smctsaction!(mdp, smcts, collect(actualstate))
  end
  return policy
end

include("mcts_checks.jl")
include("mcts_solver.jl")


text:
using OpenBLAS_jll, OpenBLAS32_jll, MKL_jll, CompilerSupportLibraries_jll
using Pkg, Artifacts, Base.BinaryPlatforms, Libdl, Test

include("utils.jl")

# Compile `dgemm_test.c` and `sgesv_test.c` against the given BLAS/LAPACK
function run_test((test_name, test_expected_outputs, test_success), libblas_name, libdirs, interface, backing_libs)
    # We need to configure this C build a bit
    cflags = String[
        "-g",
    ]
    if interface == :ILP64
        push!(cflags, "-DILP64")
    end

   # If our GCC is 64-bit but Julia is 32-bit, pass -m32
   if needs_m32()
       push!(cflags, "-m32")
   end

    ldflags = String[
        # Teach it to find that libblas and its dependencies at build time
        ("-L$(pathesc(libdir))" for libdir in libdirs)...,
        "-l$(libblas_name)",
    ]

    if !Sys.iswindows()
        # Teach it to find that libblas and its dependencies at run time
        append!(ldflags, ("-Wl,-rpath,$(pathesc(libdir))" for libdir in libdirs))
    end

    mktempdir() do dir
        @info("Compiling `$(test_name)` against $(libblas_name) ($(backing_libs)) in $(dir)")
        srcdir = joinpath(@__DIR__, test_name)
        make_cmd = `$(make) -sC $(pathesc(srcdir)) prefix=$(pathesc(dir)) CFLAGS="$(join(cflags, " "))" LDFLAGS="$(join(ldflags, " "))"`
        p = run(ignorestatus(make_cmd))
        if !success(p)
            @error("compilation failed", srcdir, prefix=dir, cflags=join(cflags, " "), ldflags=join(ldflags, " "))
        end
        @test success(p)

        env = Dict(
            # We need to tell it how to find CSL at run-time
            LIBPATH_env => append_libpath(libdirs),
            "LBT_DEFAULT_LIBS" => backing_libs,
            "LBT_STRICT" => 1,
            "LBT_VERBOSE" => 1,
        )
        cmd = `$(dir)/$(test_name)`
        p, output = capture_output(addenv(cmd, env))

        expected_return_value = success(p) ^ test_success
        if !expected_return_value
            @error("Test failed", env)
            println(output)
        end
        @test expected_return_value

        # Test to make sure the test ran properly
        has_expected_output = all(occursin(expected, output) for expected in test_expected_outputs)
        if !has_expected_output
            # Uh-oh, we didn't get what we expected.  Time to debug!
            @error("Test failed, got output:")
            println(output)

            # If we're not on CI, launch `gdb`
            if isempty(get(ENV, "CI", ""))
                debugger = Sys.isbsd() ? "lldb" : "gdb"
                @warn("Launching $debugger")
                cmd = `$(debugger) $(cmd)`
                run(addenv(cmd, env))
            end
        end
        @test has_expected_output
    end
end

# our tests, written in C, defined in subdirectories in `test/`
dgemm =         ("dgemm_test", ("||C||^2 is:  24.3384",),                  true)
sgesv =         ("sgesv_test", ("||b||^2 is:   3.0000",),                  true)
sgesv_failure = ("sgesv_test", ("Error: no BLAS/LAPACK library loaded!",), false)
sdot  =         ("sdot_test",  ("C is:   1.9900",),                        true)
zdotc =         ("zdotc_test", (
                     "C (cblas) is:   (  1.4700,   3.8300)",
                     "C (fortran) is: (  1.4700,   3.8300)",
                ),      true)

# Helper function to run all the tests with the given arguments
function run_all_tests(args...; tests = [dgemm, sgesv, sdot, zdotc])
    for test in tests
        run_test(test, args...)
    end
end

# Build version that links against vanilla OpenBLAS
openblas_interface = :LP64
if Sys.WORD_SIZE == 64 && Sys.ARCH != :aarch64
    openblas_interface = :ILP64
end
openblas_jll_libname = splitext(basename(OpenBLAS_jll.libopenblas_path)[4:end])[1]
@testset "Vanilla OpenBLAS_jll ($(openblas_interface))" begin
    run_all_tests(openblas_jll_libname, OpenBLAS_jll.LIBPATH_list, openblas_interface, "")
end

# Build version that links against vanilla OpenBLAS32
@testset "Vanilla OpenBLAS32_jll (LP64)" begin
    # Reverse OpenBLAS32_jll's LIBPATH_list so that we get the right openblas.so
    run_all_tests("openblas", reverse(OpenBLAS32_jll.LIBPATH_list), :LP64, "")
end

# Next, build a version that links against `libblastrampoline`, and tell
# the trampoline to forwards calls to `OpenBLAS_jll`
lbt_link_name, lbt_dir = build_libblastrampoline()
lbt_dir = joinpath(lbt_dir, binlib)

@testset "LBT -> OpenBLAS_jll ($(openblas_interface))" begin
    libdirs = unique(vcat(lbt_dir, OpenBLAS_jll.LIBPATH_list..., CompilerSupportLibraries_jll.LIBPATH_list...))
    run_all_tests("blastrampoline", libdirs, openblas_interface, OpenBLAS_jll.libopenblas_path)
end

# And again, but this time with OpenBLAS32_jll
@testset "LBT -> OpenBLAS32_jll (LP64)" begin
    libdirs = unique(vcat(lbt_dir, OpenBLAS32_jll.LIBPATH_list..., CompilerSupportLibraries_jll.LIBPATH_list...))
    run_all_tests("blastrampoline", libdirs, :LP64, OpenBLAS32_jll.libopenblas_path)
end

# Test against MKL_jll using `libmkl_rt`, which is :LP64 by default
if MKL_jll.is_available()
    @testset "LBT -> MKL_jll (LP64)" begin
        libdirs = unique(vcat(lbt_dir, MKL_jll.LIBPATH_list..., CompilerSupportLibraries_jll.LIBPATH_list...))
        run_all_tests("blastrampoline", libdirs, :LP64, MKL_jll.libmkl_rt_path)
    end

    # Test that we can set MKL's interface via an environment variable to select ILP64, and LBT detects it properly
    if Sys.WORD_SIZE == 64
        @testset "LBT -> MKL_jll (ILP64, via env)" begin
            withenv("MKL_INTERFACE_LAYER" => "ILP64") do
                libdirs = unique(vcat(lbt_dir, MKL_jll.LIBPATH_list..., CompilerSupportLibraries_jll.LIBPATH_list...))
                run_all_tests("blastrampoline", libdirs, :ILP64, MKL_jll.libmkl_rt_path)
            end
        end
    end
end

# Do we have Accelerate available?
veclib_blas_path = "/System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/libBLAS.dylib"
if isfile(veclib_blas_path)
    # Test that we can run BLAS-only tests without LAPACK loaded (`sgesv` test requires LAPACK symbols)
    @testset "LBT -> vecLib/libBLAS" begin
        run_all_tests("blastrampoline", [lbt_dir], :LP64, veclib_blas_path; tests=[dgemm, sdot, zdotc])
    end

    # With LAPACK as well, run all tests
    veclib_lapack_path = "/System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/libLAPACK.dylib"
    @testset "LBT -> vecLib/libLAPACK" begin
        run_all_tests("blastrampoline", [lbt_dir], :LP64, string(veclib_blas_path, ";", veclib_lapack_path))
    end
end

# Do we have a `blas64.so` somewhere?  If so, test with that for fun
blas64 = dlopen("libblas64", throw_error=false)
if blas64 !== nothing
    # Test that we can run BLAS-only tests without LAPACK loaded (`sgesv` test requires LAPACK symbols, blas64 doesn't have CBLAS)
    @testset "LBT -> libblas64 (ILP64, BLAS)" begin
        run_all_tests("blastrampoline", [lbt_dir], :ILP64, dlpath(blas64); tests=[dgemm, sdot])
    end

    # Check if we have a `liblapack` and if we do, run again, this time including `sgesv`
    lapack = dlopen("liblapack64", throw_error=false)
    if lapack !== nothing
        @testset "LBT -> libblas64 + liblapack64 (ILP64, BLAS+LAPACK)" begin
            run_all_tests("blastrampoline", [lbt_dir], :ILP64, "$(dlpath(blas64));$(dlpath(lapack))"; tests=[dgemm, sdot, sgesv])
        end
    end
end

# Finally the super-crazy test: build a binary that links against BOTH sets of symbols!
if openblas_interface == :ILP64
    inconsolable = ("inconsolable_test", ("||C||^2 is:  24.3384", "||b||^2 is:   3.0000"), true)
    @testset "LBT -> OpenBLAS 32 + 64 (LP64 + ILP64)" begin
        libdirs = unique(vcat(OpenBLAS32_jll.LIBPATH_list..., OpenBLAS_jll.LIBPATH_list..., CompilerSupportLibraries_jll.LIBPATH_list..., lbt_dir))
        run_test(inconsolable, lbt_link_name, libdirs, :wild_sobbing, "$(OpenBLAS32_jll.libopenblas_path);$(OpenBLAS_jll.libopenblas_path)")
    end
end

# Run our "direct" tests within Julia
include("direct.jl")

max_stars_repo_name: JuliaLinearAlgebra/libblastrampoline

max_issues_repo_name: JuliaLinearAlgebra/libblastrampoline

max_forks_repo_name: JuliaLinearAlgebra/libblastrampoline


text:
const DEFAULT_PRECISION = Ref{Int}(256)

"""
    precision(<:Union{Arf, Arb, Acb, arf_struct, arb_struct, acb_struct})
    precision(<:Ptr{<:Union{arf_struct, arb_struct, acb_struct}})
    precision(x::Union{arf_struct, arb_struct, acb_struct})
    precision(x::Ptr{<:Union{arf_struct, arb_struct, acb_struct}})
Get the default precision (in bits) currently used for `Arblib` arithmetic.
"""
Base.precision(::Type{<:Union{ArbTypes,ArbStructTypes}}) = DEFAULT_PRECISION[]
Base.precision(::Type{<:Ptr{<:ArbStructTypes}}) = DEFAULT_PRECISION[]

Base.precision(x::ArbStructTypes) = DEFAULT_PRECISION[]
Base.precision(x::Ptr{<:ArbStructTypes}) = DEFAULT_PRECISION[]
Base.precision(x::ArbTypes) = x.prec
# MagLike <: ArbTypes
Base.precision(x::MagLike) = DEFAULT_PRECISION[]
# disambiguation
Base.precision(::MagOrRef) = DEFAULT_PRECISION[]

@inline _precision(x::ArbTypes) = precision(x)
@inline _precision(x::BigFloat) = precision(x)
@inline _precision(z::Complex) = max(_precision(real(z)), _precision(imag(z)))
@inline _precision(
    (a, b)::Tuple{S,T},
) where {
    S<:Union{ArbTypes,BigFloat,Complex{<:Union{ArbTypes,BigFloat}}},
    T<:Union{ArbTypes,BigFloat,Complex{<:Union{ArbTypes,BigFloat}}},
} = max(_precision(a), _precision(b))
@inline _precision(
    (a, b)::Tuple{S,T},
) where {S<:Union{ArbTypes,BigFloat,Complex{<:Union{ArbTypes,BigFloat}}},T} = _precision(a)
@inline _precision(
    (a, b)::Tuple{S,T},
) where {S,T<:Union{ArbTypes,BigFloat,Complex{<:Union{ArbTypes,BigFloat}}}} = _precision(b)
@inline _precision((a, b)::Tuple{S,T}) where {S,T} = max(_precision(a), _precision(b))
@inline _precision(@nospecialize _) = DEFAULT_PRECISION[]

"""
    setprecision(::Type{<:Union{Arf, Arb, Acb}}, precision::Int)
Set the precision (in bits) to be used for `Arblib` arithmetic.
!!! warning
    This function is not thread-safe. It will affect code running on all threads, but
    its behavior is undefined if called concurrently with computations that use the
    setting.
"""
function Base.setprecision(::Type{<:ArbTypes}, precision::Integer)
    if precision < 2
        throw(DomainError(precision, "`precision` cannot be less than 2."))
    end
    DEFAULT_PRECISION[] = precision
    return precision
end

Base.setprecision(x::T, prec::Integer) where {T<:Union{Arf,ArfRef,Arb,ArbRef,Acb,AcbRef}} =
    T(x, prec = prec)
Base.setprecision(v::T, prec::Integer) where {T<:Union{ArbVector,ArbRefVector}} =
    T(v.arb_vec, prec)
Base.setprecision(v::T, prec::Integer) where {T<:Union{AcbVector,AcbRefVector}} =
    T(v.acb_vec, prec)
Base.setprecision(A::T, prec::Integer) where {T<:Union{ArbMatrix,ArbRefMatrix}} =
    T(A.arb_mat, prec)
Base.setprecision(A::T, prec::Integer) where {T<:Union{AcbMatrix,AcbRefMatrix}} =
    T(A.acb_mat, prec)

Base.setprecision(poly::ArbPoly, prec::Integer) = ArbPoly(poly.arb_poly, prec = prec)
Base.setprecision(series::ArbSeries, prec::Integer) =
    ArbSeries(series.arb_poly, degree = degree(series), prec = prec)
Base.setprecision(poly::AcbPoly, prec::Integer) = AcbPoly(poly.acb_poly, prec = prec)
Base.setprecision(series::AcbSeries, prec::Integer) =
    AcbSeries(series.acb_poly, degree = degree(series), prec = prec)


text:
mutable struct Tracker{T<:Integer}
    best_a :: T
    best_error :: BigInt
end
function (t::Tracker)(a, e::BigInt)
    if abs(e)<abs(t.best_error)
        t.best_a=a
        t.best_error=e
    end
    nothing
end
(t::Tracker)() = (t.best_a, t.best_error)
is_error_zero(t::Tracker) = t.best_error == 0


text:
@testset "Merge time series" begin
	data1 = DataFrame(Temp=[10.,20,30,40],
	       datetime=[DateTime(2010,1,1,0),DateTime(2010,1,1,1),
	         DateTime(2010,1,1,2),DateTime(2010,1,1,4)]);
	data2 = DataFrame(grav=[400.,300,200,100],
			        datetime=[DateTime(2010,1,1,0),DateTime(2010,1,1,1),
			          DateTime(2010,1,1,2),DateTime(2010,1,1,3)]);
	data3 = DataFrame(pres=[1000.,2000,3000,4000],
	  		        datetime=[DateTime(2010,1,1,1),DateTime(2010,1,1,2),
	  		          DateTime(2010,1,1,3),DateTime(2010,1,1,4)]);
	dataout = mergetimeseries(data1,data2,data3,timecol=:datetime,kind=:outer)

	@test size(dataout) == (5,4)
	@test Dates.value.(dataout[:datetime]) == Dates.value.(
				collect(DateTime(2010,1,1,0):Dates.Hour(1):DateTime(2010,1,1,4)))

	@test dataout[:Temp][[1,2,3,5]] == data1[:Temp]
	@test isnan(dataout[:Temp][4])

	@test dataout[:grav][[1,2,3,4]] == data2[:grav]
	@test isnan(dataout[:grav][end])

	@test dataout[:pres][[2,3,4,5]] == data3[:pres]
	@test isnan(dataout[:pres][1])
end


text:
using Flux

function naive_implementation()
    W₁ = param(rand(3,5))
    b₁ = param(rand(3))
    layer₁(x) = W₁ * x .+ b₁

    W₂ = param(rand(2,3))
    b₂ = param(rand(2))
    layer₂(x) = W₂ * x .+ b₂

    model(x) = layer₂(σ.(layer₁(x)))
    @show model(rand(5))
end

struct Affine
    W
    b
end

Affine(in::Integer, out::Integer) = Affine(param(randn(out, in)), param(randn(out)))

# overload call, so that object can be used as a function
(m::Affine)(x) = m.W * x .+ m.b

function struct_implementation()
    aff = Affine(10, 5)
    @show aff(rand(10))
end

function main()
    naive_implementation()
    struct_implementation()
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Documenter, ReTest

makedocs(sitename = "ReTest.j",
         modules = [ReTest, ReTest.InlineTest])

deploydocs(
    repo = "github.com/JuliaTesting/ReTest.jl.git",
)


text: solution(str,suffix) = endswith(str,suffix)

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
module EGraphs

include("../docstrings.jl")

using DataStructures
using TermInterface
using TimerOutputs
using Parameters
using Metatheory: alwaystrue, cleanast, binarize, @log
using Metatheory.Patterns
using Metatheory.Rules
using Metatheory.EMatchCompiler
using Dates

import ThreadsX

include("intdisjointmap.jl")
export IntDisjointSet
export in_same_set

include("egraph.jl")
export AbstractENode
export ENodeLiteral
export ENodeTerm
export EClassId
export EClass
export hasdata
export getdata
export setdata!
export find
export lookup
export arity
export EGraph
export merge!
export in_same_class
export addexpr!
export rebuild!
export settermtype!
export gettermtype

include("analysis.jl")
export analyze!
export extract!
export astsize
export astsize_inv
export AbstractAnalysis
export MetadataAnalysis
export getcost!

include("ematch.jl")
export Sub

include("Schedulers.jl")
export Schedulers
using .Schedulers

include("saturation.jl")
export SaturationGoal
export EqualityGoal
export reached
export SaturationParams
export saturate!
export areequal
export @areequal
export @areequalg

end


text:
#    Problem : GROUP A
#    *********
#    A quartic problem whose Hessian is an arrow-head (downwards) with
#    diagonal central part and border-width of 1.
#
#    Origonal SIF Source:
#    A.R. Conn, N.I.M. Gould, M. Lescrenier and Ph.L. Toint,
#    "Performance of a multifrontal scheme for partially separable
#    optimization",
#    Report 88/4, Dept of Mathematics, FUNDP (Namur, B), 1988.
#
#    Implementation translated from Source:
#    http://eprints.tsu.ge/234/14/Tests%20collection-K-F.pdf
#
#    ARWHEAD.SIF classification OUR2-AN-V-0
#
#    Number of variables is variable

f = (x) -> begin
    fx = 0.0
    xn_sqr = x[lastindex(x)]^2
    for i in 1:lastindex(x)-1
        fx += (x[i]^2 + xn_sqr)^2 + 3.0 - 4.0x[i]
    end
    return fx
end

g! = (g, x) -> begin
    gn = 0.0
    cxn = 4.0*x[lastindex(x)]
    xn_sqr = cxn^2/16.0
    for i in firstindex(x):lastindex(x)-1
        α = x[i]^2 + xn_sqr
        g[i] = 4.0x[i]*α - 4.0
        gn += cxn*α
    end
    g[lastindex(x)] = gn
    return g
end

fg! = (g, x) -> begin
    gn = fx = 0.0
    cxn = 4.0*x[lastindex(x)]
    xn_sqr = cxn^2/16.0
    for i in firstindex(x):lastindex(x)-1
        α = x[i]^2 + xn_sqr
        fx += (x[i]^2 + xn_sqr)^2 + 3.0 - 4.0*x[i]
        g[i] = 4.0*x[i]*α - 4.0
        gn += cxn*α
    end
    g[lastindex(x)] = gn
    return fx, g
end

init = (n::Int=5000) -> begin
	n < 2 && @warn("ARWHEAD: number of variables must be ≥ 2")
	n = max(n, 2)

    return n, ones(n)
end

TestSet["ARWHEAD"] = UncProgram("ARWHEAD", f, g!, fg!, init)


text:
struct Mlr
    B::Matrix{Float64}
    int::Matrix{Float64}
    weights::Vector{Float64}
end

"""
    mlr(X, Y, weights = ones(size(X, 1)))
Compute a mutiple linear regression model (MLR) by using the QR algorithm.
* `X` : X-data.
* `Y` : Y-data.
* `weights` : Weights of the observations.

Safe but little slower.

`X` and `Y` are internally centered. The model is computed with an intercept.
"""
function mlr(X, Y, weights = ones(size(X, 1)))
    mlr!(copy(X), copy(Y), weights)
end

function mlr!(X, Y, weights = ones(size(X, 1)))
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    weights = mweights(weights)
    xmeans = colmeans(X, weights)
    ymeans = colmeans(Y, weights)
    center!(X, xmeans)
    center!(Y, ymeans)
    sqrtD = Diagonal(sqrt.(weights))
    B = (sqrtD * X) \ (sqrtD * Y)
    int = ymeans' .- xmeans' * B
    Mlr(B, int, weights)
end

"""
    mlrchol(X, Y, weights = ones(size(X, 1)))
Compute a mutiple linear regression model (MLR)
using the Normal equations and a Choleski factorization.
* `X` : X-data, with nb. columns >= 2 (required by function cholesky).
* `Y` : Y-data.
* `weights` : Weights of the observations.

Faster but can be less accurate (squared element X'X).

`X` and `Y` are internally centered. The model is computed with an intercept.
"""
function mlrchol(X, Y, weights = ones(size(X, 1)))
    mlrchol!(copy(X), copy(Y), weights)
end

function mlrchol!(X, Y, weights = ones(size(X, 1)))
    @assert size(X, 2) > 1 "Method only working for X with > 1 column."
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    weights = mweights(weights)
    xmeans = colmeans(X, weights)
    ymeans = colmeans(Y, weights)
    center!(X, xmeans)
    center!(Y, ymeans)
    XtD = X' * Diagonal(weights)
    B = cholesky!(Hermitian(XtD * X)) \ (XtD * Y)
    int = ymeans' .- xmeans' * B
    Mlr(B, int, weights)
end

"""
    mlrpinv(X, Y, weights = ones(size(X, 1)))
Compute a mutiple linear regression model (MLR)  by using a pseudo-inverse.
* `X` : X-data.
* `Y` : Y-data.
* `weights` : Weights of the observations.

Safe but can be slower.

`X` and `Y` are internally centered. The model is computed with an intercept.
"""
function mlrpinv(X, Y, weights = ones(size(X, 1)))
    mlrpinv!(copy(X), copy(Y), weights)
end

function mlrpinv!(X, Y, weights = ones(size(X, 1)))
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    weights = mweights(weights)
    xmeans = colmeans(X, weights)
    ymeans = colmeans(Y, weights)
    center!(X, xmeans)
    center!(Y, ymeans)
    sqrtD = Diagonal(sqrt.(weights))
    sqrtDX = sqrtD * X
    tol = sqrt(eps(real(float(one(eltype(sqrtDX))))))      # see ?pinv
    B = pinv(sqrtDX, rtol = tol) * (sqrtD * Y)
    int = ymeans' .- xmeans' * B
    Mlr(B, int, weights)
end

"""
    mlrpinv_n(X, Y, weights = ones(size(X, 1)))
Compute a mutiple linear regression model (MLR)
by using the Normal equations and a pseudo-inverse.
* `X` : X-data.
* `Y` : Y-data.
* `weights` : Weights of the observations.

Safe and fast for p not too large.

`X` and `Y` are internally centered. The model is computed with an intercept.
"""
function mlrpinv_n(X, Y, weights = ones(size(X, 1)))
    mlrpinv_n!(copy(X), copy(Y), weights)
end

function mlrpinv_n!(X, Y, weights = ones(size(X, 1)))
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    weights = mweights(weights)
    xmeans = colmeans(X, weights)
    ymeans = colmeans(Y, weights)
    center!(X, xmeans)
    center!(Y, ymeans)
    XtD = X' * Diagonal(weights)
    XtDX = XtD * X
    tol = sqrt(eps(real(float(one(eltype(XtDX))))))
    B = pinv(XtD * X, rtol = tol) * (XtD * Y)
    int = ymeans' .- xmeans' * B
    Mlr(B, int, weights)
end

"""
    mlrvec(x, Y, weights = ones(length(x)))
Compute a simple linear regression model (univariate x).
* `x` : Univariate X-data.
* `Y` : Y-data.
* `weights` : Weights of the observations.

`x` and `Y` are internally centered. The model is computed with an intercept.
"""
function mlrvec(x, Y, weights = ones(length(x)))
    mlrvec!(copy(x), copy(Y), weights)
end

function mlrvec!(x, Y, weights = ones(length(x)))
    @assert size(x, 2) == 1 "Method only working for univariate x."
    x = ensure_mat(x)
    Y = ensure_mat(Y)
    weights = mweights(weights)
    xmeans = colmeans(x, weights)
    ymeans = colmeans(Y, weights)
    center!(x, xmeans)
    center!(Y, ymeans)
    xtD = x' * Diagonal(weights)
    B = (xtD * Y) ./ (xtD * x)
    int = ymeans' .- xmeans' * B
    Mlr(B, int, weights)
end

"""
    coef(object::Mlr)
Compute the coefficients of the fitted model.
* `object` : The fitted model.
"""
function coef(object::Mlr)
    (B = object.B, int = object.int)
end

"""
    predict(object::Mlr, X)
Compute the Y-predictions from the fitted model.
* `object` : The fitted model.
* `X` : X-data for which predictions are computed.
"""
function predict(object::Mlr, X)
    z = coef(object)
    pred = z.int .+ X * z.B
    (pred = pred,)
end


text:
using StochasticDiffEq, DiffEqBase,
      DiffEqProblemLibrary, OrdinaryDiffEq
using Test, Random, Statistics

using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
import DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_2Dlinear,
       prob_sde_additivesystem, prob_sde_lorenz
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear

prob = prob_sde_2Dlinear
prob2 = EnsembleProblem(prob)
sim = solve(prob2,SRIW1(),dt=1//2^(3),trajectories=10)

@test sim[1] isa DiffEqBase.RODESolution
@test sim[1,2] isa Matrix
@test sim[1,2,1] isa Float64

# Test Deprecations
@test_logs (:warn, "num_monte has been replaced by trajectories") (:warn, "parallel_type has been deprecated. Please refer to the docs for the new dispatch-based system.") solve(prob2,SRIW1(),dt=1//2^(3),num_monte=10,parallel_type=:threads)

sim = solve(prob2,SRIW1(),EnsembleThreads(),dt=1//2^(3),trajectories=10)
err_sim = DiffEqBase.calculate_ensemble_errors(sim;weak_dense_errors=true)
@test length(sim) == 10

sim = solve(prob2,SRIW1(),EnsembleThreads(),dt=1//2^(3),trajectories=10,batch_size=2)
err_sim = DiffEqBase.calculate_ensemble_errors(sim;weak_dense_errors=true)
@test length(sim) == 10

sim = solve(prob2,SRIW1(),EnsembleThreads(),dt=1//2^(3),adaptive=false,trajectories=10)
err_sim = DiffEqBase.calculate_ensemble_errors(sim;weak_timeseries_errors=true)

sim = solve(prob2,SRIW1(),EnsembleThreads(),dt=1//2^(3),trajectories=10)
DiffEqBase.calculate_ensemble_errors(sim)
@test length(sim) == 10

sim = solve(prob2,SRIW1(),EnsembleSplitThreads(),dt=1//2^(3),trajectories=10)
DiffEqBase.calculate_ensemble_errors(sim)
@test length(sim) == 10

sim = solve(prob2,SRIW1(),EnsembleSerial(),dt=1//2^(3),trajectories=10)
DiffEqBase.calculate_ensemble_errors(sim)
@test length(sim) == 10

prob = prob_sde_additivesystem
prob2 = EnsembleProblem(prob)
sim = solve(prob2,SRA1(),dt=1//2^(3),trajectories=10)
DiffEqBase.calculate_ensemble_errors(sim)

output_func = function (sol,i)
  last(last(sol))^2,false
end
prob2 = EnsembleProblem(prob,output_func=output_func)
sim = solve(prob2,SRA1(),dt=1//2^(3),trajectories=10)

prob = prob_sde_lorenz
prob2 = EnsembleProblem(prob)
sim = solve(prob2,SRIW1(),dt=1//2^(3),trajectories=10)

output_func = function (sol,i)
  last(sol),false
end

prob = prob_ode_linear
prob_func = function (prob,i,repeat)
  ODEProblem(prob.f,rand()*prob.u0,prob.tspan,1.01)
end

Random.seed!(100)
reduction = function (u,batch,I)
  u = append!(u,batch)
  u,((var(u)/sqrt(last(I)))/mean(u)<0.5) ? true : false
end

prob2 = EnsembleProblem(prob,prob_func=prob_func,output_func=output_func,reduction=reduction,u_init=Vector{Float64}())
sim = solve(prob2,Tsit5(),trajectories=10000,batch_size=20)
@test sim.converged == true

Random.seed!(100)
reduction = function (u,batch,I)
  u = append!(u,batch)
  u,false
end

prob2 = EnsembleProblem(prob,prob_func=prob_func,output_func=output_func,reduction=reduction,u_init=Vector{Float64}())
sim = solve(prob2,Tsit5(),trajectories=100,batch_size=20)
@test sim.converged == false

Random.seed!(100)
reduction = function (u,batch,I)
  u+sum(batch),false
end
prob2 = EnsembleProblem(prob,prob_func=prob_func,output_func=output_func,reduction=reduction,u_init=0.0)
sim2 = solve(prob2,Tsit5(),trajectories=100,batch_size=20)
@test sim2.converged == false
@test mean(sim.u) ≈ sim2.u/100

struct SomeUserType end
output_func = function (sol,i)
    (SomeUserType(),false)
end
prob2 = EnsembleProblem(prob,prob_func=prob_func,output_func=output_func)
sim2 = solve(prob2,Tsit5(),trajectories=2)
@test !sim2.converged && typeof(sim2.u) == Vector{SomeUserType}


text:
using Documenter, Domains

makedocs(sitename="Domains")

deploydocs(
    repo = "github.com/Hirlam/Domains.jl.git",
)


text:
const sample_stats_types = Dict(
    "mean_tree_accept" => Float64,
    "energy" => Float64,
    "energy_error" => Float64,
    "max_energy_error" => Float64,
    "step_size" => Float64,
    "step_size_bar" => Float64,
    "tree_size" => Int,
    "depth" => Int,
    "tune" => Bool,
    "diverging" => Bool,
)

compare(args...; kwargs...) = arviz.compare(args...; kwargs...) |> Pandas.DataFrame

Base.Docs.getdoc(::typeof(compare)) = Base.Docs.getdoc(arviz.compare)

@forwardfun hpd

loo(args...; kwargs...) = arviz.loo(args...; kwargs...) |> Pandas.Series

Base.Docs.getdoc(::typeof(loo)) = Base.Docs.getdoc(arviz.loo)

@forwardfun loo_pit

@forwardfun psislw

r2_score(args...; kwargs...) = arviz.r2_score(args...; kwargs...) |> Pandas.Series

Base.Docs.getdoc(::typeof(r2_score)) = Base.Docs.getdoc(arviz.r2_score)

waic(args...; kwargs...) = arviz.waic(args...; kwargs...) |> Pandas.Series

Base.Docs.getdoc(::typeof(waic)) = Base.Docs.getdoc(arviz.waic)

Base.summary(data::InferenceData) = arviz.summary(data) |> Pandas.DataFrame

Base.Docs.getdoc(::typeof(summary)) = Base.Docs.getdoc(arviz.summary)


text:
@testset "Prefix Serial Function Tests" begin

  @test isdefined(SetAndForget, :prefix_serial) == true

end


text:
module DifferentiationWrapper

using JavaCall

const CalculusFieldElement = @jimport org.hipparchus.CalculusFieldElement
const Class = @jimport java.lang.Class
const DSCompiler = @jimport org.hipparchus.analysis.differentiation.DSCompiler
const DSFactory = @jimport org.hipparchus.analysis.differentiation.DSFactory
const Derivative = @jimport org.hipparchus.analysis.differentiation.Derivative
const DerivativeStructure = @jimport org.hipparchus.analysis.differentiation.DerivativeStructure
const ExtendedUnivariateDifferentiableFunction = @jimport org.hipparchus.analysis.differentiation.ExtendedUnivariateDifferentiableFunction
const ExtendedUnivariateDifferentiableMatrixFunction = @jimport org.hipparchus.analysis.differentiation.ExtendedUnivariateDifferentiableMatrixFunction
const ExtendedUnivariateDifferentiableVectorFunction = @jimport org.hipparchus.analysis.differentiation.ExtendedUnivariateDifferentiableVectorFunction
const FDSFactory = @jimport org.hipparchus.analysis.differentiation.FDSFactory
const Field = @jimport org.hipparchus.Field
const FieldDerivative = @jimport org.hipparchus.analysis.differentiation.FieldDerivative
const FieldDerivativeStructure = @jimport org.hipparchus.analysis.differentiation.FieldDerivativeStructure
const FieldElement = @jimport org.hipparchus.FieldElement
const FieldGradient = @jimport org.hipparchus.analysis.differentiation.FieldGradient
const FieldGradientField = @jimport org.hipparchus.analysis.differentiation.FieldGradientField
const FieldSinCos = @jimport org.hipparchus.util.FieldSinCos
const FieldUnivariateDerivative = @jimport org.hipparchus.analysis.differentiation.FieldUnivariateDerivative
const FieldUnivariateDerivative1 = @jimport org.hipparchus.analysis.differentiation.FieldUnivariateDerivative1
const FieldUnivariateDerivative1Field = @jimport org.hipparchus.analysis.differentiation.FieldUnivariateDerivative1Field
const FieldUnivariateDerivative2 = @jimport org.hipparchus.analysis.differentiation.FieldUnivariateDerivative2
const FieldUnivariateDerivative2Field = @jimport org.hipparchus.analysis.differentiation.FieldUnivariateDerivative2Field
const FiniteDifferencesDifferentiator = @jimport org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator
const Gradient = @jimport org.hipparchus.analysis.differentiation.Gradient
const GradientField = @jimport org.hipparchus.analysis.differentiation.GradientField
const GradientFunction = @jimport org.hipparchus.analysis.differentiation.GradientFunction
const JString = @jimport java.lang.String
const JacobianFunction = @jimport org.hipparchus.analysis.differentiation.JacobianFunction
const MultivariateDifferentiableFunction = @jimport org.hipparchus.analysis.differentiation.MultivariateDifferentiableFunction
const MultivariateDifferentiableVectorFunction = @jimport org.hipparchus.analysis.differentiation.MultivariateDifferentiableVectorFunction
const MultivariateFunction = @jimport org.hipparchus.analysis.MultivariateFunction
const MultivariateVectorFunction = @jimport org.hipparchus.analysis.MultivariateVectorFunction
const Object = @jimport java.lang.Object
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const SparseGradient = @jimport org.hipparchus.analysis.differentiation.SparseGradient
const UnivariateDerivative = @jimport org.hipparchus.analysis.differentiation.UnivariateDerivative
const UnivariateDerivative1 = @jimport org.hipparchus.analysis.differentiation.UnivariateDerivative1
const UnivariateDerivative1Field = @jimport org.hipparchus.analysis.differentiation.UnivariateDerivative1Field
const UnivariateDerivative2 = @jimport org.hipparchus.analysis.differentiation.UnivariateDerivative2
const UnivariateDerivative2Field = @jimport org.hipparchus.analysis.differentiation.UnivariateDerivative2Field
const UnivariateDifferentiableFunction = @jimport org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction
const UnivariateDifferentiableMatrixFunction = @jimport org.hipparchus.analysis.differentiation.UnivariateDifferentiableMatrixFunction
const UnivariateDifferentiableVectorFunction = @jimport org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction
const UnivariateFunction = @jimport org.hipparchus.analysis.UnivariateFunction
const UnivariateFunctionDifferentiator = @jimport org.hipparchus.analysis.differentiation.UnivariateFunctionDifferentiator
const UnivariateMatrixFunction = @jimport org.hipparchus.analysis.UnivariateMatrixFunction
const UnivariateMatrixFunctionDifferentiator = @jimport org.hipparchus.analysis.differentiation.UnivariateMatrixFunctionDifferentiator
const UnivariateVectorFunction = @jimport org.hipparchus.analysis.UnivariateVectorFunction
const UnivariateVectorFunctionDifferentiator = @jimport org.hipparchus.analysis.differentiation.UnivariateVectorFunctionDifferentiator

export DSCompiler
export DSFactory
export Derivative
export DerivativeStructure
export ExtendedUnivariateDifferentiableFunction
export ExtendedUnivariateDifferentiableMatrixFunction
export ExtendedUnivariateDifferentiableVectorFunction
export FDSFactory
export FieldDerivative
export FieldDerivativeStructure
export FieldGradient
export FieldGradientField
export FieldUnivariateDerivative
export FieldUnivariateDerivative1
export FieldUnivariateDerivative1Field
export FieldUnivariateDerivative2
export FieldUnivariateDerivative2Field
export FiniteDifferencesDifferentiator
export Gradient
export GradientField
export GradientFunction
export JacobianFunction
export MultivariateDifferentiableFunction
export MultivariateDifferentiableVectorFunction
export SparseGradient
export UnivariateDerivative
export UnivariateDerivative1
export UnivariateDerivative1Field
export UnivariateDerivative2
export UnivariateDerivative2Field
export UnivariateDifferentiableFunction
export UnivariateDifferentiableMatrixFunction
export UnivariateDifferentiableVectorFunction
export UnivariateFunctionDifferentiator
export UnivariateMatrixFunctionDifferentiator
export UnivariateVectorFunctionDifferentiator
export abs
export acos
export acosh
export add
export add_in_place
export asin
export asinh
export atan
export atan2
export atanh
export build
export cbrt
export ceil
export check_compatibility
export compose
export constant
export copy_sign
export cos
export cosh
export create_constant
export create_variable
export differentiate
export divide
export equals
export exp
export expm1
export floor
export get_all_derivatives
export get_class
export get_compiler
export get_derivative
export get_derivative_field
export get_exponent
export get_factory
export get_field
export get_first_derivative
export get_free_parameters
export get_gradient
export get_instance
export get_nb_points
export get_one
export get_order
export get_partial_derivative
export get_partial_derivative_index
export get_partial_derivative_orders
export get_real
export get_runtime_class
export get_second_derivative
export get_size
export get_step_size
export get_univariate_derivative1_field
export get_univariate_derivative2_field
export get_value
export get_value_field
export get_zero
export hash_code
export hypot
export is_infinite
export is_na_n
export linear_combination
export log
export log10
export log1p
export multiply
export multiply_in_place
export negate
export new_instance
export notify
export notify_all
export num_vars
export pow
export reciprocal
export remainder
export rint
export root_n
export round
export scalb
export signum
export sin
export sin_cos
export sinh
export sqrt
export subtract
export tan
export tanh
export taylor
export to_degrees
export to_derivative_structure
export to_radians
export to_string
export value
export variable
export wait

include("derivative.jl")
include("derivative_structure.jl")
include("ds_compiler.jl")
include("ds_factory.jl")
include("extended_univariate_differentiable_function.jl")
include("extended_univariate_differentiable_matrix_function.jl")
include("extended_univariate_differentiable_vector_function.jl")
include("fds_factory.jl")
include("field_derivative.jl")
include("field_derivative_structure.jl")
include("field_gradient.jl")
include("field_gradient_field.jl")
include("field_univariate_derivative.jl")
include("field_univariate_derivative1.jl")
include("field_univariate_derivative1_field.jl")
include("field_univariate_derivative2.jl")
include("field_univariate_derivative2_field.jl")
include("finite_differences_differentiator.jl")
include("gradient.jl")
include("gradient_field.jl")
include("gradient_function.jl")
include("jacobian_function.jl")
include("multivariate_differentiable_function.jl")
include("multivariate_differentiable_vector_function.jl")
include("sparse_gradient.jl")
include("univariate_derivative.jl")
include("univariate_derivative1.jl")
include("univariate_derivative1_field.jl")
include("univariate_derivative2.jl")
include("univariate_derivative2_field.jl")
include("univariate_differentiable_function.jl")
include("univariate_differentiable_matrix_function.jl")
include("univariate_differentiable_vector_function.jl")
include("univariate_function_differentiator.jl")
include("univariate_matrix_function_differentiator.jl")
include("univariate_vector_function_differentiator.jl")

end


text:
#!/usr/bin/env julia

#=
 This is a modified version of controller_MPC.jl from the barc project ROS codebase.
 Modified by Vijay Govindarajan.

 Licensing Information: You are free to use or extend these projects for
 education or research purposes provided that (1) you retain this notice
 and (2) you provide clear attribution to UC Berkeley, including a link
 to http://barc-project.com

 Attibution Information: The barc project ROS code-base was developed
 at UC Berkeley in the Model Predictive Control (MPC) lab by Jon Gonzales
 (jon.gonzales@berkeley.edu). The cloud services integation with ROS was developed
 by Kiet Lam  (kiet.lam@berkeley.edu). The web-server app Dator was
 based on an open source project by Bruce Wootton
=#

# This is an implementation of MPC using the kinematic bicycle model found here:
# http://www.me.berkeley.edu/~frborrel/pdfpub/IV_KinematicMPC_jason.pdf

__precompile__()

module GPSKinMPCPathFollowerDelay

    using JuMP
    using Ipopt

    #### (1) Initialize model and model parameters and MPC gains ####
	dt_control = 0.05		# control period, ts (s)
    mdl = Model(solver = IpoptSolver(print_level=0, max_cpu_time = dt_control))

	L_a     = 1.5213 		# dist from CoG to front axle (m)
	L_b     = 1.4987 		# dist from CoG to rear axle (m)
	dt      = 0.20			# model discretization time, td (s)
	N       = 8				# horizon

	v_ref = 1.0							    # target velocity for path follower
	x_ref   = v_ref*collect(0.0:dt:N*dt)	# target x(0.0 s), ..., x(N*dt s)
	y_ref   = zeros(N+1)					# target y(0.0 s), ..., y(N*dt s)
    psi_ref = zeros(N+1)					# target psi(0.0 s), ..., psi(N*dt s)

    steer_max = 0.5 		# tire angle (d_f) bound, rad
    steer_dmax = 0.5		# tire angle (d_f) rate bound, rad/s

	a_min = -3.0			# acceleration and deceleration bounds, m/s^2
	a_max = 2.0
	a_dmax = 1.5			# jerk bound, m/s^3

	v_min = 0.0				# vel bounds (m/s)
	v_max = 20.0

	time_constant_acc = 0.4 # s
	alpha = dt_control/(time_constant_acc + dt_control) # LPF gain from https://en.wikipedia.org/wiki/Low-pass_filter

    # Cost function gains.
	@NLparameter(mdl, C_x    == 9.0) # longitudinal deviation
	@NLparameter(mdl, C_y    == 9.0) # lateral deviation
	@NLparameter(mdl, C_psi  == 10.0) # heading deviation
	@NLparameter(mdl, C_v    == 0.0) # target velocity deviation

	@NLparameter(mdl, C_dacc == 100.0) # jerk
	@NLparameter(mdl, C_ddf  == 1000.0) # slew rate
	@NLparameter(mdl, C_acc  == 0.0)   # acceleration
	@NLparameter(mdl, C_df   == 0.0)   # tire angle

	#### (2) Define State/Input Variables and Constraints ####
	# states: position (x,y), velocity (v), heading (psi)
	# inputs: tire angle (d_f) and acceleration (acc).
	#println("Creating kinematic bicycle model ....")
	@variable( mdl, x[1:(N+1)], start=0.0)
	@variable( mdl, y[1:(N+1)], start=0.0)
	#@variable( mdl, v_min <= v[1:(N+1)] <= v_max, start=0.0)
	@variable( mdl, v[1:(N+1)], start=0.0)
	@variable( mdl, psi[1:(N+1)], start=0.0)
	@variable( mdl, acc_lp[1:(N+1)], start=0.0) # low pass filtered acceleration.

	# Input Constraints
	@variable( mdl, a_min <= acc[1:N] <= a_max, start=0.0)
	@variable( mdl, -steer_max <= d_f[1:N] <= steer_max, start=0.0)

	# Input Steering Rate Constraints (with slack)
	@variable( mdl, sl_df[1:N] >= 0.0, start=0.0) # slack on steering
	@NLparameter(mdl, d_f_current == 0.0) # Current tire angle is a model parameter.
	#@NLconstraint(mdl, -steer_dmax*dt_control - sl_df[1] <= d_f[1]  - d_f_current <= steer_dmax*dt_control + sl_df[1])
	@NLconstraint(mdl, d_f[1]  - d_f_current + sl_df[1] >= -steer_dmax*dt_control)
	@NLconstraint(mdl, d_f[1]  - d_f_current - sl_df[1] <= steer_dmax*dt_control)
    for i in 1:(N-1)
        #@constraint(mdl, -steer_dmax*dt - sl_df[i+1] <= d_f[i+1] - d_f[i] <= steer_dmax*dt + sl_df[i+1])
		@constraint(mdl, d_f[i+1] - d_f[i] + sl_df[i+1] >= -steer_dmax*dt)
		@constraint(mdl, d_f[i+1] - d_f[i] - sl_df[i+1] <= steer_dmax*dt)
    end

	# Input Acceleration Rate Constraints (with slack)
	@variable( mdl, sl_acc[1:N] >= 0.0, start=0.0) # slack on acceleration
	@NLparameter(mdl, acc_current == 0.0) # Current acceleration is a model parameter.
	#@NLconstraint(mdl, -a_dmax*dt_control - sl_acc[1] <= acc[1]  - acc_current <= a_dmax*dt_control + sl_acc[1])
	@NLconstraint(mdl, acc[1]  - acc_current + sl_acc[1] >= -a_dmax*dt_control)
	@NLconstraint(mdl, acc[1]  - acc_current - sl_acc[1] <= a_dmax*dt_control)

	for i in 1:(N-1)
        #@constraint(mdl, -a_dmax*dt - sl_acc[i+1] <= acc[i+1] - acc[i] <= a_dmax*dt + sl_acc[i+1])
		@constraint(mdl, acc[i+1] - acc[i] + sl_acc[i+1] >= -a_dmax*dt)
		@constraint(mdl, acc[i+1] - acc[i] - sl_acc[i+1] <= a_dmax*dt)
    end

	#### (3) Define Objective ####

	# Reference trajectory is encoded as model parameters.
	@NLparameter(mdl, x_r[i=1:(N+1)] == x_ref[i])
	@NLparameter(mdl, y_r[i=1:(N+1)] == y_ref[i])
	@NLparameter(mdl, psi_r[i=1:(N+1)] == psi_ref[i])
	@NLparameter(mdl, v_target == v_ref)

	# Cost function.
    @NLobjective(mdl, Min, sum{ C_x*(x[i] - x_r[i])^2 + C_y*(y[i] - y_r[i])^2 + C_psi*(psi[i] - psi_r[i])^2 , i=2:(N+1)} +
                           C_v *sum{ (v[i] - v_target)^2, i = 2:N} +
                           C_acc*sum{(acc[i])^2, i=1:N} +
                           C_df*sum{(d_f[i])^2, i=1:N} +
						   C_dacc*sum{(acc[i+1] - acc[i])^2, i=1:(N-1)} +
                           C_ddf*sum{(d_f[i+1] - d_f[i])^2, i=1:(N-1)} +
						   sum{ sl_df[i] + sl_acc[i], i=1:N}
				)

	#### (4) Define System Dynamics Constraints ####
	# Reference: R.Rajamani, Vehicle Dynamics and Control, set. Mechanical Engineering Series,
	#               Spring, 2011, page 26

    # Initial condition is a model parameter.
	@NLparameter(mdl, x0     == 0.0); @NLconstraint(mdl, x[1]     == x0);
	@NLparameter(mdl, y0     == 0.0); @NLconstraint(mdl, y[1]     == y0);
	@NLparameter(mdl, psi0   == 0.0); @NLconstraint(mdl, psi[1]   == psi0);
    @NLparameter(mdl, v0     == 0.0); @NLconstraint(mdl, v[1]     == v0);

	@NLexpression(mdl, bta[i = 1:N], atan( L_b / (L_a + L_b) * tan(d_f[i]) ) )

	@NLconstraint(mdl, acc_lp[1] == acc_current)
	for i in 1:N
        # equations of motion wrt CoG
		@NLconstraint(mdl, x[i+1]   == x[i]   + dt*( v[i]*cos(psi[i] + bta[i]) ) )
		@NLconstraint(mdl, y[i+1]   == y[i]   + dt*( v[i]*sin(psi[i] + bta[i]) ) )
		@NLconstraint(mdl, psi[i+1] == psi[i] + dt*( v[i]/L_b*sin(bta[i]) ) )
        @NLconstraint(mdl, v[i+1]   == v[i]   + dt*( acc_lp[i] ) )
        @NLconstraint(mdl, acc_lp[i+1] == (1.0-alpha) * acc_lp[i] + alpha * acc[i] )
	end

    #### (5) Initialize Solver ####
	#println("MPC: Initial solve ...")
	status = solve(mdl)
	#println("MPC: Finished initial solve: ", status)

	#################################
	##### State Update Function #####
    function update_init_cond(x::Float64, y::Float64, psi::Float64, vel::Float64)
        # update mpc initial condition
        setvalue(x0,    x)
        setvalue(y0,    y)
        setvalue(psi0,  psi)
        setvalue(v0,    vel)
    end

	#####################################
	##### Reference Update Function #####
    function update_reference(x_ref::Array{Float64,1}, y_ref::Array{Float64,1}, psi_ref::Array{Float64,1}, v_ref::Float64)
    	# update the reference trajectory for MPC.
    	setvalue(x_r[i=1:(N+1)], x_ref)
    	setvalue(y_r[i=1:(N+1)], y_ref)
    	setvalue(psi_r[i=1:(N+1)], psi_ref)
    	setvalue(v_target, v_ref)
    end

	#########################################
	##### Input Update Function #####
	function update_current_input(c_swa::Float64, c_acc::Float64)
		setvalue(d_f_current, c_swa)
		setvalue(acc_current, c_acc)
	end

	#########################################
	##### Cost Update Function #####
	function update_cost(cx::Float64, cy::Float64, cp::Float64, cv::Float64,
						 cda::Float64, cdd::Float64, ca::Float64, cd::Float64)
		setvalue(C_x,   cx)
		setvalue(C_y,   cy)
		setvalue(C_psi, cp)
		setvalue(C_v,   cv)

		setvalue(C_dacc, cda)
		setvalue(C_ddf,  cdd)
		setvalue(C_acc,  ca)
		setvalue(C_df,   cd)
	end

	#################################
	##### Model Solve Function #####
    function solve_model()
        # Solve the model, assuming relevant update functions have been called by the user.
        tic()
        status = solve(mdl)
        solv_time = toq()

        # get optimal solutions
        d_f_opt = getvalue(d_f[1:N])
        acc_opt = getvalue(acc[1:N])
        return acc_opt[1], d_f_opt[1], status, solv_time
    end

	#################################
	##### Diagnostics Function ######
	function get_solver_results()
		# This function gets all relevant solver variables and returns it to the client.
		# Handy for debugging or logging full results.

		# State Variables and Reference
		x_mpc   = getvalue(x[1:(N+1)])
		y_mpc   = getvalue(y[1:(N+1)])
		v_mpc   = getvalue(v[1:(N+1)])
		psi_mpc = getvalue(psi[1:(N+1)])

		x_ref   = getvalue(x_r[1:(N+1)])
		y_ref   = getvalue(y_r[1:(N+1)])
		psi_ref = getvalue(psi_r[1:(N+1)])
		v_ref   = getvalue(v_target)

		# Optimal Solution
        d_f_opt = getvalue(d_f[1:N])
        acc_opt = getvalue(acc[1:N])

		return x_mpc, y_mpc, v_mpc, psi_mpc, x_ref, y_ref, v_ref, psi_ref, d_f_opt, acc_opt
	end
end


text:
using Pkg
tempdir = mktempdir()
Pkg.activate(tempdir)
Pkg.develop(PackageSpec(path=joinpath(@__DIR__, "..")))
Pkg.add(["BenchmarkTools", "PkgBenchmark", "Random"])
Pkg.resolve()

using YAXArrays
using BenchmarkTools

const SUITE = BenchmarkGroup()
SUITE["mapCube"] = include("bench_mapcube.jl")