text:
using PtFEM

data = Dict(
  # Plane(ndim, nst, nxe, nye, nip, direction, finite_element(nod, nodof), axisymmetric)
  :struc_el => Plane(2, 3, 20, 2, 4, :y, Quadrilateral(8, 2), false),
  :properties => [1.0e5 0.0 1.0e-5 1.0e-5;],
  :x_coords => collect(range(0.0, stop=1.0, length=21)),
  :y_coords => collect(range(0.05, stop=-0.05, length=3)),
  :temp => [],
  :support => [
    (1, [0 0]),
    (2, [0 0]),
    (3, [0 0]),
    (4, [0 0]),
    (5, [0 0])
  ],
  :dtemp => [
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50, 0.50, 0.00, -0.50,
    0.50, 0.25, 0.00, -0.25, -0.50
  ],
  :nspr => [
    (163, 2, 50.0)
  ]
)

data |> display
println()

@time m, dis_df, fm_df = p55(data)
println()

dis_df |> display
println()

println("\nThe integration point (nip = $(data[:struc_el].nip)) stresses are:\n")

fm_df |> display
println()

ext: jl


text:
using MethodInSrc
using Test

@testset "auxilliary" begin
    @test ! MethodInSrc.issubdir("/a/b/c/", "/a/g/c/d")
end

# MethodInSrc includes a type `AMatrix` just for testing.
# `sum(::AMatrix)` extends the function `Base.sum`
# `prod(::AMatrix)` is a new function defined in `MethodsInSrc`.
@testset "@isinsrc" begin
    N = 4
    m = MethodInSrc.AMatrix{Int}(N)

    @test @isinsrc sum(m)    # The extension to `Base.sum` is found in "../src/"
    @test sum(m) == N^2      # It gives the correct result.

    @test ! @isinsrc prod(m) # the generic method for `prod(::AbstractMethod)` is invoked.
    @test prod(m) == 1       # The generic method gives the correct result.
    @test @isinsrc MethodInSrc.prod(m) # This finds the function with one method defined in "../src/".
    @test MethodInSrc.prod(m) == 1     # The function in the module also gives the correct result.
end

@testset "@insrc" begin
    N = 3
    m = MethodInSrc.AMatrix{Int}(N)

    # Note that `@insrc` takes the next expression as an argument.
    # We need to test just the function call.
    @test_throws ErrorException 1 == @insrc prod(m)   # This finds the method.
    @test_throws ErrorException @insrc(prod(m)) == 1  # This does too.

    # The following methods are found in "../src", so the expressions are evaluated
    @test N^2 == @insrc sum(m)
    @test (@insrc sum(m)) == N^2
    @test @insrc(MethodInSrc.prod(m)) == 1
end

@testset "@ninsrc" begin
    N = 3
    m = MethodInSrc.AMatrix{Int}(N)
    @test 1 == @ninsrc prod(m)
end

@testset "methods in subdirectory of src" begin
    N = 4
    m = MethodInSrc.AMatrix{Int}(N)
    @test @isinsrc isone(m)
    @test ! @isinsrc iszero(m)
    @test  @isinsrc MethodInSrc.iszero(m)
    @test ! @insrc isone(m)
    @test ! @insrc MethodInSrc.iszero(m)
end


text:
# Note: Any future concrete formulation requires the definition of

# construct_device!(
#     ::OptimizationContainer,
#     ::PSY.System,
#     ::DeviceModel{<:PSY.ACBranch, MyNewFormulation},
#     ::Union{Type{CopperPlatePowerModel}, Type{AreaBalancePowerModel}},
# ) = nothing

#

# Not implemented yet
# struct TapControl <: AbstractBranchFormulation end
# struct PhaseControl <: AbstractBranchFormulation end

#################################### Branch Variables ##################################################
# Because of the way we integrate with PowerModels, most of the time PowerSimulations will create variables
# for the branch flows either in AC or DC.

get_initial_conditions_device_model(
    ::OperationModel,
    ::DeviceModel{T, <:AbstractBranchFormulation},
) where {T <: PSY.ACBranch} = DeviceModel(T, StaticBranch)

get_initial_conditions_device_model(
    ::OperationModel,
    ::DeviceModel{T, <:AbstractBranchFormulation},
) where {T <: PSY.MonitoredLine} = DeviceModel(T, StaticBranchUnbounded)

add_variables!(
    container::OptimizationContainer,
    ::Type{<:AbstractPTDFModel},
    devices::IS.FlattenIteratorWrapper{<:PSY.ACBranch},
    formulation::AbstractBranchFormulation,
) = add_variable!(container, FlowActivePowerVariable(), devices, formulation)

get_variable_binary(
    ::FlowActivePowerVariable,
    ::Type{<:PSY.ACBranch},
    ::AbstractBranchFormulation,
) = false

get_variable_multiplier(_, ::Type{<:PSY.ACBranch}, _) = NaN

function get_default_time_series_names(
    ::Type{U},
    ::Type{V},
) where {U <: PSY.ACBranch, V <: AbstractBranchFormulation}
    return Dict{Type{<:TimeSeriesParameter}, String}()
end

function get_default_attributes(
    ::Type{U},
    ::Type{V},
) where {U <: PSY.ACBranch, V <: AbstractBranchFormulation}
    return Dict{String, Any}()
end
#################################### Flow Variable Bounds ##################################################

function branch_rate_bounds!(
    container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{B},
    ::DeviceModel{B, <:AbstractBranchFormulation},
    ::Type{<:PM.AbstractDCPModel},
) where {B <: PSY.ACBranch}
    var = get_variable(container, FlowActivePowerVariable(), B)
    for d in devices
        name = PSY.get_name(d)
        for t in get_time_steps(container)
            JuMP.set_upper_bound(var[name, t], PSY.get_rate(d))
            JuMP.set_lower_bound(var[name, t], -1.0 * PSY.get_rate(d))
        end
    end
    return
end

function branch_rate_bounds!(
    container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{B},
    ::DeviceModel{B, <:AbstractBranchFormulation},
    ::Type{<:PM.AbstractPowerModel},
) where {B <: PSY.ACBranch}
    vars = [
        get_variable(container, FlowActivePowerFromToVariable(), B),
        get_variable(container, FlowActivePowerToFromVariable(), B),
    ]
    for d in devices
        name = PSY.get_name(d)
        for t in get_time_steps(container), var in vars
            JuMP.set_upper_bound(var[name, t], PSY.get_rate(d))
            JuMP.set_lower_bound(var[name, t], -1.0 * PSY.get_rate(d))
        end
    end
    return
end

################################## Rate Limits constraint_infos ############################

"""
Min and max limits for Abstract Branch Formulation
"""
function get_min_max_limits(
    device,
    ::Type{<:ConstraintType},
    ::Type{<:AbstractBranchFormulation},
) #  -> Union{Nothing, NamedTuple{(:min, :max), Tuple{Float64, Float64}}}
    return (min=-1 * PSY.get_rate(device), max=PSY.get_rate(device))
end

"""
Add branch rate limit constraints for ACBranch with AbstractActivePowerModel
"""
function add_constraints!(
    container::OptimizationContainer,
    cons_type::Type{RateLimitConstraint},
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, U},
    X::Type{<:PM.AbstractActivePowerModel},
) where {T <: PSY.ACBranch, U <: AbstractBranchFormulation}
    add_range_constraints!(container, cons_type, FlowActivePowerVariable, devices, model, X)
    return
end

"""
Add rate limit from to constraints for ACBranch with AbstractPowerModel
"""
function add_constraints!(
    container::OptimizationContainer,
    cons_type::Type{RateLimitConstraintFromTo},
    devices::IS.FlattenIteratorWrapper{B},
    model::DeviceModel{B, <:AbstractBranchFormulation},
    ::Type{T},
) where {B <: PSY.ACBranch, T <: PM.AbstractPowerModel}
    rating_data = [(PSY.get_name(h), PSY.get_rate(h)) for h in devices]

    time_steps = get_time_steps(container)
    var1 = get_variable(container, FlowActivePowerFromToVariable(), B)
    var2 = get_variable(container, FlowReactivePowerFromToVariable(), B)
    add_constraints_container!(
        container,
        cons_type(),
        B,
        [r[1] for r in rating_data],
        time_steps,
    )
    constraint = get_constraint(container, cons_type(), B)

    for r in rating_data
        for t in time_steps
            constraint[r[1], t] = JuMP.@constraint(
                container.JuMPmodel,
                var1[r[1], t]^2 + var2[r[1], t]^2 <= r[2]^2
            )
        end
    end
end

"""
Add rate limit to from constraints for ACBranch with AbstractPowerModel
"""
function add_constraints!(
    container::OptimizationContainer,
    cons_type::Type{RateLimitConstraintToFrom},
    devices::IS.FlattenIteratorWrapper{B},
    model::DeviceModel{B, <:AbstractBranchFormulation},
    ::Type{T},
) where {B <: PSY.ACBranch, T <: PM.AbstractPowerModel}
    rating_data = [(PSY.get_name(h), PSY.get_rate(h)) for h in devices]

    time_steps = get_time_steps(container)
    var1 = get_variable(container, FlowActivePowerToFromVariable(), B)
    var2 = get_variable(container, FlowReactivePowerToFromVariable(), B)
    add_constraints_container!(
        container,
        cons_type(),
        B,
        [r[1] for r in rating_data],
        time_steps,
    )
    constraint = get_constraint(container, cons_type(), B)

    for r in rating_data
        for t in time_steps
            constraint[r[1], t] = JuMP.@constraint(
                container.JuMPmodel,
                var1[r[1], t]^2 + var2[r[1], t]^2 <= r[2]^2
            )
        end
    end
end

"""
Add network flow constraints for ACBranch and NetworkModel with StandardPTDFModel
"""
function add_constraints!(
    container::OptimizationContainer,
    cons_type::Type{NetworkFlowConstraint},
    devices::IS.FlattenIteratorWrapper{B},
    model::DeviceModel{B, <:AbstractBranchFormulation},
    network_model::NetworkModel{S},
) where {B <: PSY.ACBranch, S <: StandardPTDFModel}
    ptdf = get_PTDF(network_model)
    branches = PSY.get_name.(devices)
    time_steps = get_time_steps(container)
    branch_flow = add_constraints_container!(
        container,
        NetworkFlowConstraint(),
        B,
        branches,
        time_steps,
    )
    nodal_balance_expressions = get_expression(container, ActivePowerBalance(), S)
    flow_variables = get_variable(container, FlowActivePowerVariable(), B)
    jump_model = get_jump_model(container)
    for br in devices
        name = PSY.get_name(br)
        ptdf_col = ptdf[name, :]
        flow_variables_ = flow_variables[name, :]
        for t in time_steps
            branch_flow[name, t] = JuMP.@constraint(
                jump_model,
                sum(
                    ptdf_col[i] * nodal_balance_expressions.data[i, t] for
                    i in 1:length(ptdf_col)
                ) - flow_variables_[t] == 0.0
            )
        end
    end
end

"""
Min and max limits for monitored line
"""
function get_min_max_limits(
    device::PSY.MonitoredLine,
    ::Type{<:ConstraintType},
    ::Type{<:AbstractBranchFormulation},
)
    if PSY.get_flow_limits(device).to_from != PSY.get_flow_limits(device).from_to
        @warn(
            "Flow limits in Line $(PSY.get_name(device)) aren't equal. The minimum will be used in formulation $(T)"
        )
    end
    limit = min(
        PSY.get_rate(device),
        PSY.get_flow_limits(device).to_from,
        PSY.get_flow_limits(device).from_to,
    )
    minmax = (min=-1 * limit, max=limit)
    return minmax
end

############################## Flow Limits Constraints #####################################
# TODO: Write tests for these functions
"""
Add branch flow constraints for monitored lines with DC Power Model
"""
function branch_flow_constraints!(
    container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, U},
    X::Type{<:PM.AbstractDCPModel},
) where {T <: PSY.MonitoredLine, U <: AbstractBranchFormulation}
    add_range_constraints!(
        container,
        FlowLimitConstraint,
        FlowActivePowerVariable,
        devices,
        model,
        X,
    )
    return
end

"""
Don't add branch flow constraints for monitored lines if formulation is StaticBranchUnbounded
"""
function branch_flow_constraints!(
    container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, U},
    ::Type{<:PM.AbstractDCPModel},
) where {T <: PSY.MonitoredLine, U <: StaticBranchUnbounded}
    return
end

"""
Min and max limits for flow limit from-to constraint
"""
function get_min_max_limits(
    device::PSY.MonitoredLine,
    ::Type{FlowLimitFromToConstraint},
    ::Type{<:AbstractBranchFormulation},
)
    if PSY.get_flow_limits(device).to_from != PSY.get_flow_limits(device).from_to
        @warn(
            "Flow limits in Line $(PSY.get_name(device)) aren't equal. The minimum will be used in formulation $(T)"
        )
    end
    return (
        min=-1 * PSY.get_flow_limits(device).from_to,
        max=PSY.get_flow_limits(device).from_to,
    )
end

"""
Min and max limits for flow limit to-from constraint
"""
function get_min_max_limits(
    device::PSY.MonitoredLine,
    ::Type{FlowLimitToFromConstraint},
    ::Type{<:AbstractBranchFormulation},
)
    if PSY.get_flow_limits(device).to_from != PSY.get_flow_limits(device).from_to
        @warn(
            "Flow limits in Line $(PSY.get_name(device)) aren't equal. The minimum will be used in formulation $(T)"
        )
    end
    return (
        min=-1 * PSY.get_flow_limits(device).to_from,
        max=PSY.get_flow_limits(device).to_from,
    )
end

"""
Add branch flow constraints for monitored lines
"""
function branch_flow_constraints!(
    container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, U},
    ::Type{<:PM.AbstractPowerModel},
) where {T <: PSY.MonitoredLine, U <: AbstractBranchFormulation}
    add_range_constraints!(
        container,
        FlowLimitFromToConstraint,
        FlowActivePowerFromToVariable,
        devices,
        model,
        X,
    )
    add_range_constraints!(
        container,
        FlowLimitToFromConstraint,
        FlowActivePowerToFromVariable,
        devices,
        model,
        X,
    )
    return
end

"""
Don't add branch flow constraints for monitored lines if formulation is StaticBranchUnbounded
"""
function branch_flow_constraints!(
    container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, U},
    ::Type{<:PM.AbstractPowerModel},
) where {T <: PSY.MonitoredLine, U <: StaticBranchUnbounded}
    return
end


text:
mutable struct Point
    x::Int
    y::Int
end

myParse(x) = parse.(Int, x)

Base.isequal(p1::Point, p2::Point) = p1.x == p2.x && p1.y == p2.y
Base.hash(p1::Point) = (53 + hash(p1.x))*53 + hash(p1.y)

function incrementCount!(d::Dict{Point, Int}, elem::Point)
    d[elem] = get(d, elem, 0) + 1
end

mutable struct Line
    p1::Point
    p2::Point
end

function slope(line::Line)
    return (line.p2.y - line.p1.y)/(line.p2.x - line.p1.x)
end

function countCoveredRange(line::Line, counterDict::Dict{Point, Int})
    m = slope(line)

    x1, y1 = line.p1.x, line.p1.y
    x2, y2 = line.p2.x, line.p2.y

    if m == Inf || m == -Inf
        for point in [Point(x1, y) for y in min(y1, y2):max(y1, y2)]
            incrementCount!(counterDict, point)
        end
    else
        for point in [Point(x, y1 + m*(x-x1)) for x in min(x1, x2):max(x1, x2)]
            incrementCount!(counterDict, point)
        end
    end
end

function countPoints(lines::Vector{Line}, lineFilter = line::Line -> true)
    pointDangerCounter = Dict{Point, Int}()

    foreach(line -> ( if lineFilter(line) countCoveredRange(line, pointDangerCounter) end), lines)

    count = 0

    for (_, coverage) in pointDangerCounter
        if coverage >= 2
            count += 1;
        end
    end

    return count

end

lines = [Line(Point(x[1:2]...), Point(x[3:4]...)) for x in myParse.(split.(readlines("input.txt"), r",|(->)|\s+", keepempty=false))]

#! Part 1

countPoints(lines, line::Line -> line.p1.x == line.p2.x || line.p1.y == line.p2.y)

#! Part 2

countPoints(lines)

max_stars_repo_path: Day 5/main.jl

max_issues_repo_path: Day 5/main.jl

max_forks_repo_path: Day 5/main.jl


text:
include(joinpath(@__DIR__, "..", "src/dynamics", "hopper_3D", "visuals.jl"))
T = Float64
vis = Visualizer()
open(vis)

# get hopper model
s = get_simulation("hopper_3D", "flat_3D_lc", "flat")
model = s.model
env = s.env

nq = model.dim.q
nu = model.dim.u
nc = model.dim.c
nb = model.dim.c * friction_dim(env)
nd = nq + nc + nb
nr = nq + nu + nc + nb + nd

# get trajectory
ref_traj = deepcopy(get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/hopper_3D/gaits/gait_in_place.jld2"),
    load_type = :joint_traj))

H = ref_traj.H
h = ref_traj.h
κ = 1.0e-8

# Objective function
obj = TrackingObjective(model, env, H,
    q = [[Diagonal(1.0e-4 * ones(nq)) for t = 1:H-2]; fill([Diagonal(1.0e-0 * ones(nq))],2)...],
    u = [Diagonal(1.0e-1 * [0.1, 0.1, 1.0]) for t = 1:H],
    γ = [Diagonal(1.0e-100 * ones(nc)) for t = 1:H],
    b = [Diagonal(1.0e-100 * ones(nb)) for t = 1:H])
n_opts = NewtonOptions(r_tol=3e-8, max_iter=100)
im_traj = ImplicitTraj(ref_traj, s, κ=κ)

core = Newton(s, H, h, ref_traj, im_traj, obj = obj, opts = n_opts)

β = 0.10/sqrt(2)
q0_dist = deepcopy(ref_traj.q[1] + [-β,+β,0,0,0,0,0])
q1_dist = deepcopy(ref_traj.q[2] + [-β,+β,0,0,0,0,0])
newton_solve!(core, s, im_traj, ref_traj, q0=q0_dist, q1=q1_dist)

visualize_robot!(vis, model, core.traj.q)

plot(hcat(Vector.(ref_traj.u)...)')
plot(hcat(Vector.(core.traj.u)...)')
# Check ~perfect loop
@show round.(core.traj.q[end-2] - ref_traj.q[1], digits=3)
@show round.(core.traj.q[end-1] - ref_traj.q[2], digits=3)

# Save trajectory
traj = deepcopy(core.traj)
gait_path = joinpath(@__DIR__, "..", "dynamics", "hopper_3D", "gaits", "gait_forward.jld2")
@save gait_path traj

# Reload trajectory
res = JLD2.jldopen(gait_path)
loaded_traj = res["traj"]

traj = deepcopy(get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/hopper_3D/gaits/gait_forward.jld2"),
    load_type = :joint_traj))
plot(hcat(Vector.(traj.q)...)')


text:
#include("SerialDevices/SerialDevices.jl")
include("Utils/UtilDevices.jl")
include("Amplifier/Amplifier.jl")
include("DAQ/DAQ.jl")
include("GaussMeter/GaussMeter.jl")
include("Robots/Robots.jl")
include("SurveillanceUnit/SurveillanceUnit.jl")
include("TemperatureSensor/TemperatureSensor.jl")
include("Virtual/Virtual.jl")

# List our own enums to avoid accidentally converting a different enum
# Did not list enums like LakeShoreF71GaussMeterConnectionModes atm, because their convert function uses specific strings
# and not the enum name
for enum in [:RPTriggerMode]
  @eval begin
    T = $enum
    function Base.convert(::Type{T}, x::String)
      try
        return stringToEnum(x, T)
      catch ex
        throw(ScannerConfigurationError(ex.msg))
      end
    end
  end
end


text:
using Flux
using Flux: train!, mse, params
using HTTP
using JSON

include("./utils.jl")
include("./CTOH.jl")
include("./Data.jl")
include("./SS.jl")

const black = "\e[30;47;1m"
const api_url = "https://coincheck.com/api/ticker"

function get_data()
    r = HTTP.get(api_url)
    data = JSON.parse(String(r.body))
    return collect((data["last"], data["bid"], data["ask"], data["high"], data["low"], data["volume"]))
end

function predict(N, data)
    t = time_ns()
    p = N(data)
    t = Int64(time_ns() - t)
    p2 = (p[1]>=0.5) ? "\e[32;47;1mup" : "\e[31;47;1mdown"
    return p2, t
end

function main()
    println("This software predict the price of BTC.")
    println(black, "time, price, prediction, result(previous)\e[0m")
    d1, d2, d3, d4 = Dense(6, 10, relu), Dense(10, 20), Dense(20, 10, tanh), Dense(10, 1, sigmoid)
    N = Chain(d1, d2, d3, d4)
    PN = Chain(d1, PD(d2.weight, d2.bias), PD(d3.weight, d3.bias), d4)
    DATA = DataKeeper()
    loss(x, y) = mse(N(x), y)
    opt = ADAM()
    LAST = get_data()[1]
    while true
        temp = get_data()
        t = (temp[1] >= LAST) ? 1 : 0
        push!(DATA, (temp, t))
        train!(loss, params(N), DATA.data, opt)
        t = (t==1) ? "\e[32;47;1mup" : "\e[31;47;1mdown"
        p, ti = predict(PN, temp)
        println(black, ti, ", ", temp[1], ", ", p, ", ", t, "\e[0m")
        sleep(300)
    end
end

main()


text:
using Distributions

const DE_DefaultOptions = chain(DEX_DefaultOptions, ParamsDict(
    :SamplerRadius => 8,
))

mutable struct DiffEvoOpt{P<:Population,
                          S<:IndividualsSelector,
                          X<:GeneticOperator,
                          E<:EmbeddingOperator} <: EvolutionaryOptimizer
    # TODO when sampler and bound would be parameterized, name is no longer required -- as everything is seen in the type name
    name::String
    population::P
    so::S        # selection operator
    xo::X        # genetic operator
    eo::E        # embedding operator

    function DiffEvoOpt(name, pop::P, so::S, xo::X, eo::E) where {P,S,X,E}
        return new{P,S,X,E}(name, pop, so, xo, eo)
    end
end

ask(de::DiffEvoOpt) = evolve(de, de.xo)

function evolve(de::DiffEvoOpt, xo::GeneticOperatorsMixture)
    sel_xo, tag = next(xo)
    candidates = evolve(de, sel_xo) # use the selected operator of the mixture
    # override what was set by sel_op so that adjust!() gets called for the mixture op
    for candi in candidates
        candi.extra = xo
        candi.tag = tag
    end
    return candidates
end

function evolve(de::DiffEvoOpt, mutate::MutationOperator)
    target_index = select(de.so, population(de), 1)[1]
    target = acquire_candi(population(de), target_index)
    trial = acquire_candi(population(de), target)
    apply!(mutate, trial.params, target_index)
    return post_process(de, target, trial, mutate)
end

function evolve(de::DiffEvoOpt, crossover::CrossoverOperator)
    # Sample parents and target
    indices = select(de.so, population(de), 1 + numparents(crossover))
    parent_indices = indices[1:numparents(crossover)]
    #println("parent_indices = $(parent_indices)")
    target_index = indices[end]
    target = acquire_candi(population(de), target_index)
    trial = acquire_candi(population(de), target)
    #println("target[$(target_index)] = $(target)")

    # Crossover parents and target
    @assert numchildren(crossover)==1
    apply!(crossover, trial.params, target_index,
           population(de), parent_indices)
    return post_process(de, target, trial, crossover)
end

# Now we can create specific DE optimizers that are commonly used in the
# literature.

function diffevo(problem::OptimizationProblem, options::Parameters, name::String,
                 so::IndividualsSelector = SimpleSelector(),
                 co::DiffEvoCrossoverOperator = DiffEvoRandBin1(chain(DE_DefaultOptions, options)))
    opts = chain(DE_DefaultOptions, options)
    pop = population(problem, opts)
    DiffEvoOpt(name, pop, so, co,
               RandomBound(search_space(problem)))
end

"""
The most used `DE/rand/1/bin` variant of differential evolution.
"""
de_rand_1_bin(problem::OptimizationProblem,
              options::Parameters = EMPTY_PARAMS,
              name = "DE/rand/1/bin") = diffevo(problem, options, name)

de_rand_2_bin(problem::OptimizationProblem,
              options::Parameters = EMPTY_PARAMS,
              name = "DE/rand/2/bin") = diffevo(problem, options, name,
                                                SimpleSelector(),
                                                DiffEvoRandBin2(chain(DE_DefaultOptions, options)))

"""
The most used `DE/rand/1/bin` variant with "local geography" via radius-limited sampling.
"""
de_rand_1_bin_radiuslimited(problem::OptimizationProblem,
                            options::Parameters = EMPTY_PARAMS,
                            name = "DE/rand/1/bin/radiuslimited") =
    diffevo(problem, options, name,
            RadiusLimitedSelector(chain(DE_DefaultOptions, options)[:SamplerRadius]))

de_rand_2_bin_radiuslimited(problem::OptimizationProblem,
                            options::Parameters = EMPTY_PARAMS,
                            name = "DE/rand/2/bin/radiuslimited") =
    diffevo(problem, options, name,
            RadiusLimitedSelector(chain(DE_DefaultOptions, options)[:SamplerRadius]),
            DiffEvoRandBin2(chain(DE_DefaultOptions, options)))


text:
#Módulos para cargar los datos.
using ZipFile, CSV, DataFrames, Dates

#El directorio de las bases de datos.
cd("Open_data/COVID-19")

#Fijamos la fecha de comparación:
fecha_análisis = today()

#Los datos se encuentran en formato .csv al interior del archivo .zip:
datos_zip = ZipFile.Reader(Dates.format(fecha_análisis, "yyyymm")*"/"*Dates.format(fecha_análisis, "yyyymmdd")*".zip")
#Leemos el zip directamente, con un fix para que el tipo de la columna fechas sea String:
datos = CSV.read(datos_zip.files[1], header = 1, types = Dict(:FECHA_DEF => String))

#Generamos el dataframe con los pacientes con resultado positivo (Cf. diccionario/Catálogos.)
P = filter(row -> row[:RESULTADO] == 1, datos)
#Filtramos por entidad y por ubicación de la unidad médica y de residencia.
P_CMX_UM = filter(row -> row[:ENTIDAD_UM] == 9, P)
P_CMX_RES = filter(row -> row[:ENTIDAD_RES] == 9, P)
P_MEX_UM = filter(row -> row[:ENTIDAD_UM] == 15, P)
P_MEX_RES = filter(row -> row[:ENTIDAD_RES] == 15, P)

#Vemos los tamaños, el número de filas es el primer número:
size.([P_CMX_UM, P_CMX_RES, P_MEX_UM, P_MEX_RES])


text:
@everywhere begin
    include("smote_exs.jl")
    include("ub_smote.jl")
    include("Boosting.jl")
    include("simtools.jl")
end

@everywhere function sim(n)
    # running multiple simulations

    # n = 200
    xi = 0.6
    mu_err = 2
    pct_minority = 0.1

    p = 20

    simresults = run_simulations(10, n, p, mu_err, pct_minority, ξ = xi)
    # summarise_sims(simresults)
    writecsv("simres2/$(n)_$(mu_err)_$(xi)_$(p).csv", simresults)
end

# pmap(sim, [parse(x) for x in ARGS])
pmap(sim, [300, 500, 1000])

max_stars_repo_name: paulstey/boosting

max_issues_repo_name: paulstey/boosting

max_forks_repo_name: paulstey/boosting


text:
import PMachine
using Base.Test

r1 = PMachine.Register{ UInt8 }( 1 )

r2 = PMachine.Register{ UInt8 }()

@test ( typeof( r1 ) == typeof( r2 ) )

@test get( r2 ) == 0x00

PMachine.set( r2, UInt8( 1 ) )

@test get( r1 ) == get( r2 )

m = PMachine.Memory{ UInt8, UInt8 }( 10 )

@test PMachine.getAddr( m ) == 0

@test PMachine.read( m ) == 0

PMachine.write( m, UInt8(99) )

@test PMachine.read( m ) == 99

PMachine.setAddr( m, UInt8( 2 ) )
PMachine.reset( m )

@test PMachine.read( m ) == 0
@test PMachine.getAddr( m ) == 0

ma = PMachine.Machine()

@test PMachine.codesize( ma ) == 100
@test PMachine.storesize( ma ) == 100


text:
using Plots,DelimitedFiles,TriplotRecipes

f(x,y) = exp(0.1*sin(5.1*x + -6.2*y) + 0.3*cos(4.3*x + 3.4*y))

x,y = eachcol(readdlm(joinpath(@__DIR__,"dolphin.xy")))
t = readdlm(joinpath(@__DIR__,"dolphin.t"),Int)' .+ 1
z = f.(x,y)
z_dg = similar(t, eltype(z))

for it=1:size(t,2)
    z_dg[:,it] .= z[t[:,it]] .+ 0.05*(rand() - 0.5)
end

plot(aspect_ratio=:equal,size=(800,720))
dgtripcolor!(x,y,z_dg,t,color=:magma)
trimesh!(x,y,t,fillalpha=0.0,linecolor=:white)
savefig("dolphin_dg.png")


text:
function evaluate_function(func::Union{NamedArray,Function},x_val::Real,func_name::Symbol,x_name::Symbol)
    new_dd = DataFrame()
    new_dd[!,x_name] = [x_val]
    feval = isa(func, NamedArray) ? func : func(x_val)
    if isa(feval, NamedArray)
        for name in names(feval)[1]
            new_dd[!,Symbol(name)] .= feval[name]
        end
    else
        error("FunctionEnvelope doesn't know how to interpret what is being returned by function.")
    end
    new_dd[!,:name] .= func_name
    return new_dd
end

function evaluate_function_and_add_to_dd(dd_eval::DataFrame, func::Union{NamedArray,Function},x_val::Real,func_name::Symbol,x_name::Symbol, output_for_envelope::Symbol) where T<:Real
    new_dd = evaluate_function(func,x_val,func_name,x_name)
    full_dd = vcat(dd_eval,new_dd, cols=:union)
    result = new_dd[1,output_for_envelope]
    return result, sort(full_dd, x_name)
end

function get_function_value(dd_eval::DataFrame, dd_function::DataFrame, x_val::Real, func_name::Symbol, x_name::Symbol, output_for_envelope::Symbol) where T<:Real
    if size(dd_eval)[1] > 0
        trying_point_in_cached_data = dd_eval[(dd_eval[!,:name] .== func_name) .& (dd_eval[!,x_name] .== x_val),:]
        if size(trying_point_in_cached_data)[1] > 1
            error("Duplication of same evaluation in dd_eval. Logical problem in code probably")
        elseif size(trying_point_in_cached_data)[1] == 1
            return trying_point_in_cached_data[1, output_for_envelope], dd_eval
        end
    end
    # In this case we don't have a saved result and so must evaluate.
    func = dd_function[dd_function[!,:name] .== func_name,:func][1]
    return evaluate_function_and_add_to_dd(dd_eval, func, x_val,func_name,x_name, output_for_envelope)
end


text:
# StridedView
"""
    struct StridedView{T,N,A<:DenseArray,F<:Union{FN,FC,FA,FT}} <: AbstractStridedView{T,N,F}
        parent::A
        size::NTuple{N,Int}
        strides::NTuple{N,Int}
        offset::Int
        op::F
    end

A concrete type of strided view.
"""
struct StridedView{T,N,A<:DenseArray,F<:Union{FN,FC,FA,FT}} <: AbstractStridedView{T,N,F}
    parent::A
    size::NTuple{N,Int}
    strides::NTuple{N,Int}
    offset::Int
    op::F
end

"""
    StridedView(parent::Array{S},
                size::NTuple{N,Int} = size(parent),
                strides::NTuple{N,Int} = strides(parent),
                offset::Int = 0,
                op::F = identity) where {S, N, F}

A constructor of StridedView instance where `parent` can be an `Array` with any kind of
elements `S`.
"""
function StridedView(parent::Array{S},
                        size::NTuple{N,Int} = size(parent),
                        strides::NTuple{N,Int} = strides(parent),
                        offset::Int = 0,
                        op::F = identity) where {S, N, F}

    T = Base.promote_op(op, S)
    # reshape array to vector in order to reduce number of element types
    StridedView{T,N,Vector{S},F}(reshape(parent, length(parent)), size,
                                    _normalizestrides(size, strides), offset, op)
end

"""
    StridedView(parent::A,
                size::NTuple{N,Int} = size(parent),
                strides::NTuple{N,Int} = strides(parent),
                offset::Int = 0,
                op::F = identity) where {A<:DenseArray, N, F}

A constructor of StridedView instance where `parent` is a general `DenseArray`, but not in
the form of `Array{S}`.
"""
function StridedView(parent::A,
                        size::NTuple{N,Int} = size(parent),
                        strides::NTuple{N,Int} = strides(parent),
                        offset::Int = 0,
                        op::F = identity) where {A<:DenseArray, N, F}

    T = Base.promote_op(op, eltype(parent))
    StridedView{T,N,A,F}(parent, size, strides, offset, op)
end

StridedView(a::StridedView) = a
StridedView(a::Adjoint) = StridedView(a')'
StridedView(a::Transpose) = transpose(StridedView(transpose(a)))
StridedView(a::Base.SubArray) = sview(StridedView(a.parent), a.indices...)
StridedView(a::Base.ReshapedArray) = sreshape(StridedView(a.parent), a.dims)

Base.parent(a::StridedView) = a.parent
Base.size(a::StridedView) = a.size
Base.strides(a::StridedView) = a.strides
offset(a::StridedView) = a.offset

"""
    Base.stride(a::StridedView{<:Any, N}, n::Int) where N

Return the distance in memory (in number of elements) between adjacent elements in
dimension `n`.
"""
Base.stride(a::StridedView{<:Any, 0}, n::Int) = 1
Base.stride(a::StridedView{<:Any, N}, n::Int) where N =
    (n <= N) ? a.strides[n] : a.strides[N]*a.size[N]

Base.IndexStyle(::Type{<:StridedView}) = Base.IndexCartesian()

Base.similar(a::StridedView, ::Type{T}, dims::NTuple{N,Int}) where {N,T} =
    StridedView(similar(a.parent, T, dims))

Base.copy(a::StridedView) = copy!(similar(a), a)

"Return the native address of an array."
Base.unsafe_convert(::Type{Ptr{T}}, a::StridedView{T}) where T =
    pointer(a.parent, a.offset+1)

Base.dataids(a::StridedView) = Base.dataids(a.parent)

@inline function sreshape(a::StridedView, newsize::Dims)
    if any(isequal(0), newsize)
        any(isequal(0), size(a)) || throw(DimensionMismatch())
        newstrides = one.(newsize)
    else
        newstrides = _computereshapestrides(newsize, _simplify(size(a), strides(a))...)
    end
    StridedView(a.parent, newsize, newstrides, a.offset, a.op)
end

"""
    conj(a::StridedView)

Extend `Base.conj`.
Leave the strided view `a` invariant if its elements are real.
Return a new strided view which replace `a.op` with `_conj(a.op)` in general.
"""
Base.conj(a::StridedView{<:Real}) = a
Base.conj(a::StridedView) = StridedView(a.parent, a.size, a.strides, a.offset, _conj(a.op))

"""
    permutedims(a::StridedView{<:Any,N}, p) where {N}

Extend `Base.permutedims`.
Return a new strided view with permuted size and strides according to permutation `p`.
"""
@inline function Base.permutedims(a::StridedView{<:Any,N}, p) where {N}
    (length(p) == N && TupleLabXD.isperm(p)) ||
        throw(ArgumentError("Invalid permutation of length $N: $p"))
    newsize = TupleLabXD._permute(a.size, p)
    newstrides = TupleLabXD._permute(a.strides, p)
    return StridedView(a.parent, newsize, newstrides, a.offset, a.op)
end

"""
    transpose(a::StridedView{<:Number,2})
    transpose(a::StridedView{<:Any,2})

Extend `LinearAlgebra.transpose`.

If the elements of `a.parent` are just numbers, it returns a new strided view with
transposed size and strides.

Otherwise, in general, it returns a new strided view with transposed size and strides
and also replace `a.op` with `_transpose(a.op)`.

Note that the reason for these two behaviors is similar to the difference between
`Base.permutedims` and `LinearAlgebra.transpose` in the original Julia codes.
In the general case, the elements of `a.parent` could be `Array`s themselves, therefore,
the mark `_transpose(a.op)` plays the role which is known as "acting recursively".
"""
LinearAlgebra.transpose(a::StridedView{<:Number,2}) = permutedims(a, (2,1))
LinearAlgebra.transpose(a::StridedView{<:Any,2}) = # act recursively, like Base
    permutedims(StridedView(a.parent, a.size, a.strides, a.offset, _transpose(a.op)), (2,1))

"""
    adjoint(a::StridedView{<:Number,2})
    adjoint(a::StridedView{<:Any,2})

Extend `LinearAlgebra.adjoint`.

If the elements of `a.parent` are just numbers, it returns a new strided view with
transposed size and strides, and also replace `a.op` with `_conj(a.op)`.

Otherwise, in general, it returns a new strided view with transposed size and strides
and also replace `a.op` with `_adjoint(a.op)`.
"""
LinearAlgebra.adjoint(a::StridedView{<:Number,2}) = permutedims(conj(a), (2,1))
LinearAlgebra.adjoint(a::StridedView{<:Any,2}) = # act recursively, like Base
    permutedims(StridedView(a.parent, a.size, a.strides, a.offset, _adjoint(a.op)), (2,1))

"""
    map(::FC, a::StridedView{<:Real})
    map(::FT, a::StridedView{<:Number})
    map(::FA, a::StridedView{<:Number})
    map(::FC, a::StridedView)
    map(::FT, a::StridedView)
    map(::FA, a::StridedView)

Extend `Base.map`. Return a new strided view with the correct `a.op` marked by the type.
"""
Base.map(::FC, a::StridedView{<:Real}) = a
Base.map(::FT, a::StridedView{<:Number}) = a
Base.map(::FA, a::StridedView{<:Number}) = conj(a)
Base.map(::FC, a::StridedView) =
    StridedView(a.parent, a.size, a.strides, a.offset, _conj(a.op))
Base.map(::FT, a::StridedView) =
    StridedView(a.parent, a.size, a.strides, a.offset, _transpose(a.op))
Base.map(::FA, a::StridedView) =
    StridedView(a.parent, a.size, a.strides, a.offset, _adjoint(a.op))

@inline function Base.getindex(a::StridedView{<:Any,N}, I::Vararg{Int,N}) where N
    @boundscheck checkbounds(a, I...)
    @inbounds r = a.op(a.parent[a.offset+_computeind(I, a.strides)])
    return r
end
@inline function Base.setindex!(a::StridedView{<:Any,N}, v, I::Vararg{Int,N}) where N
    @boundscheck checkbounds(a, I...)
    @inbounds a.parent[a.offset+_computeind(I, a.strides)] = a.op(v)
    return a
end

@inline function Base.getindex(a::StridedView{<:Any,N},
        I::Vararg{SliceIndex,N}) where N
    StridedView(a.parent, _computeviewsize(a.size, I), _computeviewstrides(a.strides, I),
                a.offset + _computeviewoffset(a.strides, I), a.op)
end

@propagate_inbounds Base.getindex(a::StridedView, I::ParentIndex) =
    a.op(getindex(a.parent, I.i))
@propagate_inbounds Base.setindex!(a::StridedView, v, I::ParentIndex) =
    (setindex!(a.parent, a.op(v), I.i); return a)

"""
    blasstrides(a::StridedView{<:Any,2})

If `size(a, 2) == 1` and `stride(a, 1) == 1`, return a new stride view with the strides
`(1,size(a,1))` to be compatible with the BLAS function.

Otherwise, return the original `a`.

NOTE: If `size(a,2) == 1`, `stride(a,2)` does not really have a meaning, as you are not
taking any steps in the 2nd dimension. Nonetheless, BLAS still expects that
`stride(a,2) == size(a,1)`. But sometimes in `StridedTensorXD.jl`, if you get this a from
permuting another array, it might be that `stride(a,2)` has some other value.
"""
function blasstrides(a::StridedView{<:Any,2})
    # canonialize strides to make compatible with gemm
    if size(a, 2) == 1 && stride(a, 1) == 1
        return StridedView(a.parent, a.size, (1, size(a,1)), a.offset, a.op)
    else
        return a
    end
end


text:
"""
    cubed_sphere_warp(::EquiangularSphereWarp, a, b, c)

Given points `(a, b, c)` on the surface of a cube, warp the points out to a
spherical shell of radius `R` based on the equiangular gnomonic grid proposed by
[Ronchi1996]

Returns a tuple of (x1, x2, x3), where x1, x2, x3 are the corresponding coordinates on a sphere
of radius R.

The "Cubed Sphere": A New Method for the Solution of Partial Differential Equations in Spherical Geometry
C. RONCHI,*,1 R. IACONO,* AND P. S. PAOLUCCI†
JOURNAL OF COMPUTATIONAL PHYSICS 124, 93–114 (1996)

https://github.com/CliMA/ClimateMachine.jl/blob/master/src/Numerics/Mesh/Topologies.jl
"""
function cubed_sphere_warp(
    ::EquiangularSphereWarp,
    a::FT,
    b::FT,
    c::FT,
) where {FT <: AbstractFloat}

    R = max(abs(a), abs(b), abs(c))
    function f(sR, ξ, η)
        X, Y = tan(π * ξ / 4), tan(π * η / 4)
        ζ1 = sR / sqrt(X^2 + Y^2 + 1)
        ζ2, ζ3 = X * ζ1, Y * ζ1
        ζ1, ζ2, ζ3
    end

    fdim = argmax(abs.((a, b, c)))
    if fdim == 1 && a < 0
        # (-R, *, *) : formulas for Face I from Ronchi, Iacono, Paolucci (1996)
        #              but for us face IV of the developed net of the cube
        x1, x2, x3 = f(-R, b / a, c / a)
    elseif fdim == 2 && b < 0
        # ( *,-R, *) : formulas for Face II from Ronchi, Iacono, Paolucci (1996)
        #              but for us face V of the developed net of the cube
        x2, x1, x3 = f(-R, a / b, c / b)
    elseif fdim == 1 && a > 0
        # ( R, *, *) : formulas for Face III from Ronchi, Iacono, Paolucci (1996)
        #              but for us face II of the developed net of the cube
        x1, x2, x3 = f(R, b / a, c / a)
    elseif fdim == 2 && b > 0
        # ( *, R, *) : formulas for Face IV from Ronchi, Iacono, Paolucci (1996)
        #              but for us face III of the developed net of the cube
        x2, x1, x3 = f(R, a / b, c / b)
    elseif fdim == 3 && c > 0
        # ( *, *, R) : formulas for Face V from Ronchi, Iacono, Paolucci (1996)
        #              but for us face VI of the developed net of the cube
        x3, x2, x1 = f(R, b / c, a / c)
    elseif fdim == 3 && c < 0
        # ( *, *,-R) : formulas for Face VI from Ronchi, Iacono, Paolucci (1996)
        #              but for us face I of the developed net of the cube
        x3, x2, x1 = f(-R, b / c, a / c)
    else
        error(
            "invalid case for cubed_sphere_warp(::EquiangularCubedSphere): $a, $b, $c",
        )
    end
    return x1, x2, x3
end

"""
    cubed_sphere_warp(::EquidistantSphereWarp, a, b, c)

Given points `(a, b, c)` on the surface of a cube, warp the points out to a
spherical shell of radius `R` based on the equidistant gnomonic grid outlined in
[Rancic1996] and [Nair2005]

Returns a tuple of (x1, x2, x3), where x1, x2, x3 are the corresponding coordinates on a sphere
of radius R.

The "Cubed Sphere": A New Method for the Solution of Partial Differential Equations in Spherical Geometry
C. RONCHI,*,1 R. IACONO,* AND P. S. PAOLUCCI†
JOURNAL OF COMPUTATIONAL PHYSICS 124, 93–114 (1996)

A Discontinuous Galerkin Transport Scheme on the Cubed Sphere
Ramachandran D. Nair1, Stephen J. Thomas1, and Richard D. Loft1
https://doi.org/10.1175/MWR2890.1

https://github.com/CliMA/ClimateMachine.jl/blob/master/src/Numerics/Mesh/Topologies.jl
"""
function cubed_sphere_warp(
    ::EquidistantSphereWarp,
    a::FT,
    b::FT,
    c::FT,
) where {FT <: AbstractFloat}
    R = max(abs(a), abs(b), abs(c))
    r = hypot(a, b, c)
    return R * a / r, R * b / r, R * c / r
end


text:
using MsgPack
using Documenter

makedocs(modules=[MsgPack],
         sitename="MsgPack",
         authors="Jarrett Revels and other contributors",
         pages=["API Documentation" => "index.md"])

deploydocs(repo="github.com/JuliaIO/MsgPack.jl.git")


text:
mutable struct OrthogonalSubspaceXMLElement <: MyXMLElement
    el::XMLOrNothing
    loadings::MatrixParameter
    prior::MatrixShrinkageLikelihoods

    function OrthogonalSubspaceXMLElement(loadings::MatrixParameter,
                                          prior::MyXMLElement)
        return new(nothing, loadings, prior)
    end
end


text:
abstract type NumberFold <: AbstractAnalysis end

using TermInterface

function EGraphs.make(an::Type{NumberFold}, g::EGraph, n::ENodeLiteral)
    n.value
end

# This should be auto-generated by a macro
function EGraphs.make(an::Type{NumberFold}, g::EGraph, n::ENodeTerm)
    if exprhead(n) == :call && arity(n) == 2
        op = operation(n)
        args = arguments(n)
        l = geteclass(g, args[1])
        r = geteclass(g, args[2])
        ldata = getdata(l, an, nothing)
        rdata = getdata(r, an, nothing)

        # @show ldata rdata

        if ldata isa Number && rdata isa Number
            if op == :*
                return ldata * rdata
            elseif op == :+
                return ldata + rdata
            end
        end
    end

    return nothing
end

function EGraphs.join(an::Type{NumberFold}, from, to)
    # println("joining!")
    if from isa Number
        if to isa Number
            @assert from == to
        else return from
        end
    end
    return to
end

function EGraphs.modify!(an::Type{NumberFold}, g::EGraph, id::Int64)
    # !haskey(an, id) && return nothing
    eclass = g.classes[id]
    d = getdata(eclass, an, nothing)
    if d isa Number
        newclass, _ = addexpr!(g, d)
        merge!(g, newclass.id, id)
    end
end

EGraphs.islazy(x::Type{NumberFold}) = false


text:
# Author: Mathias Louboutin, mlouboutin3@gatech.edu
# Date: July 2020

export setup_model, parse_commandline, setup_geom
export example_src_geometry, example_rec_geometry, example_model, example_info

"""
Simple 2D model setup used for the tests.
"""

function smooth(v, sigma=3)
    v0 = 1f0 .* v
    for i=-sigma:sigma
        i != 0 && (v0[:, sigma+1:end-sigma] .+= v[:, sigma+i+1:end-sigma+i])
    end
    v0[:, sigma+1:end-sigma] .*= 1/(2*sigma+1)
    return v0
end

"""
Sets up a simple 2D layered model for the wave equation operators tests
"""
function setup_model(tti=false, nlayer=2; n=(301, 151), d=(10., 10.))
    ## Set up model structure
    o = (0., 0.)
    lw = n[2] ÷ nlayer

    v = ones(Float32,n) .* 1.5f0
    vp_i = range(1.5f0, 3.5f0, length=nlayer)
    for i in range(2, nlayer, step=1)
        v[:, (i-1)*lw+ 1:end] .= vp_i[i]  # Bottom velocity
    end

    v0 = smooth(v, 7)
    rho0 = (v .+ .5f0) ./ 2
    # Slowness squared [s^2/km^2]
    m = (1f0 ./ v).^2
    m0 = (1f0 ./ v0).^2

    # Setup model structure
    if tti
        println("TTI Model")
        epsilon = smooth((v0[:, :] .- 1.5f0)/12f0, 3)
        delta =  smooth((v0[:, :] .- 1.5f0)/14f0, 3)
        theta =  smooth((v0[:, :] .- 1.5f0)/4, 3)
        model0 = Model(n, d, o, m0; epsilon=epsilon, delta=delta, theta=theta)
        model = Model(n, d, o, m; epsilon=epsilon, delta=delta, theta=theta)
    else
        model = Model(n,d,o,m, rho0)
        model0 = Model(n,d,o,m0,rho0)
        @test Model(n,d,o,m0; rho=rho0).rho == model0.rho
    end
    dm = model.m - model0.m
    return model, model0, dm
end

"""
Sets up a simple 2D acquisition for the wave equation operators tests
"""

function setup_geom(model; nsrc=1, tn=1500f0)
    ## Set up receiver geometry
    nxrec = model.n[1] - 2
    xrec = collect(range(model.d[1], stop=(model.n[1]-2)*model.d[1], length=nxrec))
    yrec = collect(range(0f0, stop=0f0, length=nxrec))
    zrec = collect(range(50f0, stop=50f0, length=nxrec))

    # receiver sampling and recording time
    T = tn   # receiver recording time [ms]
    dt = .75f0    # receiver sampling interval [ms]

    # Set up receiver structure
    recGeometry = Geometry(xrec, yrec, zrec; dt=dt, t=T, nsrc=nsrc)

    ## Set up source geometry (cell array with source locations for each shot)
    ex =  (model.n[1] - 1) * model.d[1]
    nsrc > 1 ? xsrc = range(.25f0 * ex, stop=.75f0 * ex, length=nsrc) : xsrc = .5f0 * ex
    xsrc = convertToCell(xsrc)
    ysrc = convertToCell(range(0f0, stop=0f0, length=nsrc))
    zsrc = convertToCell(range(2*model.d[2], stop=2*model.d[2], length=nsrc))

    # Set up source structure
    srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dt, t=T)

    # setup wavelet
    f0 = 0.015f0     # MHz
    wavelet = ricker_wavelet(T, dt, f0)
    q = judiVector(srcGeometry, wavelet)

    ntComp = get_computational_nt(srcGeometry, recGeometry, model; dt=dt)
    info = Info(prod(model.n), nsrc, ntComp)

    return q, srcGeometry, recGeometry, info
end

### Process command line args
function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--tti"
            help = "TTI, default False"
            action = :store_true
        "--fs"
            help = "Free surface, default False"
            action = :store_true
        "--nlayer", "-n"
            help = "Number of layers"
            arg_type = Int
            default = 2
        "--parallel", "-p"
            help = "Number of workers"
            arg_type = Int
            default = 1
    end
    return parse_args(s)
end

# Example structures
example_info(; n=(120,100), nsrc=2, ntComp=1000) = Info(prod(n), nsrc, ntComp)
example_model(; n=(120,100), d=(10f0, 10f0), o=(0f0, 0f0), m=randn(Float32, n)) = Model(n, d, o, m)

function example_rec_geometry(; nsrc=2, nrec=120, cut=false)
    startr = cut ? 150f0 : 50f0
    endr = cut ? 1050f0 : 1150f0
    xrec = range(startr, stop=endr, length=nrec)
    yrec = 0f0
    zrec = range(50f0, stop=50f0, length=nrec)
    return Geometry(xrec, yrec, zrec; dt=4f0, t=1000f0, nsrc=nsrc)
end

function example_src_geometry(; nsrc=2)
    xsrc = nsrc == 1 ? 500f0 : range(300f0, stop=900f0, length=nsrc)
    ysrc = nsrc == 1 ? 0f0 : zeros(Float32, nsrc)
    zsrc = range(50f0, stop=50f0, length=nsrc)
    return Geometry(convertToCell(xsrc), convertToCell(ysrc), convertToCell(zsrc); dt=4f0, t=1000f0)
end

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


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

    # Root constructor
    BinaryNode{T}(data) where T = new{T}(data)
    # Child node constructor
    BinaryNode{T}(data, parent::BinaryNode{T}) where T = new{T}(data, parent)
end
BinaryNode(data) = BinaryNode{typeof(data)}(data)

function leftchild(data, parent::BinaryNode)
    !isdefined(parent, :left) || error("left child is already assigned")
    node = typeof(parent)(data, parent)
    parent.left = node
end

function rightchild(data, parent::BinaryNode)
    !isdefined(parent, :right) || error("right child is already assigned")
    node = typeof(parent)(data, parent)
    parent.right = node
end

## Things we need to define
function AbstractTrees.children(node::BinaryNode)
    if isdefined(node, :left)
        if isdefined(node, :right)
            return (node.left, node.right)
        end
        return (node.left,)
    end
    isdefined(node, :right) && return (node.right,)
    return ()
end

## Things that make printing prettier
AbstractTrees.printnode(io::IO, node::BinaryNode) = print(io, node.data)

function module_tree(population)
    ancestory = [moduletracker.parentid for moduletracker in population]
    tree = BinaryNode{Int64}(1)
    for (child, parent) in enumerate(ancestory)
        if child == 1 && parent != 0
            error("no root")
        elseif child != 1
            for leaf in Leaves(tree)
                if leaf.data == parent
                    leftchild(parent, leaf)
                    rightchild(child, leaf)
                    break
                end
            end
        end
    end
    return tree
end


text:
using Turing

function Turing.randclass(urn::PolyaUrn)
  counts = get(current_task(), urn.counts)
  weights = counts ./ sum(counts)
  c = rand(Categorical(weights))
  if c == length(counts)
    urn.counts[length(counts)] = 1
    push!(urn.counts, urn.alpha)
  elseif c < length(counts) && c > 0
    urn.counts[c] = urn.counts[c] + 1
  else
    println(weights)
  end
  return Int(c)::Int
end

function f()
  urn = PolyaUrn(1.72)
  classes = tzeros(Int, 50)
  for i in 1:50
    classes[i]  = randclass(urn)
    u = unique(classes[1:i])
    Base.@assert maximum(u) == length(u)
    # println("[$(current_task())] classes: ", classes[1:i], "; urn:", urn.counts)
    produce(classes[i])
  end
end

t = Task(f)

consume(t);
a = [copy(t) for i = 1:10];

for i =1:20
  consume(t);
  map((x)->consume(x),a)
end


text:
module SparseFEM

using Compat

if VERSION >= v"0.7.0-DEV.3389"
    using SparseArrays
end
if VERSION >= v"0.7.0-DEV.3449"
    using LinearAlgebra
end

# assemble the finite-difference laplacian
function fdlaplacian(N)
    # create a 1D laplacian and a sparse identity
    fdl1 = spdiagm(-1 => ones(N-1), 0 => -2*ones(N), 1 => ones(N-1))
    # laplace operator on the full grid
    I_N = sparse(1.0I, N, N)
    return kron(I_N, fdl1) + kron(fdl1, I_N)
end

# get the list of boundary dof-indices
function get_free(N)
    L = zeros(Int, N, N)
    L[2:N-1, 2:N-1] .= 1
    return findall(!iszero, L)
end

# timing of assembly, slice and solve
function perf_sparse_fem(N)
    Ifree = [LinearIndices((N, N))[i] for i in get_free(N)]
    # assembly
    A = fdlaplacian(N)
    # boundary condition
    B = A[Ifree, Ifree]
    # solver
    @static if VERSION <= v"0.7.0-DEV.5211"
        return lufact(B)
    else
        return lu(B)
    end
end

end # module


text:
function _guess_stkbuf_offset()
    script = joinpath(@__DIR__, "guess_task_layout.jl")
    cmd = `$(Base.julia_cmd()) $script`
    return parse(Int, read(pipeline(cmd; stderr = stderr), String))
end

include_dependency("guess_task_layout.jl")
const STKBUF_OFFSET = _guess_stkbuf_offset()
const BUFSZ_OFFSET = STKBUF_OFFSET + sizeof(UInt)

@inline function getstkbuf(task)
    GC.@preserve task begin
        ptr = pointer_from_objref(task) + STKBUF_OFFSET
        stkbuf = unsafe_load(Ptr{Ptr{Cvoid}}(ptr))
    end
    return stkbuf
end

@inline function getbufsz(task)
    GC.@preserve task begin
        ptr = pointer_from_objref(task) + BUFSZ_OFFSET
        bufsz = unsafe_load(Ptr{Csize_t}(ptr))
    end
    return bufsz
end

function verify_getbufsz()
    default = Int(getbufsz(Task(nothing)))
    @assert getbufsz(Task(nothing, 2 * default)) == 2 * default
    @assert getbufsz(Task(nothing, 4 * default)) == 4 * default
    @assert getbufsz(Task(nothing, 8 * default)) == 8 * default
end

verify_getbufsz()

# Ref: https://github.com/tkf/InferableTasks.jl/blob/master/src/InferableTasks.jl
const HAS_INVOKE_IN_WORLD = isdefined(Base, :invoke_in_world)

const LAST_PTR = Ref(UInt(0))

"""
    SimpleSegmentedStacks.call(f)

Call `f()`.  Allocate the new call stack if the current call stack already
consumed too much stack space.

Error thrown in `f` may be wrapped in `TaskFailedException`.
"""
function SimpleSegmentedStacks.call(f)
    task = current_task()
    stkbuf = getstkbuf(task)
    bufsz = getbufsz(task)
    ptr = sample_stackpointer()
    LAST_PTR[] = UInt(ptr)
    if (
        task === Base.roottask ||  # bufsz of root task is not accurate
        (stkbuf !== C_NULL && ptr > stkbuf && ptr - stkbuf < bufsz ÷ 32)
    )
        iserror = Ref(true)
        parent = current_task()
        world = @static if HAS_INVOKE_IN_WORLD
            ccall(:jl_get_tls_world_age, UInt, ())
        end
        task = @task try
            @static if HAS_INVOKE_IN_WORLD
                local y = Base.invoke_in_world(world, f)
            else
                local y = f()
            end
            iserror[] = false
            yieldto(parent, y)
        catch
            yieldto(parent, nothing)
            rethrow()
        end
        y = yieldto(task)
        if iserror[]
            yield(task)
            wait(task)::Union{}
        else
            @static if HAS_INVOKE_IN_WORLD
                return y::Core.Compiler.return_type(f, Tuple{})
            else
                return y
            end
        end
    else
        return f()
    end
end


text:
using TicTacToe
using Documenter

DocMeta.setdocmeta!(TicTacToe, :DocTestSetup, :(using TicTacToe); recursive=true)

makedocs(;
    modules=[TicTacToe],
    authors="andreas <andreas@lha66.dk> and contributors",
    repo="https://github.com/0708andreas/TicTacToe.jl/blob/{commit}{path}#{line}",
    sitename="TicTacToe.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://0708andreas.github.io/TicTacToe.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/0708andreas/TicTacToe.jl",
)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Run the distributed test outside of the main driver since it needs its own
# set of dedicated workers.
include(joinpath(Sys.BINDIR, "..", "share", "julia", "test", "testenv.jl"))
disttestfile = joinpath(@__DIR__, "distributed_exec.jl")

cmd = `$test_exename $test_exeflags $disttestfile`

if !success(pipeline(cmd; stdout=STDOUT, stderr=STDERR)) && ccall(:jl_running_on_valgrind,Cint,()) == 0
    error("Distributed test failed, cmd : $cmd")
end

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
module TestCases

using ConicHulls, ConicHulls.RefHull

function test(generators::Vector)
    eval(ConicHulls.Hulls,:(numfacets=0))
    eval(ConicHulls.Primitives,:(numgenerators=0))
    NC = length(generators[1])
    test(create_hull(NC), generators)
end
function test(hull::ConicHull, generators::Vector)
    NC, F, G = nconic(hull), ftype(hull), gtype(hull)

    gs = init_hull!(hull, eye(Int, NC))

    verify(hull)
    verify_hull(hull, gs)
    for x in generators
        g = G(x)
        add!(hull, g)
        push!(gs, g)

        verify(hull)
        verify_hull(hull, gs)
    end
end

# Starting from a facet containing the primary generator in find_dominated_facet
test(Vector[[3,0,-4]])

# Starting from a facet with the primary generator in its plane in find_dominated_facet
for k=1:100; test(Vector[ [2,-1,0], [-1,2,0] ]); end

# Mapping g.x => g using a Dict instead of an ObjectIdDict in verify_hull
test(Vector[[1,0,0]])

# Errors in coplanar facter filtering
test(Vector[[1,1,0,-1]])

# Test negative orientation of initial simplex points
test(create_hull(3, -eye(Int, 3)), Vector[[-2,1,0]])

# Test elimination of identical (proportional) generator
test(Vector[[1,0]])
test(Vector[[2,0]])
test(Vector[[1,0,0]])
test(Vector[[2,0,0]])
test(Vector[[1,0,0,0]])
test(Vector[[2,0,0,0]])

# Test elimination of colinear generator
test(Vector[[2,-1]])
test(Vector[[2,0,-1]])
test(Vector[[2,0,0,-1]])

# Test elimination of coplanar generator
test(Vector[[-1,-1,2]])
test(Vector[[-1,-1,0,2]])

# Test four coplanar and conically independent generators
test(Vector[[1,1,-1]])
test(Vector[[1,1,-1,0]])

end # module


text:
Tables.istable(::Type{TableTestSet}) = true
Tables.rowaccess(::Type{TableTestSet}) = true

function Tables.rows(ts::TableTestSet)
    table = Vector{NamedTuple{(:testset, :pass, :fail, :error, :broken, :total, :depth), Tuple{String, Int, Int, Int, Int, Int, Int}}}()
    populate_table!(table, ts, 0)
    table
end

# Recursive function that prints out the results at each level of
# the tree of test sets
function populate_table!(table, ts::TableTestSet, depth)
    # Count results by each type at this level, and recursively
    # through any child test sets
    passes, fails, errors, broken, c_passes, c_fails, c_errors, c_broken = get_test_counts(ts)
    subtotal = passes + fails + errors + broken + c_passes + c_fails + c_errors + c_broken
    np = passes + c_passes
    nf = fails + c_fails
    ne = errors + c_errors
    nb = broken + c_broken
    push!(table, (testset=ts.description, pass=np, fail=nf, error=ne, broken=nb, total=subtotal, depth=depth))
    # Only print results at lower levels if we had failures
    if np + nb != subtotal
        for t in ts.results
            if isa(t, TableTestSet)
                populate_table!(table, t, depth + 1)
            end
        end
    end
end


text:
using NeuroPlots
using Test

@testset "NeuroPlots.jl" begin
    # Write your tests here.
end


text:
#!/usr/bin/env julia

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "-M"
    help = "Print out the list of files this program will create."
    action = :store_true
    "config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using DiskTracer.model
using DiskTracer.constants
using DiskTracer.geometry
using SpecialFunctions

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]
pars = convert_dict(config["parameters"], config["model"])

mol = mols[species*transition]

nr=128
rmax = 500 * AU
nz=128
zmax = 200 * AU
# Get the interpolator object for this particular setup.
interp = DiskTracer.model.get_grids(pars, mol, mol.nu_0, nr, rmax, nz, zmax)

using LaTeXStrings
using Plots
pyplot()

# pick a pixel (xprime, yprime), and trace a full ray through the disk

function plot_ray(xprime, yprime)

    v0 = 1.0 # km/s

    # Get a decent estimate of DeltaVmax
    # Calculate DeltaVmax for this ray
    rcyl_min = abs(xprime)
    DeltaVmax = sqrt(2 * kB * temperature(rcyl_min, pars)/mol.mol_weight + (pars.ksi * 1e5)^2) * 1e-5 # km/s

    # DeltaVmax = 0.5 # km/s
    rmax = 600 * AU
    nzp = 512
    zps = linspace(-600 * AU, 600 * AU, nzp)

    rcyls = Array{Float64}(nzp)
    zs = Array{Float64}(nzp)
    Ss = Array{Float64}(nzp)
    Deltavs = Array{Float64}(nzp)
    tildeDeltavs = Array{Float64}(nzp)
    DeltaVs = Array{Float64}(nzp)
    Us = Array{Float64}(nzp)
    tildeUs = Array{Float64}(nzp)
    Deltas = Array{Float64}(nzp)
    alphas = Array{Float64}(nzp)
    tildealphas = Array{Float64}(nzp)
    tildealphasG = Array{Float64}(nzp)

    z0s = DiskTracer.geometry.get_zps(xprime, yprime, pars, v0)
    println("z0s ", z0s)

    # Get the minimum rcyl
    rcyl_min = xprime

    # Occurs at what value of zprime?
    zprime_rcyl_min = cosd(pars.incl)/sind(pars.incl) * yprime

    # Actually, we want the zprime that corresponds to the midplane crossing
    zprime_midplane = -tand(pars.incl) * yprime
    # Calculate the rcyl_mid here

    # Get amplitude at midplane crossing
    # Corresponds to r_cyl_mid, 0.0
    r_cyl_mid = sqrt(xprime^2 + (cosd(pars.incl) * yprime + sind(pars.incl) * zprime_midplane)^2)

    # println("rcyl min ", rcyl_min/AU, " zprime_rcyl_min ", zprime_rcyl_min/AU, " zprime_midplane ", zprime_midplane/AU, " rcyl_mid ", r_cyl_mid/AU)

    # Get scale height at rcyl_min (will be an underestimate)
    H_rcyl_min = DiskTracer.model.Hp(rcyl_min, pars)
    H_rcyl_mid = DiskTracer.model.Hp(r_cyl_mid, pars)

    DeltaV = sqrt(2 * kB * DiskTracer.model.temperature(rcyl_min, pars)/mol.mol_weight + (pars.ksi * 1e5)^2)

    sigma_Upsilon = 10 * H_rcyl_mid / cosd(pars.incl)
    a_Upsilon = Upsilon_nu(r_cyl_mid, 0.0, mol.nu_0, pars, mol)
    mu_Upsilon = zprime_midplane

    DeltaV = sqrt(2 * kB * DiskTracer.model.temperature(rcyl_min, pars)/mol.mol_weight + (pars.ksi * 1e5)^2)

    mu_upsilon1, mu_upsilon2 = z0s
    sigma_upsilon = DeltaV * sqrt(2)/3 * (xprime * sqrt(G * pars.M_star * M_sun) * sind(pars.incl))^(4./3) / (mu_upsilon1 * (v0 * 1.e5)^(7/3))
    # a_upsilon = sqrt(2pi) * sigma_upsilon

    sigma_alpha2 = 1/(sigma_Upsilon^(-2) + sigma_upsilon^(-2)) # actually sigma^2, hence the 2 in the name
    mu_alpha1 = sigma_alpha2 * (mu_Upsilon/sigma_Upsilon^2 + mu_upsilon1/sigma_upsilon^2)
    mu_alpha2 = sigma_alpha2 * (mu_Upsilon/sigma_Upsilon^2 + mu_upsilon2/sigma_upsilon^2)

    a_alpha1 = a_Upsilon * exp(- (mu_Upsilon - mu_upsilon1)^2/(2 * (sigma_Upsilon^2 + sigma_upsilon^2)))
    a_alpha2 = a_Upsilon * exp(- (mu_Upsilon - mu_upsilon2)^2/(2 * (sigma_Upsilon^2 + sigma_upsilon^2)))

    for i=1:nzp
        zprime = zps[i]
        rcyl = DiskTracer.geometry.get_r_cyl(xprime, yprime, zprime, pars.incl)
        z = DiskTracer.geometry.get_z(xprime, yprime, zprime, pars.incl)

        rcyls[i] = rcyl
        zs[i] = z

        Deltav = v0 * 1e5 - DiskTracer.geometry.get_vlos(xprime, yprime, zprime, pars)
        Deltavs[i] = Deltav

        if zprime > 0
            tildeDeltav = DiskTracer.geometry.get_Delta_v_approx(xprime, yprime, zprime, pars, v0, z0s[1])
        else
            tildeDeltav = DiskTracer.geometry.get_Delta_v_approx(xprime, yprime, zprime, pars, v0, z0s[2])
        end
        tildeDeltavs[i] = tildeDeltav * 1e-5 #km/s

        DeltaV = sqrt(2 * kB * DiskTracer.model.temperature(rcyl, pars)/mol.mol_weight + (pars.ksi * 1e5)^2)
        DeltaVs[i] = DeltaV * 1e-5 # km/s

        Deltas[i] = exp(-Deltav^2/DeltaV^2)

        Us[i] = DiskTracer.model.Upsilon_nu(rcyl, z, mol.nu_0, pars, mol)
        tildeUs[i] = a_Upsilon * exp(-0.5 * (zprime - zprime_midplane)^2/sigma_Upsilon^2)

        alphas[i] = Us[i] * exp(-Deltav^2/DeltaV^2)
        # tildealphas[i] = tildeUs[i] * exp(-tildeDeltavs[i]^2 / DeltaVs[i]^2)

        if zprime > 0
            tildealphasG[i] = a_alpha1 * exp(-0.5 * (zprime - mu_alpha1)^2/sigma_alpha2)
        else
            tildealphasG[i] = a_alpha2 * exp(-0.5 * (zprime - mu_alpha2)^2/sigma_alpha2)
        end
    end

    # get bounding zps
    zbounds = DiskTracer.geometry.get_bounding_zps(xprime, yprime, pars, v0, DeltaVmax, rmax)
    if length(zbounds) == 2
        zstart, zend = zbounds
        println("H bound verified merged: ", DiskTracer.geometry.verify_h_bound(xprime, yprime, zstart, zend, pars, 0.5))
    else
        z1start, z1end, z2start, z2end = zbounds
        println("H bound verified front-side: ", DiskTracer.geometry.verify_h_bound(xprime, yprime, z1start, z1end, pars, 0.5))
        println("H bound verified rear-side: ", DiskTracer.geometry.verify_h_bound(xprime, yprime, z2start, z2end, pars, 0.5))
    end

    # Call the routine that actually traces the ray, and returns zp, tau, and I quantities
    sol = DiskTracer.trace.trace_pixel(xprime, yprime, v0, mol, pars, interp)
    zpts = sol.t

    # Also call our custom integrator to get tau
    zpts_c, tau_c = DiskTracer.trace_simple.trace_pixel(xprime, yprime, v0, mol, pars, interp)

    # tildetau = Array{Float64}(length(zpts))
    # for i=1:length(zpts)
    #     z = zpts[i]
    #     tildetau[i] = a_alpha1 * erfc((z - mu_c1)/(sqrt(2 * sigma_c2)))/(2 * sqrt(pi))
    # end
    #
    #
    # "Given tilde_tau, invert to find zprime"
    # function tilde_tau_inv(tilde_tau)
    #     zprime = sqrt(2 * sigma_c2) * erfcinv(tilde_tau * 2 * sqrt(pi)/a_alpha1) + mu_c1
    # end
    #
    # zpts_inv = Array{Float64}(length(zpts))
    # for i=length(zpts)
    #     zpts_inv[i] = tilde_tau_inv(tildetau[i])
    # end

    p1 = plot(zps./AU, Us, xlabel=L"$z^\prime$", ylabel=L"$\Upsilon$") #, yaxis=(scale=:log10))

    p2 = plot(zps./AU, tildeUs, ylabel=L"$\tilde{\Upsilon}$") #, yaxis=(scale=:log10))
    vline!((zbounds./AU)', line=([:red :green :blue :cyan]))

    p3 = plot(zps./AU, alphas, xlabel=L"$z^\prime$", ylabel=L"$\alpha$")
    vline!((zbounds./AU)', line=([:red :green :blue :cyan]))
    # plot!(zps./AU, tildealphas, ylabel=L"$\tilde{\alpha}$")
    plot!(zps./AU, tildealphasG, ylabel=L"$\tilde{\alpha}$")

    p4 = plot(zps./AU, alphas, xlabel=L"$z^\prime$", ylabel=L"$\alpha$", yaxis=(scale=:log10))
    vline!((zbounds./AU)', line=([:red :green :blue :cyan]))
    # plot!(zps./AU, tildealphas, ylabel=L"$\tilde{\alpha}$")
    plot!(zps./AU, tildealphasG, ylabel=L"$\tilde{\alpha}$", yaxis=(scale=:log10))

    # p5 = plot(zpts./AU, sol[1,:], xlim=(-600, 600), ylabel=L"$\tau$")
    p55 = plot(zpts_c./AU, tau_c, xlim=(-600, 600), ylabel=L"$\tau$")
    # plot!(zpts./AU, tildetau, xlim=(-600, 600))

    # p6 = plot(zpts./AU, sol[2,:], xlim=(-600, 600), ylabel=L"$I$")
    # p6 = plot(zpts./AU, zpts_inv./AU, xlim=(-600, 600))

    # println("final ", sol[2,end])

    p = plot(p1,p2,p3,p4,p55, layout=(7,1), legend=false, size=(800,1200))

    savefig(p, "approx_ray_trace.png")

    # println(DiskTracer.trace.trace_pixel_tilde(xprime, yprime, v0, mol, pars, interp))

end

plot_ray(100.*AU, -100.0*AU)


text:
import GitHub
import Dates
import Plots

if length(ARGS) <= 2
    return
end

owner = ARGS[1]
auth = GitHub.authenticate(ENV["GITHUB_TOKEN"])

params = Dict("type" => "public", "sort" => "created")
repos, next = GitHub.repos(owner; auth = auth, params = params)

range = Dict()
date = floor(Dates.now() + Dates.Month(1), Dates.Month)
for repo in repos
    if date > repo.created_at
       global date -= Dates.Month(1)
    end
    if haskey(range, date)
        push!(range[date], repo)
    else
        range[date] = [repo]
    end
end

arr = collect(range)
sort!(arr, by = x -> x.first)

count = map(x -> Pair(x.first, length(x.second)), arr)
println(count)

x = map(x -> x.first, count)
y = map(x -> x.second, count)

plot = Plots.plot(x, y)
bar = Plots.bar(rand(50))

Plots.png(plot, "plot.png")
Plots.png(bar, "bar.png")

max_stars_repo_name: CORDEA/grepon

max_issues_repo_name: CORDEA/grepon

max_forks_repo_name: CORDEA/grepon


text:
metric(x :: AbstractArray, y :: AbstractArray) =
    mapreduce((x, y) -> (x - y)^2, +, x, y)

function sumcosts(tracker1 :: AbstractArray{T, N},
                  tracker2 :: AbstractArray{T, N},
                  costfn   :: Function) where {T, N}
    @assert tracked_data(tracker1) == tracked_data(tracker2)
    return mapreduce(data -> costfn(tracker1 |> data,
                                    tracker2 |> data),
                     +, tracked_data(tracker1))
end

function sumcosts(tracker1 :: AbstractArray{T, N},
                  tracker2 :: AbstractArray{T, N},
                  costfn   :: Function,
                  weights  :: Dict{AbstractTracker{T}, Float64}) where {T, N}
    @assert tracked_data(tracker1) == tracked_data(tracker2)
    return mapreduce(data -> costfn(tracker1 |> data,
                                    tracker2 |> data) / weights[data],
                     +, tracked_data(tracker1))
end

# All cost functions accept two systems as their arguments: and return
# a distance between those systems. We will then try to minimize that
# distance.

"""
    euclid_mean(data1 :: CorrelationData, data2 :: CorrelationData)
    euclid_mean(data1 :: AbstractArray, data2 :: AbstractArray)

Calculate squared euclidean distance between values of two correlation
functions (represented as `CorrelationData` object) or two systems
(represented as `AbstractArray` object). The values are averaged
along all directions before calculation.
"""
function euclid_mean end

euclid_mean(data1 :: CorrelationData,
            data2 :: CorrelationData) =
                metric(data1 |> mean, data2 |> mean)

euclid_mean(data1 :: AbstractArray{T, N},
            data2 :: AbstractArray{T, N}) where {T, N} =
                sumcosts(data1, data2, euclid_mean)

"""
    euclid_directional(data1 :: CorrelationData, data2 :: CorrelationData)
    euclid_directional(data1 :: AbstractArray, data2 :: AbstractArray)

Calculate squared euclidean distance between values of two correlation
functions (represented as `CorrelationData` object) or two systems
(represented as `AbstractArray` object). The values calculated
along different directions are not averaged and treated separately.
"""
function euclid_directional end

function euclid_directional(data1 :: CorrelationData,
                            data2 :: CorrelationData)
    @assert directions(data1) == directions(data2)
    sum(metric(data1[dir], data2[dir]) for dir in directions(data1))
end

function euclid_directional(data1 :: AbstractArray{T, N},
                            data2 :: AbstractArray{T, N}) where {T, N}
    sumcosts(data1, data2, euclid_directional)
end

"""
    euclid_mean_weighted(data1 :: AbstractArray, data2 :: AbstractArray)

Calculate squared euclidean distance between values of two correlation
functions (represented as `CorrelationData` object) or two systems
(represented as `AbstractArray` object). The values are averaged
along all directions before calculation.

Each correlation function has its own weight as described in
Kirill M. Gerke and Marina V. Karsanina 2015 EPL 111 56002
"""
function euclid_mean_weighted(data1 :: AbstractArray{T, N},
                              data2 :: AbstractArray{T, N}) where {T, N}
    @assert tracked_data(data1) == tracked_data(data2)
    # At first, calculate weights for all correlation functions
    weights = Dict(data => euclid_mean(data1 |> data,
                                       data2 |> data)
                   for data in tracked_data(data1))
    f(:: Any, :: Any) = sumcosts(data1, data2, euclid_mean, weights)
    return f
end

"""
    euclid_directional_weighted(data1 :: AbstractArray, data2 :: AbstractArray)

Calculate squared euclidean distance between values of two correlation
functions (represented as `CorrelationData` object) or two systems
(represented as `AbstractArray` object). The values calculated
along different directions are not averaged and treated separately.

Each correlation function has its own weight as described in
Kirill M. Gerke and Marina V. Karsanina 2015 EPL 111 56002
"""
function euclid_directional_weighted(data1 :: AbstractArray{T, N},
                                     data2 :: AbstractArray{T, N}) where {T, N}
    @assert tracked_data(data1) == tracked_data(data2)
    # At first, calculate weights for all correlation functions
    weights = Dict(data => euclid_directional(data1 |> data,
                                              data2 |> data)
                   for data in tracked_data(data1))
    f(:: Any, :: Any) = sumcosts(data1, data2, euclid_directional, weights)
    return f
end

"""
    generalized_čapek_cost(tracker1, tracker2, dict)

Returns a function which calculates the cost based on two-point and
lineal-path functions for solid phase and any other functions present
in `keys(dict)`. `dict` is a dictionary which includes key-value pairs
`AbstractTracker => Float64`. Values of `dict` must be in range
`[0, 1]`. They control an initial contribution of the corresponding
correlation function. The smaller these values are the smaller is
the initial contribution. Contributions of all correlation functions
equivalizes when the cost function becomes small.
"""
function generalized_čapek_cost(data1 :: AbstractArray{T, N},
                                data2 :: AbstractArray{T, N},
                                dict  :: Dict{<:AbstractTracker{T}, Float64}) where {T, N}
    # Initial S2
    s2cost_init  = euclid_directional(Directional.s2(data1, 0),
                                      Directional.s2(data2, 0))
    # Initial L2 for solid phase
    l2scost_init = euclid_directional(Directional.l2(data1, 1),
                                      Directional.l2(data2, 1))

    # Rescale coefficients as in čapek_cost
    sdict = Dict(k => v*(s2cost_init + l2scost_init) for (k, v) in dict)

    local function f(data1 :: AbstractArray, data2 :: AbstractArray)
        # S2
        s2cost  = euclid_directional(Directional.s2(data1, 0),
                                     Directional.s2(data2, 0))
        # L2 for solid phase
        l2scost = euclid_directional(Directional.l2(data1, 1),
                                     Directional.l2(data2, 1))
        # All other cost functions
        other = sum(η/(η + s2cost + l2scost) *
                    euclid_directional(corrfn(data1),
                                       corrfn(data2))
                    for (corrfn, η) in sdict)

        return s2cost + l2scost + other
    end

    return f
end

"""
    čapek_cost(data1 :: AbstractArray, data2 :: AbstractArray, η = 0.6)

Returns a function which calculates the cost based on two-point
correlation function for void phase  and lineal-path function for
solid and void phases where contribution of the latter increases with
time.
"""
čapek_cost(data1 :: AbstractArray{T, N},
           data2 :: AbstractArray{T, N},
           η     :: Float64 = 0.6) where {T, N} =
               generalized_čapek_cost(data1, data2, Dict(T |> zero |> L2Tracker => η))


text:
export row_swap!, row_mul!, row_add!, ref!, rref!,
       row_swap,  row_mul,  row_add,  ref,  rref

# Call f on a copy of A and return the (modified) copy
function with_copy(f, A, args...; kwargs...)
  B = copy(A)
  f(B, args...; kwargs...)
  return B
end

"""
    row_swap!(A, i, j)

Swap rows `i` and `j` in matrix `A`.

Returns a copy of the two rows in the new order.

# Example

```jldoctest
julia> A = reshape(1:12, 3, 4) |> collect
3×4 Matrix{Int64}:
 1  4  7  10
 2  5  8  11
 3  6  9  12

julia> row_swap!(A, 1, 3)
2×4 Matrix{Int64}:
 3  6  9  12
 1  4  7  10

julia> A
3×4 Matrix{Int64}:
 3  6  9  12
 2  5  8  11
 1  4  7  10
```
"""
row_swap!(A, i, j) = A[[i,j], :] = A[[j,i], :]

"""
    row_swap(A, i, j)

Create a matrix equal to `A` but with rows `i` and `j` swapped.

# Example

```jldoctest
julia> A = reshape(1:12, 3, 4) |> collect
3×4 Matrix{Int64}:
 1  4  7  10
 2  5  8  11
 3  6  9  12

julia> row_swap(A, 1, 3)
3×4 Matrix{Int64}:
 3  6  9  12
 2  5  8  11
 1  4  7  10
```
"""
row_swap(A, i, j) = with_copy(row_swap!, A, i, j)

"""
    row_mul!(A, i => λ)

Multiply row `i` of matrix `A` by a factor `λ`.

Returns a copy of the new row `i`.

# Example

```jldoctest
julia> A = reshape(1:12, 3, 4) |> collect
3×4 Matrix{Int64}:
 1  4  7  10
 2  5  8  11
 3  6  9  12

julia> row_mul!(A, 2 => 100)
4-element Vector{Int64}:
  200
  500
  800
 1100

julia> A
3×4 Matrix{Int64}:
   1    4    7    10
 200  500  800  1100
   3    6    9    12
```
"""
row_mul!(A, (i, factor)) = A[i, :] *= factor

"""
    row_mul(A, i => factor)

Create a matrix equal to `A` but with row `i` multiplied by `λ`.

# Example

```jldoctest
julia> A = reshape(1:12, 3, 4) |> collect
3×4 Matrix{Int64}:
 1  4  7  10
 2  5  8  11
 3  6  9  12

julia> row_mul(A, 2 => 100)
3×4 Matrix{Int64}:
   1    4    7    10
 200  500  800  1100
   3    6    9    12
```
"""
row_mul(A, (i, factor)) = with_copy(row_mul!, A, i => factor)

"""
    row_add!(A, i => λ => j)

Add row `i` multiplied by `λ` to row `j` of matrix `A`.

Returns a copy of the new row `j`.

# Example

```jldoctest
julia> A = [i//j for i in 1:3, j in 1:4]
3×4 Matrix{Rational{Int64}}:
 1//1  1//2  1//3  1//4
 2//1  1//1  2//3  1//2
 3//1  3//2  1//1  3//4

julia> row_add!(A, 1 => -3 => 3)
4-element Vector{Rational{Int64}}:
 0//1
 0//1
 0//1
 0//1

julia> A
3×4 Matrix{Rational{Int64}}:
 1//1  1//2  1//3  1//4
 2//1  1//1  2//3  1//2
 0//1  0//1  0//1  0//1
```
"""
row_add!(A, (i, (λ, j))) = A[j, :] += λ * A[i, :]

"""
    row_add(A, i => λ => j)

Create a matrix equal to `A` but with `λ` times row `i` added to row `j`.

# Example

```jldoctest
julia> A = [i//j for i in 1:3, j in 1:4]
3×4 Matrix{Rational{Int64}}:
 1//1  1//2  1//3  1//4
 2//1  1//1  2//3  1//2
 3//1  3//2  1//1  3//4

julia> row_add(A, 1 => -3 => 3)
3×4 Matrix{Rational{Int64}}:
 1//1  1//2  1//3  1//4
 2//1  1//1  2//3  1//2
 0//1  0//1  0//1  0//1
```
"""
row_add(A, (i, (λ, j))) = with_copy(row_add!, A, i => λ => j)

"""
    ref_pass!(A)

Perform one pass of the row-echelon-form algorithm.

Returns the column index of the pivot for this pass, or `nothing` if `A` is all
zeros.
"""
function ref_pass!(A)
  # Find first non-zero column
  col = findfirst(!iszero, collect(eachcol(A)))
  isnothing(col) && return nothing

  # Move pivot to first row if not there already
  if A[1, col] == 0
    row = findfirst(!iszero, @view A[:, col])
    row_swap!(A, 1, row)
  end

  # Use row operations to bring zeros under the pivot
  for row in 2:size(A, 1)
    if A[row, col] != 0
      row_add!(A, 1 => -A[row, col]/A[1, col] => row)
    end
  end

  return col
end

"""
    log(A, txt)

Write matrix `A` with title text `txt` to standard output.
"""
function log(A, txt)
  println(txt)
  println(repeat('-', length(txt)))
  show(stdout, MIME("text/plain"), A)
  println()
  println()
end

"""
    ref!(A; show_steps=false)

Put `A` in row echelon form.

Returns the `(row, column)` indices of the pivots.

# Example

```jldoctest
julia> A = Float64[mod(i+j,3) for i in 1:4, j in 1:6]
4×6 Matrix{Float64}:
 2.0  0.0  1.0  2.0  0.0  1.0
 0.0  1.0  2.0  0.0  1.0  2.0
 1.0  2.0  0.0  1.0  2.0  0.0
 2.0  0.0  1.0  2.0  0.0  1.0

julia> ref!(A)
3-element Vector{Tuple{Int64, Int64}}:
 (1, 1)
 (2, 2)
 (3, 3)

julia> A
4×6 Matrix{Float64}:
 2.0  0.0   1.0  2.0  0.0   1.0
 0.0  1.0   2.0  0.0  1.0   2.0
 0.0  0.0  -4.5  0.0  0.0  -4.5
 0.0  0.0   0.0  0.0  0.0   0.0
```
"""
function ref!(A; show_steps=false)
  show_steps && log(A, "Start")
  pivots = Tuple{Int,Int}[]

  # Iterate on the rows, each time doing one pass and "hiding" another row.
  for startrow = 1:size(A, 1)
    col = ref_pass!(@view A[startrow:end, :])

    # No pivot found -> nothing to do
    col == nothing && break

    # Store the position of the pivot
    push!(pivots, (startrow, col))
    show_steps && log(A, "Pivot $(pivots[end])")
  end
  return pivots
end

"""
    ref(A; show_steps=false)

Compute a row echelon form of `A`.

Returns the `(row, column)` indices of the pivots.

# Example

```jldoctest
julia> A = Float64[mod(i+j,3) for i in 1:4, j in 1:6]
4×6 Matrix{Float64}:
 2.0  0.0  1.0  2.0  0.0  1.0
 0.0  1.0  2.0  0.0  1.0  2.0
 1.0  2.0  0.0  1.0  2.0  0.0
 2.0  0.0  1.0  2.0  0.0  1.0

julia> ref(A)
4×6 Matrix{Float64}:
 2.0  0.0   1.0  2.0  0.0   1.0
 0.0  1.0   2.0  0.0  1.0   2.0
 0.0  0.0  -4.5  0.0  0.0  -4.5
 0.0  0.0   0.0  0.0  0.0   0.0
```
"""
ref(A; show_steps=false) = with_copy(ref!, A; show_steps)

"""
    rref!(A; show_steps)

Put `A` in the reduced row echelon form.

Returns the `(row, column)` indices of the pivots.

# Example

```jldoctest
julia> A = Float64[mod(i+j,3) for i in 1:4, j in 1:6]
4×6 Matrix{Float64}:
 2.0  0.0  1.0  2.0  0.0  1.0
 0.0  1.0  2.0  0.0  1.0  2.0
 1.0  2.0  0.0  1.0  2.0  0.0
 2.0  0.0  1.0  2.0  0.0  1.0

julia> rref!(A)
3-element Vector{Tuple{Int64, Int64}}:
 (1, 1)
 (2, 2)
 (3, 3)

julia> A
4×6 Matrix{Float64}:
  1.0   0.0  0.0   1.0   0.0  0.0
  0.0   1.0  0.0   0.0   1.0  0.0
 -0.0  -0.0  1.0  -0.0  -0.0  1.0
  0.0   0.0  0.0   0.0   0.0  0.0
```
"""
function rref!(A; show_steps=false)
  # First calculate non-reduced form
  pivots = ref!(A; show_steps)

  # Backward phase: iterate from last to first pivot
  for (row, col) in reverse(pivots)
    # Multiply row to have pivot = 1
    pivot = A[row, col]
    if pivot != 1
        row_mul!(A, row => inv(pivot))
        show_steps && log(A, "Reduce $((row,col)) $pivot -> 1")
    end

    # Ensure all zeros above the pivot
    if !all(iszero, A[1:row-1, col])
        for i in 1:row-1
            row_add!(A, row => -A[i, col] => i)
        end
        show_steps && log(A, "Reduce above $((row,col))")
    end
  end
  return pivots
end

"""
    rref(A; show_steps)

Compute the reduced row echelon form of `A`.

# Example

```jldoctest
julia> A = Float64[mod(i+j,3) for i in 1:4, j in 1:6]
4×6 Matrix{Float64}:
 2.0  0.0  1.0  2.0  0.0  1.0
 0.0  1.0  2.0  0.0  1.0  2.0
 1.0  2.0  0.0  1.0  2.0  0.0
 2.0  0.0  1.0  2.0  0.0  1.0

julia> rref(A)
4×6 Matrix{Float64}:
  1.0   0.0  0.0   1.0   0.0  0.0
  0.0   1.0  0.0   0.0   1.0  0.0
 -0.0  -0.0  1.0  -0.0  -0.0  1.0
  0.0   0.0  0.0   0.0   0.0  0.0
```
"""
rref(A; show_steps=false) = with_copy(rref!, A; show_steps)


text:
#= test/prange.jl
=#

@testset "prange" begin
  n = 128
  M = matrixlib(:fourier, rand(n), rand(n))
  opts = LRAOptions(maxdet_tol=0., sketch_randn_niter=1)

  @time for (t, s) in ((:none,                 :none ),
                       (:RandomGaussian,       :randn),
                       (:RandomSubset,         :sub  ),
                       (:SRFT,                 :srft ),
                       (:SparseRandomGaussian, :sprn ))
    opts.sketch = s
    for T in (Float32, Float64, ComplexF32, ComplexF64)
      println("  $t/$T")

      rtol = 5*eps(real(T))
      approx_rtol = 100*rtol
      opts.rtol = rtol

      A = convert(Array{T}, T <: Real ? real(M) : M)
      nrm = norm(A)

      Q = prange(:n, A, opts)
      @test norm(A - Q*(Q'*A)     ) < approx_rtol*nrm
      Q = prange(:c, A, opts)
      @test norm(A -      (A*Q)*Q') < approx_rtol*nrm
      Q = prange(:b, A, opts)
      @test norm(A - Q*(Q'*A*Q)*Q') < approx_rtol*nrm
    end
  end
end


text:
module Tests

using Base.Test

include("../src/Markdown.jl")
include("markdown.jl")

end


text:
"""
# KM3NeT Data Definitions v1.3.1-8-gae7dfb2
https://git.km3net.de/common/km3net-dataformat
"""

module W2list_gseagen
  const W2LIST_GSEAGEN_PS = 0
  const W2LIST_GSEAGEN_EG = 1
  const W2LIST_GSEAGEN_XSEC_MEAN = 2
  const W2LIST_GSEAGEN_COLUMN_DEPTH = 3
  const W2LIST_GSEAGEN_P_EARTH = 4
  const W2LIST_GSEAGEN_WATER_INT_LEN = 5
  const W2LIST_GSEAGEN_P_SCALE = 6
  const W2LIST_GSEAGEN_BX = 7
  const W2LIST_GSEAGEN_BY = 8
  const W2LIST_GSEAGEN_ICHAN = 9
  const W2LIST_GSEAGEN_CC = 10
  const W2LIST_GSEAGEN_DISTAMAX = 11
  const W2LIST_GSEAGEN_WATERXSEC = 12
  const W2LIST_GSEAGEN_XSEC = 13
  const W2LIST_GSEAGEN_DXSEC = 14
  const W2LIST_GSEAGEN_TARGETA = 15
  const W2LIST_GSEAGEN_TARGETZ = 16
  const W2LIST_GSEAGEN_VERINCAN = 17
  const W2LIST_GSEAGEN_LEPINCAN = 18
end


text:
################################################################################
#                                    αAl₂O₃                                    #
################################################################################
export αAl₂O₃

"""
This code creates a symbolic representation of a Cauchy Equation model
for the index of refraction of amorphous alumina (αAl₂O₃).
The model fit is due to:
    West, et al. "Low-loss integrated photonics for the blue and ultraviolet regime."
    APL Photonics 4.2 (2019): 026101. [doi:10.1063/1.5052502](https://doi.org/10.1063/1.5052502)
and has been validated from 0.25-.75 μm (free space) wavelength range for films grown by
atomic layer deposition (original work) and chemical vapor deposition.
The symbolic index model and its derivatives in turn are used to generate
numerical functions for the αAl₂O₃ index, group index and GVD as a function
of temperature and wavelength.
Variable units are λ in [um] and T in [deg C]
"""

p_n_αAl₂O₃ = (
    A = 1.602,              #                           [1]
    B = 0.01193,            #                           [μm²]
    C = -0.00036,           #                           [μm⁴]
    dn_dT = 6.1e-6,         # thermo-optic coefficient  [K⁻¹]
    T₀ = 20,                # reference temperature     [°C]
)

# notes & coeffs.from original working python code for reference:
### Cauchy Equation fit coefficients for Gavin's ALD alumina films ###
# Eqn. form:
# n = A + B / lam**2 + C / lam**4
# ng = A + 3 * B / lam**2 + 5 * C / lam**4
# A_alumina = 1.602
# B_alumina = 0.01193
# C_alumina = -0.00036

function make_αAl₂O₃(;p_n=p_n_αAl₂O₃)
	@variables λ, T
	n = n_sym_cauchy( λ ; p_n...) + p_n.dn_dT  *  ( T - p_n.T₀  )
	n² = n^2
	ε 	= diagm([n², n², n²])
	ng = ng_model(n,λ)
	gvd = gvd_model(n,λ)
	models = Dict([
		:n		=>	n,
		:ng		=>	ng,
		:gvd	=>	gvd,
		:ε 		=> 	diagm([n², n², n²]),
	])
	defaults =	Dict([
		:λ		=>		0.8,		# μm
		:T		=>		p_n.T₀,		# °C

	])
	Material(models, defaults, :αAl₂O₃,colorant"royalblue1")
end

################################################################

αAl₂O₃ = make_αAl₂O₃()


text:
using NeutralLandscapes
using Test

global anyerrors = false

tests = [
   "planar gradient" => "planar.jl"
]

for test in tests
   try
      include(test.second)
      println("\033[1m\033[32m✓\033[0m\t$(test.first)")
   catch e
      global anyerrors = true
      println("\033[1m\033[31m×\033[0m\t$(test.first)")
      println("\033[1m\033[38m→\033[0m\ttest/$(test.second)")
      showerror(stdout, e, backtrace())
      println()
      break
   end
end

if anyerrors
   throw("Tests failed")
end


text:
"""
    v̄ = p_lambert_velocity(r̄ₛ, r̄ₑ, Δt, μ; dir=1)

Computes the single-revolution orbit around a body with gravity parameter `μ`, passing through positions `r̄ₛ` and `r̄ₑ`,
and separated by time interval `Δt`.

Because there are two solutions, `dir` can be used to specifiy a clockwise or anti-clockwise orbit
relative to the coordinate system.
"""
function p_lambert_velocity(r̄ₛ::AbstractVector{<:Real}, r̄ₑ::AbstractVector{<:Real}, Δt, μ; dir=1, atol=1e-12, maxit=200, warntol=atol*1000)
    rₛ, rₑ = norm(r̄ₛ), norm(r̄ₑ)

    # true anomaly change for the transfer
    Δθ = atan(norm(cross(r̄ₛ, r̄ₑ)), dot(r̄ₛ, r̄ₑ))             # angle between start and end positions
    if wrap_angle(dir*angle_in_plane(r̄ₑ, r̄ₛ, MRP(1I))) > π
        Δθ = 2π - Δθ
    end
    # p iteration constants
    k = rₛ * rₑ * (1-cos(Δθ))
    L = rₛ + rₑ
    m = rₛ * rₑ * (1+cos(Δθ))

    # bounds
    pj = k / (L + √(2m))
    pjj = k / (L - √(2m))
    if Δθ > π
        pmin = 0.
        pmax = pjj
    else
        pmin = pj
        pmax = Inf
    end

    # Newton-p-iteration
    it = 0
    err = atol + 1
    p = (pj+pjj)/2
    pnext = p
    f = g = 0.
    while (err > atol) && (it < maxit)
        it = it + 1
        p = pnext
        a = m*k*p / ((2m-L^2)*(p^2) + 2k*L*p - k^2)
        a = a==0 ? 1e-100 : a   # to avoid parabolic case ?
        f = 1 - rₑ/p * (1 - cos(Δθ))
        g = rₛ * rₑ * sin(Δθ) / √(μ*p)
        df = √(μ/p)*tan(Δθ/2)*((1-cos(Δθ))/p - 1/rₛ - 1/rₑ);
        if a > 0    # Elliptical case
            sinΔE = -rₛ*rₑ*df/√(μ*a)
            cosΔE = 1 - rₛ/a * (1-f)
            ΔE = wrap_angle(atan(sinΔE, cosΔE))
            t = g + √(a^3/μ)*(ΔE-sinΔE)
            dtdp = -g/(2p) - 1.5a*(t-g)*(k^2 + (2m-L^2)*p^2)/(m*k*p^2) + √(a^3/μ)*(2k*sinΔE)/(p*(k-L*p))
        else        # Hyperbolic case
            ΔF = acosh(1 - rₛ/a * (1-f))
            t = g + √(-a^3/μ)*(sinh(ΔF)-ΔF)
            dtdp = -g/(2p) - 1.5a*(t-g)*(k^2 + (2m-L^2)*p^2)/(m*k*p^2) - √(-a^3/μ)*(2k*sinh(ΔF))/(p*(k-L*p))
        end
        err = abs(Δt - t)/Δt
        pnext = p + (Δt - t)/dtdp
        # If the next guess is outside of allowed bounds, use bisection
        if pnext < pmin
            pnext = (p + pmin)/2
        elseif pnext > pmax
            pnext = (p + pmax)/2
        end
    end
    if err == warntol
        @warn "Lambert solver failed to converge: err = $err"
    end
    v̄ₛ = (r̄ₑ - f*r̄ₛ)/g
    return v̄ₛ
end

"""
    torb = p_lambert_velocity(sorbₛ, eorb, stime, etime, μ, s_patch_position=@SVector(zeros(3)), e_patch_position=@SVector(zeros(3)); kwargs...)

Computes the single-revolution orbit around a body with gravity parameter `μ`, passing through the positions of orbits
`sorb` and `eorb` at times `stime` and `etime`, respectively.

Positions `s_patch_position` and `e_patch_position` can be provided to modify the start and end positions, respectively.
This can be used to minimize patching error across spheres of influence.
"""
function p_lambert(sorb::Orbit, eorb::Orbit, stime, etime, s_patch_position=@SVector(zeros(3)), e_patch_position=@SVector(zeros(3)); kwargs...)
    r̄ₛ, r̄ₑ = time_orbital_position(stime,sorb) + s_patch_position, time_orbital_position(etime,eorb) + e_patch_position
    d = sorb.i<π/2 ? 1 : -1     # choose short/long way based on clockwise sense of `sorb`
    v̄ₛ = p_lambert_velocity(r̄ₛ, r̄ₑ, etime-stime, sorb.primary.μ; dir=d, kwargs...)
    return Orbit(stime, r̄ₛ, v̄ₛ, sorb.primary)
end

function plane_change_p_lambert(sorb::Orbit, eorb::Orbit, stime, etime, s_patch_position=@SVector(zeros(3)), e_patch_position=@SVector(zeros(3)); kwargs...)
    r̄ₛ, r̄ₑ = time_orbital_position(stime,sorb) + s_patch_position, time_orbital_position(etime,eorb) + e_patch_position

    # change r̄ₑ so that it is in-plane with the starting orbit, at the same angle and distance.
    r̄ₑplane = inertial_to_perifocal_bases(r̄ₑ, sorb)
    r̄ₑplane = perifocal_to_inertial_bases(norm(r̄ₑ)*normalize(@SVector([r̄ₑplane[1], r̄ₑplane[2], 0])), sorb)

    d = sorb.i<π/2 ? 1 : -1     # choose short/long way based on clockwise sense of `sorb`
    v̄ₛ = p_lambert_velocity(r̄ₛ, r̄ₑplane, etime-stime, sorb.primary.μ; dir=d, kwargs...)
    torb1 = Orbit(stime, r̄ₛ, v̄ₛ, sorb.primary)

    # identify the true anomaly at plane change
    θs = time_to_true(stime,torb1)
    pcangle = min(wrap_angle(time_to_true(etime,torb1) - θs - π/2, 0.))
    θpc = θs + pcangle
    pctime = true_to_time(θpc, torb1, stime)

    # rotate velocity vector prior to plane change burn to get post-burn velocity
    r̄pc, v̄pc1 = state_vector(θpc, torb1)
    n̂1 = normalize(cross(r̄pc, v̄pc1))
    n̂2 = normalize(cross(r̄pc, r̄ₑ))
    Ri = align_vectors(n̂1, n̂2)
    v̄pc2 = Ri * v̄pc1
    torb2 = Orbit(pctime, r̄pc, v̄pc2, torb1.primary)

    return torb1, torb2
end


text:
using Test
using FlightMechanics.Models

# translate_forces_moments
pfm1 = PointForcesMoments([0, 0 ,0], [1, 0, 0], [0, 0, 0])
trans_pfm1 = translate_forces_moments(pfm1, [0, 0, 1])
exp_pfm1 = PointForcesMoments([0, 0, 1], [1, 0, 0], [0, -1, 0])

@test isapprox(trans_pfm1.point, exp_pfm1.point)
@test isapprox(trans_pfm1.forces, exp_pfm1.forces)
@test isapprox(trans_pfm1.moments, exp_pfm1.moments)

pfm2 = PointForcesMoments([0, 0 ,0], [1, 0, 0], [1, 0, 0])
trans_pfm2 = translate_forces_moments(pfm2, [0, 0, 1])
exp_pfm2 = PointForcesMoments([0, 0, 1], [1, 0, 0], [1, -1, 0])

@test isapprox(trans_pfm2.point, exp_pfm2.point)
@test isapprox(trans_pfm2.forces, exp_pfm2.forces)
@test isapprox(trans_pfm2.moments, exp_pfm2.moments)

# isapprox
@test isapprox(pfm1, trans_pfm1)
@test ≈(pfm1, trans_pfm1)
@test isapprox(pfm2, trans_pfm2)
@test ≈(pfm2, trans_pfm2)

# +
exp_pfm_sum = PointForcesMoments([0, 0, 0], [2, 0, 0], [0, 0, 0])
@test isapprox(pfm1 + pfm1, exp_pfm_sum)
# -
exp_pfm_sub = PointForcesMoments([0, 0, 0], [0, 0, 0], [0, 0, 0])
@test isapprox(pfm2 - pfm2, exp_pfm_sub)
# *
exp_pfm_mul = PointForcesMoments([0, 0, 1], [2, 0, 0], [2, -2, 0])
@test isapprox(2 * exp_pfm2, exp_pfm_mul)

# rotate
pfm3 = PointForcesMoments([0, 0, 0], [1, -1, 1], [-1, 1, 1])
rot_pfm3 = rotate(pfm3, 0, 0, 0)
exp_rot_pfm3 = PointForcesMoments([0, 0, 0], [1, -1, 1], [-1, 1, 1])
@test isapprox(rot_pfm3, exp_rot_pfm3)

pfm3 = PointForcesMoments([0, 0, 0], [0, 1, 0], [-1, 0, 0])
rot_pfm3 = rotate(pfm3, -π/2, 0, 0)
exp_rot_pfm3 = PointForcesMoments([0, 0, 0], [1, 0, 0], [0, 1, 0])
@test isapprox(rot_pfm3, exp_rot_pfm3, atol=1e-15)

pfm3 = PointForcesMoments([0, 0, 0], [0, 1, 0], [-1, 0, 0])
rot_pfm3 = rotate(pfm3, 0, π/2, 0)
exp_rot_pfm3 = PointForcesMoments([0, 0, 0], [0, 1, 0], [0, 0, 1])
@test isapprox(rot_pfm3, exp_rot_pfm3, atol=1e-15)


text:
using ArrayInterface
using EllipsisNotation
using Test

include("basic.jl")
include("more_generic.jl")
include("cartesian.jl")


text:
module Basins

#using DifferentialEquations
using DynamicalSystems
using LinearAlgebra
using ImageFiltering
using NearestNeighbors
using Combinatorics
using Roots
using LsqFit
using RegionTrees
using DataStructures

export basin_entropy, detect_wada_grid_method, detect_wada_merge_method
export box_counting_dim
export draw_basin!, basin_map, basin_general_ds
export basin_stability, uncertainty_exponent, static_estimate

export compute_saddle,compute_basin_precise
# Write your package code here.
include("basins_attraction.jl")
include("basin_entropy.jl")
include("uncertain_dim.jl")
include("wada_detection.jl")
include("straddle.jl")

end