text:
using ModelObjectsLH, ModelParams, Random, Test

mdl = ModelParams;

function calibrated_array_test(N)
    rng = MersenneTwister(23);
    @testset "Calibrated Array $N" begin
        switches = mdl.make_test_calibrated_array_switches(N);
        # objId = ObjectId(:test);
        # sizeV = 2 : (1 + N);
        # sz = (sizeV..., );

        # T1 = Float32;
        # defaultValueM = randn(rng, T1, sz);
        # lbM = defaultValueM .- one(T1);
        # ubM = defaultValueM .+ one(T1);
        # isCalM = defaultValueM .> 0.0;

        # switches = CalibratedArraySwitches(objId, defaultValueM, lbM, ubM, isCalM);
        ca = CalibratedArray(switches);

        @test size(ca) == size(switches.defaultValueM);
        @test validate_ca(ca);

        # Change calibrated values
        pvec = get_pvector(ca);
        p = retrieve(pvec, :calValueV);
        p.value .+= 0.1;
        ca.calValueV = copy(p.value);

        @test mdl.check_fixed_values(ca);
    end
end

function all_fixed_test(N)
    @testset "All fixed" begin
        switches = mdl.make_test_calibrated_array_switches(N; allFixed = true);
        ca = CalibratedArray(switches);
        @test validate_ca(ca);
        valueM = ModelParams.values(ca);
        @test isapprox(valueM, switches.defaultValueM);
    end
end

@testset "Calibrated Array" begin
    for N = 1 : 3
        calibrated_array_test(N);
        all_fixed_test(N);
    end
end

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


text:
"""
    precompute_R2sl([;TRF_min=100e-6, TRF_max=500e-6, T2s_min=5e-6, T2s_max=21e-6, ω1_max=π/TRF_max, B1_max=1.5, greens=greens_superlorentzian])

Pre-compute and interpolate the linearized `R2sl(TRF, α, B1, T2s)` and its derivatives `dR2sldB1(TRF, α, B1, T2s)`, `R2sldT2s(TRF, α, B1, T2s)` etc. in the range specified by the arguments.

The function solves the generalized Bloch equations of an isolated semi-solid pool for values in the specified range, calulates the linearized R2sl that minimizes the error of `zs` at the end of the RF-pulse, and interpolates between the different samples.

# Optional Arguments:
- `TRF_min::Number`: lower bound of the RF-pulse duration range in seconds
- `TRF_max::Number`: upper bound of the RF-pulse duration range in seconds
- `T2s_min::Number`: lower bound of the `T2s` range in seconds
- `T2s_max::Number`: upper bound of the `T2s` range in seconds
- `ω1_max::Number`: upper bound of the Rabi frequency ω1, the default is the frequency of a 500μs long π-pulse
- `B1_max::Number`: upper bound of the B1 range, normalized so that `B1 = 1` corresponds to a perfectly calibrated RF field
- `greens=greens_superlorentzian`: Greens function in the form `G(κ) = G((t-τ)/T2s)`. This package supplies the three Greens functions `greens=greens_superlorentzian` (default), `greens=greens_lorentzian`, and `greens=greens_gaussian`

# Examples
```jldoctest
julia> R2slT = precompute_R2sl();

julia> R2sl, dR2sldB1, R2sldT2s, _ = precompute_R2sl(TRF_min=100e-6, TRF_max=500e-6, T2s_min=5e-6, T2s_max=15e-6, ω1_max=π/500e-6, B1_max=1.3, greens=greens_gaussian);

```
"""
function precompute_R2sl(;TRF_min=100e-6, TRF_max=500e-6, T2s_min=5e-6, T2s_max=21e-6, ω1_max=π/TRF_max, B1_max=1.5, greens=greens_superlorentzian)

    # interpolate super-Lorentzian Green's function for speed purposes
    if greens == greens_superlorentzian
        G = interpolate_greens_function(greens, 0, TRF_max / T2s_min)
    else
        G = greens
    end

    function calculate_z(τ, Ω)
        mfun(p, t) = [1.0]

        function hamiltonian_1D!(du, u, mfun, p::NTuple{3,Any}, t)
            ω1, T2s, g = p
            du[1] = -ω1^2 * quadgk(x -> g((t - x) / T2s) * mfun(p, x)[1], 0.0, t)[1]
        end

        z = solve(DDEProblem(hamiltonian_1D!, [1.0], mfun, (0, τ), (Ω, 1, G)), MethodOfSteps(DP8()))
        return z
    end

    function calculate_R2sl(z_fun, τ, Ω)
        z = z_fun(τ)[1]
        function f!(F, ρv)
            ρ = ρv[1]
            s = sqrt(Complex(ρ^2 - (2Ω)^2))
            x = exp(-ρ * τ / 2) * (cosh(τ * s / 2) + ρ / s * sinh(τ * s / 2))
            F[1] = real(x) - z
        end
        function j!(J, ρv)
            ρ = ρv[1]
            s = sqrt(Complex(ρ^2 - (2Ω)^2))
            x = 2 * Ω^2 * exp(-ρ * τ / 2) * (τ * s * cosh(τ * s / 2) - 2sinh(τ * s / 2)) / s^3
            J[1] = real(x)
        end

        sol = nlsolve(f!, j!, [0.1])
        return sol.zero[1]
    end

    τv = range(TRF_min / T2s_max, TRF_max / T2s_min; length=2^10)
    Ωv = range(0, B1_max * ω1_max * T2s_max; length=2^6)

    A = Matrix{Float64}(undef, length(τv), length(Ωv))
    @batch minbatch=8 for iΩ ∈ 2:length(Ωv)
        τmax = min(τv[end], TRF_max * ω1_max / Ωv[iΩ])
        z_fun = calculate_z(τmax, Ωv[iΩ])
        for iτ in eachindex(τv)
            τ = min(τv[iτ], TRF_max * ω1_max / Ωv[iΩ])
            A[iτ,iΩ] = calculate_R2sl(z_fun, τ, Ωv[iΩ])
        end
    end
    A[:,1] .= A[:,2] # extrapolation hack as the fit does not work with Ω = 0

    f = CubicSplineInterpolation((τv, Ωv), A)
    dfdτ(   τ, Ω) = Interpolations.gradient(f, τ, Ω)[1]
    dfdΩ(   τ, Ω) = Interpolations.gradient(f, τ, Ω)[2]
    d2fdτ2( τ, Ω) = Interpolations.hessian( f, τ, Ω)[1,1]
    d2fdΩ2( τ, Ω) = Interpolations.hessian( f, τ, Ω)[2,2]
    d2fdτdΩ(τ, Ω) = Interpolations.hessian( f, τ, Ω)[2,1]

    R2sl(TRF, α, B1, T2s) = f(TRF / T2s, B1 * α * T2s / TRF) / T2s

    dR2sldT2s(TRF, α, B1, T2s) = -dfdτ(TRF / T2s, B1 * α * T2s / TRF) * TRF / T2s^3 - f(TRF / T2s, B1 * α * T2s / TRF) / T2s^2 + dfdΩ(TRF / T2s, B1 * α * T2s / TRF) * B1 * α / TRF / T2s
    dR2sldB1(TRF, α, B1, T2s)  = dfdΩ(TRF / T2s, B1 * α * T2s / TRF) * α / TRF

    dR2sldω1(TRF, α, B1, T2s)  = dfdΩ(TRF / T2s, B1 * α * T2s / TRF) * B1
    dR2sldTRF(TRF, α, B1, T2s) = dfdτ(TRF / T2s, B1 * α * T2s / TRF) / T2s^2

    dR2sldT2sdω1(TRF, α, B1, T2s) = -d2fdτdΩ(TRF / T2s, B1 * α * T2s / TRF) * B1 * TRF / T2s^2 + d2fdΩ2(TRF / T2s, B1 * α * T2s / TRF) * B1^2 * α / TRF

    dR2sldB1dω1(TRF, α, B1, T2s) = dfdΩ(TRF / T2s, B1 * α * T2s / TRF) + d2fdΩ2(TRF / T2s, B1 * α * T2s / TRF) * B1 * α * T2s  / TRF
    dR2sldT2sTRF(TRF, α, B1, T2s) = -dfdτ(TRF / T2s, B1 * α * T2s / TRF) / T2s^3 - d2fdτ2(TRF / T2s, B1 * α * T2s / TRF) * TRF / T2s^4 - dfdτ(TRF / T2s, B1 * α * T2s / TRF) / T2s^3 + d2fdτdΩ(TRF / T2s, B1 * α * T2s / TRF) * B1 * α / TRF / T2s^2
    dR2sldB1dTRF(TRF, α, B1, T2s) = d2fdτdΩ(TRF / T2s, B1 * α * T2s / TRF) * α / TRF / T2s

    return (R2sl, dR2sldT2s, dR2sldB1, dR2sldω1, dR2sldTRF, dR2sldT2sdω1, dR2sldB1dω1, dR2sldT2sTRF, dR2sldB1dTRF)
end

function evaluate_R2sl_vector(α, TRF, B1, T2s, R2slT, grad_list)
    _R2s = similar(α)
    _dR2sdT2s = similar(α)
    _dR2sdB1 = similar(α)

    for i = 1:length(α)
        _R2s[i] = R2slT[1](TRF[i], α[i], B1, T2s)
    end
    if any(isa.(grad_list, grad_T2s))
        for i = 1:length(α)
            _dR2sdT2s[i] = R2slT[2](TRF[i], α[i], B1, T2s)
        end
    end
    if any(isa.(grad_list, grad_B1))
        for i = 1:length(α)
            _dR2sdB1[i] = R2slT[3](TRF[i], α[i], B1, T2s)
        end
    end
    return (_R2s, _dR2sdT2s, _dR2sdB1)
end

function evaluate_R2sl_vector_OCT(α, TRF, B1, T2s, R2slT, grad_list)
    (_R2s, _dR2sdT2s, _dR2sdB1) = evaluate_R2sl_vector(α, TRF, B1, T2s, R2slT, grad_list)

    _dR2sldω1      = similar(α)
    _dR2sldTRF     = similar(α)
    _dR2sldT2sdω1  = similar(α)
    _dR2sldB1dω1   = similar(α)
    _dR2sldT2sdTRF = similar(α)
    _dR2sldB1dTRF  = similar(α)

    for i = 1:length(α)
        _dR2sldω1[i]      = R2slT[4](TRF[i], α[i], B1, T2s)
        _dR2sldTRF[i]     = R2slT[5](TRF[i], α[i], B1, T2s)
        _dR2sldT2sdω1[i]  = R2slT[6](TRF[i], α[i], B1, T2s)
        _dR2sldB1dω1[i]   = R2slT[7](TRF[i], α[i], B1, T2s)
        _dR2sldT2sdTRF[i] = R2slT[8](TRF[i], α[i], B1, T2s)
        _dR2sldB1dTRF[i]  = R2slT[9](TRF[i], α[i], B1, T2s)
    end
    return (_R2s, _dR2sdT2s, _dR2sdB1, _dR2sldω1, _dR2sldTRF, _dR2sldT2sdω1, _dR2sldB1dω1, _dR2sldT2sdTRF, _dR2sldB1dTRF)
end

ext: jl


text:
x = randn(7,4,5)
k = randn(5,1,4)

## imfiltermtx

for border in ("replicate", "circular", "reflect", "symmetric", "value", "inner")
  @test_approx_eq  imfilter(x,k,border)  imfiltermtx(size(x),k,border)*vec(x)
end
@test_approx_eq  imfilter(x,k)  imfiltermtx(size(x),k)*vec(x)
@test_throws ArgumentError imfiltermtx(size(x),k,"value",1)

## psf2otf

@test_approx_eq  imfilter(x,Images.reflect(k),"circular")  real(ifft( fft(x) .* psf2otf(k,size(x)) ))
@test psf2otf(k) == psf2otf(k,size(k))
@test_throws DimensionMismatch psf2otf(k,(size(k)[1]-1,size(k)[2:end]...))

## padforfilter

for border in ("replicate", "circular", "reflect", "symmetric", "value", "inner")
  @test imfilter(x,k,border) == imfilter(padforfilter(x,size(k),border),k,"inner")
end


text:
using QuantileRegressions, CSV
@testset "npqreg" begin
	mycyle = CSV.read("../examples/mcycle.csv", DataFrame)
	x = mycyle.times
	y = mycyle.accel

	for t = range(0.01,0.99,length=99)
		res = QuantileRegressions.npqreg(Array(y), Array(x), t; h=2)
		@test all(!isnan, res[2])
		@test all(!isnan, res[3])
		xrange = sort(unique(x))
		res = QuantileRegressions.npqreg(Array(y), Array(x), t; xrange=xrange, h=0.5)
	end
	res = QuantileRegressions.npqreg([y[1]], [x[1]], 0.4; xrange=[x[1]], h=2)
	@test !isnan(res[2][1])
	@test isnan(res[3][1])
	res = QuantileRegressions.npqreg([y[1]], [x[1]], 0.8; xrange=[x[1]], h=0.0001)
	@test !isnan(res[2][1])
	@test isnan(res[3][1])
	res = QuantileRegressions.npqreg([y[1]], [x[1]], 0.8; xrange=[x[1].+1], h=0.1)
	@test isnan(res[2][1])
	@test isnan(res[3][1])

	@test all(QuantileRegressions.npqreg([1.0,2.0,4.0,5.0], [0.0,0.0,6.0,12.0], 0.3; h=0.001, xrange=[0.0,6.0,12.0])[2] .≈ [1.0, 4.0, 5.0])
end


text:
@symbol_func function P_F(cur_reactor::AbstractReactor)
  cur_P = K_F(cur_reactor)

  cur_P *= cur_reactor.R_0 ^ 3

  cur_P *= cur_reactor.n_bar ^ 2

  cur_P *= cur_reactor.sigma_v

  cur_P
end


text:
export RCAMspg

function RCAMspg(b, alt, r, eta, Lprev, Rprev, L::Tx, weight, iter)where {ETx<:Number, Tx<:AbstractMatrix{ETx}}

    m,n = size(b);

     #change for different rows and cols 2021.01.28
    #indi = 1:1:(maximum[m,n]);
    indin = 1:1:n;
    indim = 1:1:m;

    R = zeros(Complex{Float32},size(Rprev));

        Wr = 1.;
        Ql = 1.;

    for k = 1:alt
       time1 = @elapsed R = minFrobNormspg(L,b,1, eta,indin,Lprev,weight, iter, Wr);
       println("Time1 = ", string(time1))
       time2 = @elapsed L = minFrobNormspg(R,b',1, eta,indim,Rprev,weight, iter, Ql);
       println("Time1 = ", string(time2))
    end
    return L,R
end


text:
# Draw main plots for http calls per triple pattern per query with partial replication
# Author: Thomas Minier
using Compose
using Gadfly
using RDatasets
using CalliPlots

custom_theme = Theme(
    key_position = :top,
    bar_spacing = 7px,
    key_title_font_size=18px,
    key_label_font_size=15px,
    major_label_font_size=18px,
    minor_label_font_size=12px,
    major_label_font="Computer Modern",
    major_label_color=colorant"#000000",
    minor_label_color=colorant"#000000",
    key_label_color=colorant"#000000",
    bar_highlight=colorant"#000000"
)

Gadfly.push_theme(custom_theme)

function colors()
 return Scale.color_discrete_manual(colorant"#990000", colorant"#FF5733", colorant"#FFC300", colorant"#DAF7A6")
end

serversNames = Dict(
    "http://curiosiphi:4000/watDiv_100" => "S1",
    "http://curiosiphi:4001/watDiv_100" => "S2",
    "http://curiosiphi:4002/watDiv_100" => "S3",
    "http://curiosiphi:4003/watDiv_100" => "S4",
    "http://curiosiphi:4004/watDiv_100" => "S5",
    "http://curiosiphi:3000/watDiv_100" => "S1",
    "http://curiosiphi:3001/watDiv_100" => "S2",
    "http://curiosiphi:3002/watDiv_100" => "S3",
    "http://curiosiphi:3003/watDiv_100" => "S4",
    "http://curiosiphi:3004/watDiv_100" => "S5",
)

function loadCalls(files...)
    dfs = map(df -> by(df, [:target, :strpattern], nrow), map(readtable, files))
    sumCalls = Dict()
    dfs = map(dfs) do df
        patternPerServer = Dict()
        # patterns names are per query
        patternNames = Dict()
        patterns = Set()
        patternCpt = 0
        # rename servers & pattern
        for row in eachrow(df)
            row[:target] = serversNames[row[:target]]
            # register pattern per server
            if (! haskey(patternPerServer, row[:target]))
                patternPerServer[row[:target]] = Set()
            end

            # rename pattern
            if (! haskey(patternNames, row[:strpattern]))
                patternCpt += 1
                patternNames[row[:strpattern]] = "tp$patternCpt"
                patterns = union(patterns, Set(["tp$patternCpt"]))
            end
            row[:strpattern] = patternNames[row[:strpattern]]
            patternPerServer[row[:target]] = union(patternPerServer[row[:target]], Set([row[:strpattern]]))
            # register nb calls per pattern
            if (! haskey(sumCalls, row[:strpattern]))
                sumCalls[row[:strpattern]] = 0
            end
            sumCalls[row[:strpattern]] += row[:x1]
        end
        # fill missing values
        for (s, p) in patternPerServer
            missings = setdiff(patterns, p)
            for pattern in p
                push!(df, [s, pattern, 0])
            end
        end
        return df
    end
    return (meanDF(:x1, dfs...), sumCalls)
end

# Build dataframes

(calls_q90, sumCalls90) = loadCalls(
    "scripts/http_calls/expensive_queries/q90_http_calls_4servers_partial.csv"
)
calls_q90[:query] = "Q90"

(calls_q48, sumCalls48) = loadCalls(
    "scripts/http_calls/expensive_queries/q48_http_calls_4servers_partial.csv"
)
calls_q48[:query] = "Q48"

(calls_q17, sumCalls17) = loadCalls(
    "scripts/http_calls/expensive_queries/q17_http_calls_4servers_partial.csv"
)
calls_q17[:query] = "Q17"

(calls_q54, sumCalls54) = loadCalls(
    "scripts/http_calls/expensive_queries/q54_http_calls_4servers_partial.csv"
)
calls_q54[:query] = "Q54"

(calls_q52, sumCalls52) = loadCalls(
    "scripts/http_calls/expensive_queries/q52_http_calls_4servers_partial.csv"
)
calls_q52[:query] = "Q52"

all = [calls_q90;calls_q48;calls_q17;calls_q54;calls_q52]
group1 = [calls_q90;calls_q48;calls_q52]
group2 = [calls_q17;calls_q54]

# Build plots

plot1 = plot(group1, xgroup="query", y="x1", x="strpattern", color="target",
Geom.subplot_grid(Geom.bar(position=:stack), free_x_axis=true, Guide.xticks(orientation=:vertical)),
Guide.xlabel("Triple pattern by query"), Guide.ylabel("# HTTP requests"),
Guide.colorkey("TPF servers"), colors())

plot2 = plot(group2, xgroup="query", y="x1", x="strpattern", color="target",
Geom.subplot_grid(Geom.bar(position=:stack), free_x_axis=true, Guide.xticks(orientation=:vertical)),
Guide.xlabel("Triple pattern by query"), Guide.ylabel("# HTTP requests"),
Guide.colorkey("TPF servers"), colors())

draw(PDF("scripts/curio/http_calls_partial.pdf", 10inch, 4inch), hstack(plot1, plot2))
draw(PNG("scripts/curio/http_calls_partial.png", 10inch, 4inch), hstack(plot1, plot2))

max_stars_repo_name: Callidon/ulysses-tpf

max_issues_repo_name: Callidon/ulysses-tpf

max_forks_repo_name: Callidon/ulysses-tpf


text:
"""
Bridges `CP.BinPacking{CP.VARIABLE_CAPACITY_BINPACKING, T}` to
`CP.BinPacking{CP.NO_CAPACITY_BINPACKING, T}` by adding constraints on the
capacity variables.
"""
struct VariableCapacityBinPacking2BinPackingBridge{T} <: MOIBC.AbstractBridge
    bp::MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.BinPacking{CP.NO_CAPACITY_BINPACKING, T}}
    capa::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}}
end

function MOIBC.bridge_constraint(
    ::Type{VariableCapacityBinPacking2BinPackingBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.BinPacking{CP.VARIABLE_CAPACITY_BINPACKING, T},
) where {T}
    return MOIBC.bridge_constraint(
        VariableCapacityBinPacking2BinPackingBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{VariableCapacityBinPacking2BinPackingBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.BinPacking{CP.VARIABLE_CAPACITY_BINPACKING, T},
) where {T}
    f_scalars = MOIU.scalarize(f)

    # Create the simpler BinPacking set.
    bp_set = CP.BinPacking{CP.NO_CAPACITY_BINPACKING}(s.n_bins, s.n_items, s.weights)
    new_f = MOIU.vectorize([f_scalars[1:s.n_bins]..., f_scalars[(2 * s.n_bins + 1):end]...])
    bp = MOI.add_constraint(model, new_f, bp_set)

    # Add the capacity constraints.
    capa = [
        MOI.add_constraint(model, f_scalars[i] - f_scalars[s.n_bins + i], MOI.LessThan(zero(T)))
        for i in 1:s.n_bins
    ]

    return VariableCapacityBinPacking2BinPackingBridge(bp, capa)
end

function MOI.supports_constraint(
    ::Type{VariableCapacityBinPacking2BinPackingBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.BinPacking{CP.VARIABLE_CAPACITY_BINPACKING, T}},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{<:VariableCapacityBinPacking2BinPackingBridge})
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(::Type{VariableCapacityBinPacking2BinPackingBridge{T}}) where {T}
    return [
        (MOI.VectorAffineFunction{T}, CP.BinPacking{CP.NO_CAPACITY_BINPACKING, T}),
        (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
    ]
end

MOI.get(b::VariableCapacityBinPacking2BinPackingBridge, ::MOI.NumberOfVariables) = 0

function MOI.get(
    ::VariableCapacityBinPacking2BinPackingBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T},
        CP.BinPacking{CP.NO_CAPACITY_BINPACKING, T},
    },
) where {T}
    return 1
end

function MOI.get(
    b::VariableCapacityBinPacking2BinPackingBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T},
        MOI.LessThan{T},
    },
) where {T}
    return length(b.capa)
end

function MOI.get(
    b::VariableCapacityBinPacking2BinPackingBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T},
        CP.BinPacking{CP.NO_CAPACITY_BINPACKING, T},
    },
) where {T}
    return [b.bp]
end

function MOI.get(
    b::VariableCapacityBinPacking2BinPackingBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T},
        MOI.LessThan{T},
    },
) where {T}
    return copy(b.capa)
end


text:
using Pkg
Pkg.activate(".")
using PlotlyJS
using Dash

plot_width = 500
plot_height = 500

sd = 1.1
layout = Layout(xaxis_title="x", yaxis_title="y",
    xaxis_range=[-sd,sd],yaxis_range=[-sd,sd], title="Lissajous Curve",
    width=plot_width, height=plot_height)

l = 1
t = LinRange(0.0,l*2*pi,500)

function myplot(ω1,ω2)
    x = sin.(ω1*t)
    y = cos.(ω2*t)
    plot_data =[scatter(;x=x, y=y, mode="lines", name="Solution")]
    plot(plot_data, layout)
end

mathjax = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.4/MathJax.js?config=TeX-MML-AM_CHTML"
app = dash(external_scripts=[mathjax])
app.title = "DashLissajous"

curve = myplot(1,1)

curve_graph = dcc_graph(id="graph",figure=curve)

slider1 = dcc_slider(id="slider1",min=1, max=10,step=1,value=1, marks=Dict([i => ("$(i)") for i = 1:10]))

slider2 = dcc_slider(id="slider2",min=1,max=10, step=1,value=1, marks=Dict([i => ("$(i)") for i = 1:10]))

# App Layout
app.layout = html_div(
    [
        html_center(html_h1("Lissajous Curves")),
        html_center(curve_graph),
        html_center(html_div(style=Dict("width"=>"420px"),[slider1,html_p(raw"\(\omega_1\)"),slider2, html_p(raw"\(\omega_2\)")]))
    ]
)

callback!(app,Output("graph","figure"), Input("slider1","value"), Input("slider2","value")) do ω1, ω2
    global ω1_old, ω2_old
    myplot(ω1,ω2)
end

run_server(app, "0.0.0.0", debug=true)


text:
export SymSemiseparable

struct SymSemiseparable <: SymSemiseparableMatrix
    n::Int64
    p::Int64
    U::AbstractArray
    V::AbstractArray
end

# Constuctors
function SymSemiseparable(U::AbstractArray, V::AbstractArray)
    if size(U,1) == size(V,1) && size(U,2) == size(V,2)
        return SymSemiseparable(size(U,1),size(U,2),U,V);
    else
        error("Dimension mismatch between generators U and V")
    end
end

# Mappings
mul!(y::AbstractArray, L::SymSemiseparable, x::AbstractArray) =
    ss_mul_mat!(y, L.U, L.V, x);
mul!(y::AbstractArray, L::AdjointOperator{SymSemiseparable}, x::AbstractArray) =
    ss_mul_mat!(y, L.A.U, L.A.V, x);
function inv!(y, K::SymSemiseparable, b::AbstractArray)
	L = SymSemiseparableChol(K);
	y[:,:] = L'\(L\b);
end
function inv!(y, K::AdjointOperator{SymSemiseparable}, b::AbstractArray)
	L = SymSemiseparableChol(K.A);
	y[:,:] = L'\(L\b);
end

#####################################################################
#### Extended generator representable {p}-semiseperable matrices ####
#####################################################################

#### Matrix-matrix product ####
function ss_mul_mat!(Y::AbstractArray, U::AbstractArray, V::AbstractArray, X::AbstractArray)
    n, m = size(U);
    mx = size(X,2);
    Vbar = zeros(m,mx);
    Ubar = U'*X;
    for i = 1:n
        tmpV = V[i,:];
        tmpU = U[i,:];
        Ubar -= tmpU .* X[i:i,:];
        Vbar += tmpV .* X[i:i,:];
        Y[i,:] = Vbar'*tmpU + Ubar'*tmpV;
    end
end

#### Matrix vector product ####
function ss_create_v(U::AbstractArray, W::AbstractArray)
    n,m = size(U)
    V = zeros(n,m)
    P = zeros(m,m);
    for i = 1:n
        tmpW = W[i,:]
        tmpU = U[i,:]
        tmpV = tmpW*(tmpU'*tmpW)
        V[i,:] = tmpV + P*tmpU
        P += tmpW*tmpW';
    end
    return V
end


text:
include("weights_and_ideal.jl")

mutable struct MOEAD_DE <: AbstractParameters
    nobjectives::Int
    N::Int
    F::Float64
    CR::Float64
    λ::Array{Vector{Float64}}
    η::Float64
    p_m::Float64
    T::Int
    δ::Float64
    n_r::Float64
    z::Vector{Float64}
    B::Array{Vector{Int}}
    s1::Float64
    s2::Float64
    τ::Float64
end

"""
    MOEAD_DE(weights ;
        F = 0.5,
        CR = 1.0,
        λ = Array{Vector{Float64}}[], # ref. points
        η = 20,
        p_m = -1.0,
        T = round(Int, 0.2*length(weights)),
        δ = 0.9,
        n_r = round(Int, 0.02*length(weights)),
        z = zeros(0),
        B = Array{Int}[],
        s1 = 0.01,
        s2 = 20.0,
        information = Information(),
        options = Options())

`MOEAD_DE` implements the original version of MOEA/D-DE. It uses the contraint handling method
based on the sum of violations (for constrained optimizaton):
`g(x, λ, z) = max(λ .* abs.(fx - z)) + sum(max.(0, gx)) + sum(abs.(hx))`

To use MOEAD_DE, the output from the objective function should be a 3-touple
`(f::Vector, g::Vector, h::Vector)`, where `f` contains the objective functions,
`g` and `h` are the equality and inequality constraints respectively.

A feasible solution is such that `g_i(x) ≤ 0 and h_j(x) = 0`.

Ref. Multiobjective Optimization Problems With Complicated Pareto Sets,
MOEA/D and NSGA-II; Hui Li and Qingfu Zhang.

# Example

Assume you want to solve the following optimizaton problem:

Minimize:

`f(x) = (x_1, x_2)`

subject to:

`g(x) = x_1^2 + x_2^2 - 1 ≤ 0`

`x_1, x_2 ∈ [-1, 1]`

A solution can be:

```julia

# Dimension
D = 2

# Objective function
f(x) = ( x, [sum(x.^2) - 1], [0.0] )

# bounds
bounds = [-1 -1;
           1  1.0
        ]

nobjectives = 2
npartitions = 100

# reference points (Das and Dennis's method)
weights = gen_ref_dirs(nobjectives, npartitions)

# define the parameters
moead_de = MOEAD_DE(weights, options=Options(debug=false, iterations = 250))

# optimize
status_moead = optimize(f, bounds, moead_de)

# show results
display(status_moead)
```

"""
function MOEAD_DE(weights ;
    F = 0.5,
    CR = 1.0,
    λ = Array{Vector{Float64}}[],
    η = 20,
    p_m = -1.0,
    T = round(Int, 0.2*length(weights)),
    δ = 0.9,
    n_r = round(Int, 0.02*length(weights)),
    z::Vector{Float64} = zeros(0),
    B = Array{Int}[],
    s1 = 0.01,
    s2 = 20.0,
    information = Information(),
    options = Options(),
)

    if isempty(weights)
        error("Provide weights points")
    end

    nobjectives = length(weights[1])
    N = length(weights)

    if isempty(z)
        z = fill(Inf, nobjectives)
    end

    parameters = MOEAD_DE(nobjectives, N, F, CR, weights, η, p_m, T, δ, n_r, z, B, s1, s2, 0.0)

    initialize_closest_weight_vectors!(parameters)

    alg = Algorithm(
        parameters,
        information = information,
        options = options,
    )

    alg

end

function initialize!(
    status,
    parameters::MOEAD_DE,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
)

    if options.iterations == 0
        options.iterations = 500
    end

    if options.f_calls_limit == 0
        options.f_calls_limit = options.iterations * parameters.N + 1
    end

    status = gen_initial_state(problem,parameters,information,options,status)
    D = size(problem.bounds, 2)
    parameters.nobjectives = length(status.population[1].f)
    parameters.p_m = parameters.p_m < 0.0 ? 1.0 / D : parameters.p_m

    update_reference_point!(parameters.z, status.population)
    return status

end

function update_state!(
    status::State,
    parameters::MOEAD_DE,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
)

    F = parameters.F
    CR = parameters.CR

    D = size(problem.bounds, 2)

    N = parameters.N

    la = problem.bounds[1, :]
    lb = problem.bounds[2, :]

    population = status.population

    for i = 1:N
        if rand() < parameters.δ
            P_idx = copy(parameters.B[i])
        else
            P_idx = collect(1:N)
        end

        Vmin = minimum(sum_violations.(population[P_idx]))
        Vmax = maximum(sum_violations.(population[P_idx]))
        parameters.τ = Vmin + 0.3*(Vmax - Vmin)

        # select participats
        r1 = i
        r2 = rand(P_idx)
        while r1 == r2
            r2 = rand(P_idx)
        end

        r3 = rand(P_idx)
        while r3 == r1 || r3 == r2
            r3 = rand(P_idx)
        end

        a = get_position(population[r1])# population[r1].x
        b = get_position(population[r2])# population[r2].x
        c = get_position(population[r3])# population[r3].x

        # binomial crossover
        v = zeros(D)

        # binomial crossover
        for j = 1:D
            # binomial crossover
            if rand() < CR
                v[j] = a[j] + F * (b[j] - c[j])
            else
                v[j] = a[j]
            end
            # polynomial mutation

            if rand() < parameters.p_m
                r = rand()
                if r < 0.5
                    σ_k = (2.0 * r)^(1.0 / (parameters.η + 1)) - 1
                else
                    σ_k = 1 - (2.0 - 2.0 * r)^(1.0 / (parameters.η + 1))
                end
                v[j] = v[j] + σ_k * (lb[j] - la[j])
            end
        end

        v = replace_with_random_in_bounds!(v, problem.bounds)

        # instance child
        h = create_solution(v, problem)
        status.f_calls += 1

        update_reference_point!(parameters.z, h)

        c = 0

        z = parameters.z
        shuffle!(P_idx)
        while c < parameters.n_r && !isempty(P_idx)
            j = pop!(P_idx)
            g1 = g(h.f, parameters.λ[j], z)
            g2 = g(population[j].f, parameters.λ[j], z)
            if is_better_constrained_MOEAD_DE(g1, g2, h, population[j], parameters)
                population[j] = h
                c += 1
            end

        end

        stop_criteria!(status, parameters, problem, information, options)
        if status.stop
            break
        end
    end

end

g_te_ap(gx, V, τ, s1, s2) = V < τ ? gx + s1*V^2 : gx + s1*τ^2 + s2*(V - τ)

function is_better_constrained_MOEAD_DE(g1, g2, sol1, sol2, parameters)
    V1 = (sol1.sum_violations)
    V2 = (sol2.sum_violations)

    if V1 == 0.0 && V2 == 0.0
        return g1 <= g2
    end

    τ = parameters.τ
    s1 = parameters.s1
    s2 = parameters.s2

    V1 = (sol1.sum_violations)
    V2 = (sol2.sum_violations)

    return g_te_ap(g1, V1, τ, s1, s2) <= g_te_ap(g2, V2, τ, s1, s2)
end

function stop_criteria_moead_de(
    status::State,
    parameters::MOEAD_DE,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
    )
    return status.iteration > options.iterations
end
function final_stage!(
    status::State,
    parameters::MOEAD_DE,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
    )
    status.final_time = time()
    # status.best_sol = get_pareto_front(status.population, is_better_eca)
    # @show length(status.best_sol)
end


text:
using TropicalTensors
using TropicalTensors: potts_randpm_J, SquareLattice
using Test, Random

function exact_solve_potts3(lt::SquareLattice, jd::Dict)
    eng_list = Float64[]
    for ci in CartesianIndices((fill(3, lt.Nx*lt.Ny)...,))
        push!(eng_list, potts3_energy(lt, jd, ci.I))
    end
    max_energy = maximum(eng_list)
    degeneracy = count(==(max_energy), eng_list)
    return max_energy, degeneracy
end

function potts3_energy(lt::SquareLattice, jd::Dict, cfg)
    LI = LinearIndices((lt.Nx, lt.Ny))
    eng = 0.0
    for i=1:lt.Nx
        for j=1:lt.Ny
            i!=lt.Nx && if cfg[LI[i,j]] == cfg[LI[i+1,j]]
                eng += jd[(i,j)=>(i+1,j)]
            else
                eng -= 0.5*jd[(i,j)=>(i+1,j)]
            end
            j!=lt.Ny && if cfg[LI[i,j]] == cfg[LI[i,j+1]]
                eng += jd[(i,j)=>(i,j+1)]
            else
                eng -= 0.5*jd[(i,j)=>(i,j+1)]
            end
        end
    end
    return eng
end

@testset "check potts" begin
    L = 9
    lt = SquareLattice(L, L)
    CI = CartesianIndices(lt)
    for i=1:10
        Random.seed!(i)
        dj = potts_randpm_J(lt)
        sj = Float64[]
        for b in sgbonds(lt)
            push!(sj, dj[CI[b[1]].I => CI[b[2]].I])
        end
        res1 = solve_potts(CountingTropicalF64, Val(2), lt, dj; usecuda=false)
        res2 = solve(CountingTropicalF64, Spinglass(lt, sj, zeros(L*L)); usecuda=false)
        @test res1.n ≈ res2.n
        @test res1.c ≈ res2.c
    end
end

@testset "compare exact" begin
    L = 3
    lt = SquareLattice(L, L)
    jd = potts_randpm_J(lt)
    res1 = solve_potts(CountingTropicalF64, Val(3), lt, jd; usecuda=false)
    eng, deg = exact_solve_potts3(lt, jd)
    @test res1.n == eng
    @test res1.c == deg
end


text:
################################################################################
# Distance matrix
################################################################################
"""
    distancematrix(x [, y = x], metric = "euclidean")

Create a matrix with the distances between each pair of points of the
time series `x` and `y` using `metric`.

The time series `x` and `y` can be `AbstractDataset`s or vectors or matrices with data points
in rows.
The data point dimensions (or number of columns) must be the same for `x` and `y`.
The returned value is a `n×m` matrix, with `n` being the length (or number of rows)
of `x`, and `m` the length of `y`.

The metric can be identified by a string, or any of the `Metric`s defined in
the [`Distances` package](https://github.com/JuliaStats/Distances.jl).
The list of strings available to define the metric are:

* `"max"` or `"inf"` for the maximum or L∞ norm
  (`Chebyshev()` in the `Distances` package).
* `"euclidean"` for the L2 or Euclidean norm, used by default
  (`Euclidean()` in `Distances`).
* `"manhattan"`, `"cityblock"`, `"taxicab"` or `"min"` for the Manhattan or L1 norm
  (`Cityblock()` in `Distances`).
"""
distancematrix(x, metric::Union{Metric,String}=DEFAULT_METRIC, parallel = size(x)[1] > 500) = _distancematrix(x, getmetric(metric), Val(parallel))

# For 1-dimensional arrays (vectors), the distance does not depend on the metric
distancematrix(x::Vector, y::Vector, metric=DEFAULT_METRIC, parallel = size(x)[1] > 500) = abs.(x .- y')

# If the metric is supplied as a string, get the corresponding Metric from Distances
distancematrix(x, y, metric::String, parallel = size(x)[1] > 500) = distancematrix(x, y, getmetric(metric), parallel)

const MAXDIM = 9
function distancematrix(x::Tx, y::Ty, metric::Metric=DEFAULT_METRIC, parallel = size(x)[1] > 500) where
         {Tx<:Union{AbstractMatrix, AbstractDataset}} where {Ty<:Union{AbstractMatrix, AbstractDataset}}
    sx, sy = size(x), size(y)
    @assert sx[2] == sy[2] """
        The dimensions of the data points in `x` and `y` must be equal!
        Found dim(x)=$(sx[2]), dim(y)=$(sy[2]).
        """

    if sx[2] ≥ MAXDIM && typeof(metric) == Euclidean # Blas optimization
        return _distancematrix(Matrix(x), Matrix(y), metric, Val(parallel))
    elseif Tx <: Matrix && Ty <: Matrix && metric == Chebyshev()
        return _distancematrix(x, y, metric, Val(parallel))
    else
        return _distancematrix(Dataset(x), Dataset(y), metric, Val(parallel))
    end
end

# Core function for Matrices (wrapper of `pairwise` from the Distances package)
_distancematrix(x::AbstractMatrix, y::AbstractMatrix, metric::Metric, ::Val{false}) = pairwise(metric, x', y', dims=2)

# First we define the serial versions of the functions.
# Core function for AbstractDatasets
function _distancematrix(x::AbstractDataset{S,Tx}, y::AbstractDataset{S,Ty},
    metric::Metric, ::Val{false}) where {S, Tx, Ty}

    x = x.data
    y = y.data
    d = zeros(promote_type(Tx,Ty), length(x), length(y))
    for j in eachindex(y), i in eachindex(x)
        @inbounds d[i,j] = evaluate(metric, x[i], y[j])
    end
    return d
end

# Now, we define the parallel versions.

function _distancematrix(x::AbstractDataset{S,Tx}, y::AbstractDataset{S,Ty},
    metric::Metric, ::Val{true}) where {S, Tx, Ty}

    x = x.data
    y = y.data
    d = zeros(promote_type(Tx,Ty), length(x), length(y))
    Threads.@threads for j in eachindex(y)
        for i in eachindex(x)
            @inbounds d[i,j] = evaluate(metric, x[i], y[j])
        end
    end
    return d
end

function _distancematrix(x::Matrix{Tx}, y::Matrix{Ty},
    metric::Metric, ::Val{true}) where {Tx, Ty}

    x = x.data
    y = y.data
    d = zeros(promote_type(Tx,Ty), length(x), length(y))
    Threads.@threads for j in eachindex(y)
        for i in eachindex(x)
            @inbounds d[i,j] = evaluate(metric, x[i, :], y[j, :])
        end
    end
    return d
end

# Now, we define methods for the single trajectory.
# We can speed this up by only calculating the lower triangle,
# which halves the number of computations needed.

# Again, we'll define the serial version first:
function _distancematrix(x::Vector{T}, metric::Metric, ::Val{false}) where T
    d = zeros(T, length(x), length(x))

    for j in eachindex(x)
        for i in 1:j
            @inbounds d[i, j] = abs(x[i] - x[j])
        end
    end

    return Symmetric(d, :U)

end

function _distancematrix(x::AbstractDataset{S, T}, metric::Metric, ::Val{false}) where T where S
    d = zeros(T, length(x), length(x))

    for j in eachindex(x)
        for i in 1:j
            @inbounds d[i, j] = evaluate(metric, x[i], x[j])
        end
    end

    return Symmetric(d, :U)

end

# Now, we define the parallel version.  There's a twist, though.

# The single trajectory case is a little tricky.  If you try it naively,
# using the method we used for the serial computation above, the points are
# partioned out unequally between threads.  This means that performance actually
# **decreases** compared to the full parallel version.  To mitigate this, we need
# to partition the computation equally among all threads.
# The function I've written below does essentially this - given a length,
# it calculates the appropriate partitioning, and returns a list of indices,
# by utilizing the fact that the area is proportional to the square of the height.
# It partitions the "triangle" which needs to be computed into "trapezoids",
# which all have an equal area.
function partition_indices(len)
    indices = Vector{UnitRange{Int}}(undef, Threads.nthreads())
    length = len
    offset = 1

    # Here, we have to iterate in reverse, to get an equal area every time
    for n in Threads.nthreads():-1:1
        partition = round(Int, length / sqrt(n)) # length varies as square root of area
        indices[n] = offset:(partition .+ offset .- 1)
        length -= partition
        offset += partition
    end

    return indices
end

# And now for the actual implementation:
function _distancematrix(x::Vector{T}, metric::Metric, ::Val{true}) where T
    d = zeros(T, length(x), length(x))

    Threads.@threads for k in partition_indices(length(x))
        for j in k
            for i in 1:j
                @inbounds d[i, j] = abs(x[i] - x[j])
            end
        end
    end

    return Symmetric(d, :U)

end

function _distancematrix(x::AbstractDataset{S, T}, metric::Metric, ::Val{true}) where T where S
    d = zeros(T, length(x), length(x))

    Threads.@threads for k in partition_indices(length(x))
        for j in k
            for i in 1:j
                @inbounds d[i, j] = evaluate(metric, x[i], x[j])
            end
        end
    end

    return Symmetric(d, :U)

end


text:
using Aqua
using DirectedAcyclicGraphs
using Test

@testset "Aqua tests" begin
    Aqua.test_all(DirectedAcyclicGraphs,
                    ambiguities = false)
end


text:
# Even Fibonacci numbers

LIMIT = 4_000_000

function solve()
    a, b = 1, 2
    answer = 0
    while a < LIMIT
        if a % 2 == 0
            answer += a
        end
        a, b = b, a + b
    end
    return answer
end

print(Int(solve()))

max_stars_repo_name: wephy/project-euler

max_issues_repo_name: wephy/project-euler

max_forks_repo_name: wephy/project-euler


text:
abstract type AbstractNode{T} end

mutable struct LinkNode{T} <: AbstractNode{T}
    next::LinkNode{T}
    data::T
    LinkNode{T}(next, data) where T = new(next, data)
    LinkNode{T}() where T = (n = new(); n.next = n; n)
    LinkNode{T}(data) where T = (n = new(); n.data = data; n.next = n; n)
end

#LinkNode(next, data::T) where T = LinkNode{T}(next, data)
LinkNode(data::T) where T = LinkNode{T}(data)
LinkNode(next::LinkNode{T}, data) where T = LinkNode{T}(next, data)

mutable struct LinkedList{T} <: AbstractList{T}
    head::LinkNode{T}
    tail::LinkNode{T}
    curr::LinkNode{T}
    cnt::Int
    function LinkedList{T}() where T
        head = tail = curr = LinkNode{T}()
        cnt = 0
        new(head, tail, curr, cnt)
    end
end

isempty(l::LinkedList) = l.cnt == 0
length(l::LinkedList) = l.cnt

function insert!(l::LinkedList, item)
    l.curr.next= LinkNode(l.curr.next, item)
    if l.tail === l.curr
        l.tail = l.curr.next
    end
    l.cnt += 1
end

function append!(l::LinkedList, item)
    l.tail = l.tail.next = LinkNode{eltype(l)}(item)
    l.cnt += 1
end

function remove!(l::LinkedList)
    (l.curr === l.tail) && throw(ArgumentError("no current element"))
    item = l.curr.next.data
    if l.curr.next === l.tail
        l.tail = l.curr
    end
    l.curr.next = l.curr.next.next
    l.cnt -= 1
    return item
end

function movetostart!(l::LinkedList)
    l.curr = l.head
end

function movetoend!(l::LinkedList)
    l.curr = l.tail
end

function movetoposition!(l::LinkedList, pos)
    (1 <= pos <= l.cnt+1) || throw(BoundsError(l, pos))
    l.curr = l.head
    for i = 1:pos-1
       l.curr = l.curr.next
    end
end

function prev!(l::LinkedList)
    (l.curr === l.head) && return
    tmp = l.head
    while l.curr !== tmp.next
        tmp = tmp.next
    end
    l.curr = tmp
end

function next!(l::LinkedList)
    if l.curr !== l.tail
        l.curr = l.curr.next
    end
end

function currentposition(l::LinkedList)
    tmp = l.head
    i = 1
    while l.curr !== tmp
        tmp = tmp.next
        i += 1
    end
    return i
end

function getvalue(l::LinkedList)
    (l.curr === l.tail) && throw(ArgumentError("no current element"))
    return l.curr.next.data
end


text:
module Components
import JSON2
export Component, ComponentContainer, <|

abstract type AbstractComponent end

struct Component <: AbstractComponent
    type ::String
    namespace ::String
    props ::Dict{Symbol, Any}
    available_props ::Set{Any}
end

JSON2.@format Component begin
    #=type => (name = "type")
    namespace => (name = "namespace")
    props => (name = "props")=#
    available_props => (exclude = true,)
end

function <|(cont::Component, value::Any)
    cont.props["children"] = value
    return cont
end

function is_valid_idprop(comp::Component, idprop::Tuple{Symbol, Symbol})::Bool
    if haskey(comp.props, :id) && Symbol(comp.props[:id]) == idprop[1]
        return idprop[2] in comp.available_props
    end

    if haskey(comp.props, :children)
        if comp.props[:children] isa Vector || comp.props[:children] isa Tuple
            for child in comp.props[:children]
                if child isa Component && is_valid_idprop(child, idprop)
                    return true
                end
            end
        elseif comp.props[:children] isa Component
            return is_valid_idprop(comp.props[:children], idprop)
        end
    end
    return false
end
end

max_stars_repo_name: EricForgy/Dashboards

max_issues_repo_name: EricForgy/Dashboards

max_forks_repo_name: EricForgy/Dashboards


text:
module fuzzynum

############ Fuzzy Bool and Fuzzy Float ############

import Base: +, *, zero

primitive type FuzzyBool <: Integer 8 end

FuzzyBool(x :: Bool) = reinterpret(FuzzyBool, x)
Bool(x :: FuzzyBool) = reinterpret(Bool, x)
Base.show(io :: IO, x :: FuzzyBool) = print(io, Bool(x))
+(a :: FuzzyBool, b :: FuzzyBool) = FuzzyBool(Bool(a) || Bool(b))
*(a :: FuzzyBool, b :: FuzzyBool) = FuzzyBool(Bool(a) && Bool(b))
zero(a :: FuzzyBool) = reinterpret(FuzzyBool, false)
zero(a :: Type{FuzzyBool}) = reinterpret(FuzzyBool, false)

primitive type FuzzyFloat <: AbstractFloat 64 end

FuzzyFloat(x :: Float64) = reinterpret(FuzzyFloat, x)
Float64(x :: FuzzyFloat) = reinterpret(Float64, x)
Base.show(io :: IO, x :: FuzzyFloat) = print(io, Float64(x))
+(a :: FuzzyFloat, b :: FuzzyFloat) = FuzzyFloat(max(Float64(a), Float64(b)))
*(a :: FuzzyFloat, b :: FuzzyFloat) = FuzzyFloat(min(Float64(a), Float64(b)))
zero(a :: FuzzyFloat) = reinterpret(FuzzyFloat, 0.0)
zero(a :: Type{FuzzyFloat}) = reinterpret(FuzzyFloat, 0.0)

export FuzzyBool, FuzzyFloat, +, *, zero

############ Array Utilities ############

function setToList(s)
    dupS = copy(s)
    n = length(dupS)
    xs = []
    for i in 1:n
        push!(xs, pop!(dupS))
    end
    xs
end

function getIndexInArr(xs, k)
    len = length(xs)
    for i in 1:len
        if k == xs[i]
            return i
        end
    end
    return -1
end

# https://stackoverflow.com/a/36518149
function fastuniq(v)
  v1 = Vector{eltype(v)}()
  if length(v)>0
    laste = v[1]
    push!(v1,laste)
    for e in v
      if e != laste
        laste = e
        push!(v1,laste)
      end
    end
  end
  return v1
end

export setToList, getIndexInArr, fastuniq

############ Ghehlien (Plain Binary Relation) ############

function getTransitiveClosure(FuzzyArr)
    while true
        arr_new = FuzzyArr * FuzzyArr
        if arr_new == FuzzyArr
            return arr_new
        end
        FuzzyArr = arr_new
    end
end

function getCluster(arr, n)
    arr_new = map(x -> typeof(x) == FuzzyFloat ? Float64(x) : Bool(x), arr)

    ss = []
    for i in 1:n
        if all(xs -> !in(i, xs), ss)
            s = map(x -> x[2], Iterators.filter(x -> x[1] != zero(x[1]), zip(view(arr_new, :, i), collect(1:n))))
            if length(s) != 0
                push!(ss, s)
            end
        end
    end
    ss
end

function ghehlienAsArray(xs)
    xsLen = length(xs)
    xsUniqueSet = Set(map(x -> x[1], xs))
    union!(xsUniqueSet, Set(map(x -> x[2], xs)))
    xsUniqueLen = length(xsUniqueSet)
    xsArr = eye(Bool, xsUniqueLen)
    xsUniqueList = setToList(xsUniqueSet)
    for i in 1:xsLen
        i1 = getIndexInArr(xsUniqueList, xs[i][1])
        i2 = getIndexInArr(xsUniqueList, xs[i][2])
        xsArr[i1, i2] = true
        xsArr[i2, i1] = true
    end
    fuzzyArr = map(x -> FuzzyBool(x), xsArr)
    fuzzyTransitiveClosure = getTransitiveClosure(fuzzyArr)
    ss = getCluster(fuzzyTransitiveClosure, xsUniqueLen)
    res = []
    for i in map(m -> join(map(n -> xsUniqueList[n], m)), ss)
        push!(res, i)
    end
    res
end

function ghehlien(xs)
    for i in ghehlienAsArray(xs)
        println(i)
    end
end

export ghehlienAsArray, ghehlien

end

max_stars_repo_name: gaozhiyan/Ghehlien

max_issues_repo_name: gaozhiyan/Ghehlien

max_forks_repo_name: gaozhiyan/Ghehlien


text:
Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    csvfileAbsname =joinpath(@__DIR__ ,"csvfile.csv" )
    println(csvfileAbsname )
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/subsolver_config/config.jl
# Contains subroutines used to set default configuration for select supported
# solvers along with routines needed to adjust tolerances to mirror tolerance
# adjustments in the global solve.
#############################################################################

function set_default_config_udf!(ext::ExtensionType, m::MOI.AbstractOptimizer) end

include("config/cbc.jl")
include("config/clp.jl")
include("config/cosmo.jl")
include("config/cplex.jl")
include("config/ecos.jl")
include("config/glpk.jl")
include("config/gurobi.jl")
include("config/ipopt.jl")
include("config/knitro.jl")
include("config/mosek.jl")
include("config/scs.jl")
include("config/tulip.jl")

function set_default_config!(ext::DefaultExt, m::T) where {T <: MOI.AbstractOptimizer}
    solver_name = MOI.get(m, MOI.SolverName())
    if solver_name == "Ipopt"
        set_default_config!(Val{:ipopt}(), ext, m)
    elseif solver_name == "GLPK"
        set_default_config!(Val{:glpk}(), ext, m)
    elseif solver_name == "CPLEX"
        set_default_config!(Val{:cplex}(), ext, m)
    elseif solver_name == "Clp"
        set_default_config!(Val{:clp}(), ext, m)
    elseif solver_name == "COIN Branch-and-Cut (Cbc)"
        set_default_config!(Val{:cbc}(), ext, m)
    elseif solver_name == "SCS"
        set_default_config!(Val{:scs}(), ext, m)
    elseif solver_name == "Mosek"
        set_default_config!(Val{:mosek}(), ext, m)
    elseif solver_name == "ECOS"
        set_default_config!(Val{:ecos}(), ext, m)
    elseif solver_name == "COSMO"
        set_default_config!(Val{:cosmo}(), ext, m)
    elseif solver_name == "Tulip"
        set_default_config!(Val{:tulip}(), ext, m)
    elseif solver_name == "Gurobi"
        set_default_config!(Val{:gurobi}(), ext, m)
    elseif solver_name == "Knitro"
        set_default_config!(Val{:knitro}(), ext, m)
    else
        set_default_config_udf!(ext, m)
    end
    return nothing
end

function set_default_config!(ext::ExtensionType, m::T) where {T <: MOI.AbstractOptimizer}
    set_default_config!(DefaultExt(), m::T)
end


text:
@testset "EmbedOnehot" begin
    using Flux: data
    multi(x, n) = collect(Iterators.repeated(x, n))

    after_enc = [2,2,4,5,6,1,1]
    multi_after_enc = hcat(multi(after_enc, 5)...)

    e = Embed(10, 6)

    @test data(e(after_enc)) == data(hcat(map(i->e.embedding[:,i], after_enc)...))
    @test eltype(data(e(after_enc, 0.5))) == Float32

    @test size(e) == (10, 6)
    @test size(e, 1) == 10
    @test size(e, 2) == 6
end


text:
import .ExpressionExplorer
import .ExpressionExplorer: join_funcname_parts

"Update the cell's caches, i.e. parse code and collect metadata."
function update_caches!(notebook::Notebook, cells)
    for cell in cells
        if cell.parsedcode === nothing
            cell.parsedcode = parse_custom(notebook, cell)
			cell.module_usings = ExpressionExplorer.compute_usings(cell.parsedcode)
			cell.rootassignee = ExpressionExplorer.get_rootassignee(cell.parsedcode)
        end
	end
end

"Return a copy of `old_topology`, but with recomputed results from `cells` taken into account."
function updated_topology(old_topology::NotebookTopology, notebook::Notebook, cells)
	updated_symstates = Dict(cell => ExpressionExplorer.try_compute_symbolreferences(cell.parsedcode) for cell in cells)
	new_symstates = merge(old_topology.symstates, updated_symstates)

	# Update the combined collection of function definitions, where multiple specialisations of a function are combined into a single `SymbolsState`.
	new_funcdefs = union((symstate.funcdefs for (k, symstate) in new_symstates)...)
	new_topology = NotebookTopology(new_symstates, new_funcdefs)

	for cell in cells
		finish_cache!(new_topology, cell)
	end

	new_topology
end

"Account for globals referenced in function calls by including `SymbolsState`s from called functions in the cell itself."
function finish_cache!(topology::NotebookTopology, cell::Cell)
	calls = all_indirect_calls(topology, topology[cell])
	calls = union!(calls, keys(topology[cell].funcdefs)) # _assume_ that all defined functions are called inside the cell to trigger eager reactivity.
	filter!(in(keys(topology.combined_funcdefs)), calls)

	union!(topology[cell].references, (topology.combined_funcdefs[func].references for func in calls)...)
	union!(topology[cell].assignments, (topology.combined_funcdefs[func].assignments for func in calls)...)

	add_funcnames!(topology, cell, calls)
end

"""Add method calls and definitions as symbol references and definition, resp.

Will add `Module.func` (stored as `Symbol[:Module, :func]`) as Symbol("Module.func") (which is not the same as the expression `:(Module.func)`)."""
function add_funcnames!(topology::NotebookTopology, cell::Cell, calls::Set{Vector{Symbol}})
	push!(topology[cell].references, (topology[cell].funccalls .|> join_funcname_parts)...)
	push!(topology[cell].assignments, (keys(topology[cell].funcdefs) .|> join_funcname_parts)...)

	union!(topology[cell].references, (topology.combined_funcdefs[func].funccalls .|> join_funcname_parts for func in calls)...)
	union!(topology[cell].assignments, (keys(topology.combined_funcdefs[func].funcdefs) .|> join_funcname_parts for func in calls)...)
end


text:
"""
    setline(linewidth)

Set the line width and multiply it with the current multiplier
which is i.e. set by [`appear`](@ref) and [`disappear`](@ref).

Normal behavior without any animation is the same as `Luxor.setline`.

# Example
```
setline(10)
line(O, Point(10, 10))
```

# Arguments:
- `linewidth`: the line width in pixel
"""
function setline(linewidth)
    cs = get_current_setting()
    cs.line_width = linewidth
    current_line_width = cs.line_width * cs.mul_line_width
    Luxor.setline(current_line_width)
end

"""
    setopacity(opacity)

Set the opacity and multiply it with the current multiplier
which is i.e. set by [`appear`](@ref) and [`disappear`](@ref).

Normal behavior without any animation is the same as `Luxor.setopacity`.

# Example
```
setopacity(0.5)
circle(O, 20, :fill)
```

# Arguments:
- `opacity`: the opacity between 0.0 and 1.0
"""
function setopacity(opacity)
    cs = get_current_setting()
    cs.opacity = opacity
    current_opacity = cs.opacity * cs.mul_opacity
    Luxor.setopacity(current_opacity)
end

"""
    fontsize(fsize)

Same as `Luxor.fontsize`: Sets the current font size.

# Example
```
fontsize(12)
text("Hello World!")
```

# Arguments:
- `fsize`: the new font size
"""
function fontsize(fsize)
    cs = get_current_setting()
    cs.fontsize = fsize
    Luxor.fontsize(fsize)
end

"""
    get_fontsize(fsize)

Same as `Luxor.get_fontsize` but works with every version of Luxor that is supported by Javis.

# Example
```
fontsize(12)
fsize = get_fontsize()
text("Hello World! \$fsize")
```

# Returns
- `Float64`: the current font size
"""
function get_fontsize()
    cs = get_current_setting()
    return cs.fontsize
end

"""
    scale(scl)

Set the scale and multiply it with the current multiplier
which is i.e. set by [`appear`](@ref) and [`disappear`](@ref).

Normal behavior without any animation is the same as `Luxor.scale`.

# Example
```
scale(0.5)
circle(O, 20, :fill) # the radius would be 10 because of the scaling
```

# Arguments:
- `scl`: the new default scale
"""
function scale(scl::Number)
    scale(scl, scl)
end

"""
    scale(scl_x, scl_y)

Same as [`scale`](@ref) but the x scale and y scale can be changed independently.

# Arguments:
- `scl_x`: scale in x direction
- `scl_y`: scale in y direction
"""
function scale(scl_x, scl_y)
    cs = get_current_setting()
    cs.desired_scale = (scl_x, scl_y)
    current_scale = cs.desired_scale .* cs.mul_scale
    Luxor.scale(current_scale...)
    cs.current_scale = cs.current_scale .* current_scale
    # println("cs.current_scale: $(cs.current_scale)")
end

"""
    scaleto(x, y)

Scale to a specific scaling instead of multiplying it with the current scale.
For scaling on top of the current scale have a look at [`scale`](@ref).
"""
function scaleto(x, y)
    cs = get_current_setting()
    cs.desired_scale = (x, y)
    scaling = (x, y) ./ cs.current_scale
    # we divided by 0 but clearly we want to scale to 0
    # -> we want scaling to be 0 not Inf
    if x ≈ 0
        scaling = (0.0, scaling[2])
    end
    if y ≈ 0
        scaling = (scaling[1], 0.0)
    end
    Luxor.scale(scaling...)
    cs.current_scale = (x, y)
end


text:
using LLRBTrees
using FactCheck

include("randomscenarios.jl")

facts("Inflate and deflate 5 elements") do

    tree = LLRBTree{Int, ASCIIString}()
    push!(tree, 1, "a")
    push!(tree, 2, "b")
    push!(tree, 3, "c")
    push!(tree, 4, "d")
    push!(tree, 5, "e")

    context("Integrity") do
        @fact tree --> check_height "Not balanced"
        @fact getmaxdepth(tree) == 3 --> true "Max depth does not match"
        @fact maximum(tree).key.value == 5 --> true "Maximun does not match"

        listed = orderedpairs(tree)
        i=1
        flag = true
        for pair = listed
            flag = pair[1] == i ? flag : false
            i += 1
        end
        @fact flag --> true "Didn't order correctly"
    end

    context("Structure") do

        flag = true
        i = 1
        for letter = ["a","b","c","d","e"]
            flag = letter == tree[i].value ? flag : false
            i+=1
        end
        @fact flag --> true "Values do not correspond to keys"

        tree[1]="z"
        @fact tree[1].value == "z" --> true "Didn't assign new value to key 1"
        tree[1]="a"
        @fact tree[1].value == "a" --> true "Didn't assign new value to key 1"

        flag = true
        for i in 1:5
            flag = haskey(tree, i) ? flag : false
        end
        @fact flag --> true "All keys are present"

        @fact tree.root.left.key.value == 2 --> true "Tree structure differs from expected in inflation"
        @fact tree.root.right.key.value == 5 --> true "Tree structure differs from expected in inflation"
        @fact tree.root.left.right.key.value == 3 --> true "Tree structure differs from expected in inflation"
        @fact tree.root.left.left.key.value == 1 --> true "Tree structure differs from expected in inflation"
    end

    context("Deflation") do
        delete!(tree,2)
        @fact tree.root.left.left.key.value == 1 --> true "Tree structure differs from expected in deflation"

        delete!(tree, 4)
        @fact tree.root.isRed --> false "Tree structure differs from expected in deflation"

        delete!(tree, 1)
        delete!(tree, 3)
        delete!(tree, 5)
        @fact isleaf(tree.root) --> true "Tree wasn't empty at the end of deflation"
    end
end

randomscenarios()

FactCheck.exitstatus()


text:
# Fix of kron to work with Convex.jl in seesawAlice
# This is a fix of the bug in the kron extension
# of the latest registered version of Convex for Julia v0.5.
# The bug does not allow to apply kron to a variable and fixed matrix.
# To solve that we define a new function "simplekron" which deals with this case.

using Convex

export simplekron

function simplekron(a::Union{AbstractArray, Convex.Constant, Convex.AbstractExpr}, b::Union{AbstractArray, Convex.Constant, Convex.AbstractExpr})
	Rs = AbstractExpr[]
	for i in 1:size(a)[1]
		Vs = Convex.AbstractExpr[]
		for j in 1:size(a)[2]
			push!(Vs, a[i, j] * b)
		end
		push!(Rs, foldl(hcat, Vs))
	end
	return foldl(vcat, Rs)
end

function simplekron(a::AbstractMatrix, b::Convex.AbstractExpr)
	return simplekron(Constant(a), b)
end

function simplekron(a::AbstractExpr, b::AbstractMatrix)
	return simplekron(a, Constant(b))
end


text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = test_page_model()
    include("../src/components/LGcycle.jl")

    add_comp!(m, LGcycle)

    set_param!(m, :LGcycle, :e_globalLGemissions, readpagedata(m, "test/validationdata/$valdir/e_globalLGemissions.csv"))
    set_param!(m, :LGcycle, :y_year, [2020.,2030.,2040.,2050.,2075.,2100.,2150.,2200.,2250.,2300.]) # real value
    set_param!(m, :LGcycle, :y_year_0, 2015.) # real value
    set_param!(m, :LGcycle, :rtl_g_landtemperature, readpagedata(m, "test/validationdata/$valdir/rtl_g_landtemperature.csv"))

    p = load_parameters(m)
    set_leftover_params!(m, p) # important for setting left over component values

    # run Model
    run(m)

    conc = m[:LGcycle,  :c_LGconcentration]
    conc_compare = readpagedata(m, "test/validationdata/$valdir/c_LGconcentration.csv")

    @test conc ≈ conc_compare rtol = 1e-4
end


text:
using DrWatson
@quickactivate "FractalDimension" # uses DynamicalSystems, PyPlot
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot
using ProgressMeter, Statistics
using JLD2

include("judd.jl")

sample_size = 100
realizations = 100
deg_a = 1
F = 8.

datasets[3] = :lorenz96_chaotic
labels[3] = "Lorenz96 (\$D\$=8)"

datasets[1] = :henon_chaotic
labels[1] = "Hénon map"

datasets[2] = :lorenz_chaotic
labels[2] = "Lorenz system"

function run_exs(sample_size, realizations, F, deg_a)

    name = savename("judd_sensitivity", @dict(sample_size, realizations,deg_a,F),"jld2")

    isfile(datadir("judd", name)) && return

    systems = ["Henon", "Lorenz", "L96-D=5","L96-D=6","L96-D=7","L96-D=8","L96-D=9"]
    systs = [Systems.henon(), Systems.lorenz(),Systems.lorenz96(5, F=F),Systems.lorenz96(6, F=F),Systems.lorenz96(7, F=F),Systems.lorenz96(8, F=F),Systems.lorenz96(9, F=F)]

    Deltas_ic = []

    println("Estimation of sensitivity on initial conditions with fixed sample")
    for (system, ds) in zip(systems, systs)

        println("System: $system")

        ts = regularize(trajectory(ds,Ttr = 100, 5000))
        x = interpoint_distances(ts, sample_size)
        εs, bs, ε₀ = bin_judd(x)

        Deltas_initial_cond = zeros(realizations)
        @showprogress for k in 1:realizations
            d₀, a₀ = rand()*dimension(ds)*1.5, rand(deg_a+1)*500
            d, a = nested_optimizer(εs, bs; ε₀, d₀, a₀, η = 1e-4)
            Deltas_initial_cond[k] = d
        end

        println("For 100 different random initial conditions, the fractal dimension of $system is estimated as $(mean(Deltas_initial_cond))±$(std(Deltas_initial_cond))")
        push!(Deltas_ic, Deltas_initial_cond)
    end

    Deltas_s = []

    println("Estimation of sensitivity on sample with fixed initial conditions")

    for (system, ds) in zip(systems, systs)

        ts = regularize(trajectory(ds,Ttr = 100, 5000))

        d₀, a₀ = rand()*dimension(ds)*1.5, rand(deg_a+1)*500
        Deltas_samples = zeros(realizations)
        @showprogress for k in 1:realizations
            x = interpoint_distances(ts, sample_size)
            εs, bs, ε₀ = bin_judd(x)
            d, a = nested_optimizer(εs, bs; ε₀, d₀, a₀, η = 1e-4)
            Deltas_samples[k] = d
        end

        println("For 100 different random samples, the fractal dimension of $system is estimated as $(mean(Deltas_samples))±$(std(Deltas_samples))")
        push!(Deltas_s, Deltas_samples)
    end

    result = @strdict systems Deltas_ic Deltas_s
    @tagsave(datadir("judd", name), result)
end

function plot_exs(sample_size,realizations,F,deg_a)
    name = savename("judd_sensitivity", @dict(sample_size, realizations,deg_a,F),"jld2")
    plotname = savename("judd_sensitivity", @dict(sample_size, realizations,deg_a,F),"png")

    # isfile(plotsdir("judd", plotname)) && return

    @load datadir("judd", name) systems Deltas_ic Deltas_s

    popat!(Deltas_s,findall(x->x=="Lorenz",systems)[1])
    popat!(systems,findall(x->x=="Lorenz",systems)[1])

    fig = figure()

    pos = collect(1:length(Deltas_s))
    fp = Dict("marker"=>"o", "markerfacecolor"=>COLORS[1], "markersize"=>6,
                  "linestyle"=>LINESTYLES[1], "markeredgecolor"=>COLORS[1])
    box1 = boxplot(Deltas_s, notch = true, positions = pos, flierprops=fp,patch_artist=true)
    for item in ["boxes", "whiskers", "fliers", "medians", "caps"]
        setp(box1[item], color=COLORS[1])
    end

    xticks(pos, systems, rotation=45)
    xlabel("System")
    ylabel(L"$ \Delta^{(J)}$")
    tight_layout()
    wsave(plotsdir("judd", plotname), fig)
end

run_exs(sample_size,realizations,F,deg_a)

plot_exs(sample_size,realizations,F,deg_a)

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
using Random

function makesim(sts::Union{STS{:Barkley}, STS{:bk}})
    @unpack N, T, Δt, S, ic = sts
    Nx, Ny = typeof(N) <: Tuple ? (N[1], N[2]) : (N, N)
    FType = get(sts.p, :FType, Float32)
    # Assertions
    @assert Δt/BK_INTEG_DT == round(Int, Δt/BK_INTEG_DT) "Δt must be multiple of the integration step (0.01)"
    every = round(Int, Δt/BK_INTEG_DT)
    @assert Nx ≥ 40
    @assert Ny ≥ 40
    @assert Nx%10==0
    @assert Ny%10==0
    # prepare initial condition...
    if ic isa Union{Nothing, Integer}
        rng = MersenneTwister(isnothing(ic) ? 42 : ic)
        αα, ββ = 10, 10
        init = rand(rng,αα,ββ)
        u = zeros(FType,Nx,Ny)
        v = zeros(FType,Nx,Ny)
        v .= u .= repeat(init, inner=(Nx÷αα,Ny÷ββ))
        v .= v .> 0.2
    elseif ic isa Tuple
        u, v = ic
    else
        error(ICERROR)
    end
    t, U, V = barkley(u, v, S, T, Δt; sts.p...)
    return @strdict t sts U V
end

"""
    barkley(U0, V0, S, T, Δt; kwargs...) → t, U, V
Simulate the nonlinear Barkley model with initial conditions `U0, V0` for
the fields `U, V`.

## Keywords
```
periodic = true, FType = Float32,
skip=100, a=0.75, b=0.06, ε=0.08, D=1/50, h=0.1
```
The simulation uses `periodic` boundary conditions by default, or absorbing otherwise.
`skip` is an integer, multiple of `Δt`, which notes how many
`Δt` frames to skip before starting to save data. `FType` is the
number type of the output fields.
`a,b,ε,D` are formal parameters of the system while `h` is the spatial sampling
(the evolution algorithm is the one from the Scholarpedia entry:
http://www.scholarpedia.org/article/Barkley_model )
(notice that the integration time is fixed to 0.01, irrespectively of `Δt`).

`a` can also be a `Matrix` (of size `Nx×Ny`), to e.g. represent inhomogeneity.
"""
function barkley(u::Matrix, v::Matrix, S, T, Δt; periodic=true,
                 skip=100, a=0.75, b=0.06, ε=0.08, D=1/50, h=0.1,
                 )

    FType = eltype(u)
    bk_integ_dt = FType(BK_INTEG_DT)
    total = Int(T÷Δt) # how many steps to save
    Nx, Ny = size(u)
    @assert size(u) == size(v)
    u, v = copy(u), copy(v)
    U = Matrix{FType}[]
    V = Matrix{FType}[]
    sizehint!(U, total); sizehint!(V, total)

    Σ = zeros(FType, Nx, Ny, 2)
    r, s = 1, 2
    every = Δt÷bk_integ_dt

    function F(u, uth)
        if u < uth
            u/(1-(bk_integ_dt/ε)*(1-u)*(u-uth))
        else
            (u + (bk_integ_dt/ε)*u*(u-uth))/(1+(bk_integ_dt/ε)*u*(u-uth))
        end
    end

    function periodic_step!(u,v,Σ,a,b,Nx,Ny,s,r,D,bk_integ_dt,h)
        @inbounds for i=1:Nx, j=1:Ny
            uth = (v[i,j] + b)/_get_a(a, i, j)
            v[i,j] = v[i,j] + bk_integ_dt*(u[i,j]^3 - v[i,j])
            u[i,j] = F(u[i,j], uth) + D*bk_integ_dt/h^2 *Σ[i,j,r]
            Σ[i,j,s] -= 4u[i,j]
            Σ[  mod(i-1-1,Nx)+1,j,s] += u[i,j]
            Σ[  mod(i+1-1,Nx)+1,j,s] += u[i,j]
            Σ[i,mod(j-1-1,Ny)+1,  s] += u[i,j]
            Σ[i,mod(j+1-1,Ny)+1,  s] += u[i,j]
            Σ[i,j,r] = 0
        end
    end
    function constant_step!(u,v,Σ,a,b,Nx,Ny,s,r,D,bk_integ_dt,h)
        @inbounds for i=1:Nx, j=1:Ny
            uth = (v[i,j] + b)/_get_a(a, i, j)
            v[i,j] = v[i,j] + bk_integ_dt*(u[i,j]^3 - v[i,j])
            u[i,j] = F(u[i,j], uth) + D*bk_integ_dt/h^2 *Σ[i,j,r]
            Σ[i,j,s] -= 4*u[i,j]
            i > 1  && (Σ[i-1,j,s] += u[i,j])
            i < Nx && (Σ[i+1,j,s] += u[i,j])
            j > 1  && (Σ[i,j-1,s] += u[i,j])
            j < Ny && (Σ[i,j+1,s] += u[i,j])
            Σ[i,j,r] = 0
        end
    end

    # evolve for skip time
    for j in 1:S÷bk_integ_dt
        periodic ? periodic_step!(u,v,Σ,a,b,Nx,Ny,s,r,D,bk_integ_dt,h) :
                   constant_step!(u,v,Σ,a,b,Nx,Ny,s,r,D,bk_integ_dt,h)
        r,s = s,r
    end
    push!(U, copy(u)); push!(V, copy(v)) # actual initial condition
    for j in 1:total
        for i in 1:every
            periodic ? periodic_step!(u,v,Σ,a,b,Nx,Ny,s,r,D,bk_integ_dt,h) :
                       constant_step!(u,v,Σ,a,b,Nx,Ny,s,r,D,bk_integ_dt,h)
            r,s = s,r
        end
        push!(U, copy(u)); push!(V, copy(v))
    end
    return 0:Δt:T, U, V
end

_get_a(a::Number, i, j) = a
_get_a(a::AbstractMatrix, i, j) = @inbounds a[i,j]
const BK_INTEG_DT = 0.01


text:
using JLD

i = readdlm("i")
q = readdlm("q")
u = readdlm("u")
v = readdlm("v")
pa = readdlm("pa")

save(ARGS[1]*".jld", "i", i, "q", q, "u", u, "v", v, "pa", pa)

max_stars_repo_name: KajWiik/PolarimetryCalibration

max_issues_repo_name: KajWiik/PolarimetryCalibration

max_forks_repo_name: KajWiik/PolarimetryCalibration


text:
struct Lennnon2000Air <: ThermoState.ThermoModel end
const TAU_MAX_EXP_87 = 0.4207493606569795
const lemmon2000_air_R = 8.314510

const lemmon2000_air_T_reducing = 132.6312
const lemmon2000_air_P_reducing = 3.78502E6
const lemmon2000_air_rho_reducing = 10447.7
const lemmon2000_air_rho_reducing_inv = 1.0/lemmon2000_air_rho_reducing

const lemmon2000_air_MW = 28.9586
const lemmon2000_air_P_max = 2000E6
const lemmon2000_air_T_max = 2000.

molecular_weight(::Lennnon2000Air) = lemmon2000_air_MW

function _f0(::Lennnon2000Air, delta,tau)
    tau_inv = one(tau)/tau
    A0 =  (-0.00019536342*tau*sqrt(tau) + 17.275266575*tau + tau_inv*(tau_inv*(6.057194e-8*tau_inv
            - 0.0000210274769) - 0.000158860716) + log(delta) + 2.490888032*log(tau)

            # These two logs both fail for tau < 1e-18, can be truncated but should not be necessary.
            + 0.791309509*log(1.0 - exp(-25.36365*tau)) + 0.212236768*log(1.0 - exp(-16.90741*tau))
            - 13.841928076)

            if tau < TAU_MAX_EXP_87
        A0 -= 0.197938904*log(exp(87.31279*tau) + (2.0/3.0))
    else
        A0 -= 17.282597957782162*tau # 17.282... = 87.31279*0.197938904
    return A0
    end
end

function _fr(::Lennnon2000Air,delta,tau)

    delta2 = delta*delta
    delta3 = delta*delta2
    delta4 = delta2*delta2
    delta5 = delta*delta4
    delta6 = delta2*delta4
    taurt2 = sqrt(tau)
    taurt4 = sqrt(taurt2)
    tau2 = tau*tau
    tau3 = tau*tau2
    tau6 = tau3*tau3
    tau12 = tau6*tau6
    tau_100 = tau^0.01
    tau2_100 = tau_100*tau_100
    tau4_100 = tau2_100*tau2_100
    tau5_100 = tau_100*tau4_100
    tau10_100 = tau5_100*tau5_100
    tau15_100 = tau5_100*tau10_100
    tau8_100 = tau4_100*tau4_100
    tau16_100 = tau8_100*tau8_100
    tau20_100 = tau4_100*tau16_100
    tau32_100 = tau16_100*tau16_100
    tau33_100 = tau_100*tau32_100
    tau64_100 = tau32_100*tau32_100
    tau80_100 = tau16_100*tau64_100
    tau40_100 = tau20_100*tau20_100
    tau97_100 = tau33_100*tau64_100
    tau45_100 = tau5_100*tau40_100
    tau90_100 = tau45_100*tau45_100
    tau160_100 = tau80_100*tau80_100
    tau320_100 = tau160_100*tau160_100
    x0 = exp(-delta)
    x1 = exp(-delta2)
    x2 = tau3*exp(-delta3)
    return (-0.101365037911999994*delta*tau160_100*x0 + 0.713116392079000017*delta*tau33_100
            - 0.146629609712999986*delta*tau40_100*x1*tau320_100
            - 1.61824192067000006*delta*tau4_100*tau97_100 + 0.0148287891978000005*delta*taurt2*x2
            + 0.118160747228999996*delta + 0.0714140178971000017*delta2 + 0.134211176704000013*delta3*tau15_100
            - 0.031605587982100003*delta3*tau6*x1 - 0.17381369096999999*delta3*tau80_100*x0
            - 0.00938782884667000057*delta3*x2*tau12 - 0.0865421396646000041*delta3 - 0.042053322884200002*delta4*tau20_100
            + 0.0349008431981999989*delta4*tau2_100*tau33_100 + 0.0112626704218000001*delta4
            - 0.0472103183731000034*delta5*tau15_100*x0*tau80_100 + 0.000233594806141999996*delta5*tau3*taurt4*x1*delta6
            - 0.0122523554252999996*delta6*tau*taurt4*x0 + 0.000164957183186000006*delta6*tau45_100*tau90_100)
end

function αR_impl(mt::SingleVT,::Lennnon2000Air, _rho, T)
    R = lemmon2000_air_R
    delta = rho*lemmon2000_air_rho_reducing_inv
    tau = lemmon2000_air_T_reducing/T
    return _fr(model, delta, tau)
end

function α0_impl(mt::SingleVT,::Lennnon2000Air, _rho, T)
    R = lemmon2000_air_R
    delta = rho*lemmon2000_air_rho_reducing_inv
    tau = lemmon2000_air_T_reducing/T
    return _f0(model, delta, tau)
end

function mol_helmholtzR_impl(mt::SingleVT,::Lennnon2000Air, v, t)
    rho = 1.0e-3 / v
    R = lemmon2000_air_R
    delta = lemmon2000_air_rho_reducing_inv/v
    tau = lemmon2000_air_T_reducing/T
    return R*t*_fr(model, delta, tau)
end

function mol_helmholtz0_impl(mt::SingleVT,::Lennnon2000Air, v, t)
    rho = 1.0e-3 / v
    R = lemmon2000_air_R
    delta = lemmon2000_air_rho_reducing_inv/v
    tau = lemmon2000_air_T_reducing/T
    return R*t*_f0(model, delta, tau)
end

function mol_helmholtz_impl(mt::SingleVT,::Lennnon2000Air, v, t)
    rho = 1.0e-3 / v
    R = lemmon2000_air_R
    delta = lemmon2000_air_rho_reducing_inv/v
    tau = lemmon2000_air_T_reducing/T
    return R*t*(_f0(model, delta, tau)+_fr(model, delta, tau))
end


text:
function MV(path, forest; rounded = false, bg = "#ffffff", fc = "#000000", label = "Tree", fs = "18", nc="#ffffff") # Model Visualization
    if path[length(path) - 3 : length(path)] != ".dot"
        throw("The file path you passed is not prefixed with'dot'.")
    end
    open(path, "w") do file
        println(file, "digraph Tree_model{")
        println(file, "graph[")
        println(file, "label = \""*label*"\",")
        println(file, "bgcolor = \""*bg*"\",")
        println(file, "fontcolor = \""*fc*"\",")
        println(file, "fontsize = \""*fs*"\",")
        println(file, "style = \"filled\",")
        println(file, "margin = 0.2\n];")
        println(file, "node[")
        println(file, "shape = box,")
        if rounded
            println(file, "style = \"rounded, filled\",")
        end
        println(file, "fillcolor = \""*nc*"\"\n];")
        i = 1
        edges = []
        function tree_write(tree, i; return_i = false)
            s = []
            push!(s, (nothing, tree))
            while length(s) != 0
                old_name, v = pop!(s)
                name = "node"*string(i)
                println(file, name*" [")
                st = ""
                if v["left"] != nothing
                    st = "X["*string(v["feature_id"])*"]<"*string(v["threshold"])*"\n"
                end
                samples = sum(v["class_count"])
                st *= "samples="*string(samples)*"\n"
                st *= "value="*string(v["class_count"])
                println(file, "label = \""*st*"\"")
                println(file, "];")
                if old_name == nothing
                elseif old_name != "node1"
                    push!(edges, old_name*"->"*name*";")
                elseif name == "node2"
                    push!(edges, "node1->node2[label = \"True\"];")
                else
                    push!(edges, "node1->"*name*"[label = \"False\"];")
                end
                i += 1
                if v["left"] != nothing
                    push!(s, (name, v["right"]))
                    push!(s, (name, v["left"]))
                end
                if return_i
                    return i
                end
            end
        end
        if typeof(forest) == Dict{String, Any}
           tree_write(forest, i)
        else
            for tree in forest
                i = tree_write(tree, i, return_i = true)
            end
        end
        for edge in edges
            println(file, edge)
        end
        print(file, "}")
    end
end


text:
using Test
import ONNXRunTime as OX
import Ghost: play!, Tape, Input
import ONNX: ONNXCtx, push_call!, from_nnlib, from_onnx, save, load

function ort_run(path, ort_args...)
    model = OX.load_inference(OX.testdatapath(path))
    ort_inputs = Dict([OX.input_names(model)[i] => ort_args[i] for i=1:length(ort_args)])
    return model(ort_inputs)
end

function ort_test(tape::Tape, args...)
    mktemp() do path, _
        r1 = play!(tape, args...)
        save(path, tape)
        r2 = ort_run(path, from_nnlib.(args)...) |> values |> first |> from_onnx
        tape2 = load(path, args...; exec=true)
        r3 = tape2[tape2.result].val
        @test isapprox(r1, r2)
        @test isapprox(r1, r3)
        # for more flexibility we return the tape before saving and after loading
        return tape, tape2
    end
end

function ort_test(fn::Function, args...; kwargs...)
    tape = Tape(ONNXCtx())
    inp = [push!(tape, Input(arg)) for arg in args]
    res = push_call!(tape, fn, inp...; kwargs...)
    tape.result = res
    return ort_test(tape, args...)
end


text:
using VideoIO, Dates

"""
videowritelistener(;compression=0,overwrite=false, options=``)
Write video to file

Based on https://discourse.julialang.org/t/creating-a-video-from-a-stack-of-images/646/8
Reference for H.264: https://trac.ffmpeg.org/wiki/Encode/H.264#LosslessH.264
The range of the compression (CRF) scale is 0–51, where 0 is lossless, 23 is the default, and 51 is worst quality possible
"""
function videowritelistener(;compression=0,overwrite=false,threads=0)

    ow = overwrite ? `-y` : `-n`
    preset = (compression == 0) ? "ultrafast" : "medium"

    global sessionStat
    global camImageFrameBuffer

    while !sessionStat.terminate
        if length(camImageFrameBuffer) > 0
            sessionStat.savedframes = 0
            fname = string(Dates.format(now(), "yyyy-mm-dd_HHMMSS"),".mp4")
            h, w = size(camImageFrameBuffer[1])
            if (h % 2 != 0) || (w % 2 != 0)
                @info "Image dims for H264 video encoding need to be even ($h x $w)"
                sessionStat.recording = false
                camImageFrameBuffer = Vector{Array{UInt8}}(undef,0)
                sessionStat.bufferedframes = 0
            else
                fps = camSettings.acquisitionFramerate
                withenv("PATH" => VideoIO.libpath,
                    "LD_LIBRARY_PATH" => VideoIO.libpath,
                    "DYLD_LIBRARY_PATH" => VideoIO.libpath) do
                    open(`$(VideoIO.ffmpeg)
                        -loglevel warning
                        -threads $threads
                        $ow -f rawvideo -pix_fmt gray -s:v $(h)x$(w)
                        -r $fps -i pipe:0 -c:v libx264 -preset: $preset
                        -crf $compression -pix_fmt yuv422p
                        $fname`, "w") do out
                        while (length(camImageFrameBuffer) > 0 || sessionStat.recording) && !sessionStat.terminate
                            while length(camImageFrameBuffer) == 0 && sessionStat.recording
                                wait(Timer(0.01)) #wait without blocking during recording time
                            end
                            if length(camImageFrameBuffer) >0
                                write(out, camImageFrameBuffer[1])
                                deleteat!(camImageFrameBuffer, 1)
                                sessionStat.bufferedframes = length(camImageFrameBuffer)
                                sessionStat.savedframes += 1
                            end
                            yield()
                        end
                        filepath = joinpath(pwd(),fname)
                        @info "Video saved [$(sessionStat.savedframes) frames, $(fps) fps]: $(filepath)"
                    end
                end
            end
        end
        yield()
    end
end


text:
function _precompile_()
#    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(searchsortedfirst, (Array{Int64, 1}, Base.StepRange{Int64, Int64}, Int64, Int64,))
    precompile(searchsortedfirst, (Array{Int64, 1}, Int64, Int64, Int64,))
    precompile(tail, (RLEVector{Int64, Int64}, Int64,))
    precompile(searchsortedfirst, (Array{Int64, 1}, Array{Int64, 1}, Int64, Int64,))
    precompile(length, (RLEVector{Int64, Int64},))
    precompile(numruns, (Array{Int64, 1}, Array{Int64, 1},))
    precompile(getindex, (RLEVector{Int64, Int64}, Base.UnitRange{Int64},))
    precompile(ind2run, (RLEVector{Int64, Int64}, Base.UnitRange{Int64},))
    precompile(numruns, (Array{Int64, 1},))
    precompile(ree, (Array{Int64, 1}, Array{Int64, 1}, Int64,))
    precompile(ind2runcontext, (RLEVector{Int64, Int64}, Int64,))
    precompile(getindex, (RLEVector{Int64, Int64}, Base.StepRange{Int64, Int64},))
    precompile(inverse_ree, (Array{Int64, 1}, Array{Int64, 1},))
    precompile(ree, (Array{Int64, 1},))
    precompile(ind2runcontext, (RLEVector{Int64, Int64}, Base.UnitRange{Int64},))
    precompile(setindex!, (RLEVector{Int64, Int64}, Int64, Int64,))
    precompile(setindex!, (RLEVector{Int64, Int64}, Int64, Base.UnitRange{Int64},))
    precompile(Type, (Type{RLEVector}, Array{Int64, 1}, Array{Int64, 1},))
    precompile(setindex!, (RLEVector{Int64, Int64}, Array{Int64, 1}, Array{Bool, 1},))
    precompile(ind2run, (RLEVector{Int64, Int64}, Array{Int64, 1},))
    precompile(ind2run, (RLEVector{Int64, Int64}, Int64,))
    precompile(getindex, (RLEVector{Int64, Int64}, Array{Int64, 1},))
    precompile(ind2run, (RLEVector{Int64, Int64}, Base.StepRange{Int64, Int64},))
    precompile(Type, (Type{RLEVector}, Array{Int64, 1},))
end


text:
# process a julia file and add comments with the output of the file (STDOUT)
# after the expression the generated that output
using ArgParse
OUTPUT_COMMENTS = "#> "

function removeold(code::String)
	# get rid of old output comments
	# we have to be particularly careful to correct the
	# end of the file so new output comments dont go on the same
	# line as the last command
	# endsnewline = endswith(code, "\n")
	code=join(filter(x-> !beginswith(x, OUTPUT_COMMENTS), split(code, ['\n'])), "\n")
	return !endswith(code, "\n") ? string(code, "\n") : code
end

function getoutput(code::String)
	# evaluate code expression-by-expression and collect STDOUT
	# together with where to put it back.
	# variables from eval contaminate the function
	# so we use the hack __var to avoid variable mixup
	pos = 1
	const STDOUT_OLD = STDOUT
	output = (Int, String)[]
	ns = Module()
	begin
		while true
			statement, newpos = parse(code, pos, raise=false)
			statement.head == :error ? break: (pos = newpos)
			stdout, stdin = redirect_stdout()
			eval(ns, statement)
			redirect_stdout(STDOUT_OLD)
			close(stdin)
			result = readall(stdout)
			close(stdout)
			if length(result) > 0
				push!(output, (pos, result))
			end
		end
	end
	return output
end

function insert_output(code::String, output::Array{(Int64,String), 1})
	# insert the output from above back into the code after the
	# expression which created it
	function tocomment(text::String)
		stripped = strip(text, [' ', '\n', '\t'])
		return length(stripped) == 0 ?
			"\n" :
			string(OUTPUT_COMMENTS, replace(stripped, "\n", "\n$OUTPUT_COMMENTS"), "\n")
	end
	reverse!(output)
	for (pos, text) in output
		if pos > length(code)
			start = code
			extra = ""
		else
			start = code[1:(pos - 1)]
			extra = code[(pos):end]
		end
		comment = tocomment(text)
		code = string(start, comment, extra)
	end
	return code
end

function convertfile(fname::String)
	if !(endswith(fname, ".jl") || endswith(fname, ".JL"))
		println("## ERROR: the file \"$fname\" does not appear to be a julia file (i.e., it doesn't end with .jl)")
		return
	end
	const code = removeold(open(readall, fname, "r"))
	output = getoutput(code)
	newcode = insert_output(code, output)
	open(fname, "w") do f
		write(f, newcode)
	end
	println("--------------\nOutput comments added to \"$fname\"\n--------------")
end

function parseargs()
	if length(ARGS) == 0
		println("## ERROR: File name required as command line argument.")
		return
	end
	for fname in ARGS
		convertfile(fname)
	end
end

parseargs()

max_stars_repo_name: inkydragon/JuliaByExample

max_issues_repo_name: inkydragon/JuliaByExample

max_forks_repo_name: inkydragon/JuliaByExample


text:
using SafeTestsets
@safetestset "my f tests" begin include("f_tests.jl") end
@safetestset "my Operators tests" begin include("Operators_tests.jl") end
@safetestset "my Operators tests" begin include("Tools_tests.jl") end


text:
Base.strides(n::NamedArray) = Base.strides(n.array)

Base.unsafe_convert(p::Type{Ptr{T}}, n::NamedArray) where T = Base.unsafe_convert(p, n.array)

# # Based in LinearAlgebra/src/matmul.jl: *(A::AbstractArray{T,2}, B::AbstractArray{T,2})
# function Base.:*(A::NamedArray{TA, 2, ATA, NTA},
#                  B::NamedArray{TB, 2, ATB, NTB}) where {TA, TB, ATA, ATB, NTA, NTB}
#     TS = Base.promote_op(LinearAlgebra.matprod, TA, TB)
#     out_size = (size(A,1), size(B,2))
#     dump(out_size)
#     out = NamedArray(Array{TS}(undef, out_size),
#                      (A.dicts[1], B.dicts[2]),
#                      (A.dimnames[1], B.dimnames[2]))
#     mul!(out, A.array, B.array)
# end
#
# LinearAlgebra.qr(n::NamedArray, args...) = qr(n.array, args...)


text:
using CodecZlib,DataFrames, uCSV
using CSV
cd("C:\\Users\\Clinton Tepper\\Dropbox\\Projects\\Endowment Project\\data\\NCCS")
p="nccs_core_2015_pc_full.csv.gz"

b = open(GzipDecompressorStream, p) do f
  b = IOBuffer("$(readline(f))\n$(readline(f))")
  (b)
end

df = CSV.read(b,# strings=:raw,
  categorical=false,  silencewarnings=true, strict=false)
nCols = size(df,2)
dfNames = (s::Symbol->Symbol(lowercase(string(s)))).(names(df))

# NOTE: Below throws an error
@time begin
  b = open(GzipDecompressorStream, p) do iStream
    b = IOBuffer(read(iStream))
    (b)
  end

  df = CSV.File(b) |> DataFrame
end

#NOTE: Below throws an error
CSV.File("nccs_core_2015_pc_full.csv") |> DataFrame

#limit=209_715

#=local df::DataFrame

try
  df = CSV.File(b) |> DataFrame
catch
  chunksize = 100_000

  df = CSV.File(b, limit=100_000)
  rowsread = 100_000

  oldsize = -1
  newsize = size(df,1)
  while oldsize ≠ newsize
    oldsize=newsize

  end
end=#

types = Dict(i=>Union{String,Missing} for i ∈ 1:nCols)

#=types = Dict(i=>Vector for i ∈ 1:nCols)

b = open(GzipDecompressorStream, p) do iStream
  b = IOBuffer(read(iStream))
  (b)
end
uCSV.read(b, coltypes=types, allowmissing=true, escape='\\', quotes='\"',
  encondings=Dict(""=>missing))=#

#=@time begin
  df = DataFrame(load(File(format"CSV", p)))
  #DataFrame(load(File(format"CSV", p), colparsers=types))
end=#

#=open("nccs_core_2013_pc_full.csv.gz") do io
    load(Stream(format"CSV", GzipDecompressorStream(io))) |> DataFrame
end=#

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
@testset "NO_FIELDS" begin
    # Following doesn't work because of some deprecate_binding weirdness with not printing
    # @test (@test_deprecated NO_FIELDS) isa NoTangent
    # So just test it gives the old behavour
    @test NO_FIELDS isa NoTangent
end

@testset "extern" begin
    @test (@test_deprecated extern(@thunk(3))) == 3
    @test (@test_deprecated extern(@thunk(@thunk(3)))) == 3

    @test (@test_deprecated extern(Tangent{Foo}(x=2.0))) == (;x=2.0)
    @test (@test_deprecated extern(Tangent{Tuple{Float64,}}(2.0))) == (2.0,)
    @test (@test_deprecated extern(Tangent{Dict}(Dict(4 => 3)))) == Dict(4 => 3)

    # with differentials on the inside
    @test (@test_deprecated extern(Tangent{Foo}(x=@thunk(0+2.0)))) == (;x=2.0)
    @test (@test_deprecated extern(Tangent{Tuple{Float64,}}(@thunk(0+2.0)))) == (2.0,)
    @test (@test_deprecated extern(Tangent{Dict}(Dict(4 => @thunk(3))))) == Dict(4 => 3)

    z = ZeroTangent()
    @test (@test_deprecated extern(z)) === false

    # @test_throws doesn't play nice with `@test_deprecated` so have to be loud
    dne = NoTangent()
    @test_throws Exception extern(dne)
    ni = @not_implemented("no")
    @test_throws ChainRulesCore.NotImplementedException extern(ni)
end

@testset "Deprecated: calling thunks should call inner function" begin
    @test (@test_deprecated (@thunk(3))()) == 3
    @test (@test_deprecated (@thunk(@thunk(3)))()) isa Thunk
end

@testset "Deprecated: Inplacable Thunk argument order" begin
    @test (@test_deprecated InplaceableThunk(@thunk([1]), x->x.+=1)) isa InplaceableThunk
end

@testset "Deprecated: convert from Tangent" begin
    @test (@test_deprecated convert(NamedTuple, Tangent{Foo}(x=2.5))) == (; x=2.5)
    @test (@test_deprecated convert(Tuple, Tangent{Tuple{Float64,}}(2.0))) == (2.0,)
    @test (@test_deprecated convert(Dict, Tangent{Dict}(Dict(4 => 3)))) == Dict(4 => 3)
end


text:
# This file was generated, do not modify it. # hide
fp = fitted_params(mach_uni)
@show fp.coefs
@show fp.intercept


text:
mystring = read("file1", String)
open(io->write(io, mystring), "file2", "w")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
abstract type NormalizeMethod end

struct BoxCox <: NormalizeMethod end
struct YeoJohnson <: NormalizeMethod end


text:
module Sabertooth

using SerialPorts
# Motor commands
    const forward_m1 = 0x00
    const reverse_m1 = 0x01
    const forward_m2 = 0x04
    const reverse_m2 = 0x05
    const forward_mixed = 0x08
    const reverse_mixed = 0x09
    const right_mixed = 0x0A
    const left_mixed = 0x0B
    const ramp = 0x10

# Modify the serial port configs as needed
    const portname = "/dev/ttyACM0"
    const baudrate = 115200 #Default baudrate for Sabertooth (115200)
    const address = 128 # The address of the Sabertooth. By default, this is 128.
    const mask = 127

s = SerialPort(portname, baudrate) # Instantiate the serial port

# sendCommand(command, data): Core motor controller function, this sends an update out to the motor controller hardware, all other functions should call this function instead of rewriting.
    function sendCommand(command, data)
        checksum = ((address + command + data) & 127) # Calculate the required checksum with the bitwise operator for the bytesize (mask)
        newUpdate = [UInt8(address), UInt8(command), UInt8(data), UInt8(checksum)] # Create an array buffer of bytes to send to the motor controller hardware

        write(s, newUpdate) # Write to the motor controller hardware
    end

# drive(motor, speed): This controls one motor at a time
     function drive(motor, speed)
        # Reverse commands are equal to forward + 1
        cmds = [forward_m1, forward_m2]

        try
            command = cmds[motor]
        catch
            @warn "Drive error, invalid motor number: $motor"
        end

        if speed < 0
            speed = -speed
            command = command + 1
        end
        if speed > 100
            @warn "Drive error, invalid speed: $speed"
        end

        calculatedspeed = round(Int, *(127, speed) / 100)
        sendCommand(command, calculatedspeed)
    end
# stop(): Stops both motors
    function stop()
        drive(1, 0)
        drive(2, 0)
    end
# driveBoth(speed_m1, speed_m2):
    function driveBoth(speed_m1, speed_m2)
        drive(1, speed_m1)
        drive(2, speed_m2)
    end
# driveTank(speed)
    function driveTank(speed)
        drive(1, speed)
        drive(2, speed)
    end
end


text:
using Distributions

struct MPointvariate <: VariateForm
end

struct TVectorvariate <: VariateForm
end

"""
    MPointSupport(manifold::Manifold)

Value support for manifold-valued distributions (values from given
manifold `M`).
"""
struct MPointSupport{TM<:Manifold} <: ValueSupport
    manifold::TM
end

"""
    MPointDistribution{TM<:Manifold}

An abstract distribution for points on manifold of type `TM`.
"""
abstract type MPointDistribution{TM<:Manifold} <: Distribution{MPointvariate, MPointSupport{TM}}
end

"""
    TVectorSupport(manifold::Manifold, x)

Value support for tangent vector-valued distributions (values from tangent
space at point `x` from given manifold).
"""
struct TVectorSupport{TM<:Manifold, T} <: ValueSupport
    manifold::TM
    x::T
end

"""
    TVectorDistribution{TM<:Manifold, T}

An abstract distribution for tangent vectors at point of type `T`
from manifold of type `TM`.
"""
abstract type TVectorDistribution{TM<:Manifold, T} <: Distribution{TVectorvariate, TVectorSupport{TM, T}}
end

"""
    support(d::TVectorDistribution)

Get the object of type `TVectorSupport` for the distribution `d`.
"""
function Distributions.support(::T) where T<:TVectorDistribution
    error("support not implemented for type $T")
end


text:
include("Core.jl")
include("RealSpace.jl")
include("SiteSpace.jl")


text:
struct FlatNode{NBranches,T,Values<:Union{Nothing,MutableNTuple{NBranches,T}}}
    count::NTuple{2,Threads.Atomic{Int}}
    values::NTuple{2,Values}
end

function FlatNode{NBranches,T}() where {NBranches,T}
    if Base.issingletontype(T)
        Values = Nothing
        values = (nothing, nothing)
    else
        Values = MutableNTuple{NBranches,T}
        values = (Values(), Values())
    end
    count = (Threads.Atomic{Int}(0), Threads.Atomic{Int}(0))
    return FlatNode{NBranches,T,Values}(count, values)
end

struct FlatTreeBarrier{NBranches,T,Op,Values} <: SyncBarriers.FlatTreeBarrier{NBranches,T}
    n::Int
    flags::Matrix{OneWayObservable{Bool}}
    locals::typeof(sense_parity_states(1))
    nodes::Vector{FlatNode{NBranches,T,Values}}
    values::MutableNTuple{2,T}
    op::Op
end
# TODO: check if parity is required

FlatTreeBarrier{NBranches}(n) where {NBranches} =
    FlatTreeBarrier{NBranches,Nothing}(right, n)

function FlatTreeBarrier{NBranches,T}(op::Op, n::Integer) where {NBranches,T,Op}
    @argcheck NBranches isa Integer && NBranches > 1 && T isa Type && n > 0
    nnodes = foldl_leaf_to_root((_, x) -> x.stop, nothing, n, Val(NBranches), 1)
    nodes = [FlatNode{NBranches,T}() for _ in 1:nnodes]
    flags = [OneWayObservable{Bool}(false) for _ in 1:n, _ in 1:2]
    values = MutableNTuple{2,T}()
    return FlatTreeBarrier(
        n,
        flags,
        sense_parity_states(n),
        nodes,
        values,
        op,
    )::FlatTreeBarrier{NBranches,T,Op}
end

function SyncBarriers.reduce!(
    handle::BarrierHandle{<:FlatTreeBarrier},
    value,
    spin::Union{Integer,Nothing} = nothing,
)
    acc1 = SyncBarriers.reduce_arrive!(handle, value)
    acc2 = SyncBarriers.depart!(handle, spin)
    if acc1 isa Some
        @_assert acc2 === something(acc1)
        return something(acc1)
    else
        return acc2
    end
end

function SyncBarriers.reduce_arrive!(
    handle::BarrierHandle{<:FlatTreeBarrier{NBranches,T}},
    value,
) where {NBranches,T}
    value = convert(T, value)
    barrier = handle.barrier
    i = handle.i
    s, parity = inc_sense_parity!(view(barrier.locals, handle.i))
    acc = foldl_leaf_to_root(Some(value), barrier.n, Val(NBranches), i) do acc, x
        node = barrier.nodes[x.inode]
        cnt = node.count[parity+1]
        vals = node.values[parity+1]
        if vals !== nothing
            vals[x.iself] = something(acc)
        end
        if Threads.atomic_add!(cnt, 1) == length(x.branches) - 1
            cnt[] = 0
            if vals === nothing
                return acc  # continue
            else
                a = vals[1]
                for j in x.branches[2:end]
                    a = barrier.op(a, vals[j])
                end
                return Some(a)
            end
        else
            return Break(nothing)
        end
    end
    if acc isa Break
        return nothing
    else
        barrier.values[parity+1] = something(acc)
        for (j, flag) in pairs(@view barrier.flags[:, parity+1])
            flag[] = s
        end
        return acc::Some
    end
end

function SyncBarriers.depart!(
    handle::BarrierHandle{<:FlatTreeBarrier},
    spin::Union{Integer,Nothing} = nothing,
)
    barrier = handle.barrier
    s, parity = prev_sense_parity(view(barrier.locals, handle.i))
    flag = barrier.flags[handle.i, parity+1]
    waitif(!=(s), flag, spin)
    return barrier.values[parity+1]
end