text:
module BenefitGenerosity
#
# This is crazily over/under designed and needs re-doing
# Je n'ai fait celle-ci plus longue que parce que je n'ai pas eu le loisir de la faire plus courte.
#
using ScottishTaxBenefitModel
using .ModelHousehold: OneIndex,Person,Household
using .Definitions
using .STBIncomes
using .STBParameters: NonMeansTestedSys
using .RunSettings: Settings
using .FRSHouseholdGetter: get_household_of_person
using DataFrames, CSV
using ArgCheck

export initialise, to_set, adjust_disability_eligibility!, change_status

struct GenEntry # FIXME make a concrete type {T} ???
    cum_popn :: Real
    pid      :: BigInt
    data_year :: Int
end

const GenVec = Vector{GenEntry}

struct DisabilityChanges
    which  :: Incomes
    is_positive  :: Bool
    people :: Set{OneIndex}
end

#
# FIXME re-arrange to make this immutable
#
mutable struct EntryWrapper
    negative_candidates_aa::GenVec
    negative_candidates_pip_mob::GenVec
    positive_candidates_pip_care::GenVec
    negative_candidates_dla_children::GenVec
    positive_candidates_aa::GenVec
    positive_candidates_pip_mob::GenVec
    negative_candidates_pip_care::GenVec
    positive_candidates_dla_children::GenVec
end

const ENTRIES = EntryWrapper(
    GenVec(),
    GenVec(),
    GenVec(),
    GenVec(),
    GenVec(),
    GenVec(),
    GenVec(),
    GenVec()
)

function load_one( filename :: String ) :: GenVec
    d = CSV.File( filename ) |> DataFrame
    n = size(d)[1]
    gv = Vector{GenEntry}(undef, n)
    popn = 0.0
    for i in 1:n
        r = d[i,:]
        hh = get_household_of_person( BigInt(r.pid), r.data_year )
        if hh !== nothing
            popn += hh.weight
            gv[i] = GenEntry( popn, BigInt(r.pid), r.data_year )
        else
            println( "no such hh $(r.pid) $(r.data_year)")
        end
    end
    return gv
end

function initialise( dir :: String )
    ENTRIES.negative_candidates_aa = load_one( "$dir/negative_candidates_aa.csv" )
    ENTRIES.negative_candidates_pip_mob = load_one( "$dir/negative_candidates_pip_mob.csv" )
    ENTRIES.positive_candidates_pip_care = load_one( "$dir/positive_candidates_pip_care.csv" )
    ENTRIES.negative_candidates_dla_children = load_one( "$dir/negative_candidates_dla_children.csv" )
    ENTRIES.positive_candidates_aa = load_one( "$dir/positive_candidates_aa.csv" )
    ENTRIES.positive_candidates_pip_mob = load_one( "$dir/positive_candidates_pip_mob.csv" )
    ENTRIES.negative_candidates_pip_care = load_one( "$dir/negative_candidates_pip_care.csv" )
    ENTRIES.positive_candidates_dla_children = load_one( "$dir/positive_candidates_dla_children.csv" )
end

function make_one_set( extra_people::Number, candidates :: GenVec ) :: Set{OneIndex}
    s = Set{OneIndex}()
    for c in candidates
        if c.cum_popn >= extra_people
            break
        end
        push!(s, OneIndex( c.pid, c.data_year ))
    end
    return s
end
#
# FIXME just merge with the thing below
#
function to_set( which :: Incomes, extra_people :: Real ) :: Set{OneIndex}
    @argcheck which in SICKNESS_ILLNESS
    s = Set{OneIndex}()
    if extra_people > 0
        if which == ATTENDANCE_ALLOWANCE
            return make_one_set( extra_people, ENTRIES.positive_candidates_aa )
        elseif which == PERSONAL_INDEPENDENCE_PAYMENT_DAILY_LIVING
            return make_one_set( extra_people, ENTRIES.positive_candidates_pip_care )
        elseif which == PERSONAL_INDEPENDENCE_PAYMENT_MOBILITY
            return make_one_set( extra_people, ENTRIES.positive_candidates_pip_mob )
        elseif which == DLA_SELF_CARE || which == DLA_MOBILITY
            return make_one_set( extra_people, ENTRIES.positive_candidates_dla_children )
        end
    elseif extra_people < 0
        if which == ATTENDANCE_ALLOWANCE
            return make_one_set( -1*extra_people, ENTRIES.negative_candidates_aa )
        elseif which == PERSONAL_INDEPENDENCE_PAYMENT_DAILY_LIVING
            return make_one_set( -1*extra_people, ENTRIES.negative_candidates_pip_care )
        elseif which == PERSONAL_INDEPENDENCE_PAYMENT_MOBILITY
            return make_one_set( -1*extra_people, ENTRIES.negative_candidates_pip_mob )
        elseif which == DLA_SELF_CARE || which == DLA_MOBILITY
            return make_one_set( -1*extra_people, ENTRIES.negative_candidates_dla_children )
        end
    end
    return s
end

function adjust_disability_eligibility!( nmt_bens :: NonMeansTestedSys )
    nmt_bens.attendance_allowance.candidates = to_set(  ATTENDANCE_ALLOWANCE, nmt_bens.attendance_allowance.extra_people )
    nmt_bens.dla.candidates= to_set( DLA_SELF_CARE, nmt_bens.dla.extra_people )
    nmt_bens.pip.dl_candidates = to_set( PERSONAL_INDEPENDENCE_PAYMENT_DAILY_LIVING, nmt_bens.pip.extra_people )
    nmt_bens.pip.mobility_candidates = to_set( PERSONAL_INDEPENDENCE_PAYMENT_MOBILITY, nmt_bens.pip.extra_people )
end

function change_status(
    ;
    candidates::Set{OneIndex},
    pid::BigInt,
    change ::Real,
    choices,
    current_value,
    disqual_value )
    if change == 0
        return current_value
    elseif ! in_indexes( candidates, pid )
        return current_value
    else
        if change < 0
            return disqual_value
        else
            return rand(choices)
        end
    end
end

end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
function run!(ebm::EBM, cn::Constants, snowdepth::Vector{Float64}, SWE::Vector{Float64}, Tsurf::Vector{Float64}, in::Input)

    D = zeros(ebm.Nsmax)
    S = zeros(ebm.Nsmax)
    W = zeros(ebm.Nsmax)

    for i in 1:length(in.year)

        year = in.year[i]
        month = in.month[i]
        day = in.day[i]
        hour = in.hour[i]
        SW = in.SW[i]
        LW = in.LW[i]
        Sf = in.Sf[i]
        Rf = in.Rf[i]
        Ta = in.Ta[i]
        RH = in.RH[i]
        Ua = in.Ua[i]
        Ps = in.Ps[i]

        Qs = qsat(true, Ps, Ta, cn)
        Qa = (RH / 100) * Qs

        Ua = max(Ua, 0.1)

        surf_props(ebm, cn, Sf)

        for i in 1:6
            surf_exch(ebm, cn, Ta, Ua)
            surf_ebal(ebm, cn, Ta, Qa, Ua, Ps, SW, LW)
        end

        snowdepth[i], SWE[i] = snow(ebm, cn, Sf, Rf, Ta, D, S, W)

        soil(ebm)

        Tsurf[i] = ebm.Tsurf

    end

end

ext: jl


text:
using FindClosest
using StaticArrays, Random
using Test

@testset "FindClosest.jl" begin
    for D in 1:10, l in 2:100, i in 1:100
        Random.seed!(i)
        pts = rand(SVector{D,Float64},l)
        @test findclosest(pts) == FindClosest.bruteforce(pts)
    end
end


text:
abstract type HMCTuningAlgorithm <: MCMCTuningAlgorithm end

"""
    struct HamiltonianMC <: MCMCAlgorithm

Hamiltonian Monte Carlo sampling algorithm.

The gradient of the target density is computed via auto-differentiation.

* Note: The fields of `HamiltonianMC` are still subject to change, and not
yet part of stable public API!*

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
@with_kw struct HamiltonianMC{
    MT<:HMCMetric,
    IT<:HMCIntegrator,
    PR<:HMCProposal,
    TN<:HMCTuningAlgorithm
} <: MCMCAlgorithm
    metric::MT = DiagEuclideanMetric()
    integrator::IT = LeapfrogIntegrator()
    proposal::PR = NUTS()
    tuning::TN = StanHMCTuning()
end

export HamiltonianMC

bat_default(::Type{MCMCSampling}, ::Val{:trafo}, mcalg::HamiltonianMC) = PriorToGaussian()

bat_default(::Type{MCMCSampling}, ::Val{:nsteps}, mcalg::HamiltonianMC, trafo::AbstractDensityTransformTarget, nchains::Integer) = 10^4

bat_default(::Type{MCMCSampling}, ::Val{:init}, mcalg::HamiltonianMC, trafo::AbstractDensityTransformTarget, nchains::Integer, nsteps::Integer) =
    MCMCChainPoolInit(nsteps_init = 25) # clamp(div(nsteps, 100), 25, 250)

bat_default(::Type{MCMCSampling}, ::Val{:burnin}, mcalg::HamiltonianMC, trafo::AbstractDensityTransformTarget, nchains::Integer, nsteps::Integer) =
    MCMCMultiCycleBurnin(nsteps_per_cycle = max(div(nsteps, 10), 250), max_ncycles = 4)

get_mcmc_tuning(algorithm::HamiltonianMC) = algorithm.tuning

# MCMCIterator subtype for HamiltonianMC
mutable struct AHMCIterator{
    AL<:HamiltonianMC,
    D<:AbstractDensity,
    R<:AbstractRNG,
    PR<:RNGPartition,
    SV<:DensitySampleVector,
    HA<:AdvancedHMC.Hamiltonian,
    TR<:AdvancedHMC.Transition,
    PL<:AdvancedHMC.AbstractProposal
} <: MCMCIterator
    algorithm::AL
    density::D
    rng::R
    rngpart_cycle::PR
    info::MCMCIteratorInfo
    samples::SV
    nsamples::Int64
    stepno::Int64
    hamiltonian::HA
    transition::TR
    proposal::PL
end

function AHMCIterator(
    rng::AbstractRNG,
    algorithm::MCMCAlgorithm,
    density::AbstractDensity,
    info::MCMCIteratorInfo,
    x_init::AbstractVector{P},
) where {P<:Real}
    stepno::Int64 = 0

    npar = totalndof(density)

    params_vec = Vector{P}(undef, npar)
    params_vec .= x_init
    !(params_vec in var_bounds(density)) && throw(ArgumentError("Parameter(s) out of bounds"))

    log_posterior_value = eval_logval(density, params_vec, default_dlt())

    T = typeof(log_posterior_value)
    W = Float64 # ToDo: Support other sample weight types

    sample_info = AHMCSampleID(info.id, info.cycle, 1, CURRENT_SAMPLE, 0.0, 0, false, 0.0)
    current_sample = DensitySample(params_vec, log_posterior_value, one(W), sample_info, nothing)
    samples = DensitySampleVector{Vector{P},T,W,AHMCSampleID,Nothing}(undef, 0, npar)
    push!(samples, current_sample)

    nsamples::Int64 = 0

    rngpart_cycle = RNGPartition(rng, 0:(typemax(Int16) - 2))

    metric = ahmc_metric(algorithm.metric, npar)

    f = logdensityof(density)
    fg = fg = valgradof(f)

    init_hamiltonian = AdvancedHMC.Hamiltonian(metric, f, fg)
    hamiltonian, init_transition = AdvancedHMC.sample_init(rng, init_hamiltonian, params_vec)
    integrator = ahmc_integrator(rng, algorithm.integrator, hamiltonian, params_vec)
    proposal = ahmc_proposal(algorithm.proposal, integrator)

    # Perform a dummy step to get type-stable transition value:
    transition = AdvancedHMC.step(deepcopy(rng), deepcopy(hamiltonian), deepcopy(proposal), init_transition.z)

    chain = AHMCIterator(
        algorithm,
        density,
        rng,
        rngpart_cycle,
        info,
        samples,
        nsamples,
        stepno,
        hamiltonian,
        transition,
        proposal
    )

    reset_rng_counters!(chain)

    chain
end

function MCMCIterator(
    rng::AbstractRNG,
    algorithm::HamiltonianMC,
    density::AbstractDensity,
    chainid::Integer,
    startpos::AbstractVector{<:Real}
)
    cycle = 0
    tuned = false
    converged = false
    info = MCMCIteratorInfo(chainid, cycle, tuned, converged)
    AHMCIterator(rng, algorithm, density, info, startpos)
end

@inline _current_sample_idx(chain::AHMCIterator) = firstindex(chain.samples)
@inline _proposed_sample_idx(chain::AHMCIterator) = lastindex(chain.samples)

getalgorithm(chain::AHMCIterator) = chain.algorithm

getdensity(chain::AHMCIterator) = chain.density

getrng(chain::AHMCIterator) = chain.rng

mcmc_info(chain::AHMCIterator) = chain.info

nsteps(chain::AHMCIterator) = chain.stepno

nsamples(chain::AHMCIterator) = chain.nsamples

current_sample(chain::AHMCIterator) = chain.samples[_current_sample_idx(chain)]

sample_type(chain::AHMCIterator) = eltype(chain.samples)

function reset_rng_counters!(chain::AHMCIterator)
    set_rng!(chain.rng, chain.rngpart_cycle, chain.info.cycle)
    rngpart_step = RNGPartition(chain.rng, 0:(typemax(Int32) - 2))
    set_rng!(chain.rng, rngpart_step, chain.stepno)
    nothing
end

function samples_available(chain::AHMCIterator)
    i = _current_sample_idx(chain::AHMCIterator)
    chain.samples.info.sampletype[i] == ACCEPTED_SAMPLE
end

function get_samples!(appendable, chain::AHMCIterator, nonzero_weights::Bool)::typeof(appendable)
    if samples_available(chain)
        samples = chain.samples

        for i in eachindex(samples)
            st = samples.info.sampletype[i]
            if (
                (st == ACCEPTED_SAMPLE || st == REJECTED_SAMPLE) &&
                (samples.weight[i] > 0 || !nonzero_weights)
            )
                push!(appendable, samples[i])
            end
        end
    end
    appendable
end

function next_cycle!(chain::AHMCIterator)
    _cleanup_samples(chain)

    chain.info = MCMCIteratorInfo(chain.info, cycle = chain.info.cycle + 1)
    chain.nsamples = 0
    chain.stepno = 0

    reset_rng_counters!(chain)

    resize!(chain.samples, 1)

    i = _proposed_sample_idx(chain)
    @assert chain.samples.info[i].sampletype == CURRENT_SAMPLE
    chain.samples.weight[i] = 1

    t_stat = chain.transition.stat

    chain.samples.info[i] = AHMCSampleID(
        chain.info.id, chain.info.cycle, chain.stepno, CURRENT_SAMPLE,
        t_stat.hamiltonian_energy, t_stat.tree_depth,
        t_stat.numerical_error, t_stat.step_size
    )

    chain
end

function _cleanup_samples(chain::AHMCIterator)
    samples = chain.samples
    current = _current_sample_idx(chain)
    proposed = _proposed_sample_idx(chain)
    if (current != proposed) && samples.info.sampletype[proposed] == CURRENT_SAMPLE
        # Proposal was accepted in the last step
        @assert samples.info.sampletype[current] == ACCEPTED_SAMPLE
        samples.v[current] .= samples.v[proposed]
        samples.logd[current] = samples.logd[proposed]
        samples.weight[current] = samples.weight[proposed]
        samples.info[current] = samples.info[proposed]

        resize!(samples, 1)
    end
end

function mcmc_step!(chain::AHMCIterator)
    _cleanup_samples(chain)

    samples = chain.samples
    algorithm = getalgorithm(chain)

    chain.stepno += 1
    reset_rng_counters!(chain)

    rng = getrng(chain)
    density = getdensity(chain)

    # Grow samples vector by one:
    resize!(samples, size(samples, 1) + 1)
    samples.info[lastindex(samples)] = AHMCSampleID(
        chain.info.id, chain.info.cycle, chain.stepno, PROPOSED_SAMPLE,
        0.0, 0, false, 0.0
    )

    current = _current_sample_idx(chain)
    proposed = _proposed_sample_idx(chain)
    @assert current != proposed

    current_params = samples.v[current]
    proposed_params = samples.v[proposed]

    # Propose new variate:
    samples.weight[proposed] = 0

    chain.transition = AdvancedHMC.step(rng, chain.hamiltonian, chain.proposal, chain.transition.z)
    proposed_params[:] = chain.transition.z.θ

    current_log_posterior = samples.logd[current]
    T = typeof(current_log_posterior)

    # Evaluate prior and likelihood with proposed variate:
    proposed_log_posterior = eval_logval(density, proposed_params, default_dlt())

    samples.logd[proposed] = proposed_log_posterior

    accepted = current_params != proposed_params

    if accepted
        samples.info.sampletype[current] = ACCEPTED_SAMPLE
        samples.info.sampletype[proposed] = CURRENT_SAMPLE
        chain.nsamples += 1

        tstat = AdvancedHMC.stat(chain.transition)
        samples.info.hamiltonian_energy[proposed] = tstat.hamiltonian_energy
        # ToDo: Handle proposal-dependent tstat (only NUTS has tree_depth):
        samples.info.tree_depth[proposed] = tstat.tree_depth
        samples.info.divergent[proposed] = tstat.numerical_error
        samples.info.step_size[proposed] = tstat.step_size
    else
        samples.info.sampletype[proposed] = REJECTED_SAMPLE
    end

    delta_w_current, w_proposed = if accepted
        (0, 1)
    else
        (1, 0)
    end

    samples.weight[current] += delta_w_current
    samples.weight[proposed] = w_proposed

    nothing
end

eff_acceptance_ratio(chain::AHMCIterator) = nsamples(chain) / nsteps(chain)


text:
using SIAN

println("Setting up the problem")

ode = @ODEmodel(
  x4'(t) = - k5 * x4(t) // (k6 + x4(t)),
  x5'(t) = k5 * x4(t) // (k6 + x4(t)) - k7 * x5(t) / (k8 + x5(t) + x6(t)),
  x6'(t) = k7 * x5(t) // (k8 + x5(t) + x6(t)) - k9 * x6(t) * (k10 - x6(t)) // k10,
  x7'(t) = k9 * x6(t) * (k10 - x6(t)) // k10,
  y1(t) = x4(t),
  y2(t) = x5(t)
)

res = identifiability_ode(ode, get_parameters(ode); p=0.999, p_mod=0, nthrds=1)

println(res)

max_stars_repo_name: iliailmer/SIAN-Julia

max_issues_repo_name: iliailmer/SIAN-Julia

max_forks_repo_name: iliailmer/SIAN-Julia


text:
@testset "2 stages with $solver" for solver in lp_solvers
    include("prob5.2_data.jl")

    numScen = 2
    M1 = StructuredModel(num_scenarios=numScen)

    @variable M1 x1[1:n] >= 0
    @variable M1 v1[1:n] >= 0
    @constraints M1 begin
        x1 .== v1
    end
    @objective(M1, Min, dot(ic, v1))

    for s in 1:numScen
        M2 = StructuredModel(parent=M1, prob=p2[s], id=s)
        @variable(M2, y2[1:n, 1:m] >= 0)
        @constraints M2 begin
            demand[j=1:m], sum(y2[:,j]) == D2[j,s]
            ylim[i=1:n], sum(y2[i,:]) <= x1[i]
        end
        @objective(M2, Min, dot(C, y2 * T))
    end

    # With detectlb = true and MultiCut
    # Iteration | Status    |     LB    |      UB    | #FC | #OC |
    #     1     | Optimal   |      0.0  |     Inf    |  2  |  0  |
    #     2     | Optimal   |  22338    | 1196514.84 |  0  |  2  |
    #     3     | Optimal   |  89893.68 |  676982.10 |  0  |  2  |
    #     4     | Optimal   | 182717.07 |  353268.97 |  0  |  2  |
    #     5     | Optimal   | 303678.32 |  366193.45 |  0  |  2  |
    #     6     | Optimal   | 311405.77 |  356197.41 |  0  |  2  |
    #     7     | Optimal   | 319203.65 |  344833.83 |  0  |  2  |
    #     8     | Optimal   | 334922.82 |  368285.45 |  0  |  2  |
    #     9     | Optimal   | 339047.44 |  349576.21 |  0  |  1  |
    #    10     | Optimal   | 340315.52 |  340315.52 |  0  |  0  |
    # With detectlb = true and AveragedCut
    #     1     | Optimal   |      0.0  |     Inf    | #FC | #OC |
    #     2     | Optimal   |  22338    | 1196514.84 |  2  |  0  |
    #     3     | Optimal   |  73386.05 |  660474.47 |  0  |  1  |
    #     4     | Optimal   | 171778.69 |  356387.31 |  0  |  1  |
    #     5     | Optimal   | 297813.60 |  391523.18 |  0  |  1  |
    #     6     | Optimal   | 309298.85 |  353098.55 |  0  |  1  |
    #     7     | Optimal   | 314734.72 |  344892.19 |  0  |  1  |
    #     8     | Optimal   | 331094.09 |  376701.57 |  0  |  1  |
    #     9     | Optimal   | 336691.25 |  346276.39 |  0  |  1  |
    #    10     | Optimal   | 337198.76 |  346727.58 |  0  |  1  |
    #    11     | Optimal   | 339235.89 |  344363.19 |  0  |  1  | -> Pereira stops here
    #    12     | Optimal   | 339779.51 |  340949.37 |  0  |  1  |
    #    13     | Optimal   | 340315.52 |  340315.52 |  0  |  0  |
    function testniter(niter, K, maxncuts, cutgen, detectlb)
        if isa(cutgen, StructProg.MultiCutGenerator)
            @test niter == 10
        else
            if maxncuts == -1
                if K == -1
                    if detectlb
                        @test niter == 13
                    else
                        # 10 on Gurobi/CPLEX, 14 otherwise
                        @test niter == 10 || niter == 14
                    end
                else
                    # 9 on Gurobi/CPLEX, 11 otherwise
                    @test niter == 9 || niter == 11 # Pereira stops earlier
                end
            else
                if K == -1
                    if detectlb
                        # up to 15 on Mac OS, <= 14 on Linux
                        @test 13 <= niter <= 15
                    else
                        # 10 on Gurobi/CPLEX, between 14 and 15 otherwise
                        @test 10 <= niter <= 15
                    end
                else
                    # 9 on Gurobi/CPLEX, between 11 and 12 otherwise
                    @test 9 <= niter <= 12 # Pereira stops earlier
                end
            end
        end
    end
    fulltest(M1, 2, 340315.52, [5085,1311,3919,854], 334687.754566, 15869.996575, testniter, solver)
    # root = model2lattice(M1, 2, solver, cutmode)
    # sol = SDDP(root, 2, cutmode, :All, verbose)
    #
    # v11value = sol.sol[1:4]
    # @test sol.status == :Optimal
    # @test abs(sol.objval - 340315.52) < 0.1
    # @test norm(v11value - [5085,1311,3919,854]) < 0.1
    # SDDPclear(M1)
end


text:
using Libtask
using Test

include("tf.jl")
include("tapedtask.jl")
include("tarray.jl")
include("tref.jl")

if haskey(ENV, "BENCHMARK")
    include("benchmarks.jl")
end

max_stars_repo_name: TuringLang/Libtask

max_issues_repo_name: TuringLang/Libtask

max_forks_repo_name: TuringLang/Libtask


text:
@testset "Static Real" begin
    # `getvalues` is used to ensure that $f(sval, sval) --> sval

#    static_real_set = ([@inferred(S{B(1)}()) for (S,B) in S2B]...,)
#    base_real_set = ([B(1) for (S,B) in (S2B)]...,)

    @testset "values, eltype" begin
        for (S,B) in S2B
            @test @inferred(values(S(1))) == B(1)
            @test eltype(values(S(1))) == eltype(B(1))
        end
    end

    @testset "Addition" begin
        for (S1,B1) in S2B
            for (S2,B2) in S2B
                @test @inferred(getvalues(S1(1) + S2(1))) == B1(1) + B2(1)
                @test @inferred(values(S1(1) + B1(1))) === B1(1) + B1(1)
                @test @inferred(values(S1(1) + B2(1))) === B1(1) + B2(1)
            end
        end
    end

    @testset "Substraction" begin
        for (S1,B1) in S2B
            for (S2,B2) in S2B
                @test values(@inferred(S1(1) - S2(1))) === B1(1) - B2(1)
            end
        end
    end

    @testset "Multiplication" begin
        for (S1,B1) in S2B
            for (S2,B2) in S2B
                @test @inferred(one(S1) * one(S2)) == one(B1) * one(B2)
            end
        end
    end

    @testset "Division" begin
        for (S1,B1) in S2B
            for (S2,B2) in S2B
                @test @inferred(one(S1) / one(S2)) == one(B1) / one(B2)
            end
        end
    end

    @testset "fma" begin
        for (S1,B1) in S2B
            @test @inferred(fma(one(S1), one(S1), one(S1))) == fma(one(B1), one(B1), one(B1))
        end
    end

    # TODO: Inferrence problem
    @testset "muladd" begin
        for (S1,B1) in S2B
            @test @inferred(muladd(one(S1), one(S1), one(S1))) == muladd(one(B1), one(B1), one(B1))
        end
    end

    @testset "fld" begin
        for (S,B) in S2B
            @test @inferred(fld(S(5),S(2))) == fld(B(5),B(2))
        end
    end

    @testset "cld" begin
        for (S,B) in S2B
            @test @inferred(cld(S(5),S(2))) == cld(B(5),B(2))
        end
    end

    @testset "rem" begin
        for (S,B) in S2B
            @test @inferred(rem(S(15),S(4))) == rem(B(15),B(4))
        end
    end

    @testset "max" begin
        for (S,B) in S2B
            @test @inferred(max(one(S),S(3))) == max(one(B),B(3))
        end

    end

    @testset "min" begin
        for (S,B) in S2B
            @test @inferred(min(one(S),S(3))) == min(one(B),B(3))
        end
    end

    @testset "minmax" begin
        for (S,B) in S2B
            @test @inferred(minmax(one(S),S(3))) == minmax(one(B),B(3))
        end
    end
end

using StaticValues: sfloatrange, linspace, linspace1, srangehp
@testset "Ranges" begin
    @testset "colon" begin
        @inferred(SInt(10):SInt(1):SInt(0))

        @inferred(SInt(1):SVal(.2):SVal(2))

        @inferred(SVal(1.):SVal(.2):SVal(2.))

        @inferred(SVal(2):SVal(-.2):SVal(1))

        @inferred(SVal(0.0):SVal(-0.5))

        @inferred(SOne:SZero)
    end

    @testset "floatrange" begin
        T = Float32
        b = SVal(1)
        s = SVal(2)
        l = SVal(10)
        d = SVal(1)
        @test @inferred(sfloatrange(T, b, s, l, d)) == StepSRangeLen{Float32,SFloat64{1.0},SFloat64{2.0},SFloat64{19.0},SInt64{10},SInt64{1}}()

        b = SVal(1.0)
        s = SVal(2.0)
        d = SVal(1.0)
        @test @inferred(sfloatrange(b, s, l, d)) == StepSRangeLen{SFloat64,SFloat64{1.0},SFloat64{2.0},SFloat64{19.0},SInt64{10},SInt64{1}}()
    end

    @testset "linspace" begin
        T = Int
        b = SVal(1)
        e = SVal(10)
        l = SVal(5)
        @test @inferred(linspace(T, b, e, l)) == LinSRange{Int}(b, e, l, SInt(4))

        @inferred(linspace(Float64, b, e, l))

        b = SVal(1)
        e = SVal(10)
        l = SVal(5)
        d = SVal(1)
        @inferred(linspace(Float64, b, e, l, d))

        b = SVal(1.0)
        e = SVal(10.0)
        l = SVal(5)
        @inferred(linspace(b, e, l))

        T = Float16
        b = SVal(1.0)
        e = SVal(1.0)
        l = SVal(1)
        @test @inferred(linspace1(T, b, e, l)) == StepSRangeLen{Float16,SFloat64{1.0},SFloat64{0.0},SFloat64{1.0},SInt64{1},SInt64{1}}()

        # src/linspace.jl line 159
#        @test @inferred(linrange(T,b,e,l)) == StaticRanges.SRange{Float16,SVal{1.0,Float64},SVal{0.0,Float64},Float16(1.0),1,1}()
    end

    @testset "srangehp" begin
        b = (SInt128(1), SInt128(1))
        s = (SInt128(1), SInt128(1))
        nb = SOne
        l = SInt(2)
        f = SInt(1)
        T = Float64

        @test @inferred(srangehp(T, b, s, nb, l, f)) == StepSRangeLen{Float64,TPVal{SFloat64{1.0},SFloat64{0.0}},TPVal{SFloat64{1.0},SFloat64{0.0}},SFloat64{2.0},SInt64{2},SInt64{1}}()

        @test @inferred(srangehp(Float32, b, s, nb, l, f)) == StepSRangeLen{Float32,SFloat64{1.0},SFloat64{1.0},SFloat64{2.0},SInt64{2},SInt64{1}}()

        b = SVal(1.0)
        s = SVal(1.0)
        @test @inferred(srangehp(T, b, s, nb, l, f)) == StepSRangeLen{Float64,TPVal{SFloat64{1.0},SFloat64{0.0}},TPVal{SFloat64{1.0},SFloat64{0.0}},SFloat64{2.0},SInt64{2},SInt64{1}}()

        b = (SVal(1.0), SVal(0.0))
        s = (SVal(1.0), SVal(0.0))
        @test @inferred(srangehp(T, b, s, nb, l, f)) == StepSRangeLen{Float64,TPVal{SFloat64{1.0},SFloat64{0.0}},TPVal{SFloat64{1.0},SFloat64{0.0}},SFloat64{2.0},SInt64{2},SInt64{1}}()

        @test @inferred(srangehp(Float32, b, s, nb, l, f)) == StepSRangeLen{Float32,SFloat64{1.0},SFloat64{1.0},SFloat64{2.0},SInt64{2},SInt64{1}}()

        b = SVal(1.0)
        s = SVal(1.0)
        @test @inferred(srangehp(Float32, b, s, nb, l, f)) == StepSRangeLen{Float32,SFloat64{1.0},SFloat64{1.0},SFloat64{2.0},SInt64{2},SInt64{1}}()
    end

    #=
    @testset "steprange" begin
        # steprange_length Int
        @test length(steprange(SVal(1), SVal(1), SVal(5))) == 5
        # steprange_length general
        @test length(steprange(SVal(1.0), SVal(1.0), SVal(5.0))) == 5
        # steprange_last start == stop
        @test last(steprange(SVal(4), SVal(1), SVal(4))) == 4
        # steprange_last general
        @test last(steprange(SVal(1.0), SVal(1.0), SVal(5.0))) == 5.0
        # steprange_last_empty integer
        # steprange_last_empty general
    end
    =#
end


text:
using Pkg.TOML, StaticArrays

@testset "Receivers.jl" begin
    r = rand(SVector{3,Float64})
    t = rand(Float64)

    params = TOML.parsefile("Receivers.toml")

    xSensitivity = MPISimulations.initialize(MPISimulations.StaticHomogeneousField,"Receiver.xSensitivity",params,Float64)
    ySensitivity = MPISimulations.initialize(MPISimulations.StaticHomogeneousField,"Receiver.ySensitivity",params,Float64)
    zSensitivity = MPISimulations.initialize(MPISimulations.StaticHomogeneousField,"Receiver.zSensitivity",params,Float64)

    samplingRate = params["Receiver"]["samplingRate"]

    receiver1D = MPISimulations.ReceiveCoils(samplingRate,xSensitivity)
    @test MPISimulations.receiveCoilSensitivities(receiver1D,r,t) == SMatrix{1,3}(0.1,0.0,0.0)

    receiver2D = MPISimulations.ReceiveCoils(samplingRate,xSensitivity,ySensitivity)
    @test MPISimulations.receiveCoilSensitivities(receiver2D,r,t) == SMatrix{2,3}(0.1,0.0,0.0,0.11,0.0,0.0)

    receiver3D = MPISimulations.ReceiveCoils(samplingRate,xSensitivity,ySensitivity,zSensitivity)
    @test MPISimulations.receiveCoilSensitivities(receiver3D,r,t) == SMatrix{3,3}(0.1,0.0,0.0,0.0,0.11,0.0,0.0,0.0,0.12)
end


text:
# Functions to calculate the inital state vector, the inital masses and initial springs

function init_springs!(s::KPS4)
    l_0     = s.set.l_tether / s.set.segments
    particles = KiteUtils.get_particles(s.set.height_k, s.set.h_bridle, s.set.width, s.set.m_k)
    for j in 1:size(SPRINGS_INPUT)[1]
        # build the tether segments
        if j == 1
            for i in 1:s.set.segments
                k = s.set.e_tether * (s.set.d_tether/2000.0)^2 * pi  / l_0  # Spring stiffness for this spring [N/m]
                c = s.set.damping/l_0                                       # Damping coefficient [Ns/m]
                s.springs[i] = SP(i, i+1, l_0, k, c)
            end
        # build the bridle segments
        else
            p0, p1 = SPRINGS_INPUT[j, 1]+1, SPRINGS_INPUT[j, 2]+1 # point 0 and 1
            if SPRINGS_INPUT[j, 3] == -1
                l_0 = norm(particles[Int(p1)] - particles[Int(p0)]) * PRE_STRESS
                k = s.set.e_tether * (s.set.d_line/2000.0)^2 * pi / l_0
                p0 += s.set.segments - 1 # correct the index for the start and end particles of the bridle
                p1 += s.set.segments - 1
                c = s.set.damping/ l_0
                s.springs[j+s.set.segments-1] = SP(Int(p0), Int(p1), l_0, k, c)
            end
        end
    end
    s.springs
end

function init_masses!(s::KPS4)
    s.masses = zeros(s.set.segments+KITE_PARTICLES+1)
    l_0 = s.set.l_tether / s.set.segments
    for i in 1:s.set.segments
        s.masses[i]   += 0.5 * l_0 * s.set.rho_tether * (s.set.d_tether/2000.0)^2 * pi
        s.masses[i+1] += 0.5 * l_0 * s.set.rho_tether * (s.set.d_tether/2000.0)^2 * pi
    end
    s.masses[s.set.segments+1] += s.set.kcu_mass
    k2 = s.set.rel_top_mass * (1.0 - s.set.rel_nose_mass)
    k3 = 0.5 * (1.0 - s.set.rel_top_mass) * (1.0 - s.set.rel_nose_mass)
    k4 = 0.5 * (1.0 - s.set.rel_top_mass) * (1.0 - s.set.rel_nose_mass)
    s.masses[s.set.segments+2] += s.set.rel_nose_mass * s.set.mass
    s.masses[s.set.segments+3] += k2 * s.set.mass
    s.masses[s.set.segments+4] += k3 * s.set.mass
    s.masses[s.set.segments+5] += k4 * s.set.mass
    s.masses
end

function init_pos_vel_acc(s::KPS4, X=zeros(2 * (s.set.segments+KITE_PARTICLES)+1); old=false, delta = 0.0)
    pos = zeros(SVector{s.set.segments+1+KITE_PARTICLES, KVec3})
    vel = zeros(SVector{s.set.segments+1+KITE_PARTICLES, KVec3})
    acc = zeros(SVector{s.set.segments+1+KITE_PARTICLES, KVec3})
    pos[1] .= [0.0, delta, 0.0]
    vel[1] .= [delta, delta, delta]
    acc[1] .= [delta, delta, delta]
    sin_el, cos_el = sin(s.set.elevation / 180.0 * pi), cos(s.set.elevation / 180.0 * pi)
    for i in 1:s.set.segments
        radius = -i * (s.set.l_tether/s.set.segments)
        pos[i+1] .= [-cos_el * radius + X[i], delta, -sin_el * radius + X[s.set.segments+KITE_PARTICLES-1+i]]
        vel[i+1] .= [delta, delta, 0]
        acc[i+1] .= [delta, delta, -9.81]
    end
    vec_c = pos[s.set.segments] - pos[s.set.segments+1]
    if old
        particles = KiteUtils.get_particles(s.set.height_k, s.set.h_bridle, s.set.width, s.set.m_k)
    else
        particles = KiteUtils.get_particles(s.set.height_k, s.set.h_bridle, s.set.width, s.set.m_k, pos[s.set.segments+1], rotate_in_xz(vec_c, deg2rad(KITE_ANGLE)), s.v_apparent)
    end
    j = 1
    for i in [1,2,3] # set p8, p9, p10
        pos[s.set.segments+1+i] .= particles[i+2] + [X[s.set.segments+j], 0, X[2*s.set.segments+KITE_PARTICLES-1+j]]
        vel[s.set.segments+1+i] .= [delta, delta, delta]
        acc[s.set.segments+1+i] .= [delta, delta, -9.81]
        j +=1
    end
    acc[s.set.segments+1+4] .= [delta, delta, -9.81]
    vel[s.set.segments+1+4] .= [delta, delta, delta]
    # set p10=C and p11=D
    # x and z component of the right and left particle must be equal
    pos[s.set.segments+1+4][1] = pos[s.set.segments+1+3][1]  # D.x = C.x
    pos[s.set.segments+1+4][3] = pos[s.set.segments+1+3][3]  # D.z = C.z
    pos[s.set.segments+1+3][2] += X[end]                     # Y position of point C
    pos[s.set.segments+1+4][2] = -pos[s.set.segments+1+3][2] # Y position of point D
    for i in 1:length(pos)
        s.pos[i] .= pos[i]
    end
    pos, vel, acc
end

# Calculate the initial vectors pos and vel. Tether with the initial elevation angle
# se().elevation, particle zero fixed at origin.
# X is a vector of deviations in x and z positions, to be varied to find the inital equilibrium
function init_pos_vel(s::KPS4, X=zeros(2 * (s.set.segments+KITE_PARTICLES)))
    pos, vel, acc = init_pos_vel_acc(s, X; old=true)
    pos, vel
end

function init_inner(s::KPS4, X=zeros(2 * (s.set.segments+KITE_PARTICLES-1)+1); old=false, delta=0.0)
    pos, vel, acc = init_pos_vel_acc(s, X; old=old, delta=delta)
    vcat(pos[2:end], vel[2:end]), vcat(vel[2:end], acc[2:end])
end

# same as above, but returns a tuple of two one dimensional arrays
function init(s::KPS4, X=zeros(2 * (s.set.segments+KITE_PARTICLES-1)+1); old=false, delta=0.0)
    res1_, res2_ = init_inner(s, X; old=old, delta = delta)
    res1, res2  = vcat(reduce(vcat, res1_), [s.l_tether, 0]), vcat(reduce(vcat, res2_),[0,0])
    MVector{6*(s.set.segments+KITE_PARTICLES)+2, Float64}(res1), MVector{6*(s.set.segments+KITE_PARTICLES)+2, Float64}(res2)
end

# rotate a 3d vector around the y axis
function rotate_in_xz(vec, angle)
    result = similar(vec)
    result[1] = cos(angle) * vec[1] - sin(angle) * vec[3]
    result[2] = vec[2]
    result[3] = cos(angle) * vec[3] + sin(angle) * vec[1]
    result
end


text:
using Oceananigans.Coriolis: Ω_Earth
using Oceananigans.Coriolis: VectorInvariantEnergyConserving, VectorInvariantEnstrophyConserving

function instantiate_fplane_1(FT)
    coriolis = FPlane(FT, f=π)
    return coriolis.f == FT(π)
end

function instantiate_fplane_2(FT)
    coriolis = FPlane(FT, rotation_rate=2, latitude=30)
    return coriolis.f == FT(2)
end

function instantiate_ntfplane_1(FT)
    coriolis = NonTraditionalFPlane(FT, fz=π, fy=π)
    @test coriolis.fz == FT(π)
    @test coriolis.fy == FT(π)
end

function instantiate_ntfplane_2(FT)
    coriolis = NonTraditionalFPlane(FT, rotation_rate=2, latitude=30)
    @test coriolis.fz == FT(2)
    @test coriolis.fy == FT(2*√3)
end

function instantiate_betaplane_1(FT)
    coriolis = BetaPlane(FT, f₀=π, β=2π)
    @test coriolis.f₀ == FT(π)
    @test coriolis.β  == FT(2π)
end

function instantiate_betaplane_2(FT)
    coriolis = BetaPlane(FT, latitude=70, radius=2π, rotation_rate=3π)
    @test coriolis.f₀ == FT(6π * sind(70))
    @test coriolis.β  == FT(6π * cosd(70) / 2π)
end

function instantiate_ntbetaplane_1(FT)
    coriolis = NonTraditionalBetaPlane(FT, fz=π, fy=ℯ, β=1//7, γ=5)
    @test coriolis.fz == FT(π)
    @test coriolis.fy == FT(ℯ)
    @test coriolis.β  == FT(1//7)
    @test coriolis.γ  == FT(5)
end

function instantiate_ntbetaplane_2(FT)
    Ω, φ, R = π, 17, ℯ
    coriolis = NonTraditionalBetaPlane(FT, rotation_rate=Ω, latitude=φ, radius=R)
    @test coriolis.fz == FT(+2Ω*sind(φ))
    @test coriolis.fy == FT(+2Ω*cosd(φ))
    @test coriolis.β  == FT(+2Ω*cosd(φ)/R)
    @test coriolis.γ  == FT(-4Ω*sind(φ)/R)
end

function instantiate_hydrostatic_spherical_coriolis1(FT)
    coriolis = HydrostaticSphericalCoriolis(FT, scheme=VectorInvariantEnergyConserving())
    @test coriolis.rotation_rate == FT(Ω_Earth)
    @test coriolis.scheme isa VectorInvariantEnergyConserving

    coriolis = HydrostaticSphericalCoriolis(FT, scheme=VectorInvariantEnstrophyConserving())
    @test coriolis.rotation_rate == FT(Ω_Earth)
    @test coriolis.scheme isa VectorInvariantEnstrophyConserving
end

function instantiate_hydrostatic_spherical_coriolis2(FT)
    coriolis = HydrostaticSphericalCoriolis(FT, rotation_rate=π)
    @test coriolis.rotation_rate == FT(π)
    @test coriolis.scheme isa VectorInvariantEnergyConserving # default
end

@testset "Coriolis" begin
    @info "Testing Coriolis..."

    @testset "Coriolis" begin
        for FT in float_types
            @test instantiate_fplane_1(FT)
            @test instantiate_fplane_2(FT)

            instantiate_ntfplane_1(FT)
            instantiate_ntfplane_2(FT)
            instantiate_betaplane_1(FT)
            instantiate_betaplane_2(FT)
            instantiate_hydrostatic_spherical_coriolis1(FT)
            instantiate_hydrostatic_spherical_coriolis2(FT)

            # Test that FPlane throws an ArgumentError
            @test_throws ArgumentError FPlane(FT)
            @test_throws ArgumentError FPlane(FT, rotation_rate=7e-5)
            @test_throws ArgumentError FPlane(FT, f=1, latitude=40)
            @test_throws ArgumentError FPlane(FT, f=1, rotation_rate=7e-5, latitude=40)

            # Test that NonTraditionalFPlane throws an ArgumentError
            @test_throws ArgumentError NonTraditionalFPlane(FT)
            @test_throws ArgumentError NonTraditionalFPlane(FT, rotation_rate=7e-5)
            @test_throws ArgumentError NonTraditionalFPlane(FT, fz=1, latitude=40)
            @test_throws ArgumentError NonTraditionalFPlane(FT, fz=1, rotation_rate=7e-5, latitude=40)
            @test_throws ArgumentError NonTraditionalFPlane(FT, fy=1, latitude=40)
            @test_throws ArgumentError NonTraditionalFPlane(FT, fy=1, rotation_rate=7e-5, latitude=40)
            @test_throws ArgumentError NonTraditionalFPlane(FT, fz=1, fy=2, latitude=40)
            @test_throws ArgumentError NonTraditionalFPlane(FT, fz=1, fy=2, rotation_rate=7e-5, latitude=40)

            # Non-exhaustively test that BetaPlane throws an ArgumentError
            @test_throws ArgumentError BetaPlane(FT)
            @test_throws ArgumentError BetaPlane(FT, f₀=1)
            @test_throws ArgumentError BetaPlane(FT, β=1)
            @test_throws ArgumentError BetaPlane(FT, f₀=1e-4, β=1e-11, latitude=70)

            # Test that NonTraditionalBetaPlane throws an ArgumentError
            @test_throws ArgumentError NonTraditionalBetaPlane(FT)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, rotation_rate=7e-5)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fz=1, latitude=40)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fz=1, rotation_rate=7e-5, latitude=40)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fy=1, latitude=40)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fy=1, rotation_rate=7e-5, latitude=40)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fz=1, fy=2, latitude=40)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fz=1, fy=2, rotation_rate=7e-5, latitude=40)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fz=1, fy=2, β=3, latitude=40)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fz=1, fy=2, β=3, rotation_rate=7e-5, latitude=40)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fz=1, fy=2, β=3, γ=4, latitude=40)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fz=1, fy=2, β=3, γ=4, rotation_rate=7e-5, latitude=40)

            # Test show functions
            ✈ = FPlane(FT, latitude=45)
            show(✈); println()
            @test ✈ isa FPlane{FT}

            ✈ = NonTraditionalFPlane(FT, latitude=45)
            show(✈); println()
            @test ✈ isa NonTraditionalFPlane{FT}

            ✈ = BetaPlane(FT, latitude=45)
            show(✈); println()
            @test ✈ isa BetaPlane{FT}

            ✈ = NonTraditionalBetaPlane(FT, latitude=45)
            show(✈); println()
            @test ✈ isa NonTraditionalBetaPlane{FT}
        end
    end
end


text:
using PyPlot

"""
    polyplot(p, color)

Plot a polygon given by the point list `p` with the given `color`.
"""
function polyplot(p, color)
    n = size(p, 2)
    if p[:,n] == p[:,1]
        plot(p[1,:], p[2,:], color=color, linewidth=1)
    else
        plot([p[1,:] p[1,1]], [p[2,:] p[2,1]], color=color, linewidth=1)
    end
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
# This file contains functions for the Cost Efficiency DEA model
"""
    CostDEAModel
An data structure representing a cost DEA model.
"""
struct CostDEAModel <: AbstractCostDEAModel
    n::Int64
    m::Int64
    s::Int64
    rts::Symbol
    eff::Vector
    lambda::SparseMatrixCSC{Float64, Int64}
    techeff::Vector
    alloceff::Vector
end

"""
    deacost(X, Y, W)
Compute cost efficiency using data envelopment analysis for
inputs `X`, outputs `Y` and price of inputs `W`.

# Optional Arguments
- `rts=:VRS`: chooses variable returns to scale. For constant returns to scale choose `:CRS`.

# Examples
```jldoctest
julia> X = [5 3; 2 4; 4 2; 4 8; 7 9.0];

julia> Y = [7 4; 10 8; 8 10; 5 4; 3 6.0];

julia> W = [2 1; 2 1; 2 1; 2 1; 2 1.0];

julia> deacost(X, Y, W)
Cost DEA Model
DMUs = 5; Inputs = 2; Outputs = 2
Orientation = Input; Returns to Scale = VRS
──────────────────────────────────
       Cost  Technical  Allocative
──────────────────────────────────
1  0.615385      0.75     0.820513
2  1.0           1.0      1.0
3  1.0           1.0      1.0
4  0.5           0.5      1.0
5  0.347826      0.375    0.927536
──────────────────────────────────
```
"""
function deacost(X::Matrix, Y::Matrix, W::Matrix; rts::Symbol = :VRS)::CostDEAModel
    # Check parameters
    nx, m = size(X)
    ny, s = size(Y)

    nw, mw = size(W)

    if nx != ny
        error("number of observations is different in inputs and outputs")
    end
    if nw != nx
        error("number of observations is different in input prices and inputs")
    end
    if mw != m
        error("number of input prices and intputs is different")
    end

    # Compute efficiency for each DMU
    n = nx

    Xefficient = zeros(n,m)
    cefficiency = zeros(n)
    clambdaeff = spzeros(n, n)

    for i=1:n
        # Value of inputs and outputs to evaluate
        y0 = Y[i,:]
        w0 = W[i,:]

        # Create the optimization model
        deamodel = Model(with_optimizer(GLPK.Optimizer))
        @variable(deamodel, Xeff[1:m])
        @variable(deamodel, lambda[1:n] >= 0)

        @objective(deamodel, Min, sum(w0[j] .* Xeff[j] for j in 1:m))

        @constraint(deamodel, [j in 1:m], sum(X[t,j] * lambda[t] for t in 1:n) <= Xeff[j])
        @constraint(deamodel, [j in 1:s], sum(Y[t,j] * lambda[t] for t in 1:n) >= y0[j])

        # Add return to scale constraints
        if rts == :CRS
            # No contraint to add for constant returns to scale
        elseif rts == :VRS
            @constraint(deamodel, sum(lambda) == 1)
        else
            error("Invalid returns to scale $rts. Returns to scale should be :CRS or :VRS")
        end

        # Optimize and return results
        JuMP.optimize!(deamodel)

        Xefficient[i,:]  = JuMP.value.(Xeff)
        clambdaeff[i,:] = JuMP.value.(lambda)

    end

    # Cost, technical and allocative efficiency
    cefficiency  = vec( sum(W .* Xefficient, dims = 2) ./ sum(W .* X, dims = 2) )
    techefficiency = efficiency(dea(X, Y, orient = :Input, rts = rts, slack = false))
    allocefficiency = cefficiency ./ techefficiency
    return CostDEAModel(n, m, s, rts, cefficiency, clambdaeff, techefficiency, allocefficiency)

end

function deacost(X::Vector, Y::Matrix, W::Vector, rts::Symbol = :VRS)::CostDEAModel
    X = X[:,:]
    W = W[:,:]
    return deacost(X, Y, W, rts = rts)
end

function deacost(X::Matrix, Y::Vector, W::Matrix; rts::Symbol = :VRS)::CostDEAModel
    Y = Y[:,:]
    return deacost(X, Y, W, rts = rts)
end

function deacost(X::Vector, Y::Vector, W::Vector; rts::Symbol = :VRS)::CostDEAModel
    X = X[:,:]
    W = W[:,:]
    Y = Y[:,:]
    return deacost(X, Y, W, rts = rts)
end

function Base.show(io::IO, x::CostDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    eff = efficiency(x)
    techeff = efficiency(x, :Technical)
    alloceff = efficiency(x, :Allocative)

    if !compact
        print(io, "Cost DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Orientation = Input")
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")
        show(io, CoefTable(hcat(eff, techeff, alloceff), ["Cost", "Technical", "Allocative"], ["$i" for i in 1:n]))

    else

    end
end


text:
#=
Please read and understand the basics of perceptron learning procedure before
further looking the code here are a couple of good ways to look through resources
https://class.coursera.org/neuralnets-2012-001/lecture/10

=#
using DataArrays

function learn_perceptron(train_data_biased, target, weights)
  #=
  To evaluate the perceptron based on a calculation of x' * w (dot product)
  target is a vector containing -1 for negative examples and 1 for positive
  examples
  =#
  for i in range(1, size(train_data_biased)[1])
    activation = train_data_biased[i,:] * weights
    if !check_output(activation, target[i])
      weights = weights + target[i] * train_data_biased[i,:]'  # Update the weight vector
    end
  end
  return weights
end

function check_output(activation, real_value)
  #=
  2 correct cases:
  - activation>0 and real_value=1 ==> positive examples
  - activation<0 and real_value=-1 ==> negative examples
  Please notice that in both cases the multiplication of activation and real_value
  is positive number
  Otherwise negative number
  =#
  if activation[1] * real_value > 0
    return true
  end
    return false
end

function add_bias(train_data)
  #=
  To add bias unit to a given data set
  =#
  number_of_columns = size(train_data)[1]
  return [@data(ones(number_of_columns)) train_data]
end

function eval_perceptron(train_data_biased, weights, target)
  #=
  Count the number of errors in a given weight set and return the number of
  errors
  =#
  number_of_errors = 0
  for i in range(1, size(train_data_biased)[1])
    activation = train_data_biased[i,:] * weights
    if !check_output(activation, target[i])
      number_of_errors += 1
    end
  end
  return number_of_errors
end

max_stars_repo_name: kirnap/Machine-Learning-in-Julia

max_issues_repo_name: kirnap/Machine-Learning-in-Julia

max_forks_repo_name: kirnap/Machine-Learning-in-Julia


text:
# This file is a part of BinaryCommitteeMachineFBP.jl. License is MIT: http://github.com/carlobaldassi/BinaryCommitteeMachineFBP.jl/LICENCE.md

primitive type MagP64 <: Mag64 64 end

f2mP(a::Float64) = f2m(MagP64, a)

magformat(::Type{MagP64}) = :plain
parseinner(::Type{Val{:plain}}, s::AbstractString) = MagP64(parse(Float64, s))

convert(::Type{MagP64}, y::Float64) = f2mP(y)
convert(::Type{Float64}, y::MagP64) = m2f(y)

MagP64(y::Float64) = convert(MagP64, y)

forcedmag(::Type{MagP64}, y::Float64) = MagP64(y)

mtanh(::Type{MagP64}, x::Float64) = f2mP(tanh(x))
atanh(x::MagP64) = atanh(m2f(x))

MagP64(pp::Real, pm::Real) = MagP64((pp - pm) / (pp + pm))

isfinite(a::MagP64) = isfinite(m2f(a))

@inline function ⊗(a::MagP64, b::MagP64)
    xa = m2f(a)
    xb = m2f(b)
    # return f2mP(clamp((xa + xb) / (1 + xa * xb), -1.0, 1.0))
    # NOTE: no clamping should be required as long as the xs are in [-1,1]
    # return f2mP((xa + xb) / (1 + xa * xb))
    return f2mP((xa + xb) / muladd(xa, xb, 1.0))

end
@inline ⊘(a::MagP64, b::MagP64) = ifelse(a == b, f2mP(0.0), a ⊗ (-b))

reinforce(m0::MagP64, γ::Float64) = MagP64(tanh(atanh(m2f(m0)) * γ))

damp(newx::MagP64, oldx::MagP64, λ::Float64) = f2mP(m2f(newx) * (1 - λ) + m2f(oldx) * λ)

(*)(x::MagP64, y::MagP64) = MagP64(Float64(x) * Float64(y))

merf(::Type{MagP64}, x::Float64) = f2mP(erf(x))

function exactmix(H::MagP64, p₊::MagP64, p₋::MagP64)
    vH = m2f(H)
    pd = (m2f(p₊) + m2f(p₋)) / 2
    pz = (m2f(p₊) - m2f(p₋)) / 2

    return f2mP(pz * vH / (1 + pd * vH))
end

erfmix(H::MagP64, m₊::Float64, m₋::Float64) = MagP64(erfmix(Float64(H), m₊, m₋))

log1pxy(x::MagP64, y::MagP64) = log((1 + Float64(x) * Float64(y)) / 2)

# cross entropy with magnetizations:
#
# -(1 + x) / 2 * log((1 + y) / 2) + -(1 - x) / 2 * log((1 - y) / 2)
#
# == -x * atanh(y) - log(1 - y^2) / 2 + log(2)
function mcrossentropy(x::MagP64, y::MagP64)
    fx = m2f(x)
    fy = m2f(y)
    return -fx * atanh(fy) - log(1 - fy^2) / 2 + log(2)
end

logmag2pp(x::MagP64) = log((1 + m2f(x)) / 2)
logmag2pm(x::MagP64) = log((1 - m2f(x)) / 2)

function logZ(u0::MagP64, u::Vector{MagP64})
    zkip = logmag2pp(u0)
    zkim = logmag2pm(u0)
    for ui in u
        zkip += logmag2pp(ui)
        zkim += logmag2pm(ui)
    end
    zki = exp(zkip) + exp(zkim)
    return log(zki)
end


text:
struct CRHCCache{VN<:AbstractVector, F1<:AbstractMatrix, F2<:AbstractMatrix, V<:AbstractVector, IN<:AbstractVector}
    q::F1
    X::F1
    x::F2
    v::V
    w::V
    η::V
    u::V
    M::F1
    clusidx::IN
    clus::VN
end

function CRHCCache(X::AbstractMatrix{T1}, cl::AbstractVector{T2}) where {T1,T2}
    n, p = size(X)
    CRHCCache(similar(X), similar(X), Array{T1, 2}(undef, p, p),
             Array{T1, 1}(undef, n), Array{T1, 1}(undef, n),
             Array{T1, 1}(undef, n), Array{T1, 1}(undef, n),
             Array{T1, 2}(undef, p, p), Array{Int, 1}(undef, n),
             Array{T2, 1}(undef, n))
end

function installsortedxuw!(cache, m, k, ::Type{Val{true}})
    copyto!(cache.X, modelmatrix(m))
    copyto!(cache.u, residuals(m))
    copyto!(cache.clus, k.cl)
    if !isempty(m.rr.wts)
        cache.w .= sqrt.(smplweights(m))
        broadcast!(*, cache.u, cache.u, cache.w)
        broadcast!(*, cache.u, cache.u, cache.w)
    end

end

function installsortedxuw!(cache, m, k, ::Type{Val{false}})
    n, p = size(cache.X)
    sortperm!(cache.clusidx, k.cl)
    cidx = cache.clusidx
    u = CovarianceMatrices.residuals(m)
    w = CovarianceMatrices.smplweights(m)
    X = CovarianceMatrices.modelmatrix(m)
    uu = cache.u
    XX = cache.X
    ww = cache.w
    c  = k.cl
    cc = cache.clus
    @inbounds for i in eachindex(cidx)
        uu[i] = u[cidx[i]]
    end
    @inbounds for j in 1:p, i in eachindex(cache.clusidx)
        XX[i,j] = X[cidx[i], j]
    end
    @inbounds for i in eachindex(cache.clusidx)
        cc[i] = c[cidx[i]]
    end
    if !isempty(w)
        @inbounds for i in eachindex(cidx)
            ww[i] = sqrt(w[cidx[i]])
        end
        broadcast!(*, XX, XX, ww)
        broadcast!(*, uu, uu, ww)
    end
end

adjresid!(k::CRHC0, cache, ichol, bstarts) = nothing
adjresid!(k::CRHC1, cache, ichol, bstarts) = nothing
adjresid!(k::CRHC2, cache, ichol, bstarts) = getqii(k, cache, ichol, bstarts)
adjresid!(k::CRHC3, cache, ichol, bstarts) = getqii(k, cache, ichol, bstarts)

function dof_adjustment(cache, k::CRHC0, bstarts)
    g = length(bstarts)
    g/(g-1)
end

function dof_adjustment(cache, k::CRHC1, bstarts)
    g, (n, p) = length(bstarts), size(cache.X)
    (n-1)/(n-p) * g/(g-1)
end

dof_adjustment(cache, k::T, bstarts) where T<:Union{CRHC2} = 1.0

function dof_adjustment(cache, k::CRHC3, bstarts)
    g, (n, p) = length(bstarts), size(cache.X)
    g/(g-1)
end

function getqii(v::CRHC2, cache, A, bstarts)
    X, u  = cache.X, cache.u
    @inbounds for rnge in bstarts
        uv = view(u, rnge)
        xv = view(X, rnge, :)
        uu = copy(uv)
        xx = copy(xv)
        BB = Symmetric(I - xx*A*xx')
        uv .= cholesky!(BB).L\uu
    end
    return u
end

function getqii(v::CRHC3, cache, A, bstarts)
    X, u  = cache.X, cache.u
    @inbounds for rnge in bstarts
        uv = view(u, rnge)
        xv = view(X, rnge, :)
        uu = copy(uv)
        xx = copy(xv)
        ldiv!(cholesky!(Symmetric(I - xx*A*xx')), uu)
        uv .= uu
    end
    return u
end

function clusterize!(cache, bstarts)
    M = cache.M
    fill!(M, zero(eltype(M)))
    U = cache.q
    p, p = size(M)
    s = Array{eltype(M)}(undef, p)
    for m in bstarts
        for i = 1:p
            @inbounds s[i] = zero(Float64)
        end
        for j = 1:p, i in m
            @inbounds s[j] += U[i, j]
        end
        for j = 1:p, i = 1:p
            @inbounds M[i, j] += s[i]*s[j]
        end
    end
    return M
end


text:
# this code was closely ported from the rlpy implementaiton https://github.com/rlpy/rlpy/blob/master/rlpy/Domains/Pinball.py
# configs/*.cfg come directly from rlpy, but have been renamed making medium easy, simple_single is medium, and hard is hard.

struct PinBallConfig{T}
    start_pos::Tuple{T,T}
    target_pos::Tuple{T,T}
    target_radius::T
    ball_radius::T
    noise::T
    drag::T
    force::T

    function PinBallConfig(::Type{T}, start_pos::Tuple, target_pos::Tuple, target_radius) where {T}
        new{T}(convert.(T, start_pos), convert.(T, target_pos), T(target_radius), T(0.02), T(0.), T(0.995), T(1. / 5.))
    end
    function PinBallConfig(::Type{T}, start_pos::Tuple, target_pos::Tuple, target_radius, ball_radius) where {T}
        new{T}(convert.(T, start_pos), convert.(T, target_pos), T(target_radius), T(ball_radius), T(0.), T(0.995), T(1. / 5.))
    end
    function PinBallConfig(::Type{T}, start_pos::Tuple, target_pos::Tuple, target_radius, ball_radius, noise, drag, force) where {T}
        new{T}(convert.(T, start_pos), convert.(T, target_pos), T(target_radius), T(ball_radius), T(noise), T(drag), T(force))
    end
end

struct PinballObstacle{T} <: Any where {T}
    points::Array{Tuple{T,T}, 1}
    minx::T
    miny::T
    maxx::T
    maxy::T

    function PinballObstacle(::Type{T}, points::Array{Tuple{TX, TY}, 1}) where {T,TX,TY}
        pts = [T.(p) for p in points]
        minx, miny = min.(pts...)
        maxx, maxy = max.(pts...)
        new{T}(pts, minx, miny, maxx, maxy)
    end
end

function pinball_empty()
    return pinball_finitetime("pinball_empty.cfg", maxT=1000, stochastic_start=true, randomize=true, num_episodes=100)
end

function pinball_box()
    return pinball_finitetime("pinball_box.cfg", maxT=1000, stochastic_start=true, randomize=true, num_episodes=100)
end

function pinball_easy()
    return pinball_finitetime("pinball_easy.cfg", maxT=1000, stochastic_start=true, randomize=true, num_episodes=200)
end

function pinball_medium()
    return pinball_finitetime("pinball_medium.cfg", maxT=2000, stochastic_start=true, randomize=true, num_episodes=400, threshold=8500)
end

function pinball_hard()
    return pinball_finitetime("pinball_hard.cfg", maxT=5000, stochastic_start=true, randomize=true, num_episodes=1000, threshold=8000)
end

function pinball_finitetime(config::String; maxT=1000, stochastic_start=false, randomize=false, num_episodes=100, threshold=9000)
    X = zeros((4,2))
    X[1,:] .= [0., 1.]  # x range
	X[2,:] .= [0., 1.]  # y range
	X[3,:] .= [-2., 2.] # xdot range
    X[4,:] .= [-2., 2.] # ydot range
    S = ([0.0, maxT],
        X
    )
    A = 1:5

    obstacles, conf = read_config(Float64, config)
    if randomize
        conf = pinball_randomize(conf)
    end
    x = zeros(4)
    x .= conf.start_pos[1], conf.start_pos[2], 0.0, 0.0
    # ball = BallState(Float64, (conf.start_pos[1], conf.start_pos[2]), conf.ball_radius)
    dt = 0.05
    # x = zeros(4)

    p = (s,a)->pinball_step!(s, a, conf, obstacles, dt, stochastic_start, maxT)
    d0 = ()->pinball_d0!(conf.start_pos, stochastic_start)

    meta = Dict{Symbol,Any}()
    meta[:minreward] = -5.0
    meta[:maxreward] = 10000.0
    meta[:minreturn] = -5 * ceil(maxT / (dt * 20))  # time moves at 20*dt per step
    meta[:maxreturn] = 10000  # actually lower than this, but if you started in the goal state this would be the case.
    meta[:stochastic] = true
    meta[:minhorizon] = 40  # not sure the the true minimum is. This seems like a good lower bound
    meta[:maxhorizon] = ceil(maxT / (dt * 20))
    meta[:discounted] = false
    meta[:episodes] = num_episodes
    meta[:threshold] = threshold
    render = (state,clearplot=false)->pinballplot(state, obstacles, conf)
    m = SequentialProblem(S,X,A,p,d0,meta,render)

	return m
end

function pinball_randomize(conf::PinBallConfig)
	target_radius = conf.target_radius + rand(Uniform(-0.1, 0.1))  * conf.target_radius  # scaled target radius randomly by standard deviation of 10% of the specified radius
    ball_radius = conf.ball_radius + rand(Uniform(-0.1, 0.1)) * conf.ball_radius
    noise = rand() * 0.25    # chance for random action
    drag = 1. - exp(rand(Uniform(log(0.001), log(0.1))))  # friction coefficient for ball
    force =  rand(Uniform(0.1, 0.3))  # force applied to the ball
	conf = PinBallConfig(Float64, conf.start_pos, conf.target_pos, target_radius, ball_radius, noise, drag, force)
    return conf
end

function pinball_step!(state, action, config, obstacles, dt, stochastic_start, maxT)
    # pinball_update!($state, $action, $config, $obstacles, $dt, $maxT)
    t,x = state
    t, x, reward, γ = pinball_update!(t,x, action, config, obstacles, dt, maxT)
    if γ == 0.0
        reset_ball!(x, config.start_pos, stochastic_start)
        t = 0.0
    end

    return (t,x), x, reward, γ
    # return nothing
end

function pinball_d0!(start_pos, stochastic_start)
    x = zeros(4)
    reset_ball!(x, start_pos, stochastic_start)
    # reset_ball!(ball, start_pos, stochastic_start)
    # update_observation!(x, ball)
    t = 0.0
    return (t,x), x
end

function reset_ball!(x, start_pos::Tuple, stochastic_start::Bool)
    x[1] = start_pos[1]
    x[2] = start_pos[2]
    # ball.x = start_pos[1]
	# ball.y = start_pos[2]
	if stochastic_start
		# ball.x += 0.02 * randn()
        # ball.y += 0.02 * randn()
        x[1] += 0.02 * randn()
		x[2] += 0.02 * randn()
	end
	# ball.xDot = 0.
    # ball.yDot = 0.
    x[3] = 0.0
    x[4] = 0.0
    return nothing
end

function pinball_sim_update!(ball, dt, ball_radius, obstacles, iteration::Int)
    stepball!(ball, dt, ball_radius)

    ncollision = 0
    dxdy = (0.0, 0.0)

    for obs in obstacles
        hit, double_collision, intercept = collision(obs, ball, ball_radius)
        if hit
            dxdy = dxdy .+ collision_effect(ball, hit, double_collision, intercept)
            ncollision += 1
        end
    end
    if ncollision == 1
        # ball.xDot = dxdy[1]
        # ball.yDot = dxdy[2]
        ball[3:4] .= dxdy
        if iteration == 19
            stepball!(ball, dt, ball_radius)
        end
    elseif ncollision > 1
        ball[3] = -ball[3]
        ball[4] = -ball[4]
        # ball.xDot = -ball.xDot
        # ball.yDot = -ball.yDot
    end
    return nothing
end

function pinball_update!(t, ball, action::Int, config::PinBallConfig{T}, obstacles::Array{PinballObstacle{T},1}, dt::T, maxT) where {T}
    if action <= 0 || action > 5
        error("Action needs to be an integer in [1, 5]")
    end

	if rand() < config.noise
		action = rand(1:5)
	end
    # add action effect
    if action == 1
        add_impulse!(ball, config.force, 0.)  # Acc x
    elseif action ==2
        add_impulse!(ball, 0., -config.force) # Dec y
    elseif action ==3
        add_impulse!(ball, -config.force, 0.) # Dec x
    elseif action == 4
        add_impulse!(ball, 0., config.force) # Acc y
    else
        add_impulse!(ball, 0., 0.)  # No action
    end

    reward = 0.0
    for i in 1:20
        pinball_sim_update!(ball, dt, config.ball_radius, obstacles, i)
		t += dt
        found_goal = at_goal(ball, config)
		done = found_goal || t > maxT

        if done
            reward = (10000. * found_goal) - (1 - found_goal)

            return t,ball,reward,0.0
        end
    end

    add_drag!(ball, config.drag)
    checkbounds!(ball)

    if action == 5
        reward = -1.
    else
        reward = -5.
    end

    return t,ball,reward,1.0
end

function checkbounds!(ball)
    x = ball[1]
    y = ball[2]
    if x > 1.0
        ball[1] = 0.95
    elseif x < 0.0
        ball[1] = 0.05
    end
    if y > 1.0
        ball[2] = 0.95
    elseif y < 0.0
        ball[2] = 0.05
    end
    return nothing
end

function at_goal(ball, config::PinBallConfig)
    res = √sum(@. ((ball[1], ball[2]) - config.target_pos)^2)
    return res < config.target_radius
end

function add_impulse!(ball, Δx::T, Δy::T) where {T}
    xDot = ball[3] + Δx
    yDot = ball[4] + Δy
    ball[3] = clamp(xDot, -2., 2.)
    ball[4] = clamp(yDot, -2., 2.)
    return nothing
end

function add_drag!(ball, drag::T) where {T}
    ball[3] *= drag
    ball[4] *= drag
    return nothing
end

function stepball!(ball, dt::T, radius) where {T}
    ball[1] += ball[3] * (radius * dt)
    ball[2] += ball[4] * (radius * dt)
    return nothing
end

function read_config(::Type{T}, source) where {T}
	source = joinpath(@__DIR__, "configs", source)
    obstacles = Array{PinballObstacle{T}, 1}()
    target_pos = Tuple{T,T}((0.,0.))
    target_radius = T(0.04)
    ball_radius = T(0.02)
    start_pos = Tuple{T,T}((0.,0.))
    noise = T(0.)
    drag = T(0.995)
    force = T(1. / 5.)
    lines = readlines(source)
    for line in lines
        tokens = split(strip(line))
        if length(tokens) <= 0
            continue
        elseif tokens[1] == "polygon"
			nums = map(x->parse(T, x), tokens[2:end])
			points = [(x,y) for (x,y) in zip(nums[1:2:end], nums[2:2:end])]
            push!(obstacles, PinballObstacle(T, points))
        elseif tokens[1] == "target"
            target_pos = Tuple{T,T}(map(x->parse(T,x), tokens[2:3]))
            target_radius = parse(T,tokens[4])
        elseif tokens[1] == "start"
            start_pos = Tuple{T,T}(map(x->parse(T,x), tokens[2:3]))
        elseif tokens[1] == "ball"
            ball_radius = parse(T, tokens[2])
        end
    end

    conf = PinBallConfig(T, start_pos, target_pos, target_radius, ball_radius, noise, drag, force)

    return obstacles, conf
end

function collision(obs::PinballObstacle{T}, ball, radius) where {T}
    double_collision = false
    intercept_found = false
    intercept = ((0.0,0.0),(0.0,0.0))
    if ball[1] - radius > obs.maxx
        return intercept_found, double_collision, intercept
	end
	if ball[1] + radius < obs.minx
        return intercept_found, double_collision, intercept
	end
    if ball[2] - radius > obs.maxy
        return intercept_found, double_collision, intercept
	end
    if ball[2] + radius < obs.miny
        return intercept_found, double_collision, intercept
	end

    i = 1
    j = 2
    while i ≤ length(obs.points)
        p1, p2 = obs.points[i], obs.points[j]
        if intercept_edge(p1, p2, ball, radius)
            if intercept_found
                intercept = select_edge((p1, p2), intercept, ball)
                double_collision = true
            else
                intercept = (p1, p2)
                intercept_found = true
            end
        end
        i += 1
        j += 1
        if j > length(obs.points)
            j = 1
        end
    end
    return intercept_found, double_collision, intercept
end

function collision_effect(ball, intercept_found::Bool, double_collision::Bool, intercept::Tuple{Tuple{T,T},Tuple{T,T}})::Tuple{T,T} where {T}
    if double_collision
        return -ball[3], -ball[4]
    end

    obstacle_vector = intercept[2] .- intercept[1]
    if obstacle_vector[1] < 0.
        obstacle_vector = intercept[1] .- intercept[2]
    end

    velocity_vector = (ball[3], ball[4])
    θ = compute_angle(velocity_vector, obstacle_vector) - π
    if θ < 0.
        θ += 2π
    end

    intercept_theta = compute_angle((-1, 0), obstacle_vector)
    θ += intercept_theta

    velocity = √sum(velocity_vector.^2)
    # velocity = norm(velocity_vector)

    return velocity * cos(θ), velocity * sin(θ)
end

function compute_angle(v1, v2)
    angle_diff = atan(v1[1], v1[2]) - atan(v2[1], v2[2])
    if angle_diff < 0.
        angle_diff += 2π
    end
    return angle_diff
end

function intercept_edge(p1::Tuple{T,T}, p2::Tuple{T,T}, ball, radius) where {T}
    edge = p2 .- p1
    pball = (ball[1], ball[2])
    difference = pball .- p1

    scalar_proj = dot(difference, edge) / dot(edge, edge)
    scalar_proj = clamp(scalar_proj, 0., 1.)

    closest_pt = p1 .+ (edge .* scalar_proj)
    obstacle_to_ball = pball .- closest_pt
    distance = dot(obstacle_to_ball, obstacle_to_ball)

    if distance <= radius^2
        # collision if the ball is not moving away
        velocity = (ball[3], ball[4])
        ball_to_obstacle = closest_pt .- pball

        angle = compute_angle(ball_to_obstacle, velocity)
        if angle > π
            angle = 2π - angle
        end

        if angle > (π / 1.99)
            return false
        end
        return true
    else
        return false
    end
end

function select_edge(intersect1::Tuple{Tuple{T,T},Tuple{T,T}}, intersect2::Tuple{Tuple{T,T},Tuple{T,T}}, ball) where {T}
    velocity = (ball[3], ball[4])
    obstacle_vector1 = intersect1[2] .- intersect1[1]
    obstacle_vector2 = intersect2[2] .- intersect2[1]
    angle1 = compute_angle(velocity, obstacle_vector1)
    if angle1 > π
        angle1 -= π
    end

    angle2 = compute_angle(velocity, obstacle_vector2)
    if angle2 > π
        angle2 -= π
    end

    if abs(angle1 - π / 2.) < abs(angle2 - π / 2.)
        return intersect1
    else
        return intersect2
    end
end

@userplot PinBallPlot
@recipe function f(ap::PinBallPlot)
    state, obstacles, config = ap.args
    t, ball = state
    ballx = ball[1]
    bally = ball[2]
    bradius = config.ball_radius
	tx, ty = config.target_pos
	tr = config.target_radius

	legend := false
	xlims := (0., 1.)
	ylims := (0., 1.)
	grid := false
	ticks := nothing
	foreground_color := :white
	aspect_ratio := 1.

	# obstacles
	for ob in obstacles
		@series begin
			seriestype := :shape
			seriescolor := :blue
			xpts = [p[1] for p in ob.points]
			ypts = [p[2] for p in ob.points]

			xpts, ypts
		end
	end

	# goal
	@series begin
		seriestype := :shape
		linecolor := nothing
		seriescolor := :green
		aspect_ratio := 1.
		fillalpha := 0.4

		circle_shape(tx, ty, tr)
	end

	# ball
	@series begin
		seriestype := :shape
		linecolor := nothing
		seriescolor := :red
		aspect_ratio := 1.

		circle_shape(ballx, bally, bradius)
	end

end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl,md
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.0
#   kernelspec:
#     display_name: Julia 1.4.0
#     language: julia
#     name: julia-1.4
# ---

using PyPlot, NNlib,  ImageFiltering, Images, MEngProject, MEngProject.LamKernels, MEngProject.Laminart, MEngProject.Utils

# +
# matplotlib.rcParams["text.usetex"] = true
# -

img = convert(Array{Float64,2}, load("../input_img/Iine_100_100_gs.png"));

const σ_1 = 1
const σ_2 = 0.5
const γ = 10;
const K = 12;

# +
u_p = Laminart.I_u(img)
u_m = - u_p
v_p = fun_equ.(u_p)
v_m = fun_equ.(u_m)
v = relu.(v_p)-relu.(v_m)
V = exp(-1/8) .* imfilter(v, Kernel.gaussian(σ_2), "circular")

R = reshape(Array{eltype(V)}(undef, size(V)[1], size(V)[2]*K),size(V)[1],size(V)[2],K)
L = copy(R)

for k in 1:K
    θ = π*(k-1)/K
    R[:,:,k] = imfilter(V, relu.(LamKernels.kern_d(σ_2, θ)), "circular")
    L[:,:,k] = -imfilter(V, relu.(-LamKernels.kern_d(σ_2, θ)), "circular")
end

S_a = R .+ L
S_b = -abs.(R .- L)
S = γ .* relu.(S_a .+ S_b)
C = Laminart.fun_v_C(v_p, v_m, σ_2, K);
# -

Utils.plot_rb(img, "img", true, -1, 1, true)

Utils.plot_rb(u_p)

Utils.plot_rb(u_m)

Utils.plot_rb(v_p)

Utils.plot_rb(v_m,  "img", true, -2,2)

Utils.plot_rb(v)

Utils.plot_rb(V)

Utils.save_orientations_rb(R, "img")

# +

Utils.save_orientations_rb(L, "img")

# +

Utils.save_orientations_gs(S, "img",0,1.2,true)
# -

Utils.save_orientations_gs(C, "img", 0,2)

# +
R_ = reshape(Array{eltype(V)}(undef, size(V)[1], size(V)[2],2*K),size(V)[1],size(V)[2],2*K)
L_ = copy(R_)

for k in 1:2*K
    θ = π*(k-1)/K
    R_[:,:,k] = imfilter(V, relu.(LamKernels.kern_d(σ_2, θ)), "circular")
    L_[:,:,k] = -imfilter(V, relu.(-LamKernels.kern_d(σ_2, θ)), "circular")
end

S_a_ = R_ .+ L_
S_b_ = -abs.(R_ .- L_)
S_ = γ .* relu.((S_a_ .+ S_b_));

# +
C_ = reshape(Array{eltype(V)}(undef, size(V)[1], size(V)[2],K),size(V)[1],size(V)[2],K)

for k in 1:K
    C_[:,:,k] = S_[:,:,k] .+ S_[:,:,k+K]
end
# -

Utils.save_orientations_gs(C_, "img", 0, 2.2)

findmax(abs.(C-C_))


text:
using Test
using ModelingToolkit
import ModelingToolkit: ExtraVariablesSystemException, ExtraEquationsSystemException

include("../examples/electrical_components.jl")

function UnderdefinedConstantVoltage(;name, V = 1.0)
    val = V
    @named p = Pin()
    @named n = Pin()
    @parameters V
    eqs = [
           V ~ p.v - n.v
           # Remove equation
           # 0 ~ p.i + n.i
          ]
    ODESystem(eqs, t, [], [V], systems=[p, n], defaults=Dict(V => val), name=name)
end

function OverdefinedConstantVoltage(;name, V = 1.0, I = 1.0)
    val = V
    val2 = I
    @named p = Pin()
    @named n = Pin()
    @parameters V I
    eqs = [
           V ~ p.v - n.v
           # Overdefine p.i and n.i
           n.i ~ I
           p.i ~ I
          ]
    ODESystem(eqs, t, [], [V], systems=[p, n], defaults=Dict(V => val, I => val2), name=name)
end

R = 1.0
C = 1.0
V = 1.0
@named resistor = Resistor(R=R)
@named capacitor = Capacitor(C=C)
@named source = UnderdefinedConstantVoltage(V=V)

rc_eqs = [
          connect(source.p, resistor.p)
          connect(resistor.n, capacitor.p)
          connect(capacitor.n, source.n)
         ]

@named rc_model = ODESystem(rc_eqs, t, systems=[resistor, capacitor, source])
@test_throws ModelingToolkit.ExtraVariablesSystemException structural_simplify(rc_model)

@named source2 = OverdefinedConstantVoltage(V=V, I=V/R)
rc_eqs2 = [
          connect(source2.p, resistor.p)
          connect(resistor.n, capacitor.p)
          connect(capacitor.n, source2.n)
         ]

@named rc_model2 = ODESystem(rc_eqs2, t, systems=[resistor, capacitor, source2])
@test_throws ModelingToolkit.ExtraEquationsSystemException structural_simplify(rc_model2)

max_stars_repo_name: paulrschrater/ModelingToolsJulia

max_issues_repo_name: paulrschrater/ModelingToolsJulia

max_forks_repo_name: paulrschrater/ModelingToolsJulia


text:
#= Finite field aritmhetic implemented in unsigned integers =#
#= It is generally faster than AbstractAlgebra.GF interface =#

# returns x - y*a modulo ch
@inline function umultsubmod(
        x::T, y::T, a::T, ch::T) where {T<:Unsigned}
    ya = y*a % ch
    (x + ch - ya) % ch
end

# returns x + y*a modulo ch
@inline function umultsummod(
        x::T, y::T, a::T, ch::T) where {T<:Unsigned}
    ya = y*a % ch
    (x + ya) % ch
end

# returns x⁻¹ modulo ch
@inline function uinvmod(x::T, ch::T) where {T<:Unsigned}
    # seems like we dont need % ch
    invmod(x, ch) % ch
end

# returns -x modulo ch
@inline function ucompmod(x::T, ch::T) where {T<:Unsigned}
    (ch - x) % ch
end

# returns x*y modulo ch
@inline function umultmod(x::T, y::T, ch::T) where {T<:Unsigned}
    (x*y) % ch
end

# returns x + y modulo ch
@inline function usummod(x::T, y::T, ch::T) where {T<:Unsigned}
    (x + y) % ch
end


text:
using Quaternions
using LinearAlgebra
using QuaternionArgs
using Test

const N=100
@testset "QuaternionArgs Tests" begin
    @testset "getter" begin
        for typ in [QuaternionArgF16, QuaternionArgF32, QuaternionArgF64]
            a = rand(typ)
            @test amp(a)==a.q
            @test phase1(a)==a.phi
            @test phase2(a)==a.theta
            @test phase3(a)==a.psi
        end
    end

    for typ in [QuaternionF16, QuaternionF32, QuaternionF64]
        @testset "$typ" begin
            @testset "conversion" begin
                for i in 1:N
                    a = QuaternionArg(randn(typ))
                    @test isapprox(a,Quaternion(a))
                end
            end

            @testset "four arithmetic" begin
                for i in 1:N
                    a = randn(typ)
                    b = randn(typ)
                    ag = QuaternionArg(a)
                    bg = QuaternionArg(b)

                    anorm = normalize(ag)
                    @test amp(anorm)==one(anorm.q)
                    @test isapprox(a+b, Quaternion(ag+bg))
                    @test isapprox(a-b, Quaternion(ag-bg))
                    @test isapprox(a*b, Quaternion(ag*bg))
                    @test isapprox(a/b, Quaternion(ag/bg))
                end
            end
        end
    end
end