text:
using Random
using Distributions
using StatsBase
using LinearAlgebra
using SpecialFunctions

using DelimitedFiles

using Printf
using ProgressMeter

include("./utils.jl")
include("./load_data.jl")

###############################################################################
###############################################################################
### disease model

abstract type AbstractDiseaseModel end

"""
    DiseaseModel

define disease model that characeterises per-agent disease progression starting
from infection at time 0, to being infectious after the latent period, to begin
no longer infectious after the infectious period.
"""
struct UnivariateDiseaseModel{T1,T2} <: AbstractDiseaseModel
    P_latent::UnivariateDistribution{T1}
    P_infectious::UnivariateDistribution{T2}
end

struct DeltaDiseaseModel{T<:Real} <: AbstractDiseaseModel
    P_latent::Normal{T}
    P_infectious::Normal{T}

    function DeltaDiseaseModel{T}(mean_latent, mean_infectious) where T
        new{T}(Normal(mean_latent,0), Normal(mean_infectious, 0))
    end
end
DeltaDiseaseModel(mean_latent, mean_infectious) = DeltaDiseaseModel{Float64}(mean_latent, mean_infectious)

# default example is delta
DiseaseDefault() = DeltaDiseaseModel(seconds_from_days(3.0), seconds_from_days(5.0))

# generically distributed periods (cf. [LLoyed-Smith et al, Nature, 2001] for motivation of Gamma)
GammaDiseaseModel(mean_latent, mean_infectious, k) = UnivariateDiseaseModel(Gamma(k, mean_latent/k), Gamma(k, mean_infectious/k))

# literature based guesses for pandemic diseases
Covid19() =  DeltaDiseaseModel(NaN,NaN)
H1N1()    =  DeltaDiseaseModel(NaN,NaN)

###############################################################################
###############################################################################
### Viral load

"""
    viral_load(model, time_range; samples=Int(1e6), seed=1000)

evaluates for a disease `model` the ensemble probability to be infectious along
a `time_range` (for each step).

# Arguments
    * `model`: DiseaseModel
    * `time_range`: range over which the probability is evaluated.

# Keyword argeuments
    * `samples`: the number of disease progressions realized (ensemble size)
    * `seed`: seed of MersenneTwister random number generator
"""
function viral_load(model::AbstractDiseaseModel, time_range; samples::Int=Int(1e6), seed::Int=1000)
    print("... evaluate viral load\n")
    rng = MersenneTwister(seed)
    load = Histogram(time_range)

    # define local function
    time_to_index(time) = StatsBase.binindex(load,time)

    p = Progress(samples, 1, "Viral Load: ", offset=1)
    for i in 1:samples
        next!(p)
        index_start, index_end = time_to_index.(infectious_interval(model, rng))
        if index_end > length(load.weights)
            load.weights[index_start:index_end-1] .+= 1
        else
            load.weights[index_start:index_end] .+= 1
        end
    end

    load = float(load)
    load.weights ./= float(samples)

    return load
end

function infectious_interval(model::AbstractDiseaseModel, rng::AbstractRNG)
    time_infectious_start = rand(rng, model.P_latent)
    time_infectious_end   = time_infectious_start + rand(rng, model.P_infectious)
    return time_infectious_start, time_infectious_end
end

function infectious_interval(model::DeltaDiseaseModel)
    time_infectious_start = mean(model.P_latent)
    time_infectious_end   = time_infectious_start + mean(model.P_infectious)
    return time_infectious_start, time_infectious_end
end

###############################################################################
###############################################################################
### evaluate/sample infectious encounter

"""
generate samples only with start times that are below time_max.

This is a control.
"""
function samples_infectious_encounter_with_infection_in_early_interval(
        model::AbstractDiseaseModel,
        ets::encounter_trains{I,T},
        time_max::Number;
        #optional
        samples::Int=Int(1e6),
        rng::AbstractRNG=MersenneTwister(1000)
    ) where {I,T}
    # manipulate the cumulative sum of valid encounters to only draw from early
    # interval
    times_onset = similar(trains(ets));
    for (i,train) in enumerate(trains(ets))
        times_onset[i] = trains(ets)[i][trains(ets)[i] .< time_max];
    end
    ets_onset = encounter_trains(times_onset, ids(ets), duration(ets), timestep(ets));

    if sum(length.(ets_onset)) > 0
        data = [Int64[] for i in 1:length(ets)]
        p = Progress(samples, 1, "Sample Infectious Encounter: ", offset=1)
        for i in 1:samples
            next!(p)
            index_train, index_time = random_encoutner(rng, ets_onset)
            times = trains(ets)[index_train]
            time = times[index_time]
            interval_count = time .+ infectious_interval(model, rng)
            encounter = count_times_in_interval(times[index_time+1:end], interval_count)
            push!(data[index_train], encounter)
        end
        return data
    else
        error("The chosen interval leaves no valid disease start points")
    end
end

function samples_infectious_encounter_random_onset(
        model::AbstractDiseaseModel,
        ets::encounter_trains{I,T};
        #optional
        weighted_train::Bool=false,
        weighted_time::Bool=false,
        samples::Int=Int(1e6),
        rng::AbstractRNG=MersenneTwister(1000)
    ) where {I,T}
    support_random_time = 0:timestep(ets):duration(ets)
    if weighted_train
        weight_random_train = ProbabilityWeights(length.(ets))
    end
    if weighted_time
        encounter_rate = rate(ets, 0:timestep(ets):seconds_from_days(7))
        weight_random_time = ProbabilityWeights(vcat(encounter_rate.weights*4))
    end

    data = [Int64[] for i in 1:length(ets)]
    p = Progress(samples, 1, "Sample Infectious Encounter: ", offset=1)
    for i in 1:samples
        next!(p)
        # try several times to sample
        count_bad_tries = 0
        # first draw and fix realization of diases to not introduce a bias when
        # long diseasess are not fitting into data set
        interval_infectious = infectious_interval(model, rng)
        while true
            # draw random starting time unitl disases fits into remaining duration
            if weighted_train
                rand_train = sample(rng, 1:length(ets), weight_random_train)
            else
                rand_train = rand(rng, 1:length(ets))
            end
            if weighted_time
                rand_time = sample(rng, support_random_time, weight_random_time)
            else
                rand_time  = rand(rng, support_random_time)
            end
            train = trains(ets)[rand_train]
            interval_count = rand_time .+ interval_infectious
            if last(interval_count) < duration(ets)
                encounter = count_times_in_interval(train[searchsortedlast(train, rand_time)+1:end], interval_count)
                push!(data[rand_train], encounter)
                break
            end
            count_bad_tries += 1
            if count_bad_tries == 1000
                throw(BoundsError(ets, "infectious interval does not fit into contact trains"))
            end
        end
    end
    return data
    return samples_infectious_encounter_random_onset(model, ets, support_random_time = support_random_time, weight_random_time = weight_random_time)
end

"""
    samples_infectious_encounter(model, ets [, samples, rng])

returns a data object (A vector of length trains, each containing a vector with
samples of number of infectious encounter for a particular train). A particular
sample is obtained with `sample_infectious_encounter`.
"""
function samples_infectious_encounter(
        model::AbstractDiseaseModel,
        ets::encounter_trains{I,T};
        #optional
        samples::Int=Int(1e6),
        rng::AbstractRNG=MersenneTwister(1000)
    ) where {I,T}
    data = [Int64[] for i in 1:length(ets)]
    p = Progress(samples, 1, "Sample Infectious Encounter: ", offset=1)
    for i in 1:samples
        next!(p)
        # try several times to sample
        num_encounter, index_train = sample_infectious_encounter(ets, model, rng)
        push!(data[index_train], num_encounter)
    end
    return data
end

"""
   sample_infectious_encounter(source, infectious_interval, rng)
"""
function sample_infectious_encounter(
        ets::encounter_trains{I,T},
        model::AbstractDiseaseModel,
        rng::AbstractRNG
    ) where {I,T}
    count_bad_tries = 0
    # first draw and fix realization of diases to not introduce a bias when
    # long diseasess are not fitting into data set
    interval_infectious = infectious_interval(model, rng)
    while true
        # draw random starting time unitl disases fits into remaining duration
        index_train, index_time = random_encoutner(rng, ets)
        times = trains(ets)[index_train]
        time = times[index_time]
        interval_count = time .+ interval_infectious
        if last(interval_count) < duration(ets)
            encounter = count_times_in_interval(times[index_time+1:end], interval_count)
            return encounter, index_train
        end
        count_bad_tries += 1
        if count_bad_tries == 1000
            throw(BoundsError(ets, "infectious interval does not fit into contact trains"))
        end
    end
end

#special case for deltaDisease
function samples_infectious_encounter(
        model::DeltaDiseaseModel,
        ets::encounter_trains{I,T};
        #optional
    ) where {I,T}
    interval_infectious  = infectious_interval(model)
    train_duration = duration(ets)
    data = [Int64[] for j in 1:length(ets)]
    for (j,train) in enumerate(trains(ets))
        for (i, time_contact) in enumerate(train)
            interval_infectious_contact = time_contact .+ interval_infectious
            # counter infectious counters if remaining train fits into infectious window
            if last(interval_infectious_contact) < train_duration
                encounter = count_times_in_interval(train[i+1:end], interval_infectious_contact)
                push!(data[j], encounter)
            else # break the for loop because all disease will not fit into remaining duration for later times either
                break
            end
        end
    end
    return data
end

"""
    mean_infectious_encounter(model, contact; samples, rng)

estimatee the mean of infectious encounters for a univariate disease
`model` and a source of `contacts` (can be contact_trains object or something
else)
"""
function mean_infectious_encounter(
        model::AbstractDiseaseModel,
        ets::encounter_trains{I,T};
        #optional
        samples::Int=Int(1e6),
        rng::AbstractRNG=MersenneTwister(1000)
    ) where {I,T}
    data = samples_infectious_encounter(model, ets, samples=samples, rng=rng)
    return mean(vcat(data...))
end
function mean_infectious_encounter(
        model::DeltaDiseaseModel,
        ets::encounter_trains{I,T};
        #optional
    ) where {I,T}
    data = samples_infectious_encounter(model, ets)
    return mean_from_samples_infectious_encounter(data)
end
mean_from_samples_infectious_encounter(data) = mean(vcat(data...))

"""
    distribution_infectious_encounter(model, contact; samples, seed)

sample the distribution of infectious encounters for a univariate disease
`model` and a source of `contacts` (can be contact_trains object or something
else)
"""
function distribution_infectious_encounter(
        model::AbstractDiseaseModel,
        ets::encounter_trains{I,T};
        ### optional
        edges=missing,
        samples::Int=Int(1e6),
        rng::AbstractRNG=MersenneTwister(1000)
    ) where {I,T}
    data = samples_infectious_encounter(model, ets)
    return distribution_from_samples_infectious_encounter(data, edges=edges)
end
function distribution_infectious_encounter(
        model::DeltaDiseaseModel,
        ets::encounter_trains{I,T};
        ### optional
        edges=missing,
    ) where {I,T}
    data = samples_infectious_encounter(model, ets)
    return distribution_from_samples_infectious_encounter(data, edges=edges)
end
function distribution_from_samples_infectious_encounter(
        data;
        ### optional
        edges=missing,
    )
    if ismissing(edges)
        list_encounter = vcat(data...)
        dist = fit(Histogram{Float64}, list_encounter, 0:1:maximum(list_encounter)+2)
    else
        dist = fit(Histogram{Float64}, vcat(data...), edges)
    end
    normalize!(dist)
    return dist
end

function count_times_in_interval(times::Vector{T1}, interval_count::Tuple{T2,T2}) where {T1<:Number,T2<:Number}
    encounter = 0
    for time in times
        encounter += ininterval(time, interval_count)
    end
    return encounter
end
function ininterval(time::T1, interval_count::Tuple{T2,T2}) where {T1<:Number,T2<:Number}
    return (getindex(interval_count,1) <= time < getindex(interval_count,2))
end


text:
using Plots
pyplot()
#using Random
using Distributions
using Optim

taskNr = 3;

####### Code used for several tasks #######
function myNormal(x,m,s)
    # PDF of Gaussian (watch out for elementwise . !)
    y = 1/sqrt(2*pi*s^2) * exp.(-(x.-m).^2 ./ (2*s^2))

    return y
end

# Parameters:
N = 10000;
mu = 3.4;
sigma = sqrt(3);

# Draw samples from built-in normal distribution
r = rand(Normal(mu,sigma), N);
###################################

if taskNr == 1
    # Draw a histogram, normalised
    histogram(r, normalize = :pdf,
        label = string(N, " Samples"))

    # Draw the PDF
    lo = minimum(r);
    hi = maximum(r);
    symRange = maximum([abs(lo)-mu, abs(hi)-mu]);
    x = range(mu-symRange, mu+symRange, length = 100);
    plot!(x, myNormal(x,mu,sigma),
        w = 2,
        label = "PDF")

    savefig("TaskA-2-a.eps")

elseif taskNr == 2
    # Determine granularity and create x and y:
    dx = .1;
    x = collect(-10:dx:20);
    p = myNormal(x,mu,sigma);

    # Riemann sums of
    # Expected value (weigh elementwise by x):
    E_Riemann = sum(p .* x) * dx;
    # Variance (weigh square-distance elementwise by x):
    variance_Riemann = sum(p .* (x.-mu).^2) * dx;

    # Check closeness:
    println(string("Distance of calculated E from true value: ",
        mu - E_Riemann))
    println(string("Distance of calculated E from sample value: ",
        mean(r) - E_Riemann))
    println(string("Distance of calculated Variance from true value ",
        sigma^2 - variance_Riemann))
    println(string("Distance of calculated Variance from sample value ",
        var(r) - variance_Riemann))
    # Note: var(r) automatically scales with n-1 instead of n

elseif taskNr == 3
    # Parameters:
    N = 10000;
    mu = 0;
    sigma = 1;

    # Draw samples from built-in normal distribution
    x = rand(Normal(mu,sigma), N);

    # Make a log normally distributed vector
    y = exp.(x);

    # Draw a histogram, normalised
    histogram(y,
        normalize = :pdf,
        bins = 0:.1:maximum(y),
        w = .1,     # line width of bin edges
        label = string(N, " Samples"))

    # transform the normal pdf to log-normal pdf:
    pLog(x,m,s) = myNormal(log.(x),m,s) ./ x

    # Draw the PDF
    lo = minimum(y);
    hi = maximum(y);
    y_pdf = collect(lo:.1:hi);
    z = pLog(y_pdf,mu,sigma)
    plot!(y_pdf, z,
        w = 1.5,
        label = "PDF")

    savefig("TaskA-2-d-i.eps")

    # Modes using the points in the PDF:
    @time begin
        pdfMax = maximum(z);
        maxIdx = findfirst(z .== pdfMax);
    end

    modeManual = y_pdf[maxIdx]
    println(string("The manual mode is @ y = ", modeManual))

    # Modes using optimisation:
    myObjective(a) = -pLog(a[1],mu,sigma)    # fix the mean and variance
    @time begin
        optim = optimize(myObjective,[0.0, 0.0],NelderMead());
    end

    modeOptim = optim.minimizer[1];
    println(string("The optimised mode is @ y = ", modeOptim))

end

ext: jl


text:
# Test simple evaluations:
options = Options(;
    default_params..., binary_operators=(+, *, /, -), unary_operators=(cos, sin)
)

# Here, we unittest the fast function evaluation scheme
# We need to trigger all possible fused functions, with all their logic.
# These are as follows:

## We fuse (and compile) the following:
##  - op(op2(x, y)), where x, y, z are constants or variables.
##  - op(op2(x)), where x is a constant or variable.
##  - op(x), for any x.
## We fuse (and compile) the following:
##  - op(x, y), where x, y are constants or variables.
##  - op(x, y), where x is a constant or variable but y is not.
##  - op(x, y), where y is a constant or variable but x is not.
##  - op(x, y), for any x or y
for fnc in [
    # deg2_l0_r0_eval
    (x1, x2, x3) -> x1 * x2,
    (x1, x2, x3) -> x1 * 3.0f0,
    (x1, x2, x3) -> 3.0f0 * x2,
    (((x1, x2, x3) -> 3.0f0 * 6.0f0), ((x1, x2, x3) -> Node(3.0f0) * 6.0f0)),
    # deg2_l0_eval
    (x1, x2, x3) -> x1 * sin(x2),
    (x1, x2, x3) -> 3.0f0 * sin(x2),

    # deg2_r0_eval
    (x1, x2, x3) -> sin(x1) * x2,
    (x1, x2, x3) -> sin(x1) * 3.0f0,

    # deg1_l2_ll0_lr0_eval
    (x1, x2, x3) -> cos(x1 * x2),
    (x1, x2, x3) -> cos(x1 * 3.0f0),
    (x1, x2, x3) -> cos(3.0f0 * x2),
    (((x1, x2, x3) -> cos(3.0f0 * -0.5f0)), ((x1, x2, x3) -> cos(Node(3.0f0) * -0.5f0))),

    # deg1_l1_ll0_eval
    (x1, x2, x3) -> cos(sin(x1)),
    (((x1, x2, x3) -> cos(sin(3.0f0))), ((x1, x2, x3) -> cos(sin(Node(3.0f0))))),

    # everything else:
    (x1, x2, x3) -> (sin(cos(sin(cos(x1) * x3) * 3.0f0) * -0.5f0) + 2.0f0) * 5.0f0,
]

    # check if fnc is tuple
    if typeof(fnc) <: Tuple
        realfnc = fnc[1]
        nodefnc = fnc[2]
    else
        realfnc = fnc
        nodefnc = fnc
    end

    global tree = nodefnc(Node("x1"), Node("x2"), Node("x3"))

    N = 100
    nfeatures = 3
    X = randn(MersenneTwister(0), Float32, nfeatures, N)

    test_y = eval_tree_array(tree, X, options)[1]
    true_y = realfnc.(X[1, :], X[2, :], X[3, :])

    zero_tolerance = 1e-6
    @test all(abs.(test_y .- true_y) / N .< zero_tolerance)
end


text:
module MyExample

greet()=println("This is my example")
greet2()=println("greet2")
greet3()=println("greet3")

include("extrafile.jl")

export my_f

end

max_stars_repo_name: giandopal/MyExample

max_issues_repo_name: giandopal/MyExample

max_forks_repo_name: giandopal/MyExample


text:
import Win32API
using Base.Test

@test 1 == Win32API.MessageBoxA(C_NULL, "world!", "Hello", 0)
@test 1 == Win32API.MessageBoxW(C_NULL, [0x3044, 0x0000], [0x3042, 0x0000], 0)


text:
# Note to devs. Use GLPK or Cbc for models with linear constraints and linear cost functions
# Use OSQP for models with quadratic cost function and linear constraints and ipopt otherwise
@testset "All PowerModels models construction" begin
    networks = [
        (PM.ACPPowerModel, fast_ipopt_optimizer),
        (PM.ACRPowerModel, fast_ipopt_optimizer),
        (PM.ACTPowerModel, fast_ipopt_optimizer),
        #(PM.IVRPowerModel, fast_ipopt_optimizer), #instantiate_ivp_expr_model not implemented
        (PM.DCPPowerModel, fast_ipopt_optimizer),
        (PM.DCMPPowerModel, fast_ipopt_optimizer),
        (PM.NFAPowerModel, fast_ipopt_optimizer),
        (PM.DCPLLPowerModel, fast_ipopt_optimizer),
        (PM.LPACCPowerModel, fast_ipopt_optimizer),
        (PM.SOCWRPowerModel, fast_ipopt_optimizer),
        (PM.SOCWRConicPowerModel, scs_solver),
        (PM.QCRMPowerModel, fast_ipopt_optimizer),
        (PM.QCLSPowerModel, fast_ipopt_optimizer),
        #(PM.SOCBFPowerModel, fast_ipopt_optimizer), # not implemented
        (PM.BFAPowerModel, fast_ipopt_optimizer),
        #(PM.SOCBFConicPowerModel, fast_ipopt_optimizer), # not implemented
        (PM.SDPWRMPowerModel, scs_solver),
        #(PM.SparseSDPWRMPowerModel, scs_solver), # bug in PM: https://github.com/lanl-ansi/PowerModels.jl/issues/769, uncomment when v0.18.1 is tagged
        (PTDFPowerModel, fast_ipopt_optimizer),
    ]
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    for (network, solver) in networks
        template = get_thermal_dispatch_template_network(network)
        ps_model = OperationsProblem(template, c_sys5; optimizer = solver)
        @test build!(ps_model; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT
        @test !isnothing(ps_model.internal.optimization_container.pm)
        @test :nodal_balance_active in
              keys(ps_model.internal.optimization_container.expressions)
    end
end

@testset "Network Copper Plate" begin
    template = get_thermal_dispatch_template_network(CopperPlatePowerModel)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    c_sys14_dc = PSB.build_system(PSITestSystems, "c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    parameters = [true, false]
    test_results = IdDict{System, Vector{Int}}(
        c_sys5 => [120, 0, 120, 120, 24],
        c_sys14 => [120, 0, 120, 120, 24],
        c_sys14_dc => [120, 0, 120, 120, 24],
    )
    constraint_names = [:CopperPlateBalance]
    objfuncs = [GAEVF, GQEVF, GQEVF]
    test_obj_values = IdDict{System, Float64}(
        c_sys5 => 240000.0,
        c_sys14 => 142000.0,
        c_sys14_dc => 142000.0,
    )

    for (ix, sys) in enumerate(systems), p in parameters
        ps_model =
            OperationsProblem(template, sys; optimizer = OSQP_optimizer, use_parameters = p)

        @test build!(ps_model; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT
        psi_constraint_test(ps_model, constraint_names)
        moi_tests(
            ps_model,
            p,
            test_results[sys][1],
            test_results[sys][2],
            test_results[sys][3],
            test_results[sys][4],
            test_results[sys][5],
            false,
        )
        psi_checkobjfun_test(ps_model, objfuncs[ix])
        psi_checksolve_test(ps_model, [MOI.OPTIMAL], test_obj_values[sys], 10000)
    end
    ps_model_re = OperationsProblem(
        template,
        PSB.build_system(PSITestSystems, "c_sys5_re");
        optimizer = GLPK_optimizer,
        use_parameters = true,
        balance_slack_variables = true,
    )
    @test build!(ps_model_re; output_dir = mktempdir(cleanup = true)) ==
          PSI.BuildStatus.BUILT
    psi_checksolve_test(ps_model_re, [MOI.OPTIMAL], 240000.0, 10000)
end

@testset "Network DC-PF with PTDF Model" begin
    template = get_thermal_dispatch_template_network(StandardPTDFModel)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    c_sys14_dc = PSB.build_system(PSITestSystems, "c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    objfuncs = [GAEVF, GQEVF, GQEVF]
    constraint_names =
        [:RateLimit_lb__Line, :RateLimit_ub__Line, :CopperPlateBalance, :network_flow__Line]
    parameters = [true, false]
    PTDF_ref = IdDict{System, PTDF}(
        c_sys5 => PTDF(c_sys5),
        c_sys14 => PTDF(c_sys14),
        c_sys14_dc => PTDF(c_sys14_dc),
    )
    test_results = IdDict{System, Vector{Int}}(
        c_sys5 => [264, 0, 264, 264, 168],
        c_sys14 => [600, 0, 600, 600, 504],
        c_sys14_dc => [600, 48, 552, 552, 456],
    )
    test_obj_values = IdDict{System, Float64}(
        c_sys5 => 340000.0,
        c_sys14 => 142000.0,
        c_sys14_dc => 142000.0,
    )
    for (ix, sys) in enumerate(systems), p in parameters
        ps_model = OperationsProblem(
            template,
            sys;
            optimizer = OSQP_optimizer,
            use_parameters = p,
            PTDF = PTDF_ref[sys],
        )

        @test build!(ps_model; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT
        psi_constraint_test(ps_model, constraint_names)
        moi_tests(
            ps_model,
            p,
            test_results[sys][1],
            test_results[sys][2],
            test_results[sys][3],
            test_results[sys][4],
            test_results[sys][5],
            false,
        )
        psi_checkobjfun_test(ps_model, objfuncs[ix])
        psi_checksolve_test(
            ps_model,
            [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL],
            test_obj_values[sys],
            10000,
        )
    end
    # PTDF input Error testing
    ps_model = OperationsProblem(template, c_sys5; optimizer = GLPK_optimizer)
    @test build!(
        ps_model;
        console_level = Logging.AboveMaxLevel,  # Ignore expected errors.
        output_dir = mktempdir(cleanup = true),
    ) == PSI.BuildStatus.FAILED
end

@testset "Sparse Network DC-PF with PTDFPowerModel" begin
    template = get_thermal_dispatch_template_network(PTDFPowerModel)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    c_sys14_dc = PSB.build_system(PSITestSystems, "c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    objfuncs = [GAEVF, GQEVF, GQEVF]
    constraint_names =
        [:RateLimit_lb__Line, :RateLimit_ub__Line, :CopperPlateBalance, :network_flow__Line]
    parameters = [true, false]
    test_results = IdDict{System, Vector{Int}}(
        c_sys5 => [264, 0, 264, 264, 168],
        c_sys14 => [600, 0, 600, 600, 504],
        c_sys14_dc => [600, 48, 552, 552, 456],
    )
    test_obj_values = IdDict{System, Float64}(
        c_sys5 => 340000.0,
        c_sys14 => 142000.0,
        c_sys14_dc => 142000.0,
    )
    for (ix, sys) in enumerate(systems), p in parameters
        ps_model =
            OperationsProblem(template, sys; optimizer = OSQP_optimizer, use_parameters = p)

        @test build!(ps_model; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT
        psi_constraint_test(ps_model, constraint_names)
        moi_tests(
            ps_model,
            p,
            test_results[sys][1],
            test_results[sys][2],
            test_results[sys][3],
            test_results[sys][4],
            test_results[sys][5],
            false,
        )
        psi_checkobjfun_test(ps_model, objfuncs[ix])
        psi_checksolve_test(
            ps_model,
            [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL],
            test_obj_values[sys],
            10000,
        )
    end
end

@testset "Network DC lossless -PF network with PowerModels DCPlosslessForm" begin
    template = get_thermal_dispatch_template_network(DCPPowerModel)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    c_sys14_dc = PSB.build_system(PSITestSystems, "c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    objfuncs = [GAEVF, GQEVF, GQEVF]
    constraint_names = [
        :RateLimit_ub__Line,
        :RateLimit_lb__Line,
        PSI.make_constraint_name(PSI.NODAL_BALANCE_ACTIVE, PSY.Bus),
    ]
    parameters = [true, false]
    test_results = IdDict{System, Vector{Int}}(
        c_sys5 => [384, 0, 408, 408, 288],
        c_sys14 => [936, 0, 1080, 1080, 840],
        c_sys14_dc => [984, 48, 984, 984, 840],
    )
    test_obj_values = IdDict{System, Float64}(
        c_sys5 => 342000.0,
        c_sys14 => 142000.0,
        c_sys14_dc => 142000.0,
    )
    for (ix, sys) in enumerate(systems), p in parameters
        ps_model =
            OperationsProblem(template, sys; optimizer = OSQP_optimizer, use_parameters = p)
        @test build!(ps_model; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT
        psi_constraint_test(ps_model, constraint_names)
        moi_tests(
            ps_model,
            p,
            test_results[sys][1],
            test_results[sys][2],
            test_results[sys][3],
            test_results[sys][4],
            test_results[sys][5],
            false,
        )
        psi_checkobjfun_test(ps_model, objfuncs[ix])
        psi_checksolve_test(
            ps_model,
            [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL],
            test_obj_values[sys],
            1000,
        )
    end
end

@testset "Network Solve AC-PF PowerModels StandardACPModel" begin
    template = get_thermal_dispatch_template_network(ACPPowerModel)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    c_sys14_dc = PSB.build_system(PSITestSystems, "c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    objfuncs = [GAEVF, GQEVF, GQEVF]
    # Check for voltage and angle constraints
    constraint_names = [
        :RateLimitFT__Line,
        :RateLimitTF__Line,
        PSI.make_constraint_name(PSI.NODAL_BALANCE_ACTIVE, PSY.Bus),
        PSI.make_constraint_name(PSI.NODAL_BALANCE_REACTIVE, PSY.Bus),
    ]
    parameters = [true, false]
    test_results = IdDict{System, Vector{Int}}(
        c_sys5 => [1056, 0, 384, 384, 264],
        c_sys14 => [2832, 0, 720, 720, 696],
        c_sys14_dc => [2832, 96, 672, 672, 744],
    )
    test_obj_values = IdDict{System, Float64}(
        c_sys5 => 340000.0,
        c_sys14 => 142000.0,
        c_sys14_dc => 142000.0,
    )
    for (ix, sys) in enumerate(systems), p in parameters
        ps_model = OperationsProblem(
            template,
            sys;
            optimizer = ipopt_optimizer,
            use_parameters = p,
        )
        @test build!(ps_model; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT
        psi_constraint_test(ps_model, constraint_names)
        moi_tests(
            ps_model,
            p,
            test_results[sys][1],
            test_results[sys][2],
            test_results[sys][3],
            test_results[sys][4],
            test_results[sys][5],
            false,
        )
        psi_checkobjfun_test(ps_model, objfuncs[ix])
        psi_checksolve_test(
            ps_model,
            [MOI.TIME_LIMIT, MOI.OPTIMAL, MOI.LOCALLY_SOLVED],
            test_obj_values[sys],
            10000,
        )
    end
end

# TODO: Add constraint tests for these models, other is redundant with first test
@testset "Network Solve AC-PF PowerModels NFAPowerModel" begin
    template = get_thermal_dispatch_template_network(NFAPowerModel)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    c_sys14_dc = PSB.build_system(PSITestSystems, "c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    objfuncs = [GAEVF, GQEVF, GQEVF]
    constraint_names = [PSI.make_constraint_name(PSI.NODAL_BALANCE_ACTIVE, PSY.Bus)]
    parameters = [true, false]
    test_results = Dict{System, Vector{Int}}(
        c_sys5 => [264, 0, 264, 264, 120],
        c_sys14 => [600, 0, 600, 600, 336],
        c_sys14_dc => [648, 48, 552, 552, 384],
    )
    test_obj_values = IdDict{System, Float64}(
        c_sys5 => 300000.0,
        c_sys14 => 142000.0,
        c_sys14_dc => 142000.0,
    )
    for (ix, sys) in enumerate(systems), p in parameters
        ps_model =
            OperationsProblem(template, sys; optimizer = OSQP_optimizer, use_parameters = p)
        @test build!(ps_model; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT
        psi_constraint_test(ps_model, constraint_names)
        moi_tests(
            ps_model,
            p,
            test_results[sys][1],
            test_results[sys][2],
            test_results[sys][3],
            test_results[sys][4],
            test_results[sys][5],
            false,
        )
        psi_checkobjfun_test(ps_model, objfuncs[ix])
        psi_checksolve_test(
            ps_model,
            [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL],
            test_obj_values[sys],
            10000,
        )
    end
end

# TODO: Add constraint tests for these models, other is redundant with first test
@testset "Other Network AC PowerModels models" begin
    networks = [#ACPPowerModel, Already tested
        ACRPowerModel,
        ACTPowerModel,
    ]
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    c_sys14_dc = PSB.build_system(PSITestSystems, "c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    parameters = [true, false]
    # TODO: add model specific constraints to this list. Voltages, etc.
    constraint_names = [
        PSI.make_constraint_name(PSI.NODAL_BALANCE_ACTIVE, PSY.Bus),
        PSI.make_constraint_name(PSI.NODAL_BALANCE_REACTIVE, PSY.Bus),
    ]
    ACR_test_results = Dict{System, Vector{Int}}(
        c_sys5 => [1056, 0, 240, 240, 264],
        c_sys14 => [2832, 0, 240, 240, 696],
        c_sys14_dc => [2832, 96, 240, 240, 744],
    )
    ACT_test_results = Dict{System, Vector{Int}}(
        c_sys5 => [1344, 0, 384, 384, 840],
        c_sys14 => [3792, 0, 720, 720, 2616],
        c_sys14_dc => [3696, 96, 672, 672, 2472],
    )
    test_results = Dict(zip(networks, [ACR_test_results, ACT_test_results]))
    for network in networks, sys in systems, p in parameters
        template = get_thermal_dispatch_template_network(network)
        ps_model = OperationsProblem(
            template,
            sys;
            optimizer = fast_ipopt_optimizer,
            use_parameters = p,
        )
        @test build!(ps_model; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT
        psi_constraint_test(ps_model, constraint_names)
        moi_tests(
            ps_model,
            p,
            test_results[network][sys][1],
            test_results[network][sys][2],
            test_results[network][sys][3],
            test_results[network][sys][4],
            test_results[network][sys][5],
            false,
        )
        @test !isnothing(ps_model.internal.optimization_container.pm)
    end
end

# TODO: Add constraint tests for these models, other is redundant with first test
@testset "Network DC-PF PowerModels quadratic loss approximations models" begin
    networks = [DCPLLPowerModel, LPACCPowerModel]
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    c_sys14_dc = PSB.build_system(PSITestSystems, "c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    parameters = [true, false]
    # TODO: add model specific constraints to this list. Bi-direccional flows etc
    constraint_names = [PSI.make_constraint_name(PSI.NODAL_BALANCE_ACTIVE, PSY.Bus)]
    test_obj_values = IdDict{System, Float64}(
        c_sys5 => 340000.0,
        c_sys14 => 142000.0,
        c_sys14_dc => 142000.0,
    )
    DCPLL_test_results = Dict{System, Vector{Int}}(
        c_sys5 => [528, 0, 408, 408, 288],
        c_sys14 => [1416, 0, 1080, 1080, 840],
        c_sys14_dc => [1416, 48, 984, 984, 840],
    )
    LPACC_test_results = Dict{System, Vector{Int}}(
        c_sys5 => [1200, 0, 384, 384, 840],
        c_sys14 => [3312, 0, 720, 720, 2616],
        c_sys14_dc => [3264, 96, 672, 672, 2472],
    )
    test_results = Dict(zip(networks, [DCPLL_test_results, LPACC_test_results]))
    for network in networks, sys in systems, p in parameters
        template = get_thermal_dispatch_template_network(network)
        ps_model = OperationsProblem(
            template,
            sys;
            optimizer = ipopt_optimizer,
            use_parameters = p,
        )
        @test build!(ps_model; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT
        psi_constraint_test(ps_model, constraint_names)
        moi_tests(
            ps_model,
            p,
            test_results[network][sys][1],
            test_results[network][sys][2],
            test_results[network][sys][3],
            test_results[network][sys][4],
            test_results[network][sys][5],
            false,
        )
        @test !isnothing(ps_model.internal.optimization_container.pm)
        psi_checksolve_test(
            ps_model,
            [MOI.OPTIMAL, MOI.LOCALLY_SOLVED],
            test_obj_values[sys],
            10000,
        )
    end
end

@testset "Network Unsupported Power Model Formulations" begin
    for network in PSI.UNSUPPORTED_POWERMODELS
        template = get_thermal_dispatch_template_network(network)
        ps_model = OperationsProblem(
            template,
            PSB.build_system(PSITestSystems, "c_sys5");
            optimizer = ipopt_optimizer,
        )
        @test build!(
            ps_model;
            console_level = Logging.AboveMaxLevel,  # Ignore expected errors.
            output_dir = mktempdir(cleanup = true),
        ) == PSI.BuildStatus.FAILED
    end
end


text:
module FinancialBlotter

# documentation
using Docile

if VERSION < v"0.4-"
  using Dates
else
  using Base.Dates
end
using Reexport
@reexport using TimeSeries, FinancialSeries

import Base: show, push! #getindex, add!

@document

export Blotter, Portfolio, Account, Trade, tradearray,
       blottercolnames, blotterticker,
       push!,#add!,
       merge, parsedatetime, parsedatetime1, parsedatetime_from_TOS, makedatetime, datetolastsecond, discretesignal,
       startedtrades, closedtrades

# Phase 2 methods and const
#        OrderBook,
#        orderbookcolnames, orderbookbidvalues, orderbookoffervalues, orderbooksellvalues, orderbookcovervalues, orderbookticker,
#        fillorderbook, fillblotter

include("blotter.jl")
include("show.jl")

### TODO
#include("trades.jl")
#include("portfolio.jl")
#include("account.jl")
#include("getindex.jl")
#include("utilities.jl")

# include("Phase2/blotter1.jl")
# include("Phase2/orderbook.jl")
# include("Phase2/trades.jl")
# include("Phase2/statemachine.jl")

end


text:
# Not in FiniteDiff because `u` -> scalar isn't used anywhere else,
# but could be upstreamed.
mutable struct UGradientWrapper{fType,tType,P} <: Function
  f::fType
  t::tType
  p::P
end

(ff::UGradientWrapper)(uprev) = ff.f(uprev,ff.p,ff.t)

mutable struct ParamGradientWrapper{fType,tType,uType} <: Function
  f::fType
  t::tType
  u::uType
end

(ff::ParamGradientWrapper)(p) = ff.f(ff.u,p,ff.t)

# the next four definitions are only needed in case of non-diagonal SDEs

mutable struct ParamNonDiagNoiseGradientWrapper{fType,tType,uType} <: Function
  f::fType
  t::tType
  u::uType
end

(ff::ParamNonDiagNoiseGradientWrapper)(p) = vec(ff.f(ff.u,p,ff.t))

mutable struct ParamNonDiagNoiseJacobianWrapper{fType,tType,uType,duType} <: Function
  f::fType
  t::tType
  u::uType
  du::duType
end

function (ff::ParamNonDiagNoiseJacobianWrapper)(p)
  du1 = similar(p, size(ff.du))
  ff.f(du1,ff.u,p,ff.t)
  return vec(du1)
end

function (ff::ParamNonDiagNoiseJacobianWrapper)(du1,p)
  ff.f(du1,ff.u,p,ff.t)
  return vec(du1)
end

mutable struct UNonDiagNoiseGradientWrapper{fType,tType,P} <: Function
  f::fType
  t::tType
  p::P
end

(ff::UNonDiagNoiseGradientWrapper)(uprev) = vec(ff.f(uprev,ff.p,ff.t))

mutable struct UNonDiagNoiseJacobianWrapper{fType,tType,P,duType} <: Function
  f::fType
  t::tType
  p::P
  du::duType
end

(ff::UNonDiagNoiseJacobianWrapper)(uprev) = (du1 = similar(ff.du); ff.f(du1,uprev,ff.p,ff.t); vec(du1))

function (ff::UNonDiagNoiseJacobianWrapper)(du1,uprev)
  ff.f(du1,uprev,ff.p,ff.t)
  return vec(du1)
end

Base.@pure function determine_chunksize(u,alg::DiffEqBase.AbstractSensitivityAlgorithm)
  determine_chunksize(u,get_chunksize(alg))
end

Base.@pure function determine_chunksize(u,CS)
  if CS != 0
    return CS
  else
    return ForwardDiff.pickchunksize(length(u))
  end
end

function jacobian(f, x::AbstractArray{<:Number}, alg::DiffEqBase.AbstractSensitivityAlgorithm)
  if alg_autodiff(alg)
    J = ForwardDiff.jacobian(f, x)
  else
    J = FiniteDiff.finite_difference_jacobian(f, x)
  end
  return J
end

function jacobian!(J::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number},
                   fx::Union{Nothing,AbstractArray{<:Number}}, alg::DiffEqBase.AbstractSensitivityAlgorithm, jac_config)
  if alg_autodiff(alg)
    if fx === nothing
      ForwardDiff.jacobian!(J, f, x)
    else
      ForwardDiff.jacobian!(J, f, fx, x, jac_config)
    end
  else
    FiniteDiff.finite_difference_jacobian!(J, f, x, jac_config)
  end
  nothing
end

function gradient!(df::AbstractArray{<:Number}, f,
                   x::Union{Number,AbstractArray{<:Number}},
                   alg::DiffEqBase.AbstractSensitivityAlgorithm, grad_config)
    if alg_autodiff(alg)
        ForwardDiff.gradient!(df, f, x, grad_config)
    else
        FiniteDiff.finite_difference_gradient!(df, f, x, grad_config)
    end
    nothing
end

"""
  jacobianvec!(Jv, f, x, v, alg, (buffer, seed)) -> nothing

``Jv <- J(f(x))v``
"""
function jacobianvec!(Jv::AbstractArray{<:Number}, f, x::AbstractArray{<:Number},
                      v, alg::DiffEqBase.AbstractSensitivityAlgorithm, config)
  if alg_autodiff(alg)
    buffer, seed = config
    TD = typeof(first(seed))
    T  = typeof(first(seed).partials)
    @. seed = TD(x, T(tuple(v)))
    f(buffer, seed)
    Jv .= ForwardDiff.partials.(buffer, 1)
  else
    buffer1, buffer2 = config
    f(buffer1,x)
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    @. x += ϵ*v
    f(buffer2,x)
    @. x -= ϵ*v
    @. du = (buffer2 - buffer1)/ϵ
  end
  nothing
end

function vecjacobian!(dλ, y, λ, p, t, S::SensitivityFunction;
                      dgrad=nothing, dy=nothing)
  _vecjacobian!(dλ, y, λ, p, t, S, S.sensealg.autojacvec, dgrad, dy)
  return
end

function _vecjacobian!(dλ, y, λ, p, t, S::SensitivityFunction, isautojacvec::Bool, dgrad, dy)
  @unpack sensealg, f = S
  prob = getprob(S)
  if isautojacvec isa Bool && !isautojacvec
    @unpack J, uf, f_cache, jac_config = S.diffcache
    if !(prob isa DiffEqBase.SteadyStateProblem)
      if DiffEqBase.has_jac(f)
        f.jac(J,y,p,t) # Calculate the Jacobian into J
      else
        uf.t = t
        uf.p = p
        jacobian!(J, uf, y, f_cache, sensealg, jac_config)
      end
      mul!(dλ',λ',J)
    end
    if dgrad !== nothing
      @unpack pJ, pf, paramjac_config = S.diffcache
      if DiffEqBase.has_paramjac(f)
        # Calculate the parameter Jacobian into pJ
        f.paramjac(pJ,y,prob.p,t)
      else
        pf.t = t
        pf.u = y
        if DiffEqBase.isinplace(prob)
          jacobian!(pJ, pf, prob.p, f_cache, sensealg, paramjac_config)
        else
          temp = jacobian(pf, prob.p, sensealg)
          pJ .= temp
        end
      end
      mul!(dgrad',λ',pJ)
    end
    if dy !== nothing
      if DiffEqBase.isinplace(prob)
         f(dy, y, p, t)
      else
        dy[:] .= vec(f(y, p, t))
      end
    end
  elseif DiffEqBase.isinplace(prob)
    _vecjacobian!(dλ, y, λ, p, t, S, ReverseDiffVJP(), dgrad, dy)
  else
    _vecjacobian!(dλ, y, λ, p, t, S, ZygoteVJP(), dgrad, dy)
  end
  return
end

function _vecjacobian!(dλ, y, λ, p, t, S::SensitivityFunction, isautojacvec::TrackerVJP, dgrad, dy)
  @unpack sensealg, f = S
  prob = getprob(S)
  isautojacvec = get_jacvec(sensealg)
  if DiffEqBase.isinplace(prob)
    _dy, back = Tracker.forward(y, prob.p) do u, p
      out_ = map(zero, u)
      f(out_, u, p, t)
      Tracker.collect(out_)
    end
    tmp1, tmp2 = Tracker.data.(back(λ))
    dλ[:] .= vec(tmp1)
    dgrad !== nothing && (dgrad[:] .= vec(tmp2))
    dy !== nothing && (dy[:] .= vec(Tracker.data(_dy)))
  else
    _dy, back = Tracker.forward(y, prob.p) do u, p
      Tracker.collect(f(u, p, t))
    end
    tmp1, tmp2 = Tracker.data.(back(λ))
    dλ[:] .= vec(tmp1)
    dgrad !== nothing && (dgrad[:] .= vec(tmp2))
    dy !== nothing && (dy[:] .= vec(Tracker.data(_dy)))
  end
  return
end

function _vecjacobian!(dλ, y, λ, p, t, S::SensitivityFunction, isautojacvec::ReverseDiffVJP, dgrad, dy)
  @unpack sensealg, f = S
  prob = getprob(S)
  isautojacvec = get_jacvec(sensealg)

  if typeof(prob) <: SteadyStateProblem || (eltype(λ) <: eltype(prob.u0) && typeof(t) <: eltype(prob.u0) && compile_tape(sensealg.autojacvec))
    tape = S.diffcache.paramjac_config

  ## These other cases happen due to autodiff in stiff ODE solvers
  elseif DiffEqBase.isinplace(prob)
    _y = eltype(y) === eltype(λ) ? y : convert.(promote_type(eltype(y),eltype(λ)),y)
    tape = ReverseDiff.GradientTape((_y, prob.p, [t])) do u,p,t
      du1 = similar(u, size(u))
      f(du1,u,p,first(t))
      return vec(du1)
     end
  else
    _y = eltype(y) === eltype(λ) ? y : convert.(promote_type(eltype(y),eltype(λ)),y)
    tape = ReverseDiff.GradientTape((_y, prob.p, [t])) do u,p,t
      vec(f(u,p,first(t)))
    end
  end

  if prob isa DiffEqBase.SteadyStateProblem
    tu, tp = ReverseDiff.input_hook(tape)
  else
    tu, tp, tt = ReverseDiff.input_hook(tape)
  end
  output = ReverseDiff.output_hook(tape)
  ReverseDiff.unseed!(tu) # clear any "leftover" derivatives from previous calls
  ReverseDiff.unseed!(tp)
  if !(prob isa DiffEqBase.SteadyStateProblem)
    ReverseDiff.unseed!(tt)
  end
  ReverseDiff.value!(tu, y)
  ReverseDiff.value!(tp, prob.p)
  if !(prob isa DiffEqBase.SteadyStateProblem)
    ReverseDiff.value!(tt, [t])
  end
  ReverseDiff.forward_pass!(tape)
  ReverseDiff.increment_deriv!(output, λ)
  ReverseDiff.reverse_pass!(tape)
  copyto!(vec(dλ), ReverseDiff.deriv(tu))
  dgrad !== nothing && copyto!(vec(dgrad), ReverseDiff.deriv(tp))
  ReverseDiff.pull_value!(output)
  dy !== nothing && copyto!(vec(dy), ReverseDiff.value(output))
  return
end

function _vecjacobian!(dλ, y, λ, p, t, S::SensitivityFunction, isautojacvec::ZygoteVJP, dgrad, dy)
  @unpack sensealg, f = S
  prob = getprob(S)

  isautojacvec = get_jacvec(sensealg)
  if DiffEqBase.isinplace(prob)
    _dy, back = Zygote.pullback(y, prob.p) do u, p
      out_ = Zygote.Buffer(similar(u))
      f(out_, u, p, t)
      vec(copy(out_))
    end
    tmp1,tmp2 = back(λ)
    dλ[:] .= vec(tmp1)
    dgrad !== nothing && tmp2 != nothing && (dgrad[:] .= vec(tmp2))
    dy !== nothing && (dy[:] .= vec(_dy))
  else
    _dy, back = Zygote.pullback(y, prob.p) do u, p
      vec(f(u, p, t))
    end
    tmp1,tmp2 = back(λ)
    dλ[:] .= vec(tmp1)

    dy !== nothing && (dy[:] .= vec(_dy))
    dgrad !== nothing && tmp2 != nothing && (dgrad[:] .= vec(tmp2))
  end
  return
end

function jacNoise!(λ, y, p, t, S::SensitivityFunction;
                      dgrad=nothing, dλ=nothing, dy=nothing)
  _jacNoise!(λ, y, p, t, S, S.sensealg.noise, dgrad, dλ, dy)
  return
end

function _jacNoise!(λ, y, p, t, S::SensitivityFunction, isnoise::Bool, dgrad, dλ, dy)
  @unpack sensealg, f = S
  prob = getprob(S)
  if isnoise isa Bool && !isnoise
    if dgrad !== nothing
      @unpack pJ, pf, f_cache, paramjac_noise_config = S.diffcache
      if DiffEqBase.has_paramjac(f)
        # Calculate the parameter Jacobian into pJ
        f.paramjac(pJ,y,prob.p,t)
      else
        pf.t = t
        pf.u = y
        if DiffEqBase.isinplace(prob)
          jacobian!(pJ, pf, prob.p, nothing, sensealg, nothing)
          #jacobian!(pJ, pf, prob.p, f_cache, sensealg, paramjac_noise_config)
        else
          temp = jacobian(pf, prob.p, sensealg)
          pJ .= temp
        end
      end

      if StochasticDiffEq.is_diagonal_noise(prob)
        pJt = transpose(λ).*transpose(pJ)
        dgrad[:] .= vec(pJt)
      else
        m = size(prob.noise_rate_prototype)[2]
        for i in 1:m
          tmp = λ'*pJ[(i-1)*m+1:i*m,:]
          dgrad[:,i] .= vec(tmp)
        end
      end
    end

    if dλ !== nothing && (isnoisemixing(sensealg) || !StochasticDiffEq.is_diagonal_noise(prob))
      @unpack J, uf, f_cache, jac_noise_config = S.diffcache
      if dy!=nothing
        if DiffEqBase.isinplace(prob)
          f(dy, y, p, t)
        else
          dy .= f(y, p, t)
        end
      end

      if DiffEqBase.has_jac(f)
        f.jac(J,y,p,t) # Calculate the Jacobian into J
      else
        if DiffEqBase.isinplace(prob)
          if dy != nothing
            ForwardDiff.jacobian!(J,uf,dy,y)
          else
            if StochasticDiffEq.is_diagonal_noise(prob)
              dy = similar(y)
            else
              dy = similar(prob.noise_rate_prototype)
              f(dy, y, p, t)
              ForwardDiff.jacobian!(J,uf,dy,y)
            end
            f(dy, y, p, t)
            ForwardDiff.jacobian!(J,uf,dy,y)
          end
        else
          tmp = ForwardDiff.jacobian(uf,y)
          J .= tmp
        end
      #  uf.t = t
      #  uf.p = p
      #  jacobian!(J, uf, y, nothing, sensealg, nothing)
      end

      if StochasticDiffEq.is_diagonal_noise(prob)
        Jt = transpose(λ).*transpose(J)
        dλ[:] .= vec(Jt)
      else
        for i in 1:m
          tmp = λ'*J[(i-1)*m+1:i*m,:]
          dλ[:,i] .= vec(tmp)
        end
      end
    end
  elseif DiffEqBase.isinplace(prob)
    _jacNoise!(λ, y, p, t, S, ReverseDiffNoise(), dgrad, dλ, dy)
  else
    _jacNoise!(λ, y, p, t, S, ZygoteNoise(), dgrad, dλ, dy)
  end
  return
end

function _jacNoise!(λ, y, p, t, S::SensitivityFunction, isnoise::ReverseDiffNoise, dgrad, dλ, dy)
  @unpack sensealg, f = S
  prob = getprob(S)

  for (i, λi) in enumerate(λ)
    tapei = S.diffcache.paramjac_noise_config[i]
    tu, tp, tt = ReverseDiff.input_hook(tapei)
    output = ReverseDiff.output_hook(tapei)
    ReverseDiff.unseed!(tu) # clear any "leftover" derivatives from previous calls
    ReverseDiff.unseed!(tp)
    ReverseDiff.unseed!(tt)
    ReverseDiff.value!(tu, y)
    ReverseDiff.value!(tp, p)
    ReverseDiff.value!(tt, [t])
    ReverseDiff.forward_pass!(tapei)
    if StochasticDiffEq.is_diagonal_noise(prob)
      ReverseDiff.increment_deriv!(output, λi)
    else
      ReverseDiff.increment_deriv!(output, λ)
    end
    ReverseDiff.reverse_pass!(tapei)

    deriv = ReverseDiff.deriv(tp)
    dgrad[:,i] .= vec(deriv)
    ReverseDiff.pull_value!(output)

    if StochasticDiffEq.is_diagonal_noise(prob)
      dλ !== nothing && (dλ[:,i] .= vec(ReverseDiff.deriv(tu)))
      dy !== nothing && (dy[i] = ReverseDiff.value(output))
    else
      dλ !== nothing && (dλ[:,i] .= vec(ReverseDiff.deriv(tu)))
      dy !== nothing && (dy[:,i] .= vec(ReverseDiff.value(output)))
    end
  end
  return
end

function _jacNoise!(λ, y, p, t, S::SensitivityFunction, isnoise::ZygoteNoise, dgrad, dλ, dy)
  @unpack sensealg, f = S
  prob = getprob(S)

  if StochasticDiffEq.is_diagonal_noise(prob)
    if DiffEqBase.isinplace(prob)
      for (i, λi) in enumerate(λ)
        _dy, back = Zygote.pullback(y, prob.p) do u, p
          out_ = Zygote.Buffer(similar(u))
          f(out_, u, p, t)
          copy(out_[i])
        end
        tmp1,tmp2 = back(λi) #issue: tmp2 = zeros(p)
        dgrad[:,i] .= vec(tmp2)
        dλ !== nothing && (dλ[:,i] .= vec(tmp1))
        dy !== nothing && (dy[i] = _dy)
      end
    else
      for (i, λi) in enumerate(λ)
        _dy, back = Zygote.pullback(y, prob.p) do u, p
          f(u, p, t)[i]
        end
        tmp1,tmp2 = back(λi)
        dgrad[:,i] .= vec(tmp2)
        dλ !== nothing && (dλ[:,i] .= vec(tmp1))
        dy !== nothing && (dy[i] = _dy)
      end
    end
  else
    if DiffEqBase.isinplace(prob)
      for (i, λi) in enumerate(λ)
        _dy, back = Zygote.pullback(y, prob.p) do u, p
          out_ = Zygote.Buffer(similar(prob.noise_rate_prototype))
          f(out_, u, p, t)
          copy(out_[:,i])
        end
        tmp1,tmp2 = back(λ)#issue with Zygote.Buffer
        dgrad[:,i] .= vec(tmp2)
        dλ !== nothing && (dλ[:,i] .= vec(tmp1))
        dy !== nothing && (dy[:,i] .= vec(_dy))
      end
    else
      for (i, λi) in enumerate(λ)
         _dy, back = Zygote.pullback(y, prob.p) do u, p
          f(u, p, t)[:,i]
        end
        tmp1,tmp2 = back(λ)
        dgrad[:,i] .= vec(tmp2)
        dλ !== nothing && (dλ[:,i] .= vec(tmp1))
        dy !== nothing && (dy[:,i] .= vec(_dy))
      end
    end
  end
  return
end

function accumulate_cost!(dλ, y, p, t, S::SensitivityFunction, dgrad=nothing)
  @unpack dg, dg_val, g, g_grad_config = S.diffcache
  if dg != nothing
    if !(dg isa Tuple)
      dg(dg_val,y,p,t)
      dλ .+= vec(dg_val)
    else
      dg[1](dg_val[1],y,p,t)
      dλ .+= vec(dg_val[1])
      if dgrad !== nothing
        dg[2](dg_val[2],y,p,t)
        dgrad .-= vec(dg_val[2])
      end
    end
  else
    g.t = t
    gradient!(dg_val, g, y, S.sensealg, g_grad_config)
    dλ .+= vec(dg_val)
  end
  return nothing
end

function build_jac_config(alg,uf,u)
  if alg_autodiff(alg)
    jac_config = ForwardDiff.JacobianConfig(uf,u,u,
                 ForwardDiff.Chunk{determine_chunksize(u,alg)}())
  else
    if diff_type(alg) != Val{:complex}
      jac_config = FiniteDiff.JacobianCache(similar(u),similar(u),
                                                 similar(u),diff_type(alg))
    else
      tmp = Complex{eltype(u)}.(u)
      du1 = Complex{eltype(u)}.(du1)
      jac_config = FiniteDiff.JacobianCache(tmp,du1,nothing,diff_type(alg))
    end
  end
  jac_config
end

function build_param_jac_config(alg,uf,u,p)
  if alg_autodiff(alg)
    jac_config = ForwardDiff.JacobianConfig(uf,u,p,
                 ForwardDiff.Chunk{determine_chunksize(p,alg)}())
  else
    if diff_type(alg) != Val{:complex}
      jac_config = FiniteDiff.JacobianCache(similar(p),similar(u),
                                                 similar(u),diff_type(alg))
    else
      tmp = Complex{eltype(p)}.(p)
      du1 = Complex{eltype(u)}.(u)
      jac_config = FiniteDiff.JacobianCache(tmp,du1,nothing,diff_type(alg))
    end
  end
  jac_config
end

function build_grad_config(alg,tf,du1,t)
  if alg_autodiff(alg)
    grad_config = ForwardDiff.GradientConfig(tf,du1,
                                    ForwardDiff.Chunk{determine_chunksize(du1,alg)}())
  else
    grad_config = FiniteDiff.GradientCache(du1,t,diff_type(alg))
  end
  grad_config
end


text:
# This is the code that handles getting the JLSO file itself on to and off of the disk
# However, it does not describe how to serialize or deserialize the individual objects
# that is done lazily and the code for that is in serialization.jl
function Base.write(io::IO, jlso::JLSOFile)
    _versioncheck(jlso.version, WRITEABLE_VERSIONS)

    # Setup an IOBuffer for serializing the manifest
    project_toml = sprint(Pkg.TOML.print, jlso.project)
    # @show jlso.manifest
    manifest_toml = read(
        compress(jlso.compression, IOBuffer(sprint(Pkg.TOML.print, jlso.manifest)))
    )

    # We declare the dict types to be more explicit about the output format.
    metadata = Dict{String, Union{String, Vector{UInt8}}}(
        "version" => string(jlso.version),
        "julia" => string(jlso.julia),
        "format" => string(jlso.format),
        "compression" => string(jlso.compression),
        "image" => jlso.image,
        "project" => project_toml,
        "manifest" => manifest_toml,
    )

    if jlso.version >= v"4"
        # NOTE: We intentionally use a hypthen in the object names and nbytes to make it
        # harder to accidentally accessing them outside of the BSON context

        doc = Dict{String, Union{Dict, Vector}}(
            "metadata" => metadata,
            "object-names" => string.(keys(jlso.objects)),
            "object-nbytes" => Int64.(length.(values(jlso.objects))),
            "object-pos" => zeros(Int64, length(jlso.objects))
        )

        # Early exit condition if there are no objects
        isempty(jlso.objects) && return BSON.bson(io, doc)

        buf = IOBuffer()

        # Allocate our BSON doc so we can determine object positions in the file
        BSON.bson(buf, doc)

        # Determine each objects starting position and update the doc dict
        # Calculate an offset that can be applied to the nbytes array
        offsets = cumsum([0, doc["object-nbytes"][1:end-1]...])
        doc["object-pos"] = Int64.(position(buf) .+ offsets)

        # Now we can write the doc to our actual IO
        BSON.bson(io, doc)

        # Now write our raw data bytes to the end of the IO
        for obj in values(jlso.objects)
            write(io, obj)
        end
    else
        BSON.bson(
            io,
            Dict{String, Dict}(
                "metadata" => metadata,
                "objects" => Dict(string(k) => v for (k, v) in jlso.objects),
            )
        )
    end
end

# `read`, unlike `load` does not deserialize any of the objects within the JLSO file,
# they will be `deserialized` when they are indexed out of the returned JSLOFile object.
function Base.read(io::IO, ::Type{JLSOFile})
    parsed = BSON.load(io)
    _versioncheck(parsed["metadata"]["version"], READABLE_VERSIONS)
    d = upgrade_jlso(parsed)
    compression = Symbol(d["metadata"]["compression"])

    # Try decompressing the manifest, otherwise just return a dict with the raw data
    manifest = try
        Pkg.TOML.parse(
            read(
                decompress(compression, IOBuffer(d["metadata"]["manifest"])),
                String
            )
        )
    catch e
        warn(LOGGER, e)
        Dict("raw" => d["metadata"]["manifest"])
    end

    return JLSOFile(
        VersionNumber(d["metadata"]["version"]),
        VersionNumber(d["metadata"]["julia"]),
        Symbol(d["metadata"]["format"]),
        compression,
        d["metadata"]["image"],
        Pkg.TOML.parse(d["metadata"]["project"]),
        manifest,
        _extract_objects(io, d),
        ReentrantLock()
    )
end

"""
    save(io, data)
    save(path, data)

Creates a JLSOFile with the specified data and kwargs and writes it back to the io.
"""
save(io::IO, data::JLSOFile) = write(io, data)
save(io::IO, data; kwargs...) = save(io, JLSOFile(data; kwargs...))
save(io::IO, data::Pair...; kwargs...) = save(io, Dict(data...); kwargs...)
function save(path::Union{AbstractPath, AbstractString}, args...; kwargs...)
    return open(io -> save(io, args...; kwargs...), path, "w")
end

"""
    load(io, objects...) -> Dict{Symbol, Any}
    load(path, objects...) -> Dict{Symbol, Any}

Load the JLSOFile from the io and deserialize the specified objects.
If no object names are specified then all objects in the file are returned.
"""
load(path::Union{AbstractPath, AbstractString}, args...) = open(io -> load(io, args...), path)
function load(io::IO, objects::Symbol...)
    jlso = read(io, JLSOFile)
    objects = isempty(objects) ? names(jlso) : objects
    result = Dict{Symbol, Any}()

    @sync for o in objects
        @spawn begin
            # Note that calling getindex on the jlso triggers the deserialization of the object
            deserialized = jlso[o]
            lock(jlso.lock) do
                result[o] = deserialized
            end
        end
    end

    return result
end

# FileIO Interface
fileio_save(f, args...; kwargs...) = save(f.filename, args...; kwargs...)
fileio_load(f, args...) = load(f.filename, args...)

# Objects were initialially saved inline with the BSON document prior to v4
# To support larger objects/files we started writing objects to the end of file and only save
# names and sizes to the objects dict in the bson doc.
# Since our ideal structure is one dict with object bytes we don't solve this using the
# `upgrade` function, but rather extract both formats here.
function _extract_objects(io::IO, d)
    if VersionNumber(d["metadata"]["version"]) < v"4"
        return Dict{Symbol, Vector{UInt8}}(Symbol(k) => v for (k, v) in d["objects"])
    else
        objects = Dict{Symbol, Vector{UInt8}}()

        for (k, pos, nb) in zip(d["object-names"], d["object-pos"], d["object-nbytes"])
            seek(io, pos)
            bytes = read(io, nb)
            length(bytes) < nb && throw(EOFError())
            objects[Symbol(k)] = bytes
        end

        return objects
    end
end

max_stars_repo_name: invenia/JLSO

max_issues_repo_name: invenia/JLSO

max_forks_repo_name: invenia/JLSO


text:
# broadcasting

using Base.Broadcast: BroadcastStyle, Broadcasted

struct ROCArrayStyle{N} <: AbstractGPUArrayStyle{N} end
ROCArrayStyle(::Val{N}) where N = ROCArrayStyle{N}()
ROCArrayStyle{M}(::Val{N}) where {N,M} = ROCArrayStyle{N}()

BroadcastStyle(::Type{<:ROCArray{T,N}}) where {T,N} = ROCArrayStyle{N}()

Base.similar(bc::Broadcasted{ROCArrayStyle{N}}, ::Type{T}) where {N,T} =
    similar(ROCArray{T}, axes(bc))

Base.similar(bc::Broadcasted{ROCArrayStyle{N}}, ::Type{T}, dims...) where {N,T} =
    ROCArray{T}(undef, dims...)

# broadcasting type ctors isn't GPU compatible
Broadcast.broadcasted(::ROCArrayStyle{N}, f::Type{T}, args...) where {N, T} =
    Broadcasted{ROCArrayStyle{N}}((x...) -> T(x...), args, nothing)


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

Copyright (c) 2015 Peter Kovesi
pk@peterkovesi.com

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

The Software is provided "as is", without warranty of any kind.

----------------------------------------------------------------------------=#

import ColorTypes, Images

export applycolourmap, applydivergingcolourmap, applycycliccolourmap, ternaryimage
export applycolormap, applydivergingcolormap, applycycliccolormap

#-------------------------------------------------------------------------------
"""
applycolourmap/applycolormap: Applies colourmap to a single channel
image to obtain an RGB result

```
Usage: rgbimg = applycolourmap(img, cmap, rnge)

Arguments:  img - Single channel image to apply colourmap to.
                  ::ImageMeta{T,2} or ::Array{Float64,2}
           cmap - RGB colourmap as generated by cmap().
                  ::Array{ColorTypes.RGBA{Float64},1}
           rnge - Optional 2-vector specifying the min and max values in
                  the image to be mapped across the colour map.  Values
                  outside this range are mapped to the end points of the
                  colour map.  If range is omitted the full range
                  of image values are used.

Returns: rgbimg - RGB image of floating point values in the range 0-1.
                  NaN values in the input image are mapped to black.
                  ::ImageMeta{Float64,3} or ::Array{Float64,3}
```
Why use this function when you can simply set a colour map?

Well, actually you probably want to use the functions
applycycliccolourmap() and applydivergingcolourmap() which make use of
this function.

Many visualisation packages may automatically apply an offset and
perform some scaling of your data to normalise it to a range of, say,
0-255 before applying a colour map and rendering it on your screen.
In many cases this is useful. However, if you are wanting to render
your data with a diverging or cyclic colour map then this behaviour is
definitely not appropriate because these types of colour maps requires
that data values are honoured in some way to make any sense.

By providing a 'range' parameter this function allows you to apply a
colour map in a way that respects your data values.

See also: cmap, applycycliccolourmap, applydivergingcolourmap
"""
function applycolourmap(imgin::Array{T,2}, cmap::Array{ColorTypes.RGBA{Float64},1}, rnge::Array) where T
    # This function is also used by relief() as a base image upon which to
    # apply relief shading.

    # November  2013 - Original Matlab code
    # November  2015 - Ported to Julia

    ncolours = length(cmap)
    @assert ndims(imgin) == 2   "Image must be single channel"
    @assert rnge[1] < rnge[2] "rnge[1] must be less than rnge[2]"

    # Convert cmap to a ncolours x 3 array for convenience
    rgbmap = RGBA2FloatArray(cmap)

    img = copy(imgin)
    # Set any Nan entries to rnge[1] so that they will be mapped to the first
    # colour in the colour map without an error being thrown.  Later all NaN
    # regions are set to black.
    mask = isnan.(img)
    img[mask] .= rnge[1]

    # Convert image values to integers that can be used to index into colourmap
    rimg = round.(Int16, (img.-rnge[1])/(rnge[2]-rnge[1]) * (ncolours-1) ) .+ 1

    # Clamp any out of range values (note the values specified in range[] may
    # be well inside the range of values in the image)
    rimg[rimg .< 1] .= 1
    rimg[rimg .> ncolours] .= ncolours

    (rows,cols) = size(img)
    rgbimg = zeros(rows,cols,3)

    for r = 1:rows, c = 1:cols
        rgbimg[r,c,:] = rgbmap[rimg[r,c],:] * !mask[r,c]
    end

    return rgbimg
end

# Case when range is not supplied: use minimum and maximum of image.
# AbstractArray captures both Image and Array{Float64,2}
function applycolourmap(img::AbstractArray, cmap::AbstractArray)
    return applycolourmap(img, cmap, float([minimum(img) maximum(img)]))
end

# Case for img::ImageMeta
function applycolourmap(img::ImageMeta{T,2}, cmap, rnge) where T
    rgbimg = applycolourmap(float(Images.data(img)), cmap, rnge)
    return ImageMeta(rgbimg, colordim = 3, colorspace=Images.RGB,
                        spatialorder = img.properties["spatialorder"])
end

function applycolourmap(img::AbstractArray{T,2}, cmap::Array{Float64, 2}, rnge) where T<:Real
    return applycolourmap(img, FloatArray2RGBA(cmap), rnge)
end

# For those who spell colour without a 'u'
"""
applycolormap: Applies colourmap to a single channel image to obtain
an RGB result
```
Usage: rgbimg = applycolormap(img, cmap, rnge)

Arguments:  img - Single channel image to apply colourmap to.
                  ::ImageMeta{T,2} or ::Array{Float64,2}
           cmap - RGB colourmap as generated by cmap().
                  ::Array{ColorTypes.RGBA{Float64},1}
           rnge - Optional 2-vector specifying the min and max values in
                  the image to be mapped across the colour map.  Values
                  outside this range are mapped to the end points of the
                  colour map.  If range is omitted the full range
                  of image values are used.

Returns: rgbimg - RGB image of floating point values in the range 0-1.
                  NaN values in the input image are rendered as black.
                  ::ImageMeta{Float64,3} or ::Array{Float64,3}

For full documentation see applycolourmap
                                    ^
```
See also: cmap, applycycliccolourmap, applydivergingcolourmap
"""
function applycolormap(imgin::Array{Float64,2}, cmap::Array{ColorTypes.RGBA{Float64},1}, rnge::Array)
    applycolourmap(imgin, cmap, rnge)
end

# Case when range is not supplied: use minimum and maximum of image.
# AbstractArray captures both Image and Array{Float64,2}
function applycolormap(img::AbstractArray, cmap::AbstractArray)
    return applycolourmap(img, cmap, float([minimum(img) maximum(img)]))
end

# Case for img::ImageMeta
function applycolormap(img::ImageMeta{T,2}, cmap, rnge) where T
    rgbimg = applycolourmap(float(Images.data(img)), cmap, rnge)
    return ImageMeta(rgbimg, colordim = 3, colorspace=Images.RGB,
                     spatialorder = img.properties["spatialorder"])
end

function applycolormap(img::AbstractArray{T,2}, cmap::Array{Float64, 2}, rnge) where T<:Real
    return applycolourmap(img, FloatArray2RGBA(cmap), rnge)
end

#------------------------------------------------------------------
"""
ternaryimage:  Perceptualy uniform ternary image from 3 bands of data

This function generates a ternary image using 3 basis colours that are
closely matched in lightness, as are their secondary colours.  The
colours are not as vivid as the RGB primaries but they produce ternary
images with consistent feature salience no matter what permutation of
channel-colour assignement is used.  This is in contrast to ternary
images constructed with the RGB primaries where the channel that
happens to get encoded in green dominates the perceptual result.

Useful for Landsat imagery or radiometric images.

```
Usage: rgbimg = ternaryimage(img; bands, histcut, RGB)

Argument:
            img - Multiband image with at least 3 bands.
                  ::ImageMeta{T,3} or ::Array{T<:Real,3}

Keyword Arguments:
          bands - Array of 3 values indicating the bands, to be assigned to
                  the red, green and blue basis colour maps.  If omitted
                  bands defaults to [1, 2, 3].
        histcut - Percentage of image band histograms to clip.  It can be
                  useful to clip 1-2%. If you see lots of white in your
                  ternary image you have clipped too much. Defaults to 0.
            RGB - Boolean flag, if set to true the classical RGB primaries
                  are used to construct the ternary image rather than the
                  lightness matched primaries. Defaults to false.

Returns:
          rgbimg - RGB ternary image
                  ::ImageMeta{T,3} or ::Array{T<:Real,3}
```
For the derivation of the three primary colours see:
Peter Kovesi. Good Colour Maps: How to Design Them.
arXiv:1509.03700 [cs.GR] 2015.

See also: applycolourmap, linearrgbmap
"""
function ternaryimage(img::Array{T,3}; bands::Array = [1, 2, 3],
                      histcut::Real = 0.0, RGB::Bool=false) where T<:Real

    N = 256
    (rows, cols, chan) = size(img)

    if minimum(bands) < 1 || maximum(bands) > chan
        error("Band specification outside number of image channels")
    end

    if RGB          # Use classical RGB primaries
        R = [1 0 0]
        G = [0 1 0]
        B = [0 0 1]

    else           # Use lightness matched primaries.
                   # For their derivation see the reference above.
        R = [0.90 0.17 0.00]
        G = [0.00 0.50 0.00]
        B = [0.10 0.33 1.00]
    end

    Rmap = equalisecolourmap("rgb", linearrgbmap(R, N))
    Gmap = equalisecolourmap("rgb", linearrgbmap(G, N))
    Bmap = equalisecolourmap("rgb", linearrgbmap(B, N))

    if histcut > eps()
        rgbimg = applycolourmap(histtruncate(img[:,:,bands[1]], histcut), Rmap) +
        applycolourmap(histtruncate(img[:,:,bands[2]], histcut), Gmap) +
        applycolourmap(histtruncate(img[:,:,bands[3]], histcut), Bmap)
    else
        rgbimg = applycolourmap(img[:,:,bands[1]], Rmap) +
        applycolourmap(img[:,:,bands[2]], Gmap) +
        applycolourmap(img[:,:,bands[3]], Bmap)
    end

    return rgbimg
end

# Case for img::ImageMeta  ** Conversion of RGB4 FixedPointNumbers images does not work **
function ternaryimage(img::ImageMeta{T,3};
                      bands::Array = [1, 2, 3], histcut::Real = 0.0, RGB::Bool=false) where T

    rgbimg = ternaryimage(float(Images.data(img)), bands=bands, histcut=histcut, RGB=RGB)
    return ImageMeta(rgbimg, colordim = 3, colorspace=Images.RGB, spatialorder = img.properties["spatialorder"])
end

#------------------------------------------------------------------------------
"""
applycycliccolourmap/applycycliccolormap: Applies a cyclic colour map
to an image of angular data

For angular data to be rendered correctly it is important that the data values
are respected so that data values are correctly assigned to specific entries
in a cyclic colour map.  The assignment of values to colours also depends on
whether the data is cyclic over pi, or 2*pi.

In contrast, default display methods typically do not respect data values
directly and can perform inappropriate offsetting and normalisation of the
angular data before display and rendering with a colour map.

The rendering of the angular data with a specified colour map can be modulated
as a function of an associated image amplitude.  This allows the colour map
encoding of the angular information to be modulated to represent the
amplitude/reliability/coherence of the angular data.

```
Usage: rgbimg = applycycliccolourmap(ang, cmap)
       rgbimg = applycycliccolourmap(ang, cmap, keyword args ...)

Arguments:
           ang - Image of angular data to be rendered
                 ::ImageMeta or ::Array{Float64,2}
          cmap - Cyclic colour map to render the angular data with.

Keyword arguments:

           amp - Amplitude image used to modulate the mapped colours of the
                 angular data.  If not supplied no modulation of colours is
                 performed.
                 ::ImageMeta or ::Array{Float64,2}
   modtoblack  - Boolean flag/1 indicating whether the amplitude image is used to
                 modulate the colour mapped image values towards black,
                 or towards white.  The default is true, towards black.
   cyclelength - The cycle length of the angular data.  Use a value of pi
                 if the data represents orientations, or 2*pi if the data
                 represents phase values.  If the input data is in degrees
                 simply set cycle in degrees and the data will be
                 rendered appropriately. Default is 2*pi.

Returns: rgbim - The rendered image.
                 ::ImageMeta{Float64,3} or ::Array{Float64,3}
```

For a list of all cyclic colour maps that can be generated by cmap() use:

```
> cmap("cyclic")
```

See also: cmap, scalogram, ridgeorient, applycolourmap, applydivergingcolourmap
"""
function applycycliccolourmap(ang::Array{Float64,2}, cmap::Array{ColorTypes.RGBA{Float64},1};
                  amp::Array=Array{Float64}(undef, 0,0), cyclelength::Real=2*pi, modtoblack::Bool=true)
    # September 2014
    # November  2015 Ported to Julia

    # Apply colour map to angular data.  Some care is needed with this.  Unlike
    # normal 'linear' data one cannot apply shifts and/or rescale values to
    # normalise them.  The raw angular data values have to be respected.

    angmod = mod.(ang, cyclelength)   # Ensure data values are within range 0 - cyclelength
    rgbimg = applycolourmap(angmod, cmap, [0, cyclelength])

    if !isempty(amp)   # Display image with rgb values modulated by amplitude
        normamp = normalise(amp)  # Enforce amplitude  0 - 1

        if modtoblack  # Modulate rgb values by amplitude fading to black
            for n = 1:3
                rgbimg[:,:,n] = rgbimg[:,:,n].*normamp
            end

        else           # Modulate rgb values by amplitude fading to white
            for n = 1:3
                rgbimg[:,:,n] = 1 - (1 - rgbimg[:,:,n]).*normamp
            end
        end
    end

    return rgbimg
end

# Case for img::ImageMeta
function applycycliccolourmap(img::ImageMeta{T1,2}, cmap::Array{ColorTypes.RGBA{Float64},1};
                          amp::ImageMeta{T2,2}=ImageMeta(Array{Float64}(undef, 0,0)),
                          cyclelength::Real=2*pi, modtoblack::Bool=true) where {T1,T2}
    rgbimg = applycycliccolourmap(float(Images.data(img)), cmap, amp=float(Images.data(amp)),
                   cyclelength=cyclelength, modtoblack=modtoblack)
    return ImageMeta(rgbimg, colordim = 3, colorspace=Images.RGB, spatialorder = img.properties["spatialorder"])
end

# For those who spell colour without a 'u'
"""
applycycliccolormap:  Applies a cyclic colour map to an image of angular data
```
For full documentation see applycycliccolourmap
                                          ^
```
See also: applycolourmap, applydivergingcolourmap
"""
function applycycliccolormap(ang::Array{Float64,2}, cmap::Array{ColorTypes.RGBA{Float64},1};
                  amp::Array=Array{Float64}(undef, 0,0), cyclelength::Real=2*pi, modtoblack::Bool=true)
    applycycliccolourmap(ang, cmap, amp, cyclelength, modtoblack)
end

# Case for img::ImageMeta
function applycycliccolormap(img::ImageMeta{T1,2}, cmap::Array{ColorTypes.RGBA{Float64},1};
                             amp::ImageMeta{T2,2}=ImageMeta(Array{Float64}(undef, 0,0)),
                             cyclelength::Real=2*pi, modtoblack::Bool=true) where {T1,T2}
    applycycliccolourmap(img, cmap, amp, cyclelength, modtoblack)
end

#------------------------------------------------------------------------------
"""
applydivergingcolourmap/applydivergingcolormap - Applies a diverging
colour map to an image

For data to be displayed correctly with a diverging colour map it is
important that the data values are respected so that the reference value in
the data is correctly associated with the centre entry of a diverging
colour map.

In contrast, default display methods typically do not respect data values
directly and can perform inappropriate offsetting and normalisation of the
data before display and rendering with a colour map.

```
Usage:  rgbim = applydivergingcolourmap(img, map, refval)

Arguments:
           img - Image to be rendered.  ::ImageMeta or ::Array{Float64,2}
           map - Colour map to render the data with.
        refval - Reference value to be associated with centre point of
                 diverging colour map.  Defaults to 0.
Returns:
        rgbimg - The rendered image.
                 ::ImageMeta{Float64,3} or ::Array{Float64,3}
```
For a list of all diverging colour maps that can be generated by cmap()
use: > cmap("div")

See also: applycolourmap, applycycliccolourmap
"""
function applydivergingcolourmap(img::AbstractArray, cmap::AbstractArray, refval::Real = 0.0)

    minv = minimum(img)
    maxv = maximum(img)

    if refval < minv || refval > maxv
        @warn("Reference value is outside the range of image values")
    end

    dr = maximum([maxv - refval, refval - minv])
    rnge = [-dr dr] .+ refval

    return applycolourmap(img, cmap, rnge)
end

# For those who spell colour without a 'u'
"""
applydivergingcolormap - Applies a diverging colour map to an image

```
For full documentation see applydivergingcolourmap
                                             ^
```
See also: applycolourmap, applycycliccolourmap
"""
function applydivergingcolormap(img::AbstractArray, cmap::AbstractArray, refval::Real = 0.0)
    applydivergingcolourmap(img, cmap, refval)
end


text:
using Polyhedra

using Compat
using Compat.SparseArrays
using Compat.LinearAlgebra
using Compat.Test

import MultivariatePolynomials
MP = MultivariatePolynomials

using StaticArrays

include("utils.jl")

include("solvers.jl")

include("elements.jl")
include("comp.jl")
include("representation.jl")
include("polyhedron.jl")

include("redundancy.jl")
include("doubledescription.jl")

include("interval.jl")

include("polyhedra_to_lpqp.jl")
include("default.jl")

include("show.jl")

include("polyhedra.jl")
for (arith, T) in (("floating point", Float64), ("exact", Rational{BigInt}))
    @testset "Polyhedra tests in $arith arithmetic" begin
        polyhedratest(SimplePolyhedraLibrary{T}(lp_solver), ["board"])
    end
end


text:
# Run simulations with multivariate normal mixtures to compare the MFM with the DPM.
module SimCompareDPM

using BayesianMixtures
B = BayesianMixtures
Theta = B.MVN.Theta
can_plot = (Pkg.installed("PyPlot")!=nothing)
can_save = (Pkg.installed("JLD")!=nothing)
can_plot ? using PyPlot : warn("Skipping plots since PyPlot is not installed.")
can_save ? using JLD : warn("Saving disabled since JLD is not installed.")

# Settings
plot_runs = true # display plots of results or not
save_runs = true # save results to file or not
save_figures = true  # save the figures to file
from_file = false # true: use existing results from file, false: run sampler to generate results
reset_random = true # reset the random number generator before each run or not
use_hyperprior = false # put hyperprior on based distribution parameters or not
n_values = [50,100,250,1000,2500] #[50,100,250,1000,2500]  # sizes of data sets to simulate
model_types = ["MFM","DPM"]
data_ID = "k = 3" # label of data distribution to simulate from (k = 3 or k = 4)
t_max = 25 # guess at an upper bound on number of clusters that will occur
reps = 1:5 # 1:5  # replications to run for each model,n pair
n_burn = 5000 # 5000
n_total = 20*n_burn
n_keep = min(1000,n_total)
results_directory = "results-SimCompareDPM"

# Create directory for results if it doesn't exist
if !isdir(results_directory); mkdir(results_directory); end
println(results_directory)

# Simulate data and run the sampler for each n, rep, model_type
if !from_file
    for n in n_values
        for rep in reps

            # Simulate from the data distribution
            if reset_random; srand(n+rep); end
            if data_ID=="k = 3"
                k = 3
                d = 2
                probabilities = [0.45,0.3,0.25]
                angle = pi/4
                R = [cos(angle) -sin(angle);sin(angle) cos(angle)]
                V = diagm([2.5,0.2])
                theta = [Theta([4,4],eye(d)),
                         Theta([7,4],inv(R*V*R')),
                         Theta([6,2],diagm(1 ./ [3,0.1]))]
                x,z = B.MVN.mixrnd(n,probabilities,theta)
                true_distribution = (k,probabilities,theta,z)
                @assert(n>2*d)

            elseif data_ID=="k = 4"
                k = 4
                d = 2
                probabilities = [0.44,0.3,0.25,0.01]
                angle = pi/4
                R = [cos(angle) -sin(angle);sin(angle) cos(angle)]
                V = diagm([2.5,0.2])
                theta = [Theta([4,4],eye(d)),
                         Theta([7,4],inv(R*V*R')),
                         Theta([6,2],diagm(1 ./ [3,0.1])),
                         Theta([8,11],diagm(1 ./ [0.1,0.1]))]
                x,z = B.MVN.mixrnd(n,probabilities,theta)
                true_distribution = (k,probabilities,theta,z)
                @assert(n>2*d)
            else
                error("Unknown data_ID: $data_ID")
            end

            for model_type in model_types
                # Run sampler
                if reset_random; srand(n+rep); end
                options = B.options("MVN",model_type,x,n_total;
                    n_keep=n_keep, n_burn=n_burn, use_hyperprior=use_hyperprior, t_max=t_max)
                result = B.run_sampler(options)

                # Save result to file
                if can_save && save_runs
                    B.save_result(joinpath(results_directory,"n=$n-rep=$rep-type=$model_type.jld"),result)
                end
            end
            gc()  # (call the garbage collector to free up memory)
        end
    end
end

if can_save && can_plot && save_runs && plot_runs
    figure_number = 0
    color_list = Any[[0.5,1,0.5],[0.25,0.25,1],[0.75,0,1],"y","k"]
    marker_list = Any["s","^","v","o","*"]

    # Plot posterior on t for each n, model_type
    for model_type in model_types
        B.open_figure(figure_number+=1)
        for (i_n,n) in enumerate(n_values)
            results = [B.load_result(joinpath(results_directory,"n=$n-rep=$rep-type=$model_type.jld")) for rep in reps]
            B.plot_t_posterior_average(results; color=color_list[i_n], marker=marker_list[i_n], label="n = $n")
            B.labels("$model_type posterior on t","t (number of clusters)","p(t|data)")
            xlim(0,10); ylim(0,1)
            legend(loc="upper right", fontsize=10, numpoints=1)
        end
        if save_figures; savefig("simCompareDPM-t_posterior-$model_type.png",dpi=200); end
    end

    # Plot MFM posterior on k for each n
    for model_type in intersect(model_types,["MFM"])
        B.open_figure(figure_number+=1)
        for (i_n,n) in enumerate(n_values)
            results = [B.load_result(joinpath(results_directory,"n=$n-rep=$rep-type=$model_type.jld")) for rep in reps]
            B.plot_k_posterior_average(results; color=color_list[i_n], marker=marker_list[i_n], label="n = $n")
            B.labels("$model_type posterior on k","k (number of components)","p(k|data)")
            xlim(0,10); ylim(0,1)
            legend(loc="upper right", fontsize=10, numpoints=1)
        end
        if save_figures; savefig("simCompareDPM-k_posterior-$model_type.png",dpi=200); end
    end

    # Plot some density estimates for each model_type
    for model_type in model_types
        for (i_n,n) in enumerate(intersect(n_values,[50,250,1000]))
            B.open_figure(figure_number+=1;figure_size=(6,5))
            rep = 1
            result = B.load_result(joinpath(results_directory,"n=$n-rep=$rep-type=$model_type.jld"))
            B.plot_density_estimate(result; resolution=40)
            B.labels("$model_type density estimate, n = $n")
            if save_figures; savefig("simCompareDPM-density-$model_type-n=$n.png",dpi=200); end
        end
    end

    # Plot some sample clusterings for each model_type
    for model_type in model_types
        for (i_n,n) in enumerate(intersect(n_values,[250]))
            rep = 1
            indices = [800,900,1000]
            result = B.load_result(joinpath(results_directory,"n=$n-rep=$rep-type=$model_type.jld"))
            for index in intersect(1:n_keep,indices)
                B.open_figure(figure_number+=1;figure_size=(6,5))
                x = result.options.x
                z = result.z[:,index][:]
                B.plot_clusters(x,z; markersize=5)
                B.labels("$model_type sample clusters, n = $n, index = $index")
                if save_figures; savefig("simCompareDPM-clusters-$model_type-n=$n-index=$index.png",dpi=200); end
            end
        end
    end
end

end # module simCompareDPM


text:
import RowEchelon

#=
Homework 3
Section 1.5
Problem 11
=#

# Create the matrix to be evaluated
A = [1 -4 -2 0 -3 -5 ;
	 0 0 1 0 0 -1 ;
	 0 0 0 0 1 -4 ;
	 0 0 0 0 0 0 ]

# Row reduce matrix A
A_rref = RowEchelon.rref(A)
sizeA_rref = size(A_rref)

# Identify free variables
leadingEntryLocated = false
pivotCols = zeros(0)
pivotRows = zeros(0)
row = 1
col = 1
while row <= sizeA_rref[1]
	# Loop once for each row
	global leadingEntryLocated
	global pivotCols
	global pivotRows
	global row
	global col
	leadingEntryLocated = false
	col = 1
	while col <= sizeA_rref[2]
		# Loop once for each col
		if A_rref[row,col] != 0 && leadingEntryLocated == false
			# Leading Entry Located
			leadingEntryLocated = true
			append!(pivotCols,col)
			append!(pivotRows,row)
		end
		col += 1
	end
	row += 1
end

freeVars = zeros(0)

for i = 1:sizeA_rref[2]
	if !(i in pivotCols)
		# Every column that isnt a pivot column
		# Has a free variable
		append!(freeVars, i)
	end
end

# Print Solution
println()
print("x = ")
for freeVar in freeVars
	vector = zeros(6)
	for rows = 1:sizeA_rref[1]
		if A_rref[rows,Int64(freeVar)] != 0
			vector[rows] = -A_rref[rows,Int64(freeVar)]
		end
	end
	vector[Int64(freeVar)] = 1
	println()
	print("x_");print(Int64(freeVar))
	print("(")
	print(vector)
	print(")")
	freeVarLength = size(freeVars)[1]
	if freeVar != freeVars[freeVarLength]
		print(" + ")
	end
end

println()
println()
exit(0)


text:
#=
in this example, we just put all the data into the
network at once over and over again to optimizes the
parameters of network.
=#

# first class data and lable
using Plots
input1 = 3*randn(2,1024);
label1 = zeros(2,1024);
label1[1,:] .= 1.0;

# second class data and lable
t = 0:pi/1000:2*pi;
input2 = zeros(2,length(t));
label2 = zeros(2,length(t));
label2[2,:] .= 1.0;
for i = 1:length(t)
    r = randn() + 18.0;
    input2[1,i] = r * cos(t[i]);
    input2[2,i] = r * sin(t[i]);
end

# make model's input and label
input = hcat(input1,input2)
label = hcat(label1,label2)
x = Variable(input);
l = Variable(label);
N = size(x,2)

scatter(input2[1,:],input2[2,:],color=:blue)
scatter!(input1[1,:],input1[2,:],color=:red)

function runme(epoch, lr, opt)
    mlpmodel  = Chain(Dense(2,32),
                      Dense(32,32),
                      Dense(32,32),
                      Linear(32,2))

    params  = xparamsof(mlpmodel)
    optim   = opt(params;lr=lr)
    lossval = zeros(epoch,1)

    for i = 1:epoch
        outs = forward(mlpmodel, x)
        COST = crossEntropyCost(softmax(outs;dims=1), l)
        backward()
        update!(optim)
        zerograds!(optim)
        lossval[i] = COST/N
    end
    return lossval
end

plot(vec(log.(runme(400, 1e-4, SGD))), color=:red,label="SGD");
plot!(vec(log.(runme(400, 1e-4, Momentum))), color=:orange,label="Momentum");
plot!(vec(log.(runme(400, 1e-4, Adam))), color=:green,label="Adam");
plot!(vec(log.(runme(400, 1e-3, RMSProp))), color=:blue,label="RMSProp");
plot!(vec(log.(runme(400, 1e-3, AdaGrad))), color=:cyan,label="AdaGrad")

L1 = 1.0
function runmeL1(epoch, lr, opt, L1)
    mlpmodel  = Chain(Dense(2,32),
                      Dense(32,32),
                      Dense(32,32),
                      Linear(32,2))

    params  = xparamsof(mlpmodel)
    optim   = opt(params;lr=lr,L1decay=L1)
    lossval = zeros(epoch,1)

    for i = 1:epoch
        outs = forward(mlpmodel, x)
        COST = crossEntropyCost(softmax(outs;dims=1), l)
        backward()
        update!(optim)
        zerograds!(optim)
        lossval[i] = COST/N
    end
    return lossval
end
plot(vec(log.(runmeL1(400, 1e-4, SGDL1, L1))), color=:red,label="SGDL1")
plot!(vec(log.(runmeL1(400, 1e-4, MomentumL1, L1))), color=:orange,label="MomentumL1")
plot!(vec(log.(runmeL1(400, 1e-4, AdamL1, L1))), color=:green,label="AdamL1")
plot!(vec(log.(runmeL1(400, 1e-3, RMSPropL1, L1))), color=:blue,label="RMSPropL1")
plot!(vec(log.(runmeL1(400, 1e-3, AdaGradL1, L1))), color=:cyan,label="AdaGradL1")

L2 = 1.0
function runmeL2(epoch, lr, opt, L2)
    mlpmodel  = Chain(Dense(2,32),
                      Dense(32,32),
                      Dense(32,32),
                      Linear(32,2))

    params  = xparamsof(mlpmodel)
    optim   = opt(params;lr=lr,L2decay=L2)
    lossval = zeros(epoch,1)

    for i = 1:epoch
        outs = forward(mlpmodel, x)
        COST = crossEntropyCost(softmax(outs;dims=1), l)
        backward()
        update!(optim)
        zerograds!(optim)
        lossval[i] = COST/N
    end
    return lossval
end
plot(vec(log.(runmeL2(400, 1e-4, SGDL2, L2))), color=:red,label="SGDL2")
plot!(vec(log.(runmeL2(400, 1e-4, MomentumL2, L2))), color=:orange,label="MomentumL2")
plot!(vec(log.(runmeL2(400, 1e-4, AdamL2, L2))), color=:green,label="AdamL2")
plot!(vec(log.(runmeL2(400, 1e-3, RMSPropL2, L2))), color=:blue,label="RMSPropL2")
plot!(vec(log.(runmeL2(400, 1e-3, AdaGradL2, L2))), color=:cyan,label="AdaGradL2")

L1 = 1.0
L2 = 1.0
function runmeL1L2(epoch, lr, opt, L1, L2)
    mlpmodel  = Chain(Dense(2,32),
                      Dense(32,32),
                      Dense(32,32),
                      Linear(32,2))

    params  = xparamsof(mlpmodel)
    optim   = opt(params;lr=lr,L1decay=L1,L2decay=L2)
    lossval = zeros(epoch,1)

    for i = 1:epoch
        outs = forward(mlpmodel, x)
        COST = crossEntropyCost(softmax(outs;dims=1), l)
        backward()
        update!(optim)
        zerograds!(optim)
        lossval[i] = COST/N
    end
    return lossval
end
plot(vec(log.(runmeL1L2(400, 1e-4, SGDL1L2, L1, L2))),color=:red,label="SGDL1L2")
plot!(vec(log.(runmeL1L2(400, 1e-4, MomentumL1L2, L1, L2))),color=:orange,label="MomentumL1L2")
plot!(vec(log.(runmeL1L2(400, 1e-4, AdamL1L2, L1, L2))),color=:green,label="AdamL1L2")
plot!(vec(log.(runmeL1L2(400, 1e-5, RMSPropL1L2, L1, L2))),color=:blue,label="RMSPropL1L2")
plot!(vec(log.(runmeL1L2(400, 1e-5, AdaGradL1L2, L1, L2))),color=:cyan,label="AdaGradL1L2")

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta