text:
"""
    copras(decisionMat, weights, fs)
Apply COPRAS (COmplex PRoportional ASsesment) method for a given matrix, weights and, type of criteria.

# Arguments:
 - `decisionMat::DataFrame`: n × m matrix of objective values for n alternatives and m criteria
 - `weights::Array{Float64, 1}`: m-vector of weights that sum up to 1.0. If the sum of weights is not 1.0, it is automatically normalized.
 - `fs::Array{Function,1}`: m-vector of type of criteria. The benefit criteria shown with "maximum", and the cost criteria shown with "minimum".

 # Description
copras() applies the COPRAS method to rank n alternatives subject to m criteria and criteria type vector.

# Output
- `::COPRASResult`: COPRASResult object that holds multiple outputs including scores and best index.

# Examples
```julia-repl
julia> decmat = [2.50 240 57 45 1.10 0.333333;
       2.50 285 60 75 4.00 0.428571;
       4.50 320 100 65 7.50 1.111111;
       4.50 365 100 90 7.50 1.111111;
       5.00 400 100 90 11.00 1.111111;
       2.50 225 60 45 1.10 0.333333;
       2.50 270 57 60 4.00 0.428571;
       4.50 330 100 70 7.50 1.111111;
       4.50 365 100 80 7.50 1.111111;
       5.00 380 110 65 8.00 1.111111;
       2.50 285 65 80 4.00 0.400000;
       4.00 280 75 65 4.00 0.400000;
       4.50 365 102 95 7.50 1.111111;
       4.50 400 102 95 7.50 1.111111;
       6.00 450 110 95 11.00 1.176471;
       6.00 510 110 105 11.00 1.176471;
       6.00 330 140 110 18.50 1.395349;
       2.50 240 65 80 4.00 0.400000;
       4.00 280 75 75 4.00 0.400000;
       4.50 355 102 95 7.50 1.111111;
       4.50 385 102 90 7.50 1.111111;
       5.00 385 114 95 7.50 1.000000;
       6.00 400 110 90 11.00 1.000000;
       6.00 480 110 95 15.00 1.000000;
       6.00 440 140 100 18.50 1.200000;
       6.00 500 140 100 18.50 1.200000;
       5.00 450 125 100 15.00 1.714286;
       6.00 500 150 125 18.50 1.714286;
       6.00 515 180 140 22.00 2.307692;
       7.00 550 200 150 30.00 2.307692;
       6.00 500 180 140 15.00 2.307692;
       6.00 500 180 140 18.50 2.307692;
       6.00 500 180 140 22.00 2.307692;
       7.00 500 180 140 30.00 2.307692;
       7.00 500 200 140 37.00 2.307692;
       7.00 500 200 140 45.00 2.307692;
       7.00 500 200 140 55.00 2.307692;
       7.00 500 200 140 75.00 2.307692]
38×6 Array{Float64,2}:
 2.5  240.0   57.0   45.0   1.1  0.333333
 2.5  285.0   60.0   75.0   4.0  0.428571
 4.5  320.0  100.0   65.0   7.5  1.11111
 4.5  365.0  100.0   90.0   7.5  1.11111
 5.0  400.0  100.0   90.0  11.0  1.11111
 2.5  225.0   60.0   45.0   1.1  0.333333
 2.5  270.0   57.0   60.0   4.0  0.428571
 4.5  330.0  100.0   70.0   7.5  1.11111
 4.5  365.0  100.0   80.0   7.5  1.11111
 5.0  380.0  110.0   65.0   8.0  1.11111
 2.5  285.0   65.0   80.0   4.0  0.4
 ⋮                               ⋮
 6.0  500.0  150.0  125.0  18.5  1.71429
 6.0  515.0  180.0  140.0  22.0  2.30769
 7.0  550.0  200.0  150.0  30.0  2.30769
 6.0  500.0  180.0  140.0  15.0  2.30769
 6.0  500.0  180.0  140.0  18.5  2.30769
 6.0  500.0  180.0  140.0  22.0  2.30769
 7.0  500.0  180.0  140.0  30.0  2.30769
 7.0  500.0  200.0  140.0  37.0  2.30769
 7.0  500.0  200.0  140.0  45.0  2.30769
 7.0  500.0  200.0  140.0  55.0  2.30769
 7.0  500.0  200.0  140.0  75.0  2.30769

julia> df = makeDecisionMatrix(decmat)
38×6 DataFrame
 Row │ Crt1     Crt2     Crt3     Crt4     Crt5     Crt6
     │ Float64  Float64  Float64  Float64  Float64  Float64
─────┼───────────────────────────────────────────────────────
   1 │     2.5    240.0     57.0     45.0      1.1  0.333333
   2 │     2.5    285.0     60.0     75.0      4.0  0.428571
   3 │     4.5    320.0    100.0     65.0      7.5  1.11111
   4 │     4.5    365.0    100.0     90.0      7.5  1.11111
   5 │     5.0    400.0    100.0     90.0     11.0  1.11111
   6 │     2.5    225.0     60.0     45.0      1.1  0.333333
   7 │     2.5    270.0     57.0     60.0      4.0  0.428571
   8 │     4.5    330.0    100.0     70.0      7.5  1.11111
   9 │     4.5    365.0    100.0     80.0      7.5  1.11111
  ⋮  │    ⋮        ⋮        ⋮        ⋮        ⋮        ⋮
  30 │     7.0    550.0    200.0    150.0     30.0  2.30769
  31 │     6.0    500.0    180.0    140.0     15.0  2.30769
  32 │     6.0    500.0    180.0    140.0     18.5  2.30769
  33 │     6.0    500.0    180.0    140.0     22.0  2.30769
  34 │     7.0    500.0    180.0    140.0     30.0  2.30769
  35 │     7.0    500.0    200.0    140.0     37.0  2.30769
  36 │     7.0    500.0    200.0    140.0     45.0  2.30769
  37 │     7.0    500.0    200.0    140.0     55.0  2.30769
  38 │     7.0    500.0    200.0    140.0     75.0  2.30769
                                              20 rows omitted

julia> weights = [0.1667, 0.1667, 0.1667, 0.1667, 0.1667, 0.1667];

julia> fns = [maximum, maximum, maximum, maximum, maximum, minimum];

julia> result = copras(df, w, fs);

julia> result.bestIndex
38

julia> result.scores
38-element Array{Float64,1}:
 0.021727395411605937
 0.019814414550092637
 0.01789214190869233
 0.01624057709923278
 0.01507318798582843
 0.021837811311495522
 0.020735423365838293
 0.01754833367014725
 0.016658288043259514
 0.016258710042371068
 0.019388734250223458
 ⋮
 0.01064826476628102
 0.01000964101170639
 0.007056714858865703
 0.010627520153194962
 0.010419795157349462
 0.010212070161503961
 0.00876261130160425
 0.007733739643860005
 0.007258939653356005
 0.006665439665226005
 0.005478439688966004
```
# References
Kaklauskas, A., Zavadskas, E. K., Raslanas, S., Ginevicius, R., Komka, A., & Malinauskas, P. (2006). Selection of low-e windows in retrofit of public buildings by applying multiple criteria method COPRAS: A Lithuanian case. Energy and buildings, 38(5), 454-462.
Özdağoğlu, A. (2013). İmalat işletmeleri için eksantrik pres alternatiflerinin COPRAS yöntemi ile karşılaştırılması. Gümüşhane Üniversitesi Sosyal Bilimler Enstitüsü Elektronik Dergisi, 4(8), 1-22.
Yıldırım, B. F., Timor, M. (2019). "Bulanık ve Gri COPRAS Yöntemleri Kullanılarak Tedarikçi Seçim Modeli Geliştirilmesi". Optimum Ekonomi ve Yönetim Bilimleri Dergisi, 6 (2), 283-310.
"""
function copras(decisionMat::DataFrame, weights::Array{Float64,1}, fns::Array{Function,1})::COPRASResult

    #mat = convert(Matrix, decisionMat)
    mat = Matrix(decisionMat)

    nrows, ncols = size(mat)
    w = unitize(weights)
    normalizedMat = copy(mat)
    for col in 1:ncols
        for row in 1:nrows
            normalizedMat[row, col] = w[col] .* (mat[row, col] ./ sum(mat[:, col]))
        end
    end

    sPlus = zeros(Float64, nrows)
    sMinus = zeros(Float64, nrows)

    for row in 1:nrows
        for col in 1:ncols
            if fns[col] == maximum
                @inbounds sPlus[row] = sPlus[row] + normalizedMat[row, col]
            elseif fns[col] == minimum
                @inbounds sMinus[row] = sMinus[row] + normalizedMat[row, col]
            end
        end
    end

    Q = zeros(Float64, nrows)
    Z = sum(1 ./ sMinus)

    for row in 1:nrows
        Q[row] = sPlus[row] + (sum(sMinus) / (sMinus[row] * Z))
    end
    scores = zeros(Float64, nrows)
    scores = Q ./ maximum(Q)

    rankings = sortperm(scores)
    bestIndex = rankings |> last

    result = COPRASResult(
        decisionMat,
        w,
        scores,
        rankings,
        bestIndex
    )

    return result
end

"""
    copras(setting)
Apply COPRAS (COmplex PRoportional ASsesment) method for a given matrix, weights and, type of criteria.

# Arguments:
 - `setting::MCDMSetting`: MCDMSetting object.

 # Description
copras() applies the COPRAS method to rank n alternatives subject to m criteria and criteria type vector.

# Output
- `::COPRASResult`: COPRASResult object that holds multiple outputs including scores and best index.
"""
function copras(setting::MCDMSetting)::COPRASResult
    copras(
        setting.df,
        setting.weights,
        setting.fns
    )
end

max_stars_repo_name: drvinceknight/JMcDM

max_issues_repo_name: drvinceknight/JMcDM

max_forks_repo_name: drvinceknight/JMcDM


text:
function testSolveSession(client, context, variableLabels; maxSeconds=180)
    # allVariableLabels = ls(client, context, variableLabels)

    # do the solve
    resultId = solveSession(client,context)

    # Wait for them to be done before proceeding.
    waitForCompletion(client, Task[resultId;]; expectedStatuses=["Complete"], maxSeconds)

    # Get PPE's are there for the connected variables.
    # TODO - complete the factor graph.
    for v in variableLabels # getVariables(client, context, detail=SUMMARY)
        # First solve key (only) is the default result
        @test fetch( getVariable(client, context, v) )["ppes"][1]["solveKey"] == "default"
    end
end

function testVizHelpersApp(context)
    @show GraphVizApp(context)
    @show MapVizApp(context)

    nothing
end

function runSolveTests(client, context)
    @testset "solve-testset" begin
        @info "Running solve-testset"

        # TODO: Test the rest of the variables.
        variableLabels = ["x0", "x1"]

        testSolveSession(client, context, variableLabels)
    end

    @testset "appviz-testset" begin
        # generate viz help objects from current context
        testVizHelpersApp(context)
    end
end

ext: jl


text:
using BenchmarkTools
using AugmentedGaussianProcesses
using Distances, LinearAlgebra, DelimitedFiles
using MLDataUtils: splitobs
using Random
const AGP = AugmentedGaussianProcesses

compat = Dict{String,Dict{String,Vector{String}}}()
likelihoodnames = ["Gaussian","AugmentedStudentT","AugmentedLogistic","BayesianSVM","AugmentedLogisticSoftMax"]
inferencenames = ["AnalyticInference","StochasticAnalyticInference"]
modelnames = ["GP","VGP","SVGP"]
funcs = ["init","elbo","computematrices","updatevariational","updatehyperparam","predic","predicproba"]
compat["GP"] = Dict{String,Vector{String}}("Gaussian"=>["AnalyticInference"])
compat["VGP"] = Dict{String,Vector{String}}()
compat["VGP"]["AugmentedStudentT"] = ["AnalyticInference"]
compat["VGP"]["AugmentedLogistic"] = ["AnalyticInference"]
compat["VGP"]["BayesianSVM"] = ["AnalyticInference"]
compat["VGP"]["AugmentedLogisticSoftMax"] = ["AnalyticInference"]
compat["SVGP"] = Dict{String,Vector{String}}()
compat["SVGP"]["Gaussian"] = ["AnalyticInference","StochasticAnalyticInference"]
compat["SVGP"]["AugmentedStudentT"] = ["AnalyticInference","StochasticAnalyticInference"]
compat["SVGP"]["AugmentedLogistic"] = ["AnalyticInference","StochasticAnalyticInference"]
compat["SVGP"]["BayesianSVM"] = ["AnalyticInference","StochasticAnalyticInference"]
compat["SVGP"]["AugmentedLogisticSoftMax"] = ["AnalyticInference","StochasticAnalyticInference"]

# const SUITE = BenchmarkGroup(["Models"])
Random.seed!(1234)
nData = 200; nDim = 3;
X = rand(nData,nDim)
y = Dict("Gaussian"=>norm.(eachrow(X)),"AugmentedStudentT"=>norm.(eachrow(X)),"BayesianSVM"=>sign.(norm.(eachrow(X)).-1.0),"AugmentedLogistic"=>sign.(norm.(eachrow(X)).-1.0),"AugmentedLogisticSoftMax"=>floor.(norm.(eachrow(X.*2))))
n_ind = 50; batchsize = 50; ν = 10.0
convertl(lname::String) = lname*(lname != "BayesianSVM" ? "Likelihood" : "")*"("*(lname == "AugmentedStudentT" ? "ν" : "")*")"
converti(iname::String) = iname*"("*(iname[1:10] == "Stochastic" ? "batchsize" : "")*")"
add_ind(mname::String) = mname == "SVGP" ? ",n_ind" : ""
kernel = RBFKernel([2.0],variance=1.0,dim=nDim)
models = Dict{String,Dict{String,Dict{String,AbstractGP}}}()
SUITE["Models"] = BenchmarkGroup(modelnames)
for model in String.(keys(compat))
    SUITE["Models"][model] = BenchmarkGroup(String.(keys(compat[model])))
    models[model] = Dict{String,Dict{String,AbstractGP}}()
    for likelihood in String.(keys(compat[model]))
        SUITE["Models"][model][likelihood] = BenchmarkGroup(compat[model][likelihood])
        models[model][likelihood] = Dict{String,AbstractGP}()
        for i in compat[model][likelihood]
            SUITE["Models"][model][likelihood][i] = BenchmarkGroup(funcs)
            if model == "GP"
                # println(Meta.parse(model*"(\$X,\$y[\"$likelihood\"],\$kernel,Autotuning=true,atfrequency=1)"))
                models[model][likelihood][i] = eval(Meta.parse(model*"(X,y[\"$likelihood\"],kernel,Autotuning=true,atfrequency=1)"))
                SUITE["Models"][model][likelihood][i]["init"] = eval(Meta.parse("@benchmarkable $model(\$X,y_train,\$kernel,Autotuning=true,atfrequency=1) setup=(y_train = \$y[\"$likelihood\"])"))
            else
                # println(Meta.parse(model*"(X,y[\"$likelihood\"],kernel,$(convertl(likelihood)) ,$(converti(i))$(add_ind(model)),Autotuning=true,atfrequency=1)"))
                models[model][likelihood][i] = eval(Meta.parse(model*"(X,y[\"$likelihood\"],kernel,$(convertl(likelihood)) ,$(converti(i))$(add_ind(model)),Autotuning=true,atfrequency=1)"))
                SUITE["Models"][model][likelihood][i]["init"] = eval(Meta.parse("@benchmarkable $model(\$X,y_train,\$kernel,$(convertl(likelihood)),$(converti(i)) $(add_ind(model)),Autotuning=true,atfrequency=1) setup=(y_train = \$y[\"$likelihood\"])"))
            end
            SUITE["Models"][model][likelihood][i]["elbo"] = @benchmarkable ELBO(gpmodel) setup=(gpmodel=deepcopy($(models[model][likelihood][i])))
            SUITE["Models"][model][likelihood][i]["computematrices"] = @benchmarkable AGP.computeMatrices!(gpmodel) setup=(gpmodel=deepcopy($(models[model][likelihood][i])))
            SUITE["Models"][model][likelihood][i]["updatevariational"] = @benchmarkable AGP.variational_updates!(gpmodel) setup=(gpmodel=deepcopy($(models[model][likelihood][i])))
            SUITE["Models"][model][likelihood][i]["updatehyperparam"] = @benchmarkable AGP.update_hyperparameters!(gpmodel) setup=(gpmodel=deepcopy($(models[model][likelihood][i])))
            SUITE["Models"][model][likelihood][i]["predic"] = @benchmarkable predict_y(gpmodel,$X) setup=(gpmodel=deepcopy($(models[model][likelihood][i])))
            SUITE["Models"][model][likelihood][i]["predicproba"] = @benchmarkable proba_y(gpmodel,$X) setup=(gpmodel=deepcopy($(models[model][likelihood][i])))
        end
    end
end
#
# if isfile(paramfile)
#     loadparams!(suite,BenchmarkTools.load(paramfile))
# else
#     println("Tuning parameters")
#     tune!(suite,verbose=true)
#     BenchmarkTools.save(paramfile,params(suite))
# end
# println("Running benchmarks")
# results = run(suite,verbose=true,seconds=30)
# save_target = "results/multiclass_"*("$(now())"[1:10])
# i = 1
# while isfile(save_target*"_$(i).json")
#     global i += 1
# end
# BenchmarkTools.save(save_target*"_$(i).json",results)


text:
using Test
using BenchmarkTools
using LinearAlgebra
include("../src/Library.jl")

# Quadratic basis
@testset "quadratic" begin
    A = [1 2 3;4 5 6;7 8 10]
    B = quadratic(A,3,3)
    @test cond(A) < 1e6
    @test cond(B) < 1e6
end

# Cubic basis
@testset "cubic" begin
    A = [1 2;3 4]
    B = cubic(A,2,2)
    @test cond(B) <1e6

    A = [1 2 3;4 5 6;7 8 10]
    B = cubic(A,3,3)
    @test cond(B) < 1e6
end

# Test for a small unit matrix
@testset "ones" begin
    A = ones(3,3)
    o = PolynomialBasis()
    B = basis(A,o)
    @test cond(B) == Inf
end

# Test for an identity matrix
@testset "identity" begin
    A = [1 0 0;0 1 0;0 0 1]
    o = PolynomialBasis()
    B = basis(A,o)
    @test cond(B) < 1e6
end

# Condition number of a non square matrix
@testset "non square matrix" begin
    A = [1 0 0 0;0 1 0 0;0 0 0 1]
    o = PolynomialBasis()
    B = basis(A,o)
    @test cond(B) < 1e6
end

# Condition number of random non square matrix
@testset "random non square matrix" begin
    A = 1.0 .+ randn(3,2)
    o = PolynomialBasis()
    B = basis(A,o)
    @test cond(B) < 1e6
end

@testset "basis 2 -  correctness" begin
    A = [1 2;3 4]
    o = PolynomialBasis()
    B = basis(A,o)
    @test B[1,:] == [1,1]
    @test B[2,:] == [1,2]
    @test B[3,:] == [3,4]
end

@testset "basis 3 - correctness" begin
    A = [1 2 3;4 5 6;7 8 10]
    o = PolynomialBasis()
    B = basis(A,o)
    @test size(B,1) == 1 + 3 + 3*2 + 3*3
    @test size(B,2) == 3
    @test B[1,:] == [1,1,1]
    @test B[2,:] == [1,2,3]
    @test B[3,:] == [4,5,6]
    @test B[end,:] == [343,512,1000]
end

@testset "trig" begin
    A = [1 2;3 4]
    B = trig(A,2,2)
    B[1,:] ≈ [sin(1),sin(2)]
    B[2,:] ≈ [sin(3),sin(4)]
    B[3,:] ≈ [cos(1),cos(2)]
    B[4,:] ≈ [cos(3),cos(4)]
end

@testset "trig" begin
    A = [1 2 3; 4 5 6; 7 8 9]
    o = TrigBasis()
    B = basis(A,o)
    @test size(B,1) == 9
    @test size(B,2) == 3
    @test cond(B) < 1e6
end

@testset "PolyTrig basis" begin
    A = [1 2 3; 4 5 6; 7 8 9]
    o = PolyTrigBasis()
    B = basis(A,o)
    @test cond(B) < 1e6
end


text:
using Pkg

current_directory = @__DIR__

Pkg.activate(joinpath(current_directory, ".."))

using FV3Kernels, NCDatasets, Test

input_file_name = joinpath(current_directory, "../data/inputs/c_sw_12x24.nc")

input_file = NCDataset(input_file_name)

state = State(input_file)

@test state.isd == -2

# @show state


text:
@testset "Controller: Implicit Dynamics" begin
	# Quadruped reference trajectory
	s = ContactImplicitMPC.get_simulation("quadruped", "flat_2D_lc", "flat")
	model = s.model
	env = s.env

	ref_traj = deepcopy(ContactImplicitMPC.get_trajectory(model, env,
		joinpath(module_dir(), "src/dynamics/quadruped/gaits/gait2.jld2"),
		load_type = :split_traj_alt))

	ip_opts = InteriorPointOptions(
				κ_tol = 2.0 * 1.0e-4,
				r_tol = 1.0e-8,
				diff_sol = true,
				solver = :empty_solver)

	im_traj = ContactImplicitMPC.ImplicitTrajectory(ref_traj, s,
		κ = 1.0e-4,
		mode = :configuration,
		opts=ip_opts)

	ContactImplicitMPC.implicit_dynamics!(im_traj, ref_traj)

	@test all([norm(im_traj.dq2[t], Inf) < 1.0e-2 for t = 1:ref_traj.H-1])
end


text:
module HTbin
using Glob
using Parameters
using Dates
using Geodesy

using RemoteFiles

include("common.jl")
include("so_files.jl")
include("executive_files.jl")
include("system_bin.jl")

greet() = print("Hello World!")

end # module


text:
export libration_points, pos_to_inertial, state_to_inertial
export R3BPSystemProperties

#---------------------------#
# RESTRICTED 3-BODY SYSTEMS #
#---------------------------#

default_reference_frame(::Abstract_R3BPModel) = SynodicFrame()

Base.show(io::IO, x::Abstract_R3BPModel) = print(io, "$(nameof(typeof(x)))$(x.props)")

#-------------------#
# SYSTEM PROPERTIES #
#-------------------#
struct R3BPSystemProperties{L <: Unitful.Length, V <: Unitful.Velocity, T <: Unitful.Time}
    b1 :: Symbol        # Identifier of body 1
    b2 :: Symbol        # Identifier of body 1
    μ  :: Float64       # Mass ratio
    L  :: L             # Distance between centers
    V  :: V             # Orbital velocity of larger body
    T  :: T             # Time unit: orbital period of smaller body (per radian)
    e  :: Float64       # Eccentricity
    R1 :: SVector{3,L}  # Radius of central body
    R2 :: SVector{3,L}  # Radius of smaller body
end

@memoize function R3BPSystemProperties(a::Symbol, b::Symbol; kwargs...)
    # Extract the static data for this system
    a, b = @. Symbol(lowercase(String((a, b))))
    data = get(_3BPSystemPropertyData, (a, b),
           get(_3BPSystemPropertyData, (b, a), nothing))
    isnothing(data) && error("No system data available for $a, $b")

    # Get user-provided values
    data = merge(data, kwargs)

    # Run-time compute remaining values
    μ = get(data, :μ, SpiceUtils.mass_fraction(a, b))
    R1 = get(data, :R1, SVector{3}(bodvrd(String(a), "RADII")u"km"))
    R2 = get(data, :R2, SVector{3}(bodvrd(String(b), "RADII")u"km"))
    # Ensure the Length attribute matches the Radius type
    L = convert(eltype(R1), data[:L])

    R3BPSystemProperties(a, b, μ, L, data[:V], data[:T], data[:e], R1, R2)
end
R3BPSystemProperties(system::Abstract_DynamicalModel) = R3BPSystemProperties(primary_body(system), secondary_body(system))

# TODO: Compute these constants at run-time?
const _3BPSystemPropertyData = Dict(
    (:mars, :phobos) => Dict(
        :L => 9.376e3u"km",
        :V => 2.138u"km/s",
        :T => 0.31891023u"d" / (2π),
        :e => 0.0151
    ),
    (:earth, :moon) => Dict(
        :L => 384399u"km",
        :V => 1.022u"km/s",
        :T => 27.321661u"d" / (2π),
        :e => 0.0549
    ),
    (:sun, :mars) => Dict(
        :L => 2.279497905330276e8u"km",
        :V => 24.128831378998047u"km/s",
        :T => 686.99u"d" / (2π),
        :e => 0.0935643512
    ),
    (:sun, :earth) => Dict(
        # Values from Wikipedia
        :L => 149598023u"km",
        :V => 29.78u"km/s",
        :T => 365.256363004u"d" / (2π),
        :e => 0.0167086
    ),
    (:jupiter, :europa) => Dict(
        # Values from Wikipedia, [Pellegrini.2006]
        :L => 670900u"km",
        :V => 13.740u"km/s",
        :T => 3.551181u"d" / (2π),
        :e => 0.009,
        :μ => 2.528009215182033e-5
    )
)

Base.show(io::IO, x::R3BPSystemProperties) = print(io, (x.b1, x.b2))
ModelingToolkit.parameters(model::Abstract_R3BPModel) = [getfield(model.props, i.name) for i in ModelingToolkit.parameters(model.ode.ode_system)]

#---------#
# METHODS #
#---------#
@memoize function libration_points(props::R3BPSystemProperties)
    # Begin with some guesses for the libration points
    μ = props.μ
    L_guess = [
        [(1 - μ)^2, 0., 0.],         # L1
        [(1 - μ)^-1, 0., 0.],        # L2
        [-1., 0., 0.],               # L3
        [(1 - μ)^2, (1 - μ^2), 0.],  # L4
        [(1 - μ)^2, -(1 - μ^2), 0.]  # L5
    ]

    # Solve for gradient(libration point guess) == 0 using forward AD
    libration_points = nlsolve.(
        (pos) -> ForwardDiff.gradient(p -> centrifugal_potential(μ, p), pos),
        L_guess, autodiff=:forward)

    # Return the zero solution (i.e. the position) for each libration point
    map(l -> l.zero, libration_points)
end
libration_points(model::Abstract_R3BPModel) = libration_points(model.props)

function distance_vectors(μ, pos)
    x, y, z = pos
    ((μ + x)^2     + y^2 + z^2)^(1/2), # Distance to the central body
    ((1 - μ - x)^2 + y^2 + z^2)^(1/2)  # Distance to the smaller body
end

@doc "Compute the Hill sphere radius for two bodies"
function hill_sphere_radius(system::Abstract_R3BPModel)
    mass_1 = bodvrd(system.props.b1, "GM")[1]
    mass_2 = bodvrd(system.props.b2, "GM")[1]
    return  system.props.L * (1 - system.props.e) * cbrt(mass_2 / (3 * mass_1))
end

@doc """ Returns the position rotation matrix from a rotating 3BP frame to an inertial frame. """
pos_to_inertial(t) = [
    # [Pavlak 2013, eq. 2.47]
    cos(t)  -sin(t)  0;
    sin(t)   cos(t)  0;
    0        0       1
]
pos_from_inertial(t) = inv(pos_to_inertial(t))

@doc """ Returns the state rotation matrix from a rotating 3BP frame to an inertial frame. """
state_to_inertial(t) = [
    # [Pavlak 2013, eq. 2.50]
     cos(t)  -sin(t)  0  0        0       0;
     sin(t)   cos(t)  0  0        0       0;
     0        0       1  0        0       0;
    -sin(t)  -cos(t)  0  cos(t)  -sin(t)  0;
     cos(t)  -sin(t)  0  sin(t)   cos(t)  0;
     0        0       0  0        0       1
]
state_from_inertial(t) = inv(state_to_inertial(t))

primary_body(m::Abstract_R3BPModel) = m.props.b1
secondary_body(m::Abstract_R3BPModel) = m.props.b2

function collision(system::Abstract_R3BPModel, body, dist=bodvrd(String(body), "RADII")[1]; radii=1.)
    circ_props = R3BPSystemProperties(primary_body(system), secondary_body(system))
    diam = (radii * dist) ./ ustrip(u"km", circ_props.L)
    return ContinuousCallback((integrator) -> terminate!(integrator, :Crashed); interp_points=10) do u, t, integrator
        check_distance(u, t, system, body, diam)
    end
end

function check_distance(u, t, system::Abstract_R3BPModel, body::Symbol, distance)
    body_pos = (body == primary_body(system)) ? SVector{3}(-system.props.μ, 0., 0.) : SVector{3}(1 - system.props.μ, 0., 0.)
    return norm(@view(u[1:3]) - body_pos) - distance
end

function convert_to_frame(state::State{<:Abstract_R3BPModel,<:Abstract_ReferenceFrame}, frame::InertialFrame)
    to_inertial = state_to_inertial(state.tspan[1])
    converted_u0 = to_inertial * state.u0
    new_state = State(state.model, InertialFrame(), converted_u0, state.tspan)
    return new_state
end

#----------#
# INCLUDES #
#----------#
include("er3bp.jl")
include("cr3bp.jl")  # Depends on ER3BP


text:
push!(LOAD_PATH, "scripts/julia_modules/ArgParse2.jl/src")
using ArgParse2
using DataFrames
using DataFramesMeta
using CSV
using Statistics

function main()
    parser = ArgumentParser(prog = "",
                        description = "Evaluates ckstering accuracy based on expected gebus distribution",
                        epilog = "sizes..ajusted")
    #add_argument!(parser, "--read-counts","-c", help = "directory with Reads counts based on aligned clustering",type = String)
    add_argument!(parser, "--reference-values","-r", help = "Reference values",type = String)
    add_argument!(parser,"--read-counts","-c", metavar="COUNTS", nargs = "+", required = true)
    args = parse_args(parser)
    ref_data =  DataFrame!(CSV.File(args.reference_values))
    println("File CorvsMaximumcluster CorvsTotal CorvsNumberOfClusters NumberOfClusters DetectedReqGenus")
    ct = 0
    cor1 = 0
    cor2 = 0
    cor3 = 0
    for cntf in args.read_counts
        ct += 1
        counts = DataFrame!(CSV.File(cntf))
        counts = @transform( counts, Genus2 = getindex.(split.(:Genus,"/"),1) )
        counts = @transform( counts, Genus2 = getindex.(split.(:Genus2,"-"),1) )
        counts = @transform( counts, Genus2 = getindex.(split.(:Genus2),1) )
        counts1 = @select(counts, :Genus2, :Maximum_abundance)
        counts2 = @select(counts, :Genus2, :Total_abundance)
        counts3 = @select(counts, :Genus2, :Number_of_sequences)
        dd1 =join(counts1,ref_data, on = :Genus2 => :Genus, kind = :inner)
        dd2 =join(counts2,ref_data, on = :Genus2 => :Genus, kind = :inner)
        dd3 =join(counts3,ref_data, on = :Genus2 => :Genus, kind = :inner)
        if (nrow(dd1) > 1)
            cor1 = cor( dd1[:,:Maximum_abundance]  ,  dd1[:,Symbol("16S")] )
            cor2 = cor( dd2[:,:Total_abundance]  ,  dd2[:,Symbol("16S")] )
            cor3 = cor( dd3[:,:Number_of_sequences]  ,  dd3[:,Symbol("16S")] )
        end
        println(cntf," ",cor1," ",cor2," ",cor3," ", sum(dd3[:,:Number_of_sequences])," ",nrow(dd1))
    end
end

main()


text: # Placeholder empty entry point script file

max_stars_repo_name: tlienart/SIIP-Tutorial

max_issues_repo_name: tlienart/SIIP-Tutorial

max_forks_repo_name: tlienart/SIIP-Tutorial


text:
using DataFrames, MANOVABNPTest, StatsBase, LinearAlgebra, StaticArrays, Random
rng = MersenneTwister(1)

function test_sample_01(D::Int, h::Int, l::Int, H0::Int)
    N = 4 * [50, 150, 300]
    K = [
        0.625 1.498 1.162;
        0.819 1.806 1.359;
        1.051 1.994 1.652
    ]
    γ = γvector(4, H0)
    Σ = cholesky(collect(I(D)))
    x = 1 .+ (0:N[h]-1) .% 4
    y = rand(N[h], D) * Σ.U
    for i in 1:N[h]
        x[i] == 2 && (y[i, :] .+= K[l, 1] * γ[2])
        x[i] == 3 && (y[i, :] .*= K[l, 2] ^ γ[3])
        x[i] == 4 && (y[i, :] .+= K[l, 3] * γ[4] * rand([-1, 1]))
    end
    ȳ, S = mean_and_cov(y)
    y = (y .- ȳ) / cholesky(S).U
    y = [SVector{D}(y[i, :]) for i ∈ 1:N[h]]
    return y, x
end

function test_sample_02(D::Int, h::Int, l::Int, H0::Int)
    N = 4 * [50, 150, 300]
    K = [
        0.625 1.498 1.162;
        0.819 1.806 1.359;
        1.051 1.994 1.652
    ]
    γ = γvector(4, H0)
    Σ = cholesky(collect(I(D)))
    x = 1 .+ (0:N[h]-1) .% 4
    y = rand(N[h], D) * Σ.U
    for i in 1:N[h]
        x[i] == 2 && (y[i, :] .+= K[l, 1] * γ[2])
        x[i] == 3 && (y[i, :] .*= K[l, 2] ^ γ[3])
        x[i] == 4 && (y[i, :] .+= K[l, 3] * γ[4] * rand([-1, 1]))
    end
    return y, x
end

# Exp 1
D = 2
m = MANOVABNPTest.Model(D = D)
y, x = test_sample_01(D, 1, 1, 2)
pγ = MANOVABNPTest.fit(m, y, x; iter = 4000, refresh_rate = 50, rng = rng)
println(pγ)

# Exp 2
D = 2
m = MANOVABNPTest.Model(D = D)
y, x = test_sample_01(D, 1, 1, 3)
pγ = MANOVABNPTest.fit(m, y, x; iter = 4000, refresh_rate = 50, rng = rng)
println(pγ)

# Exp 3
D = 10
m = MANOVABNPTest.Model(D = D)
y, x = test_sample_01(D, 1, 1, 2)
pγ = MANOVABNPTest.fit(m, y, x; iter = 4000, refresh_rate = 50, rng = rng)
println(pγ)

# Exp 4
D = 10
m = MANOVABNPTest.Model(D = D)
y, x = test_sample_01(D, 1, 1, 3)
pγ = MANOVABNPTest.fit(m, y, x; iter = 4000, refresh_rate = 50, rng = rng)
println(pγ)

# Exp 5
D = 10
m = MANOVABNPTest.Model(D = D)
y, x = test_sample_01(D, 1, 2, 2)
pγ = MANOVABNPTest.fit(m, y, x; iter = 4000, refresh_rate = 50, rng = rng)
println(pγ)

# Exp 6
D = 10
m = MANOVABNPTest.Model(D = D)
y, x = test_sample_01(D, 1, 2, 3)
pγ = MANOVABNPTest.fit(m, y, x; iter = 4000, refresh_rate = 50, rng = rng)
println(pγ)

# Exp 7
D = 10
m = MANOVABNPTest.Model(D = D)
y, x = test_sample_01(D, 1, 3, 2)
pγ = MANOVABNPTest.fit(m, y, x; iter = 4000, refresh_rate = 50, rng = rng)
println(pγ)

# Exp 8
D = 10
m = MANOVABNPTest.Model(D = D)
y, x = test_sample_01(D, 1, 3, 3)
pγ = MANOVABNPTest.fit(m, y, x; iter = 4000, refresh_rate = 50, rng = rng)
println(pγ)

# grid = LinRange(-3, 3, 10) |> collect
# m = MANOVABNPTest.Model(D = 2)
# MANOVABNPTest.fit(m, y, x, grid; iter = 200, rng = rng);

# Exp 9
y, x = test_sample_01(2, 1, 1, 2)
Y = Matrix{Float64}(vcat(y'...))
println(train(Y, x))

# Exp 10
y, x = test_sample_01(2, 1, 1, 2)
Y = Matrix{Float64}(vcat(y'...))
println(train_γ(Y, x))


text:
using Arpack, SparseArrays, LinearAlgebra

# define functions used here
include("ades.jl")
include("acre.jl")
include("pascaltriangle.jl")
include("in2b.jl")
# include("b2in.jl")
include("epsilon.jl")
include("cutMsizeEne.jl")
include("Hintfunccutoff2.jl")

function Hsocfunccutoff!(indvec::Vector{Int64}, Msize0::Int64, Np::Int64, matp::Matrix{Int64}, energyijmat::Matrix{Float64}, ksoc::Float64, Omega::Float64, Hsoc::SparseMatrixCSC{ComplexF64})

    Msize = Msize0*2
    maxmatpcut = length(indvec)

    # matp = zeros(Int,Msize+1,Np+1);
    # pascaltriangle!(Msize,Np,matp) # the size is Msize+1 times Np+1
    # maxmatp = matp[Msize+1,Np+1] # the indices are m+1 and n+1 for N^m_ns

    # defines vectors and matrices
    Hsoc .= 0; #Hsoc = spzeros(ComplexF64,maxmatpcut,maxmatpcut); #spzeros(ComplexF64,maxmatp,maxmatp);

    # define a matrix for the Hamiltonian
    for nn = 1:maxmatpcut #maxmatp

        vecmbnn = spzeros(Int64,Msize+1)
        in2b!(indvec[nn],Msize,Np,matp,vecmbnn) #vecmbnn = in2b(indvec[nn],Msize,Np,matp) #in2b(nn,Msize,Np)
        # energyij = 0.

        # free terms
        for jj = 1:Msize

            vecmbnnj = spzeros(Int64,Msize+1)
            ades!(jj,vecmbnn,vecmbnnj) #vecmbnnj = ades(jj,vecmbnn)
            if vecmbnnj[Msize+1] == 0
               continue
            end

            for ii = 1:Msize

                vecmbnnij = spzeros(Int64,Msize+1)
                acre!(ii,vecmbnnj,vecmbnnij) #vecmbnnij = acre(ii,vecmbnnj)

                # energyij = epsilon(ii,jj,ksoc,Omega)
                energyij = energyijmat[ii,jj]
                if isapprox(energyij,0) #energyij == 0 # energy is Int and zero if ii==jj
                   continue #vecmbnnij[1+Msize] = 0
                end

                for mm = 1:nn

                    vecmbmm = spzeros(Int64,Msize+1)
                    in2b!(indvec[mm],Msize,Np,matp,vecmbmm) #vecmbmm = in2b(indvec[mm],Msize,Np)

                    if vecmbnnij[1:Msize] == vecmbmm[1:Msize]
                       Hsoc[mm,nn] = Hsoc[mm,nn] + energyij*sqrt(vecmbnnij[Msize+1])
                    end

                end

            end

        end

        # # Interactions
        # for ll = 1:Msize
        #
        # end

    end

    # since the Hamiltonian is helmitian
    Hsoc = Hsoc + Hsoc' - spdiagm(diag(Hsoc))

    # return Hsoc

end

function main2(gdown::Float64, gup::Float64, gdu::Float64, ksoc::Float64, Omega::Float64, Msize0::Int64, Np::Int64, Ene0minumhalf::Int64, specnum::Int64)

    # define cutoff of energy in Fock states
    Msize = Msize0*2
    matp = zeros(Int64,Msize+1,Np+1)
    pascaltriangle!(Msize,Np,matp) # the size is Msize+1 times Np+1
    # note the indices are m+1 and n+1 for N^m_n
    indvec = cutMsizeEne(Msize0,Np,matp,Ene0minumhalf)
    maxmatpcut = length(indvec)

    # single-particle Hamiltonian

    # calculate and save energyij
    energyijmat = zeros(Float64,Msize,Msize);
    for ii = 1:Msize
        for jj = 1:Msize
            energyijmat[ii,jj] = epsilon(ii,jj,ksoc,Omega)
        end
    end

    # mat0 = Hsocfunccutoff(indvec,Msize0,Np,ksoc,Omega)
    mat0 = spzeros(ComplexF64,maxmatpcut,maxmatpcut)
    Hsocfunccutoff!(indvec,Msize0,Np,matp,energyijmat,ksoc,Omega,mat0)

    # interaction Hamiltonian

    # # calculate and save Vijkl(n1,n2,n3,n4)
    # nmax = ceil(Int64,Msize/2)-1
    # matA = spzeros(Float64,nmax,nmax,nmax)
    # Vijkl(n1,n2,n3,n4)

    # mat1, mat2, mat3 = Hintfunccutoff(indvec,Msize0,Np)

    mat1 = spzeros(Float64,maxmatpcut,maxmatpcut)
    mat2 = spzeros(Float64,maxmatpcut,maxmatpcut)
    mat3 = spzeros(Float64,maxmatpcut,maxmatpcut)
    Hintfunccutoff2!(indvec,Msize0,Np,matp,mat1,mat2,mat3)

    # diagonalisation
    lambda, phi = eigs(mat0+gdown*mat1+gup*mat2+gdu*mat3,nev=specnum,which=:SR)

    return lambda

end

max_stars_repo_name: leios/SOCED

max_issues_repo_name: leios/SOCED

max_forks_repo_name: leios/SOCED


text:
"""
    Path() -> SystemPath
    Path(fp::Tuple) -> SystemPath
    Path(fp::AbstractString) -> AbstractPath
    Path(fp::P; overrides...) -> P

Responsible for creating the appropriate platform specific path
(e.g., [PosixPath](@ref) and [WindowsPath`](@ref) for Unix and
Windows systems respectively)
"""
function Path end

function Path(fp::T; overrides...) where T<:AbstractPath
    override_fields = keys(overrides)
    override_vals = values(overrides)

    T((s in override_fields ? override_vals[s] : getfield(fp, s) for s in fieldnames(T))...)
end

Path(str::AbstractString) = parse(AbstractPath, str)

"""
    @p_str -> Path

Constructs a [Path](@path) (platform specific subtype of [AbstractPath](@ref)),
such as `p"~/.juliarc.jl"`.
"""
macro p_str(fp)
    return :(Path($fp))
end

function Base.getproperty(fp::T, attr::Symbol) where T <: AbstractPath
    if isdefined(fp, attr)
        return getfield(fp, attr)
    elseif attr === :drive
        return ""
    elseif attr === :root
        return POSIX_PATH_SEPARATOR
    elseif attr === :anchor
        return fp.drive * fp.root
    elseif attr === :separator
        return POSIX_PATH_SEPARATOR
    else
        # Call getfield even though we know it'll error
        # so the message is consistent.
        return getfield(fp, attr)
    end
end

function Base.propertynames(::T, private::Bool=false) where T <: AbstractPath
    public_names = (:drive, :root, :anchor, :separator)
    return private ? Base.merge_names(public_names, fieldnames(T)) : public_names
end

#=
We only want to print the macro string syntax when compact is true and
we want print to just return the string (this allows `string` to work normally)
=#
function Base.print(io::IO, fp::AbstractPath)
    print(io, fp.anchor * join(fp.segments, fp.separator))
end

function Base.show(io::IO, fp::AbstractPath)
    get(io, :compact, false) ? print(io, fp) : print(io, "p\"$fp\"")
end

# Needed for path Cmd interpolation to work correctly.
# If this isn't defined the fact that a path is iterable will result in each segment being
# treated as a separate argument.
Base.arg_gen(fp::AbstractPath) = Base.arg_gen(string(fp))

# Default string constructors for AbstractPath types should fall back to calling `parse`.
(::Type{T})(str::AbstractString) where {T<:AbstractPath} = parse(T, str)

function Base.parse(::Type{P}, str::AbstractString; kwargs...) where P<:AbstractPath
    result = tryparse(P, str; kwargs...)
    result === nothing && throw(ArgumentError("$str cannot be parsed as $P"))
    return result
end

function Base.tryparse(::Type{AbstractPath}, str::AbstractString; debug=false)
    result = nothing
    types = Vector{eltype(PATH_TYPES)}()

    for P in PATH_TYPES
        r = tryparse(P, str)

        # If we successfully parsed the path then save that result
        # and break if we aren't in debug mode, otherwise record how many
        if r !== nothing
            # Only assign the result if it's currently `nothing`
            result = result === nothing ? r : result

            if debug
                push!(types, P)
            else
                break
            end
        end
    end

    if length(types) > 1
        @debug(
            string(
                "Found multiple path types that could parse the string specified ($types). ",
                "Please use a specific `parse` method if $(first(types)) is not the correct type."
            )
        )
    end

    return result
end

Base.convert(T::Type{<:AbstractPath}, x::AbstractString) = parse(T, x)
Base.convert(::Type{String}, x::AbstractPath) = string(x)
Base.promote_rule(::Type{String}, ::Type{<:AbstractPath}) = String
Base.isless(a::P, b::P) where P<:AbstractPath = isless(a.segments, b.segments)

"""
      cwd() -> SystemPath

Get the current working directory.

# Examples
```julia-repl
julia> cwd()
p"/home/JuliaUser"

julia> cd(p"/home/JuliaUser/Projects/julia")

julia> cwd()
p"/home/JuliaUser/Projects/julia"
```
"""
cwd() = Path(pwd())
home() = Path(homedir())
Base.expanduser(fp::AbstractPath) = fp
Base.broadcastable(fp::AbstractPath) = Ref(fp)

# components(fp::AbstractPath) = tuple(drive(fp), root(fp), path(fp)...)

"""
    isrelative(fp::AbstractPath) -> Bool

Returns true if `fp.root` is empty, indicating that it is a relative path.
"""
isrelative(fp::AbstractPath) = isempty(fp.root)

"""
    isabsolute(fp::AbstractPath) -> Bool

Returns true if `fp.root` is not empty, indicating that it is an absolute path.
"""
isabsolute(fp::AbstractPath) = !isempty(fp.root)

# Support immutable indexing API
Base.getindex(fp::AbstractPath, idx) = fp.segments[idx]
Base.firstindex(::AbstractPath) = 1
Base.lastindex(fp::AbstractPath) = length(fp)

# Support iteration protocol
Base.eltype(::AbstractPath) = String
Base.length(fp::AbstractPath) = length(fp.segments)
Base.iterate(fp::AbstractPath, state=1) = iterate(fp.segments, state)

#=
Path Modifiers
===============================================
The following are methods for working with and extracting
path components
=#
"""
    hasparent(fp::AbstractPath) -> Bool

Returns whether there is a parent directory component to the supplied path.
"""
hasparent(fp::AbstractPath) = length(fp.segments) > isrelative(fp)

"""
    parent{T<:AbstractPath}(fp::T) -> T

Returns the parent of the supplied path. If no parent exists
then either "/" or "." will be returned depending on whether the path
is absolute.

# Example
```jldoctest
julia> parent(p"~/.julia/v0.6/REQUIRE")
p"~/.julia/v0.6"

julia> parent(p"/etc")
p"/"

julia> parent(p"etc")
p"."

julia> parent(p".")
p"."
```
"""
Base.parent(fp::AbstractPath) = parents(fp)[end]

"""
    parents{T<:AbstractPath}(fp::T) -> Array{T}

Return all parents of the path. If no parent exists then either "/" or "."
will be returned depending on whether the path is absolute.

# Example
```jldoctest
julia> parents(p"~/.julia/v0.6/REQUIRE")
3-element Array{PosixPath,1}:
 p"~"
 p"~/.julia"
 p"~/.julia/v0.6"

julia> parents(p"/etc")
1-element Array{PosixPath,1}:
 p"/"

julia> parents(p"etc")
1-element Array{PosixPath,1}:
 p"."

julia> parents(p".")
1-element Array{PosixPath,1}:
 p"."
```
"""
function parents(fp::T) where {T <: AbstractPath}
    if hasparent(fp)
        # Iterate from 1:n-1 or 0:n-1 for relative and absolute paths respectively.
        # (i.e., include fp.root when applicable)
        return [Path(fp; segments=fp.segments[1:i]) for i in isrelative(fp):length(fp.segments) - 1]
    elseif fp.segments == tuple(".") || !isempty(fp.root)
        return [fp]
    else
        return [Path(fp; segments=tuple("."))]
    end
end

"""
  *(a::T, b::Union{T, AbstractString, AbstractChar}...) where {T <: AbstractPath} -> T

Concatenate paths, strings and/or characters, producing a new path.
This is equivalent to concatenating the string representations of paths and other strings
and then constructing a new path.

# Example
```jldoctest
julia> p"foo" * "bar"
p"foobar"
```
"""
function Base.:(*)(a::T, b::Union{T, AbstractString, AbstractChar}...) where T <: AbstractPath
    return parse(T, *(string(a), string.(b)...))
end

"""
  /(a::AbstractPath, b::Union{AbstractPath, AbstractString}...) -> AbstractPath

Join the path components into a new full path, equivalent to calling `joinpath`.

# Example
```jldoctest
julia> p"foo" / "bar"
p"foo/bar"

julia> p"foo" / "bar" / "baz"
p"foo/bar/baz"
```
"""
/(root::AbstractPath, pieces::Union{AbstractPath, AbstractString}...) = join(root, pieces...)

"""
    join(root::AbstractPath, pieces::Union{AbstractPath, AbstractString}...) -> AbstractPath

Joins path components into a full path.

# Example
```jldoctest; setup = :(using FilePathsBase: /, join)
julia> join(p"~/.julia/v0.6", "REQUIRE")
p"~/.julia/v0.6/REQUIRE"
```
"""
function join(prefix::AbstractPath, pieces::Union{AbstractPath, AbstractString}...)
    pre = prefix
    segments = String[pre.segments...]

    # Convert any strings to paths
    for p in (isa(x, AbstractPath) ? x : Path(x) for x in pieces)
        # Replace prefix if a piece is actually an absolute path.
        if isabsolute(p)
            pre = p
            segments = String[pre.segments...]
        else
            append!(segments, p.segments)
        end
    end

    return Path(pre; segments=tuple(segments...))
end

function Base.splitext(fp::AbstractPath)
    new_fp, ext = splitext(string(fp))
    return (Path(new_fp), ext)
end

Base.basename(fp::AbstractPath) = fp.segments[end]

Base.splitdir(path::AbstractPath) = (dirname(path), basename(path))

"""
    filename(fp::AbstractPath) -> AbstractString

Extracts the `basename` without the extension.

# Example
```jldoctest
julia> filename(p"~/repos/FilePathsBase.jl/src/FilePathsBase.jl")
"FilePathsBase"

julia> filename(p"~/Downloads/julia-1.4.0-linux-x86_64.tar.gz")
"julia-1.4.0-linux-x86_64.tar"
```
"""
filename(fp::AbstractPath) = rsplit(basename(fp), "."; limit=2)[1]

"""
    extension(fp::AbstractPath) -> AbstractString

Extracts the last extension from a filename if there any, otherwise it returns an empty string.

# Example
```jldoctest
julia> extension(p"~/repos/FilePathsBase.jl/src/FilePathsBase.jl")
"jl"
```
"""
function extension(fp::AbstractPath)
    name = basename(fp)

    tokenized = split(name, '.')
    if length(tokenized) > 1
        return tokenized[end]
    else
        return ""
    end
end

"""
    extensions(fp::AbstractPath) -> AbstractString

Extracts all extensions from a filename if there any, otherwise it returns an empty string.

# Example
```jldoctest
julia> extensions(p"~/repos/FilePathsBase.jl/src/FilePathsBase.jl.bak")
2-element Array{SubString{String},1}:
 "jl"
 "bak"
```
"""
function extensions(fp::AbstractPath)
    name = basename(fp)

    tokenized = split(name, '.')
    if length(tokenized) > 1
        return tokenized[2:end]
    else
        return []
    end
end

"""
    isempty(fp::AbstractPath) -> Bool

Returns whether or not a path is empty.

NOTE: Empty paths are usually only created by `Path()`, as `p""` and `Path("")` will
default to using the current directory (or `p"."`).
"""
Base.isempty(fp::AbstractPath) = isempty(fp.segments)

"""
    normalize(fp::AbstractPath) -> AbstractPath

normalizes a path by removing "." and ".." entries.
"""
function normalize(fp::T) where {T <: AbstractPath}
    p = fp.segments
    result = String[]
    rem = length(p)
    count = 0
    del = 0

    while count < length(p)
        str = p[end - count]

        if str == ".."
            del += 1
        elseif str != "."
            if del == 0
                push!(result, str)
            else
                del -= 1
            end
        end

        rem -= 1
        count += 1
    end

    return Path(fp; segments=tuple(fill("..", del)..., reverse(result)...))
end

"""
    absolute(fp::AbstractPath) -> AbstractPath

Creates an absolute path by adding the current working directory if necessary.
"""
function absolute(fp::AbstractPath)
    result = expanduser(fp)

    if isabsolute(result)
        return normalize(result)
    else
        return normalize(join(cwd(), result))
    end
end

"""
    relative{T<:AbstractPath}(fp::T, start::T=cwd())

Creates a relative path from either the current directory or an arbitrary start directory.
"""
function relative(fp::T, start::T=cwd()) where {T<:AbstractPath}
    curdir = "."
    pardir = ".."

    p = absolute(fp).segments
    s = absolute(start).segments

    p == s && return parse(T, curdir)

    i = 0
    while i < min(length(p), length(s))
        i += 1
        @static if Sys.iswindows()
            if lowercase(p[i]) != lowercase(s[i])
                i -= 1
                break
            end
        else
            if p[i] != s[i]
                i -= 1
                break
            end
        end
    end

    pathpart = p[(i + 1):findlast(x -> !isempty(x), p)]
    prefix_num = findlast(x -> !isempty(x), s) - i - 1
    if prefix_num >= 0
        relpath_ = isempty(pathpart) ?
            tuple(fill(pardir, prefix_num + 1)...) :
            tuple(fill(pardir, prefix_num + 1)..., pathpart...)
    else
        relpath_ = tuple(pathpart...)
    end

    if isempty(relpath_)
        return parse(T, curdir)
    else
        # Our assumption is that relative paths shouldn't have a root or drive.
        # This seems to be consistent with Filesystem on Posix and Windows paths.
        return Path(fp; segments=relpath_, drive="", root="")
    end
end

"""
    canonicalize(path::AbstractPath) -> AbstractPath

Canonicalize a path by making it absolute, `.` or `..` segments and resolving any symlinks
if applicable.

WARNING: Fallback behaviour ignores symlinks and should be extended for paths where
symlinks are permitted (e.g., `SystemPath`s).
"""
canonicalize(fp::AbstractPath) = normalize(absolute(fp))

Base.lstat(fp::AbstractPath) = stat(fp)

"""
    mode(fp::AbstractPath) -> Mode

Returns the `Mode` for the specified path.

# Example
```julia-repl
julia> mode(p"src/FilePathsBase.jl")
-rw-r--r--
```
"""
mode(fp::AbstractPath) = stat(fp).mode
Base.filesize(fp::AbstractPath) = stat(fp).size

"""
    modified(fp::AbstractPath) -> DateTime

Returns the last modified date for the `path`.

# Example
```julia-repl
julia> modified(p"src/FilePathsBase.jl")
2017-06-20T04:01:09
```
"""
modified(fp::AbstractPath) = stat(fp).mtime

"""
    created(fp::AbstractPath) -> DateTime

Returns the creation date for the `path`.

# Example
```julia-repl
julia> created(p"src/FilePathsBase.jl")
2017-06-20T04:01:09
```
"""
created(fp::AbstractPath) = stat(fp).ctime
Base.isdir(fp::AbstractPath) = isdir(mode(fp))
Base.isfile(fp::AbstractPath) = isfile(mode(fp))
Base.islink(fp::AbstractPath) = islink(lstat(fp).mode)
Base.issocket(fp::AbstractPath) = issocket(mode(fp))
Base.isfifo(fp::AbstractPath) = issocket(mode(fp))
Base.ischardev(fp::AbstractPath) = ischardev(mode(fp))
Base.isblockdev(fp::AbstractPath) = isblockdev(mode(fp))

"""
    cp(src::AbstractPath, dst::AbstractPath; force=false, follow_symlinks=false)

Copy the file or directory from `src` to `dst`. An existing `dst` will only be overwritten
if `force=true`. If the path types support symlinks then `follow_symlinks=true` will
copy the contents of the symlink to the destination.
"""
function Base.cp(src::AbstractPath, dst::AbstractPath; force=false, follow_symlinks=false)
    if exists(dst)
        if force
            rm(dst; force=force, recursive=true)
        else
            throw(ArgumentError("Destination already exists: $dst"))
        end
    end

    if !exists(src)
        throw(ArgumentError("Source path does not exist: $src"))
    elseif isdir(src)
        mkdir(dst)

        for fp in readdir(src)
            cp(src / fp, dst / fp; force=force)
        end
    elseif isfile(src)
        write(dst, read(src))
    elseif islink(src)
        follow_symlinks ? symlink(readlink(src), dst) : write(dst, read(canonicalize(src)))
    else
        throw(ArgumentError("Source path is not a file or directory: $src"))
    end

    return dst
end

"""
    mv(src::AbstractPath, dst::AbstractPath; force=false)

Move the file or director from `src` to `dst`. An exist `dst` will only be overwritten if
`force=true`.
"""
function Base.mv(src::AbstractPath, dst::AbstractPath; force=false)
    cp(src, dst; force=force)
    rm(src; recursive=true)
    return dst
end

"""
    sync([f::Function,] src::AbstractPath, dst::AbstractPath; delete=false, overwrite=true)

Recursively copy new and updated files from the source path to the destination.
If delete is true then files at the destination that don't exist at the source will be removed.
By default, source files are sent to the destination if they have different sizes or the source has newer
last modified date.

Optionally, you can specify a function `f` which will take a `src` and `dst` path and return
true if the `src` should be sent. This may be useful if you'd like to use a checksum for
comparison.
"""
function sync(src::AbstractPath, dst::AbstractPath; kwargs...)
    sync(should_sync, src, dst; kwargs...)
end

function sync(f::Function, src::AbstractPath, dst::AbstractPath; delete=false, overwrite=true)
    # Throw an error if the source path doesn't exist at all
    exists(src) || throw(ArgumentError("$src does not exist"))

    # If the top level source is just a file then try to just sync that
    # without calling walkpath
    if isfile(src)
        # If the destination exists then we should make sure it is a file and check
        # if we should copy the source over.
        if exists(dst)
            isfile(dst) || throw(ArgumentError("$dst is not a file"))
            if overwrite && f(src, dst)
                cp(src, dst; force=true)
            end
        else
            cp(src, dst)
        end
    else
        isdir(src) || throw(ArgumentError("$src is neither a file or directory."))
        if exists(dst) && !isdir(dst)
            throw(ArgumentError("$dst is not a directory while $src is"))
        end

        # Create an index of all of the source files
        src_paths = collect(walkpath(src))
        index = Dict(
            Tuple(setdiff(p.segments, src.segments)) => i for (i, p) in enumerate(src_paths)
        )

        if exists(dst)
            for p in walkpath(dst)
                k = Tuple(setdiff(p.segments, dst.segments))

                if haskey(index, k)
                    src_path = src_paths[index[k]]
                    if overwrite && f(src_path, p)
                        cp(src_path, p; force=true)
                    end

                    delete!(index, k)
                elseif delete
                    rm(p; recursive=true)
                end
            end

            # Finally, copy over files that don't exist at the destination
            # But we need to iterate through it in a way that respects the original
            # walkpath order otherwise we may end up trying to copy a file before its parents.
            index_pairs = collect(pairs(index))
            index_pairs = index_pairs[sortperm(last.(index_pairs))]
            for (seg, i) in index_pairs
                cp(src_paths[i], Path(dst, segments=tuple(dst.segments..., seg...)); force=true)
            end
        else
            cp(src, dst)
        end
    end
end

function should_sync(src::AbstractPath, dst::AbstractPath)
    src_stat = stat(src)
    dst_stat = stat(dst)

    if src_stat.size != dst_stat.size || src_stat.mtime > dst_stat.mtime
        @debug(
            "syncing: $src -> $dst, " *
            "size: $(src_stat.size) -> $(dst_stat.size), " *
            "modified_time: $(src_stat.mtime) -> $(dst_stat.mtime)"
        )
        return true
    else
        return false
    end
end

"""
    download(url::Union{AbstractString, AbstractPath}, localfile::AbstractPath)

Download a file from the remote `url` and save it to the `localfile` path.

NOTE: Not downloading into a `localfile` directory matches the base Julia behaviour.
https://github.com/rofinn/FilePathsBase.jl/issues/48
"""
function Base.download(url::AbstractString, localfile::P) where P <: AbstractPath
    mktemp(P) do fp, io
        download(url, fp)
        cp(fp, localfile; force=false)
    end
end

Base.download(url::AbstractPath, localfile::AbstractPath) = cp(url, localfile; force=true)

function Base.download(url::AbstractPath, localfile::AbstractString)
    download(url, Path(localfile))
    return localfile
end

"""
    readpath(fp::P; join=true, sort=true) where {P <: AbstractPath} -> Vector{P}
"""
function readpath(p::P; join=true, sort=true)::Vector{P} where P <: AbstractPath
    @static if VERSION < v"1.4"
        results = readdir(p)
        sort && sort!(results)
        return join ? joinpath.(p, results) : parse.(P, results)
    else
        return parse.(P, readdir(p; join=join, sort=sort))
    end
end

"""
    walkpath(fp::AbstractPath; topdown=true, follow_symlinks=false, onerror=throw)

Performs a depth first search through the directory structure
"""
function walkpath(fp::P; topdown=true, follow_symlinks=false, onerror=throw)  where P <: AbstractPath
    return Channel(ctype=P) do chnl
        for p in readpath(fp)
            topdown && put!(chnl, p)
            if isdir(p) && (follow_symlinks || !islink(p))
                # Iterate through children
                children = walkpath(
                    p; topdown=topdown, follow_symlinks=follow_symlinks, onerror=onerror
                )

                for c in children
                    put!(chnl, c)
                end
            end
            topdown || put!(chnl, p)
        end
    end
end

"""
  open(filename::AbstractPath; keywords...) -> FileBuffer
  open(filename::AbstractPath, mode="r) -> FileBuffer

Return a default FileBuffer for `open` calls to paths which only support `read` and `write`
methods. See base `open` docs for details on valid keywords.
"""
Base.open(fp::AbstractPath; kwargs...) = FileBuffer(fp; kwargs...)

function Base.open(fp::AbstractPath, mode)
    if mode == "r"
        return FileBuffer(fp; read=true, write=false)
    elseif mode == "w"
        return FileBuffer(fp; read=false, write=true, create=true, truncate=true)
    elseif mode == "a"
        return FileBuffer(fp; read=false, write=true, create=true, append=true)
    elseif mode == "r+"
        return FileBuffer(fp; read=true, writable=true)
    elseif mode == "w+"
        return FileBuffer(fp; read=true, write=true, create=true, truncate=true)
    elseif mode == "a+"
        return FileBuffer(fp; read=true, write=true, create=true, append=true)
    else
        throw(ArgumentError("$mode is not support for $(typeof(fp))"))
    end
end

# Fallback read write methods
Base.read(fp::AbstractPath, ::Type{T}) where {T} = open(io -> read(io, T), fp)
Base.write(fp::AbstractPath, x) = open(io -> write(io, x), fp, "w")

# Default `touch` will just write an empty string to a file
Base.touch(fp::AbstractPath) = write(fp, "")

Base.tempname(::Type{<:AbstractPath}) = Path(tempname())
tmpname() = tempname(SystemPath)

Base.tempdir(::Type{<:AbstractPath}) = Path(tempdir())
tmpdir() = tempdir(SystemPath)

Base.mktemp(P::Type{<:AbstractPath}) = mktemp(tempdir(P))
mktmp() = mktemp(SystemPath)

Base.mktemp(fn::Function, P::Type{<:AbstractPath}) = mktemp(fn, tempdir(P))
mktmp(fn::Function) = mktemp(fn, SystemPath)

Base.mktempdir(P::Type{<:AbstractPath}) = mktempdir(tempdir(P))
mktmpdir() = mktempdir(SystemPath)

Base.mktempdir(fn::Function, P::Type{<:AbstractPath}) = mktempdir(fn, tempdir(P))
mktmpdir(fn::Function) = mktempdir(fn, SystemPath)

function Base.mktemp(parent::AbstractPath)
    fp = parent / string(uuid4())
    # touch the file in case `open` isn't implement for the path and
    # we're buffering locally.
    touch(fp)
    io = open(fp, "w+")
    return fp, io
end

function Base.mktempdir(parent::AbstractPath)
    fp = parent / string(uuid4())
    mkdir(fp)
    return fp
end

function Base.mktemp(fn::Function, parent::AbstractPath)
    (tmp_fp, tmp_io) = mktmp(parent)
    try
        fn(tmp_fp, tmp_io)
    finally
        close(tmp_io)
        rm(tmp_fp)
    end
end

function Base.mktempdir(fn::Function, parent::AbstractPath)
    tmpdir = mktmpdir(parent)
    try
        fn(tmpdir)
    finally
        rm(tmpdir, recursive=true)
    end
end

mktmp(arg1, args...) = mktemp(arg1, args...)
mktmpdir(arg1, args...) = mktempdir(arg1, args...)

"""
	isdescendant(fp::P, asc::P) where {P <: AbstractPath} -> Bool

Returns `true` if `fp` is within the directory tree of the `asc`.
"""
isdescendant(fp::P, asc::P) where {P <: AbstractPath} = fp == asc || asc in parents(fp)

"""
	isascendant(fp::P, desc::P) where {P <: AbstractPath} -> Bool

Returns `true` if `fp` is a directory containing `desc`.
"""
isascendant(fp::P, desc::P) where {P <: AbstractPath} = isdescendant(desc, fp)

"""
    diskusage(fp::AbstractPath)

Returns the total size in bytes of the `AbstractPath`.  This is guaranteed to give the same
result as summing the `filesize` of all nodes in `walkpath(fp)` including directory
node block sizes.

This is the preferred method for computing file or directory sizes for remote file
systems as it should limit the number of remote calls where possible.
"""
diskusage(fp::AbstractPath) = isfile(fp) ? filesize(fp) : diskusage(walkpath(fp))
diskusage(itr) = mapreduce(filesize, +, itr)

Base.include(m::Module, path::AbstractPath) = Base.include(identity, m, path)
function Base.include(mapexpr::Function, m::Module, path::AbstractPath)
    tmp_file = cwd() / string(uuid4(), "-", basename(path))
    try
        cp(path, tmp_file; force = true)
        Base.include(mapexpr, m, string(tmp_file))
    finally
        rm(tmp_file; force = true)
    end
end

macro __INCLUDE__()
    return quote
        m = @__MODULE__
        m.include(path::AbstractPath) = Base.include(identity, m, path)
        m.include(mapexpr::Function, path::AbstractPath) = Base.include(mapexpr, m, path)
    end
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

import Base:reverse
"""
  `reverse(fst::WFST)`

Returns `rfst`, the reversed version of `fst`.
If `fst` transuces the sequence `x` to `y` with weight `w`, `rfst` transuces the reverse of `x` to the reverse of `y` with weight `reverse(w)`.
"""
function reverse(fst::WFST{W,A,D}) where {W,A,D}
  RW = FiniteStateTransducers.reversetype(W)
  rfst = WFST(get_isym(fst),get_osym(fst); W = RW)

  add_states!(rfst,length(fst)+1)
  initial!(rfst,D(1)) # 1 is an additional initial superstate

  for (q,arcs) in enumerate(get_states(fst))
    qr = D(q+1) # reversed state
    if isinitial(fst,q)
      w = get_weight(fst,q)
      final!(rfst,qr,reverse(w))
    end
    if isfinal(fst,q)
      w = get_weight(fst,q)
      push!(rfst[1], Arc{RW,D}(D(0),D(0),reverse(w),qr))
    end
    for arc in arcs
      n = get_nextstate(arc)+1
      w = get_weight(arc)
      ilab, olab= get_ilabel(arc), get_olabel(arc)
      push!(rfst[n], Arc{RW,D}(ilab,olab,reverse(w),qr))
    end
  end

  return rfst

end


text:
import ElectricFields: time_integral

@testset "Gaussian pulses" begin
    @testset "Equivalent parameters" begin
        @field(F1) do
            I₀ = 1.0
            T = 2.0
            τ = 3.0
            σmax = 6.0
        end
        env1 = envelope(F1)
        @test env1 isa ElectricFields.GaussianEnvelope
        σ1 = 3.0/(2*√(2log(2)))
        @test 6σ1 ≤ env1.tmax ≤ 6σ1 + 2
        @test time_integral(F1) ≈ σ1*√(2π)
        @test string(env1) == "Gaussian envelope of duration 72.5665 as (intensity FWHM; ±6.28σ)"

        @field(F2) do
            I₀ = 1.0
            T = 2.0
            σ = 3.0
            tmax = 6.0
        end
        env2 = envelope(F2)
        @test env2 isa ElectricFields.GaussianEnvelope
        τ2 = 3.0*(2*√(2log(2)))
        @test env2.τ ≈ τ2
        @test 6.0 ≤ env2.tmax ≤ 6.0 + 2.0

        @field(F3) do
            I₀ = 1.0
            T = 2.0
            σ = 3.0
            Tmax = 3.0
        end
        env3 = envelope(F3)
        @test env3 isa ElectricFields.GaussianEnvelope
        τ2 = 3.0*(2*√(2log(2)))
        @test env3.τ ≈ τ2
        @test env3.tmax == 6.0

        @test continuity(F3) == Inf
        @test span(F3) == -6.0..6.0
        @test time_integral(F3) ≈ 3.0*√(2π)
    end

    @testset "FWHM derivation" begin
        # FWHM is defined for the intensity envelope, which is given
        # by the field amplitude envelope squared, which in turn
        # depends on the derivative of the vector amplitude. Here we
        # check that the creation routine for Gaussian pulses manages
        # to find an FWHM for the vector amplitude that is consistent
        # with the requested FWHM for the intensity envelope.
        @field(F) do
            I₀ = 1.0
            T = 2.0
            τ = 3.0
            σmax = 6.0
        end
        @test intensity(F, 1.5) ≈ 0.5 rtol=1e-7

        @field(F2) do
            I₀ = 1.0
            T = 2.0
            τ = 3.0
            toff = 4.0
            tmax = 5.0
            env = :trunc_gauss
        end
        @test intensity(F2, 1.5) ≈ 0.5 rtol=1e-7
    end

    @testset "Smooth turn-off of truncated Gaussian" begin
        @field(F) do
            I₀ = 1.0
            T = 2.0
            τ = 3.0
            toff = 4.0
            tmax = 5.0
            env = :trunc_gauss
        end
        env = envelope(F)
        @test string(env) == "Truncated Gaussian envelope of duration 3.0000 jiffies = 72.5665 as (intensity FWHM; turn-off from 96.7554 as to 120.9442 as)"

        σ = 3.0/(2*√(2log(2)))
        @test env.σ ≈ σ

        α = env.α
        for t = -3.5:0.5:3.5
            @test env(t) ≈ exp(-α*t^2) rtol=1e-12
        end
        for t = 4.2:0.1:4.9
            @test 0 < env(t) < exp(-α*t^2)
        end

        @test continuity(F) == Inf

        @test span(F) == -5..5

        @test time_integral(F) ≈ σ*√(2π)
    end

    @testset "Truncated Gaussian turn-off via standard deviations" begin
        @field(F) do
            I₀ = 1.0
            T = 2.0
            τ = 3.0
            σoff = 4.0
            σmax = 5.0
            env = :trunc_gauss
        end

        σ = params(F)[:σ]
        @test params(F)[:toff] ≈ 4σ
        @test params(F)[:tmax] ≈ 5σ
    end
end


text:
using Revise
using AdFem
using PyPlot
using LinearAlgebra
using Statistics
using MAT
SMALL_SIZE = 20
MEDIUM_SIZE = 20
BIGGER_SIZE = 20

matplotlib.rc("font", size=SMALL_SIZE)          # controls default text sizes
matplotlib.rc("axes", titlesize=SMALL_SIZE)     # fontsize of the axes title
matplotlib.rc("axes", labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
matplotlib.rc("xtick", labelsize=SMALL_SIZE)    # fontsize of the tick labels
matplotlib.rc("ytick", labelsize=SMALL_SIZE)    # fontsize of the tick labels
matplotlib.rc("legend", fontsize=SMALL_SIZE)    # legend fontsize
matplotlib.rc("figure", titlesize=BIGGER_SIZE)  # fontsize of the figure title

mmesh = Mesh(joinpath(PDATA, "twoholes_large.stl"))
mmesh = Mesh(mmesh.nodes * 10, mmesh.elems, -1, 2)

function f(x, y)
    0.0001*(10/(1+x^2) + x * y + 10*y^2)
end

E = eval_f_on_gauss_pts(f, mmesh)

E2 = matread("fenics/bwd2-601.mat")["nu"]
E3 = matread("fenics/bwd3-601.mat")["nu"]
E4 = matread("fenics/bwd4-481.mat")["nu"]

loss2 = matread("fenics/bwd2-601.mat")["loss"]
loss3 = matread("fenics/bwd3-601.mat")["loss"]
loss4 = matread("fenics/bwd4-481.mat")["loss"]

close("all")
figure(figsize=(20, 5))
subplot(131)
title("Exact")
visualize_scalar_on_gauss_points(E, mmesh)
gca().set_rasterized(true)
subplot(132)
title("DNN")
visualize_scalar_on_gauss_points(E2, mmesh)
xlabel("x")
subplot(133)
title("Difference")
visualize_scalar_on_gauss_points(abs.(E-E2), mmesh)
xlabel("x")
tight_layout()
savefig("Stokes.png")
savefig("Stokes.pdf")


text:
type LatticeShortRateModelEngineCommon{T <: ShortRateTree}
  tg::TimeGrid
  lattice::T
end

function update!(eng::LatticeShortRateModelEngine)
  if length(eng.common.tg.times) > 0
    eng.common.lattice = tree(eng.model, eng.common.tg)
  end

  return eng
end


text:
cd(@__DIR__)
using Pkg
Pkg.activate("..")
using Revise, VisualSearchACTR, Statistics
import_gui()

experiment = Experiment(n_trials=2, trace=true, visible=true, speed=1.0);
stimuli = generate_stimuli(experiment)
model = run_trial!(experiment, stimuli...);


text:
using LinearAlgebra: mul!

function calWeightedQuantity(;
    top     :: Float64,
    bot     :: Float64,
    split_z :: Float64,
    zs      :: AbstractArray{Float64, 1},
    hs      :: AbstractArray{Float64, 1},
    layer   :: Integer,
)

    Δh     = ocn.hs[layer, i, j]
    Δh_top = ocn.zs[layer, i, j] - split_z
    Δh_bot = Δh - Δh_top

    return ( Δh_top * top + Δh_bot * bot ) / Δh

end

function stepOcean_prepare!(ocn::Ocean; cfgs...)

    adv_scheme = cfgs[:adv_scheme]

    if adv_scheme == :static
        return
    end

    # Transform input wind stress vector first
    DisplacedPoleCoordinate.project!(ocn.gi, ocn.in_flds.taux, ocn.in_flds.tauy, ocn.τx, ocn.τy, direction=:Forward)

    if adv_scheme == :ekman_HOOM_partition

        H_ek =  50.0
        H_rf = 250.0
        H_total = H_ek + H_rf

        bot_lay_ek = getLayerFromDepth(
            z  = - H_ek,
            zs = ocn.zs_bone,
            Nz = ocn.Nz_bone,
        )

        bot_lay_rf = getLayerFromDepth(
            z  = - H_total,
            zs = ocn.zs_bone,
            Nz = ocn.Nz_bone,
        )

        @loop_hor ocn i j let

            ΔM̃_half = (ocn.τx[i, j] + ocn.τy[i, j] * im) / (2.0 * ρ_sw * (ocn.ϵs[i, j] + ocn.fs[i, j] * im) )

            ṽ_ek =   ΔM̃_half / H_ek
            ṽ_rf = - ΔM̃_half / H_rf

            u_ek, v_ek = real(ṽ_ek), imag(ṽ_ek)
            u_rf, v_rf = real(ṽ_rf), imag(ṽ_rf)

            if bot_lay_ek == -1

                ocn.u[:, i, j] .= u_ek
                ocn.v[:, i, j] .= v_ek

            else

                ocn.u[1:bot_lay_ek, i, j] .= u_ek
                ocn.v[1:bot_lay_ek, i, j] .= v_ek

                # Mix the top of RF layer
                Δh     = ocn.hs[bot_lay_ek, i, j]
                Δh_top = H_ek + ocn.zs[bot_lay_ek, i, j]
                Δh_bot = Δh - Δh_top

                ocn.u[bot_lay_ek, i, j] = (Δh_top * u_ek + Δh_bot * u_rf) / Δh
                ocn.v[bot_lay_ek, i, j] = (Δh_top * v_ek + Δh_bot * v_rf) / Δh

                if bot_lay_ek < ocn.Nz[i, j] # Bottom layers exists
                    if bot_lay_rf == -1
                       ocn.u[bot_lay_ek+1:end, i, j] .= u_rf
                       ocn.v[bot_lay_ek+1:end, i, j] .= v_rf
                    else
                       ocn.u[bot_lay_ek+1:bot_lay_rf, i, j] .= u_rf
                       ocn.v[bot_lay_ek+1:bot_lay_rf, i, j] .= v_rf

                        # Mix the bottom of RF layer
                        Δh     = ocn.hs[bot_lay_rf, i, j]
                        Δh_top = H_total + ocn.zs[bot_lay_rf, i, j]
                        Δh_bot = Δh - Δh_top

                        ocn.u[bot_lay_rf, i, j] = Δh_top * u_rf / Δh
                        ocn.v[bot_lay_rf, i, j] = Δh_top * v_rf / Δh

                    end
                end

            end
        end

    elseif adv_scheme == :test

        ocn.u .= 0
        ocn.v .= 0

        for k=1:2
             ocn.u[k, :, :] = 1.0 * cos.(ocn.gi.c_lat) .* sin.(ocn.gi.c_lon)
        end

        for k=3:12
             ocn.u[k, :, :] = -1.0/5 * cos.(ocn.gi.c_lat) .* sin.(ocn.gi.c_lon)
        end

    elseif adv_scheme == :ekman_codron2012_partition

        H_ek =  50.0
        H_rf = 250.0   # Codron (2012) suggests 150 - 350 meters. Here I take the average.
        H_total = H_ek + H_rf

        bot_lay_ek = getLayerFromDepth(
            z  = - H_ek,
            zs = ocn.zs_bone,
            Nz = ocn.Nz_bone,
        )

        bot_lay_rf = getLayerFromDepth(
            z  = - H_total,
            zs = ocn.zs_bone,
            Nz = ocn.Nz_bone,
        )

        @loop_hor ocn i j let

            M̃ = (ocn.τx[i, j] + ocn.τy[i, j] * im) / (ρ_sw * (ocn.ϵs[i, j] + ocn.fs[i, j] * im) )

            ṽ_ek =   M̃ / H_ek
            ṽ_rf = - M̃ / H_rf

            u_ek, v_ek = real(ṽ_ek), imag(ṽ_ek)
            u_rf, v_rf = real(ṽ_rf), imag(ṽ_rf)

            if bot_lay_ek == -1

                ocn.u[:, i, j] .= u_ek
                ocn.v[:, i, j] .= v_ek

            else

                ocn.u[1:bot_lay_ek, i, j] .= u_ek
                ocn.v[1:bot_lay_ek, i, j] .= v_ek

                # Mix the top of RF layer
                Δh     = ocn.hs[bot_lay_ek, i, j]
                Δh_top = H_ek + ocn.zs[bot_lay_ek, i, j]
                Δh_bot = Δh - Δh_top

                ocn.u[bot_lay_ek, i, j] = (Δh_top * u_ek + Δh_bot * u_rf) / Δh
                ocn.v[bot_lay_ek, i, j] = (Δh_top * v_ek + Δh_bot * v_rf) / Δh

                if bot_lay_ek < ocn.Nz[i, j] # Bottom layers exists
                    if bot_lay_rf == -1
                       ocn.u[bot_lay_ek+1:end, i, j] .= u_rf
                       ocn.v[bot_lay_ek+1:end, i, j] .= v_rf
                    else
                       ocn.u[bot_lay_ek+1:bot_lay_rf, i, j] .= u_rf
                       ocn.v[bot_lay_ek+1:bot_lay_rf, i, j] .= v_rf

                        # Mix the bottom of RF layer
                        Δh     = ocn.hs[bot_lay_rf, i, j]
                        Δh_top = H_total + ocn.zs[bot_lay_rf, i, j]
                        Δh_bot = Δh - Δh_top

                        ocn.u[bot_lay_rf, i, j] = Δh_top * u_rf / Δh
                        ocn.v[bot_lay_rf, i, j] = Δh_top * v_rf / Δh

                    end
                end

            end
        end

    else
        throw(ErrorException("Unknown advection scheme: " * string(adv_scheme)))
    end

    #println("calHorVelBnd with spmtx")
    tmp = getSpace!(ocn.wksp, :T)
    mul!(view(ocn.u_bnd, :), ocn.ASUM.U_fluxmask_interp_T, view(ocn.u, :))
    mul!(view(ocn.v_bnd, :), ocn.ASUM.V_fluxmask_interp_T, view(ocn.v, :))

    mul_autoflat!(ocn.div,   ocn.ASUM.T_DIVx_U, ocn.u_bnd)
    mul_autoflat!(tmp      , ocn.ASUM.T_DIVy_V, ocn.v_bnd)
    @. ocn.div += tmp

    calVerVelBnd!(
        gi    = ocn.gi,
        Nx    = ocn.Nx,
        Ny    = ocn.Ny,
        Nz    = ocn.Nz,
        w_bnd = ocn.w_bnd,
        hs    = ocn.hs,
        div   = ocn.div,
        mask3 = ocn.mask3,
    )
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
##
using Revise
using GLMakie
import Pkg
Pkg.activate(".")
##
import CSV
using DataFrames
using StaticArrays
import Statistics
using NearestNeighbors
using ProgressMeter
##
filepath = "./Data/STORM/actin_desmin_600nm.csv"
localization_df = DataFrame(CSV.File(filepath))
##
function subset_xy_bounding_box(localization_df::DataFrame, xlim, ylim)
    selection_criteria_x = (xlim[1] .< localization_df.x .< xlim[2])
    selection_criteria_y = (ylim[1] .< localization_df.y .< ylim[2])
    selection_criteria = selection_criteria_x .& selection_criteria_y
    return localization_df[selection_criteria, :]
end

function split_dataframe_per_probe(localization_df::DataFrame)
    probes = sort!(Vector(localization_df[:, :probe]))
    probe_ids = unique(probes)
    return [localization_df[localization_df.probe .== id, :] for id in probe_ids]
end

function pointcloud_from_dataframe(localization_df::DataFrame)
    point_cloud = @inbounds @views [SVector{3, Float32}(localization_df[i, [:x, :y, :z]]...) for i = 1:size(localization_df)[1]]
    return point_cloud
end
##
xlim, ylim = (25000, 30000), (25000, 30000)
df = subset_xy_bounding_box(localization_df, xlim, ylim)
df_list =  split_dataframe_per_probe(df)
point_clouds = pointcloud_from_dataframe.(df_list)
αactinin_tree = KDTree(point_clouds[2])
##
import ThreadTools: tmap1
function RipleyK(r::T, tree::NNTree, area) where {T}
    counts = Statistics.mean([length(inrange(tree, p, r)) for p in tree.data])
    N = length(tree.data)
    return counts*area/N
end

function RipleyK(r, T, area, p::ProgressMeter.AbstractProgress)
    K = RipleyK(r, T, area)
    next!(p)
    return K
end

function RipleyK(r::AbstractVector, tree::NNTree, area)
    p = Progress(length(r))
    return tmap1(x -> RipleyK(x, tree, area, p), r)
end

RipleyK(r, points, area) = RipleyK(r, KDTree(points), area)


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Deterministic equivalent generation #
# =================================== #
function generate!(stochasticprogram::StochasticProgram{N}, structure::DeterministicEquivalent{N}) where N
    # Set the optimizer
    structure.model.moi_backend = optimizer(stochasticprogram)
    # Prepare decisions
    structure.model.ext[:decisions] = structure.decisions
    add_decision_bridges!(structure.model)
    # Generate all stages
    for stage in 1:N
        generate!(stochasticprogram, structure, stage)
    end
    return nothing
end

function generate!(stochasticprogram::StochasticProgram{N}, structure::DeterministicEquivalent{N}, stage::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    dep_model = structure.model
    if stage == 1
        # Define first-stage problem
        generator(stochasticprogram, :stage_1)(dep_model, stage_parameters(stochasticprogram, 1))
    else
        # Sanity check on scenario probabilities
        if num_scenarios(stochasticprogram, stage) > 0
            p = stage_probability(stochasticprogram, stage)
            abs(p - 1.0) <= 1e-6 || @warn "Scenario probabilities do not add up to one. The probability sum is given by $p"
        end
        # Cache current objective and sense
        dep_obj = objective_function(dep_model)
        obj_sense = objective_sense(dep_model)
        obj_sense = obj_sense == MOI.FEASIBILITY_SENSE ? MOI.MIN_SENSE : obj_sense
        # Null objective temporarily in case subproblem objectives are zero
        @objective(dep_model, obj_sense, 0)
        # Define second-stage problems, renaming variables according to scenario.
        stage_two_params = stage_parameters(stochasticprogram, 2)
        visited_objs = collect(keys(object_dictionary(dep_model)))
        # Loop through scenarios and incrementally build deterministic equivalent
        for (i, scenario) in enumerate(scenarios(stochasticprogram))
            # Generate model information for scenario i
            stage_key = Symbol(:stage_, stage)
            generator(stochasticprogram, stage_key)(dep_model, stage_parameters(stochasticprogram, stage), scenario)
            # Update objective and cache the subobjective function for scenario i
            (sub_sense, sub_obj) = get_stage_objective(dep_model, 2, i, Val{N}())
            if obj_sense == sub_sense
                dep_obj += probability(scenario) * sub_obj
            else
                dep_obj -= probability(scenario) * sub_obj
            end
            # Bookkeeping for objects added in scenario i
            for (objkey,obj) in filter(kv->kv.first ∉ visited_objs, object_dictionary(dep_model))
                newkey = if isa(obj, VariableRef) || isa(obj, DecisionRef)
                    # Update variable name to reflect stage and scenario
                    varname = if N > 2
                        varname = add_subscript(add_subscript(JuMP.name(obj), stage), i)
                    else
                        varname = add_subscript(JuMP.name(obj), i)
                    end
                    set_name(obj, varname)
                    # Return new key
                    newkey = Symbol(varname)
                elseif isa(obj, AbstractArray{<:VariableRef}) || isa(obj, AbstractArray{<:DecisionRef})
                    # Update object name to reflect stage and scenario
                    arrayname = if N > 2
                        arrayname = add_subscript(add_subscript(objkey, stage), i)
                    else
                        arrayname = add_subscript(objkey, i)
                    end
                    # Handle each individual variable in array
                    for var in obj
                        # Update variable name to reflect stage and scenario
                        splitname = split(name(var), "[")
                        varname = if N > 2
                            varname = add_subscript(splitname[1], stage)
                        else
                            varname = splitname[1]
                        end
                        set_name(var, @sprintf("%s[%s", add_subscript(varname,i), splitname[2]))
                    end
                    # Return new key
                    newkey = Symbol(arrayname)
                elseif isa(obj, JuMP.ConstraintRef)
                    # Update constraint name to reflect stage and scenario
                    conname = if N > 2
                        conname = add_subscript(add_subscript(JuMP.name(obj), stage), i)
                    else
                        conname = add_subscript(JuMP.name(obj), i)
                    end
                    set_name(obj, conname)
                    # Return new key
                    newkey = Symbol(conname)
                elseif isa(obj, AbstractArray{<:ConstraintRef}) ||
                       isa(obj, AbstractArray{<:GenericAffExpr}) ||
                       isa(obj, AbstractArray{<:DecisionAffExpr}) ||
                       isa(obj, AbstractArray{<:_Constant})
                    # Update obj name to reflect stage and scenario
                    arrayname = if N > 2
                        arrayname = add_subscript(add_subscript(objkey, stage), i)
                    else
                        arrayname = add_subscript(objkey, i)
                    end
                    # Handle each individual constraint/expression in array
                    if obj isa AbstractArray{<:ConstraintRef}
                        for con in obj
                            # Update constraint name to reflect stage and scenario
                            if name(con) == ""
                                # Do not need update unnamed constraints
                                continue
                            end
                            splitname = split(name(con), "[")
                            conname = if N > 2
                                conname = add_subscript(splitname[1], stage)
                            else
                                conname = splitname[1]
                            end
                            set_name(con, @sprintf("%s[%s", add_subscript(conname,i), splitname[2]))
                        end
                    end
                    # Return new key
                    newkey = Symbol(arrayname)
                elseif isa(obj, AbstractJuMPScalar)
                    # Update obj name to reflect stage and scenario
                    exprname = if N > 2
                        exprname = add_subscript(add_subscript(objkey, stage), i)
                    else
                        exprname = add_subscript(objkey, i)
                    end
                    # Return new key
                    newkey = Symbol(exprname)
                elseif isa(obj, _Constant)
                    # Update obj name to reflect stage and scenario
                    constantname = if N > 2
                        constantname = add_subscript(add_subscript(objkey, stage), i)
                    else
                        constantname = add_subscript(objkey, i)
                    end
                    # Return new key
                    newkey = Symbol(constantname)
                else
                    continue
                end
                # Change name of obj to newkey to avoid collisions
                dep_model.obj_dict[newkey] = obj
                delete!(dep_model.obj_dict, objkey)
                # Bookkeep newkey to avoid handling again
                push!(visited_objs, newkey)
            end
        end
        # Update objective
        set_objective_function(dep_model, dep_obj)
        set_objective_sense(dep_model, obj_sense)
    end
    return nothing
end

function clear!(dep::DeterministicEquivalent)
    # Clear decisions
    clear!(dep.decisions)
    # Clear deterministic equivalent model
    empty!(dep.model)
    return nothing
end


text:
using FeldtLib: swedish_postnr2position

@testset "postal numbers" begin

# Can't run this since you need to manually set your Google Maps API key...
@testset "test one known postal number" begin
    #lat, lng = swedish_postnr2position("417 56")
    #@test lat ~ 57.7121456
    #@test lat ~ 11.9484268
end

end


text:
module CircularMessage

# Write your package code here.
greet() = println("CircularMessage")

mutable struct circular_message
    m_size::Int
    #bins::SVector{720,Float64}
    #bins::SArray{720,Float64}
    #values::SVector{720,Float64}
    bins
    values
    #cx_values::SVector{720,Complex}
    # Create a circular message
    function circular_message(m_size::Int = 720)
        x = new()
        x.m_size = m_size
        x.bins = (1/m_size*2pi)*[0.:1.:m_size-1;]
        x.values = ones(m_size).* log(1.0/(2*pi))
    end
end

end


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

# TYPE

mutable struct IPW <: TwoStageModel

    first_stage::OneStageModel
    second_stage::OLS
    pscore::Vector{Float64}
    eweights::ProbabilityWeights{Float64, Float64, Vector{Float64}}

    IPW() = new()
end

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

# FIRST STAGE

function first_stage(::Type{IPW}, MM::Type{<:OneStageModel}, MD::Microdata; kwargs...)

    FSD                    = Microdata(MD)
    FSD.mapping[:response] = MD.mapping[:treatment]

    pop!(FSD.mapping, :treatment)

    fit(MM, FSD; kwargs...)
end

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

# ESTIMATION

function fit(
        ::Type{IPW}, MM::Type{<:OneStageModel}, MD::Microdata; novar::Bool = false, kwargs...
    )

    m = first_stage(IPW, MM, MD; novar = novar)
    fit(IPW, m, MD; novar = novar, kwargs...)
end

function fit(
        ::Type{IPW},
        MM::OneStageModel,
        MD::Microdata;
        novar::Bool = false,
        trim::AbstractFloat = 0.0,
    )

    w = getweights(MD)
    d = getvector(MD, :treatment)
    π = fitted(MM)
    v = [(1.0 - di) / (1.0 - πi) + di / πi for (di, πi) in zip(d, π)]

    v[((trim .> π) .| (1.0 - trim .< π))] .= 0.0

    SSD                   = Microdata(MD)
    SSD.mapping[:control] = asgn(MD.model, (MD.model[:treatment], InterceptTerm{true}()))

    obj              = IPW()
    obj.first_stage  = MM
    obj.second_stage = OLS(SSD)
    obj.pscore       = π
    obj.eweights     = ProbabilityWeights(v)

    _fit!(second_stage(obj), reweight(w, obj.eweights))
    novar || _vcov!(obj, getcorr(obj), w)

    return obj
end

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

# SCORE (MOMENT CONDITIONS)

score(obj::IPW) = lmul!(Diagonal(obj.eweights), score(second_stage(obj)))

# EXPECTED JACOBIAN OF SCORE × NUMBER OF OBSERVATIONS

jacobian(obj::IPW, ::UnitWeights) = jacobian(second_stage(obj), obj.eweights)

function jacobian(obj::IPW, w::AbstractWeights)
    jacobian(second_stage(obj), reweight(w, obj.eweights))
end

# EXPECTED JACOBIAN OF SCORE W.R.T. FIRST-STAGE PARAMETERS × NUMBER OF OBSERVATIONS

function crossjacobian(obj::IPW, ::UnitWeights)

    d = getvector(obj, :treatment)
    π = obj.pscore
    D = [(1.0 - di) / abs2(1.0 - πi) - di / abs2(πi) for (di, πi) in zip(d, π)]

    D[iszero.(obj.eweights)] .= 0.0

    g₁ = jacobexp(obj.first_stage)
    g₂ = score(second_stage(obj))

    return g₂' * lmul!(Diagonal(D), g₁)
end

function crossjacobian(obj::IPW, w::AbstractWeights)

    d = getvector(obj, :treatment)
    π = obj.pscore
    D = [wi * ((1.0 - di) / abs2(1.0 - πi) - di / abs2(πi)) for (di, πi, wi)
         in zip(d, π, w)]

    D[iszero.(obj.eweights)] .= 0.0

    g₁ = jacobexp(obj.first_stage)
    g₂ = score(second_stage(obj))

    return g₂' * lmul!(Diagonal(D), g₁)
end

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

# FITTED VALUES

predict(obj::IPW, MD::Microdata) = predict(second_stage(obj), MD)

# DERIVATIVE OF FITTED VALUES

jacobexp(obj::IPW) = jacobexp(second_stage(obj))

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

# UTILITIES

coefnames(obj::IPW) = coefnames(second_stage(obj))

max_stars_repo_name: lbittarello/Microeconometrics

max_issues_repo_name: lbittarello/Microeconometrics

max_forks_repo_name: lbittarello/Microeconometrics


text:
"""
Functions for finding piecewise linear approximations to concave
    functions at integer points.

    These functions solve the sparsest reduction problem from the paper.
"""

function check_cb_submodular(w::Vector{Float64},tol = 1e-10)
    """
    Given a set of weights defining a cardinality-based function,
    check whether it is submodular.

    w(i) = w_{i-1}
    """
    for i = 2:length(w)-1
        if 2*w[i] - w[i-1] - w[i+1] < -tol
            return false
        end
    end
    return true
end

function check_cb_symmetric(w::Vector{Float64})
    """
    Check whether weights defining a cardinality-based function
    define a symmetric splitting function.

    If there are k weights, this is for a function defined on
    k-1 nodes, since we include w_0.
    """
    k = length(w)
    for i = 1:round(Int64,floor(k/2))
        if w[i] != w[k-i+1]
            return false
        end
    end
    return true
end

function random_CombGad(k::Int64)
    """
    This generates a random combined cardinality-based gadget.

    Returns the parameters as well as the integer
    evaluations.
    """
    J = rand(round(Int64,k/2):k-1)
    a = rand(J)
    b = sort((k-1)*rand(J))
    z0 = 2*rand()
    zk = 2*rand()
    w = CombGad(z0,zk,k,a,b)
    return z0, zk, a, b, w
end

function random_scb_function(k::Int64,prob1=1.0,prob2=1.0)
    """
    Generate a random cardinality-based submodular function.
        k = number of nodes the function is defined on.

    This generates a random combined cardinality-based gadget.

    prob1 = probability that z0 > 0
    prob2 = probability that zk > 0
    """
    J = rand(round(Int64,k/2):k-1)
    a = round.(rand(J),digits = 2)
    while minimum(a) == 0
        a = round.(rand(J),digits = 2)
    end
    b = round.(sort((k-1)*rand(J)),digits = 2)
    while minimum(b) == 0
        b = round.(sort((k-1)*rand(J)),digits = 2)
    end

    if rand() < prob1
        z0 = 2*round.(rand(),digits = 2)
    else
        z0 = 0.0
    end
    if rand() < prob2
        zk = 2*round.(rand(),digits = 2)
    else
        zk = 0.0
    end
    w = CombGad(z0,zk,k,a,b)
    w = w./maximum(w)*3
end

function CombGad(z0,zk,k::Int64,a,b)
    """
    This is a k-node Combined Gadget function.
    Returs the integer values.
    """
    J = length(a)
    @assert(length(a) == length(b))
    @assert(is_increasing(b))

    @assert(b[J] < k)
    @assert(minimum(a) > 0)
    @assert(minimum(b) > 0)
    I = collect(0:k)
    w = z0*(k .- I) + zk*I + sum( a[j]*min.( I*(k-b[j]), (k .- I)*b[j]) for j = 1:J)

    return w
end

function CombGad_func(z0,zk,k::Int64,a,b)
    """
    This is a k-node Combined Gadget function.
    It returns a callable function for any x
    """
    J = length(a)
    if J > 0
        @assert(length(a) == length(b))
        @assert(is_increasing(b))
        @assert(b[J] < k)
        @assert(minimum(a) > 0)
        @assert(minimum(b) > 0)
        f = x->(z0*(k .- x) + zk*x + sum( a[j]*min.( x*(k-b[j]), (k .- x)*b[j]) for j = 1:J))
    else
        f = x->(z0*(k .- x) + zk*x)
    end
    return f
end

function is_increasing(b)
    for i = 1:length(b)-1
        if b[i+1] < b[i]
            return false
        end
    end
    return true
end

is_decreasing(b) = is_increasing(reverse(b))

function plf_to_CombGad(m::Vector{Float64},b::Vector{Float64},k::Int64,f0::Float64)
    """
    Find the Combined Gadget weights for a concave, piecewise linear function f.

    m = slope vector
    b = breakpoints
    k = defines interval [0,k]
    f0 = f(0)

    Output:
    z0, zk, a, b = parameters for Combined Gadget Function
    f = vector of function evaluations at breakpoints.

    More detail
    ------------
    A nonnegative, piecewise linear function f on [0,k] is uniquely determined
    by its breakpoints, slopes, and a single point (0,f(0)).

    Given the slopes m and breakpoints b, and point f(0) of such a function,
    which is additionally concave (meaning m is a decreasing vector).

    We return the parameters z0, zk, a, and b that define the CombGad function
    whose continuous extension is exactly the function f.

    Let J+1 be the number of linear pieces.

    We'll also return the function values at breakpoints: (b_i, f_i)
    """
    @assert(is_decreasing(m))
    J = length(m)-1
    z0 = f0/k
    a = zeros(J)
    f = zeros(J+2)
    f[1] = f0
    bext = [0;b;k]
    for j = 1:J
        a[j] = 1/k*(m[j]-m[j+1])
        f[j+1] = m[j]*(bext[j+1]-bext[j]) + f[j]
    end
    f[J+2] = m[J+1]*(bext[J+2]-bext[J+1])+f[J+1]
    zk = f[end]/k
    return z0, zk, a, b, f[2:end-1]
end

function breakpoints_to_slopes(b,f)
    """
    Given breakpoints and function evaluations at these breakpoints
        for a nonnegative, piecewise linear function on [0,k],
        return the slopes of the curves.

        b[i] = is the ith breakpoint
        f[i] = the piecewise linear function evaluated at b[i]

        Although not technically breakpoints, the first and last entry of b
        should be 0 and k for some integer k, so that we can define the function.
    """
    @assert(length(b) == length(f))
    @assert(minimum(b) == 0)
    @assert(is_increasing(b))
    @assert(b[1] == 0)
    k = b[end]
    @assert(k == round(k))

    J = length(b)-1
    m = zeros(J)
    for i = 1:J
        m[i] = (f[i+1]-f[i])/(b[i+1]-b[i])
    end
    return m

end

# functions for the greedy search
function line_through(x1,y1,x2,y2)
    """
    Return the line through points (x1,y1) and (x2,y2).
    Gives the slope as well as the line as a callable function.
    """
    m = (y2-y1)/(x2-x1)
    inter = y1-m*x1
    l_fun = x->(m*x .+ inter)
    return m, l_fun
end

function line_intersection(m1,m2,x1,x2,y1,y2)
    """
    Given two lines, defined by slopws m1, m2 and
    points (x1,y1), (x2,y2), find the point of intersection.
    """
    if m1 == m2
        println("Lines are parallel")
        return
    end
    b1 = y1 - m1*x1
    b2 = y2 - m2*x2
    x = (b2-b1)/(m1-m2)
    y = m1*x+b1
    return x,y
end

function Refined_SCB_to_CGF(w,epsi=1e-14)
    """
    This allows you to find the CGF with minimum number of linear pieces
    AND the smallest (or approximately the smallest) approximation guarantee
    you can get using that number of pieces.

    It refines the output from one call to the original algoirhtm
    for finding the minimum number of pieces.
    """
    z0, zk, a, b, cgf = SubCardFun_to_CGF_weights(w,epsi,false)
    Jtrue = length(a) + 1
    J = Jtrue
    best_eps = epsi
    # Now find the same number of linear pieces, but better
    # approximations, if possible
    while Jtrue == J && epsi > 1e-10
        epsi = epsi/2
        z0_t, zk_t, a_t, b_t, cgf_t = SubCardFun_to_CGF_weights(w,epsi,false)
        J = length(a_t) + 1
        if J == Jtrue
            z0 = z0_t
            zk = zk_t
            a = a_t
            b = b_t
            cgf = cgf_t
            best_eps = epsi
        end
    end
    return z0, zk, a, b, cgf, best_eps
end

function SubCardFun_to_CGF_weights(w,epsi=1e-14,verbose = false)
    """
    Given function evaluations for a submodular cardinality-based
    function w, return the parameters needed to construct a
    combination of CB-gadgets that approximately model that function.
    """
    # even if the goal is to perfectly fit the function,
    # need to fit it to within a small tolerance for numerical issues
    epsi = max(epsi,1e-14)
    k = length(w)-1
    Lines, Slopes, X, Y, Ranges = get_pwl_approx(w,epsi,verbose)
    z0, zk, a, b, f = Lines_to_CombGadget(Slopes,X,Y,k)
    cgf = CombGad_func(z0,zk,k,a, b)
    return z0, zk, a, b, cgf
end

function Lines_to_CombGadget(Slopes,X,Y,k)
    """
    Given information about a set of lines,
    (namely, the slope and a point for each line)
    return the parameters that describe the piecewise
    linear function as a Combined Gadget Function.
    """
    b1 = Y[1] - Slopes[1]*X[1]
    f0 = b1
    b = zeros(length(X)-1)
    for t = 1:length(X)-1
        x,y = line_intersection(Slopes[t],Slopes[t+1],X[t],X[t+1],Y[t],Y[t+1])
        b[t] = x
    end

    z0, zk, a, b, f = plf_to_CombGad(Slopes,b,k,f0)

    return z0, zk, a, b, f
end

function SymmetricSCB_to_Gadget(w,epsi,verbose = false)
    @assert(is_increasing(w))

    # in the symmetric case we also assume w_0 = 0
    @assert(w[1] == 0)
    r = length(w) - 1
    Lines, M, X, Y, Ranges = get_pwl_approx(w,epsi,false)
    mlast = M[end]

    # In the symmetric case, penalties are increasing
    # and we need to ensure the last linear piece is flat
    if mlast < 0
        # negative sloped lines not needed
        # replace with flat line
        M[end] = 0
        X[end] = r
        Y[end] = w[end]
        # pop!(Lines)
    elseif Ranges[end][1] == r-1 && Ranges[end-1][2] == r-1 && mlast > 0
        # Covered point r-1 twice. Get rid of it
        M[end] = 0
        X[end] = r
        Y[end] = w[end]
        # pop!(Lines)
    elseif mlast > 0
        push!(M,0)
        push!(X,r)
        push!(Y,w[end])
    end
    # push!(Lines, x->w[end]*ones(length(x)))

    b = zeros(length(X)-1)
    for t = 1:length(X)-1
        if M[t] == M[t+1]
            @show r, M
            @show X, Y
            @show mlast
        end
        x,y = line_intersection(M[t],M[t+1],X[t],X[t+1],Y[t],Y[t+1])
        b[t] = x
    end

    a = zeros(length(M)-1)
    for i = 1:length(a)
        a[i] = M[i]- M[i+1]
    end

    return a,b
end

function get_pwl_approx(w,epsi,verbose = false)
    """
    Given penalties w for a submodular cardinality based function,
    find a minimum sized set of linear pieces that provides
    an upper bounding (1+epsi) approximation to w everywhere.
    """
    k = length(w)-1
    @assert(epsi >= 0)
    Lines = Vector()
    Slopes = Vector{Float64}()
    X = Vector{Float64}()
    Y = Vector{Float64}()
    new_i = 1
    Ranges = Vector{Tuple{Int64,Int64}}()

    while new_i <= k+1
        i = new_i
        x2, w2, m, l_fun, new_i = next_line(w,i, epsi; verbose=verbose)

        if i == k || i == k+1
            push!(Ranges,(k-1,k))
        else
            push!(Ranges,(i-1,new_i-2))
        end
        # Save information about the lines
        push!(Lines,l_fun)
        push!(Slopes,m)
        push!(X,x2)
        push!(Y,w2)

    end
    # push!(Ranges,k+1)
    return Lines,Slopes, X, Y, Ranges
end

function next_line(w,i,epsi;verbose = false)
    """
    w is the submodular cardinality-based penalties we're trying to cover.
    i is the next index in w that currently isn't covered by a linear piece.
    Next index i means next integer we need to cover is x[i] = i-1
    """
    k = length(w)-1
    x = collect(0:k)

    # Check for end case 1: x[i] == k-1 or x[i] == k
    if x[i] == k-1 || x[i] == k
        # Add line covering (k-1,w_{k-1}) and (k,w_k)
        xi = x[end-1]
        wi = w[end-1]
        x2 = x[end]
        w2 = w[end]
        m, l_fun = line_through(xi,wi,x2,w2)
        if verbose println("Last line through ($xi,$wi) and ($x2,$w2)") end
        return x2, w2, m, l_fun, k+2
    end

    # Otherwise, start a new line at point (xi,li).
    # This is the worst approx at xi we will accept, maximizing reach.
    xi = x[i]
    li = (1+epsi)*w[i]
    if verbose println("Start a line at ($(xi),$li)") end

    # Visit next point to define a new line
    current = i+1
    x2 = x[current]
    w2 = w[current]
    m, l_fun = line_through(xi,li,x2,w2)
    if verbose println("Candidate line L: through ($x2,$w2)") end

    current += 1

    # Visit points one at a time, greedily finding best cover
    while true
        x_current = x[current]          # integer value we consider
        w_current = w[current]          # true penalty at that value
        l_current = l_fun(x_current)    # approximation provided by current line

        if verbose
            println("Visiting i = $x_current")
            println("(w[$x_current],L[$x_current],(1+eps)w[$x_current]) = ($(w_current),$(l_current),$((1+epsi)*w_current))")
        end

        if l_current < w_current

            if verbose println("Undershot, need new line.") end

            # These are the points defining the new line
            x2 = x_current
            w2 = w_current
            m, l_fun = line_through(xi,li,x2,w2)
            if verbose println("Candidate line L: through ($x2,$w2)") end

            # If this covers last point, we're done
            if x2 == k
                return x2, w2, m, l_fun, k+2
            else
                current += 1
            end

        elseif l_current > (1+epsi)*w_current
            if verbose println("Overshot, return") end

            # we can return the next line, slope m going through point (x2, w2)
            # and also the next point that is now uncovered

            return x2, w2, m, l_fun, current

        elseif l_current <= (1+epsi)*w_current && x_current == k

            # We're done, there's nothing left to cover.
            return x2, w2, m, l_fun, k+2

        else
            # We've covered a point, but there are more points to go.
            current += 1
            x_current = x[current]
        end
    end
end

function clique_round_bound(epsilon)
    """
    Return the upper bound on the number of linear pieces
    needed to bound the clique concave function from [1,k/2]
    for any k.
    """
    return ceil.(log2.(log2.(1 ./ epsilon))) .* ceil.(2*(1 .+ epsilon) ./ sqrt.(2*epsilon))
end

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM