text:
module DiscreteChoiceModels

using OnlineStats, Tables, Printf, Statistics, LinearAlgebra, DataFrames, Dagger, ForwardDiff, FLoops, SplittablesBase
import FunctionWrappers: FunctionWrapper

include("LogitModel.jl")
include("util.jl")
include("utilitymacro.jl")
include("mnl.jl")
include("compute_utility.jl")

export multinomial_logit, @utility, @β, @dummy_code

end


text:
using FITSIO, DSP

export Spectrum, goodfraction

unit_λ() = u"angstrom"
unit_flux() = u"erg" / u"s" / u"cm"^2
unit_lum() =  u"erg" / u"s"

unit_flux_density() = unit_flux() / unit_λ()
unit_lum_density() = unit_lum() / unit_λ()

log10_scale_λ()    = 0
log10_scale_flux() = -17
log10_scale_lum()  = 42

scale_λ()    = 10. ^log10_scale_λ()
scale_flux() = 10. ^log10_scale_flux()
scale_lum()  = 10. ^log10_scale_lum()

struct Spectrum
    label::String
    λ::Vector{Float64}
    flux::Vector{Float64}
    err::Vector{Float64}
    good::Vector{Bool}
    resolution::Float64  # σ km/s
    meta::Dict{Symbol, Any}
    function Spectrum(λ::Vector{T}, flux::Vector{T}, err::Vector{T};
                      good::Union{Nothing, Vector{Bool}}=nothing,
                      resolution=NaN,
                      label="") where T <: AbstractFloat
        if good == nothing
            good = fill(true, length(λ))
        end
        if length(err) == 0
            @warn "Uncertainties were not provided: assuming 0.05 * (median(flux) + abs(flux))"
            err = 0.05 .* (median(flux) .+ abs.(flux))
        end
        @assert length(λ) == length(flux) == length(err) == length(good)
        @assert minimum(err) > 0 "Uncertainties must be positive!"

        # Ensure wavelengths are sorted
        ii = sortperm(λ)
        λ = λ[ii]

        # Estimate sampling resolution in km/s
        sampling_res = median((λ[2:end] .- λ[1:end-1]) ./ ((λ[2:end] .+ λ[1:end-1]) ./ 2)) * 3e5
        if isnan(resolution)
            @warn "Resolution is not provided, assuming it is equal to twice the sampling resolution..."
            resolution = 2 * sampling_res
        end
        @assert sampling_res < resolution "Estimated resolution ($sampling_res) < provided resolution ($resolution)"

        return new(label, λ, flux[ii], err[ii], good, resolution,
                   Dict{Symbol, Any}(:sampling_resolution => sampling_res))
    end
end

Spectrum(λ::Vector{T1}, flux::Vector{T2}; kw...) where {T1 <: Quantity, T2 <: Quantity} = Spectrum(λ, flux, flux[[]]; kw...)
function Spectrum(λ::Vector{T1}, flux::Vector{T2}, err::Vector{T2}; kw...) where {T1 <: Quantity, T2 <: Quantity}
    Spectrum(getproperty.(uconvert.(Ref(unit_λ())           , λ   ), :val) ./ scale_λ(),
             getproperty.(uconvert.(Ref(unit_flux_density()), flux), :val) ./ scale_flux(),
             getproperty.(uconvert.(Ref(unit_flux_density()), err ), :val) ./ scale_flux();
             kw...)
end

function Spectrum(::Val{:SDSS_DR10}, file::AbstractString; ndrop=100, resolution=150.)  # TODO: Check resolution is correct
    f = FITS(file)
    λ = 10 .^read(f[2], "loglam")
    flux = float.(read(f[2], "flux"))
    ivar = float.(read(f[2], "ivar"))
    mask = read(f[2], "and_mask")
    close(f)

    ii = sortperm(λ)
    λ    =    λ[ii]
    flux = flux[ii]
    ivar = ivar[ii]
    mask = mask[ii]

    good = convert(Vector{Bool}, ((mask .== 0)  .&
                                  (ivar .> 0)   .&
                                  (flux .> 0)))
    if ndrop > 0
        good[1:ndrop] .= false
        good[end-ndrop+1:end] .= false
    end

    out = Spectrum(λ, flux, sqrt.(1 ./ ivar), good=good, label=file, resolution=resolution)
    return out
end

function Spectrum(::Val{:ASCII}, file::AbstractString; columns=[1,2,3], kw...)
    @assert length(columns) >= 2

    λ    = Vector{Float64}()
    flux = Vector{Float64}()
    unc  = Vector{Float64}()
    for l in readlines(file)
        l = strip(strip(l))
        (l[1] == '#')  &&  continue
        s = string.(split(l, keepempty=false))
        @assert length(s) >= maximum(columns)
        push!(λ   , Meta.parse(s[columns[1]]))
        push!(flux, Meta.parse(s[columns[2]]))
        if length(columns) == 3
            push!(unc, Meta.parse(s[columns[3]]))
        end
    end
    good = fill(true, length(λ))
    out = Spectrum(λ, flux, unc, label=file; kw...)
    return out
end

goodfraction(d::Spectrum) = length(findall(d.good)) / length(d.good)

function instrumental_broadening(λ, flux, σ_kms)
    # In some case the main reducer may evaluate to a single value
    # (rather than a vector)
    (length(flux) == 1)  &&  (return flux)

    return conv_gauss(λ, flux, σ_kms)
end

ext: jl


text:
using TaylorSeries
# computes all available derivatives, including value (0th order derivative)
derivatives(t::Taylor1) = derivatives!(similar(t.coeffs), t)
function derivatives!(d::AbstractVector, t::Taylor1)
    c, n = t.coeffs, t.order
    @. d = c * factorial(0:n)
end
# computes all derivatives up to mth order (including value)
function derivatives(f, x::Real, m::Int)
    d = zeros(eltype(x), m+1)
    derivatives!(d, f, x)
end
function derivatives(f, x::AbstractVector{<:Real}, m::Int)
    d = zeros(eltype(x), m+1, length(x))
    for (i, xi) in enumerate(x)
        d[:, i] .= derivatives(f, xi, m)
    end
    return d
end

function derivatives!(d::AbstractVector, f, x::Real)
    m = length(d) - 1
    @. d = 0
    d[1] = x
    d[2] = 1
    x = Taylor1(d, m) # mth order taylor series
    fx = f(x) # TODO: in place?
    derivatives!(d, fx) # overwrites coeffs vector with derivatives
end

# IDEA: could specialize for popular kernel types
# function derivatives!(d::AbstractVector, f::CovarianceFunctions.Exp, r::Real)
#     m = length(d)-1
#     @. d = f(r)
# end


text:
# =================================
# Parsing package data
# =================================

function _get_dimensions(Pdata)
    Pdim = Pdata["dimensions"]
    depth = Pdim["depth_cm"]
    height = Pdim["height_cm"]
    width = Pdim["width_cm"]
    return Dimensions(depth, height, width)
end

function _get_package_id(Pid)
    return Pid[11:end]
end

function _get_time_utc(time_utc)
    if !isa(time_utc, String) && isnan(time_utc)
        # packages do not have associated time_utc
        date = Date("")
        time = Time("")
    else
        date = Date(split(time_utc, " ")[1], DateFormat("y-m-d"))
        time = Time(split(time_utc, " ")[2])
    end
    return date, time
end

function _get_time_window(Pdata)
    start_time_utc = Pdata["time_window"]["start_time_utc"]
    start_date, start_time = _get_time_utc(start_time_utc)

    end_time_utc = Pdata["time_window"]["end_time_utc"]
    end_date, end_time = _get_time_utc(end_time_utc)

    time_window = TimeWindow(start_date, start_time, end_date, end_time)

    return time_window
end

function _get_scan_status(Pdata)
    status = 0
    if haskey(Pdata, "scan_status")
        val = Pdata["scan_status"]
        if val == "REJECTED"
            status = 1
        elseif val == "DELIVERY_ATTEMPTED"
            status = 2
        elseif val == "DELIVERED"
            status = 3
        end
    end
    return status
end

function _get_service_time(Pdata)
    return Pdata["planned_service_time_seconds"]
end

function parse_package_data(; file=joinpath("..", "data", "model_apply_inputs", "new_package_data.json"))
    package_data = read_json(file)
    return parse_package_data(package_data)
end

# receives a dictionary returned by read_json
function parse_package_data(data::Dict{String, Any})
    m = length(data)

    PT = Package
    ST = Task{BareStop, PT}
    RT = BareRoute{ST}

    routes = Vector{RT}(undef, m)

    i = 1
    for (R, D) in data  # loop over routes
        route_hash = _get_name(R)
        num_stops = length(D)
        stops = Vector{ST}(undef, num_stops)
        k = 1

        for (Sid, Sdata) in D  # loop over stops
            num_packages = length(Sdata) # number of packages in this stop
            packages = Vector{PT}(undef, num_packages)
            j = 1

            for (Pid, Pdata) in Sdata # loop over packages

                package_hash = _get_package_id(Pid)
                time_window = _get_time_window(Pdata)
                planned_service = _get_service_time(Pdata)
                package_dimensions = _get_dimensions(Pdata)
                status = _get_scan_status(Pdata)

                packages[j] = Package(package_hash, time_window, planned_service, package_dimensions, status)
                j += 1
            end
            stops[k] = Task(BareStop(Sid), packages)
            k += 1
        end
        routes[i] = BareRoute(route_hash, stops)
        i += 1
    end
    return routes
end

# =================================
# Parsing route data
# ================================

function _get_name(R)
    return R[9:end] # ignore RouteID
end

function _get_score(R)
    score = 0
    if haskey(R, "route_score")
        val = R["route_score"]
        if val == "Low"
            score = 1
        elseif val == "Medium"
            score = 2
        elseif val == "High"
            score = 3
        end
    end
    return score
end

function _get_zone(S)
    zone_str = S["zone_id"]
    if isa(zone_str, Float64) && isnan(zone_str)
        # stations do not have associated zone_id
        zone_id = Zone('0', "")
    else
        zone_id = Zone(zone_str)
    end
    return zone_id
end

function parse_route_data(; file=joinpath("..", "data", "model_apply_inputs", "new_route_data.json"))
    route_data = read_json(file)
    return parse_route_data(route_data)
end

# receives a dictionary returned by read_json
function parse_route_data(data::Dict{String, Any})
    m = length(data)
    routes = Vector{Route{Stop}}(undef, m)

    i = 1
    for (R, D) in data
        hex_hash = _get_name(R)
        station_code = D["station_code"]
        date = Date(D["date_YYYY_MM_DD"], DateFormat("y-m-d"))
        departure_time = Time(D["departure_time_utc"])
        capacity = D["executor_capacity_cm3"]

        Dstops = D["stops"]
        stops = Vector{Stop}(undef, length(Dstops))
        k = 1
        for (Sid, Sdata) in Dstops
            id = Sid
            lat = Sdata["lat"]
            lng = Sdata["lng"]
            stop_type = Sdata["type"] == "Station" ? 0 : 1
            zone_id = _get_zone(Sdata)
            stops[k] = Stop(id, lat, lng, stop_type, zone_id)
            k += 1
        end
        score = _get_score(D)
        routes[i] = Route(hex_hash, station_code, date, departure_time, capacity, stops, score)
        i += 1
    end
    return routes
end

# =================================
# Parsing travel times data
# =================================

function parse_travel_times(; file=joinpath("..", "data", "model_apply_inputs", "new_travel_times.json"))
    travel_times_data = read_json(file)
    return parse_travel_times(travel_times_data)
end

# receives a dictionary returned by read_json
function parse_travel_times(data::Dict{String, Any})
    m = length(data)
    travel_times = Vector{TravelTimes}(undef, m)

    i = 1
    for (R, O) in data # loop over routes

        hex_hash = _get_name(R)
        stop_names = collect(keys(O))
        num_stops = length(stop_names)

        M = zeros(Float64, num_stops, num_stops)
        times_matrix = NamedArray(M, (stop_names, stop_names))

        for orig in keys(O) # loop over origins
            for dest in keys(O[orig]) # loop over destinations
                times_matrix[orig, dest] = O[orig][dest]
            end
        end
        travel_times[i] = TravelTimes(hex_hash, stop_names, times_matrix)
        i += 1
    end
    return travel_times
end

# =================================
# Reading and writing JSON files
# =================================

"""
    read_json(src)

Read the JSON file in the string `src` and store it in a Julia dictionary.

### Notes

See https://discourse.julialang.org/t/dataframes-best-way-to-import-from-json-format-file/29946/6
"""
function read_json(src)
    # create variable to write the information
    d = Dict()
    open(src, "r") do f
        d = JSON.parse(f)  # parse and transform data
    end
    return d
end

function write_csv(R::Route)
    S = stops(R)

    # cada fila es (name, zone_id, latitud, longitud)
    xval = longitude.(S)
    yval = latitude.(S)
    names = name.(S)
    zones = zone.(S)
    A = hcat(names, zones, yval, xval)

    route_id = name(R)
    CSV.write(route_id * ".csv", DataFrame(A), writeheader=false)
end

# ==============================
# Methods to parse sequences
# ==============================

function parse_sequence(target="actual"; file=joinpath("..", "data", "model_apply_inputs", "new_actual_sequences.json"))
    sequence_data = read_json(file)
    return parse_sequence(sequence_data)
end

# receive a dictionary and return a vector of sequences
# use target="proposed" to parse our own sequences
function parse_sequence(data::Dict{String, Any}, target="actual")
    m = length(data)
    sequences = Vector{Sequence{BareStop}}(undef, m)

    i = 1
    for (R, D) in data
        route_hash = _get_name(R)
        actual_sequence = Vector{BareStop}(undef, length(D[target]))

        for (SeqId, SeqData) in D[target]
            # el vector de actual sequence ya queda ordenado
            actual_sequence[SeqData+1] = BareStop(SeqId)
        end

        sequences[i] = Sequence(route_hash, actual_sequence)
        i += 1
    end

    return sequences
end

function write_json_sequence(data; path=joinpath("..", "data", "model_apply_outputs", "proposed_sequences.json"))
    # dictionary to write
    dict = sequence_to_dict(data)

    # pass data as a json string (how it shall be displayed in a file)
    stringdata = JSON.json(dict)

    # write the file with the stringdata variable information
    open(path, "w") do f
        write(f, stringdata)
    end;
    return nothing
end

function sequence_to_dict(S::Sequence)
    dict = Dict{String, Any}()

    route_id = "RouteID_" * name(S)
    seq = Dict{String, Int}()
    for (i, stop) in enumerate(stops(S)) # loop over stops, assumed that the stops are sorted
        stop_id = name(stop)
        seq[stop_id] = i-1
    end
    dict[route_id] = Dict("proposed"=>seq)
    return dict
end

function sequence_to_dict(data::Vector{<:Sequence})
    dict = Dict{String, Any}()

    for S in data # loop over sequences
        route_id = "RouteID_" * name(S)
        seq = Dict{String, Int}()
        for (i, stop) in enumerate(stops(S)) # loop over stops, assumed that the stops are sorted
            stop_id = name(stop)
            seq[stop_id] = i-1
        end
        dict[route_id] = Dict("proposed"=>seq)
    end
    dict
end

function print_score(S::Sequence; src="evaluar.py")
    # convert the sequence to JSON format
    write_json_sequence(S)

    route_id = "RouteID_" * name(S)
    score = run(`python $src $route_id`)
    return score
end

get_score(S::Vector{<:Sequence}) = get_score.(S)

function get_score(S::Sequence)
    out = IOCapture.capture() do
                 print_score(S)
          end
    out_str = split(out.output[1:end-1], "\n")
    result = [parse(Float64, x) for x in out_str]
    return result[1]
end


text:
using JuMP, Gurobi, Ipopt
using DataFrames, JLD, CSVFiles
using MatpowerCases

### Include the file with input data functions
include("NetworkDataType.jl")
include("NetworkLoad.jl")
include("P2PTradeType.jl")
include("TradeFunction.jl")
include("DLMP.jl")

testsystem = "AP15busDN"
#testsystem = "case141_OPF"

buses, lines, generators, datamat = NetworkLoad(testsystem)
lineset = 1:length(lines)
busset = 1:length(buses)
genset = 1:length(generators)
buses[1].Vmax = 1
buses[1].Vmin = 1
windset = Int[]
for g in genset
    if datamat["gen"][g,2] == "Wind"
        push!(windset, generators[g].gindex)
    end
end

root = -1
for g in genset
    if generators[g].gtype=="Root"
        root = g
    end
end
if root == -1
    display("No generator at root node")
    generators[1].gtype = "Root"
    root = 1
end

temp = Generator
for g in genset
    if generators[g].gtype == "Root"
        temp = Generator(length(genset)+1, "DG", generators[g].location, generators[g].Pg, generators[g].Qg, generators[g].Qmax,
        generators[g].Qmin, generators[g].Vg, generators[g].mBase, generators[g].status, generators[g].Pmax, generators[g].Pmin,
        generators[g].cost, generators[g].SUcost, generators[g].SDcost, generators[g].RU, generators[g].RD, generators[g].UPtime, generators[g].DNtime)
        #push!(generators, temp)
    end
end

genset = 1:length(generators)
gensetU = root
gensetP = setdiff(genset, gensetU)

B_g = []
for g in genset
   if ~any(x->x==generators[g].location, B_g)
      push!(B_g, generators[g].location)
      buses[generators[g].location].Pd = 0
   end
end
B_d = setdiff(busset,B_g)

B_gn = Array{Array{Int64}}(length(buses),1)
for ii in 1:length(buses)
    B_gn[ii] = Int64[]
end
for g in genset
    push!(B_gn[generators[g].location], g)
end
############
partlevel = 0.0
gam = partlevel*ones(busset)
############
###
trade_scale = 1e-3 # [MWh]

agents = Agent[]
for b in 1:length(buses)
    if isempty(B_gn[b])
        temp_agent = Agent(b, b, 1, 0, 0, [], gam[b]*buses[b].Pd, gam[b]*buses[b].Qd)
        push!(agents,temp_agent)
    else
        for g in B_gn[b]
            if generators[g].gtype != "Root"
                temp_agent = Agent(b, b, 1, generators[g].Pmin, generators[g].Pmax, generators[g].cost, gam[b]*buses[b].Pd, gam[b]*buses[b].Qd)
                push!(agents,temp_agent)
            end
        end
    end
end
# agents = Agent[]
# for b in 1:length(buses)
#     if isempty(B_gn[b])
#         temp_agent = Agent(b, b, 1, 0, 0, [], gam[b]*buses[b].Pd, gam[b]*buses[b].Qd)
#         push!(agents,temp_agent)
#     elseif generators[B_gn[b][1]].gtype != "Root"
#         g = B_gn[b][1]
#         temp_agent = Agent(b, b, 1, generators[g].Pmin, generators[g].Pmax, generators[g].cost, gam[b]*buses[b].Pd, gam[b]*buses[b].Qd)
#         push!(agents,temp_agent)
#     end
# end

# Trade(tradeID, kind, Pes, Peb, As, Ab, costNw)
trades = generateTrade(agents, trade_scale)
trades_old = generateTrade(agents, trade_scale)

trades_mat = Array{Any}(length(agents),length(agents))
trades_mat_buy = Array{Any}(length(agents),length(agents))
trades_mat_sell = Array{Any}(length(agents),length(agents))
for ii in 1:length(agents), jj in 1:length(agents)
    trades_mat[ii,jj] = []
    trades_mat_buy[ii,jj] = []
    trades_mat_sell[ii,jj] = []
end
trades_selected_mat = Array{Any}(length(buses),length(buses))
trades_selected_mat_buy = Array{Any}(length(buses),length(buses))
trades_selected_mat_sell = Array{Any}(length(buses),length(buses))
for ii in 1:length(buses), jj in 1:length(buses)
    trades_selected_mat[ii,jj] = []
    trades_selected_mat_buy[ii,jj] = []
    trades_selected_mat_sell[ii,jj] = []
end

for w in 1:length(trades)
    push!(trades_mat[trades[w].As,trades[w].Ab],trades[w].tradeID)
end

trades_selected = Int[]
trades_rest = Int[]
trades_selected_old = Int[]
trades_rest_old = Int[]
dispatch = zeros(length(buses))
GenCon = zeros(length(buses))
Generation = zeros(length(buses))
Consumption = zeros(length(buses))
GenCon_old = zeros(length(buses))
dispatch_old = zeros(length(buses))
status = Symbol()
DLMP = Float64[]
DLMP_stack = Any[]
dispatch_stack = Any[]
DLMPInfo = DataFrame()
NodeInfo = DataFrame()
LineInfo = DataFrame()
dispatch_peerG = 0.0*zeros(genset)
param_delta = 1.0
tic()
while true
    #for _ in 1:10
    accepted = Array{Int64}[]
    rejected = Array{Int64}[]
    pg = Float64[]
    Revenue_Sell = Float64[]
    Cost_Buy = Float64[]
    Cost_Network = Float64[]
    Cost_DG = Float64[]
    trades_old[:] = trades[:]
    trades_selected_old = trades_selected
    for w in 1:length(trades)
        trades[w].Pes = 0
        trades[w].Peb = 0
    end

    while true
        accepted = Array{Int64}[]
        rejected = Array{Int64}[]
        pg = Float64[]

        gap = 1e-3
        for agentID in 1:length(agents)
            temp_accepted, temp_rejected, temp_pg, Revenue_Sell, Cost_Buy, Cost_Network, Cost_DG = selectTrade(trades, agents, agentID, gap, trade_scale)
            push!(accepted, temp_accepted)
            push!(rejected, temp_rejected)
            push!(pg, temp_pg)
        end

        deltaP = 5*trade_scale #$5/MWh = $5e-3/kWh
        for w in 1:length(trades)
            if (w in accepted[trades[w].Ab]) && ~(w in accepted[trades[w].As])
                if trades[w].Peb > trades[w].Pes
                    trades[w].Pes += deltaP
                else
                    trades[w].Peb += deltaP
                end
            end
        end
        # diff = sum((w in accepted[trades[w].Ab]) && ~(w in accepted[trades[w].As]) for w in 1:length(trades))
        # if diff == 0
        #     break
        # end

        ###diff = sum(abs(trades[w].Pes - trades_old[w].Pes) for w in 1:length(trades)) + sum(abs(trades[w].Peb - trades_old[w].Peb) for w in 1:length(trades))
        if isempty(trades)
            diff = 0
        else
            diff = sum((w in accepted[trades[w].Ab]) && ~(w in accepted[trades[w].As]) for w in 1:length(trades))
        end
        if diff == 0
            break
        end
    end

    trades_selected = Int[]
    trades_rest = Int[]
    for w in 1:length(trades)
        if (w in accepted[trades[w].Ab]) && (w in accepted[trades[w].As])
            push!(trades_selected, w)
        else
            push!(trades_rest, w)
        end
    end

    for i in 1:length(agents)
        dispatch[i] = ceil(agents[i].Pd/trade_scale) # MW to kW conversion
    end

    for w in trades_selected
        dispatch[agents[trades[w].Ab].location] -= 1
        dispatch[agents[trades[w].As].location] += 1
        GenCon[agents[trades[w].Ab].location] -= 1
        GenCon[agents[trades[w].As].location] += 1
        Consumption[agents[trades[w].Ab].location] -= 1
        Generation[agents[trades[w].As].location] += 1
    end

    GenCon = GenCon*trade_scale
    dispatch = dispatch*trade_scale # Kw to MW conversion
    Generation = Generation*trade_scale
    Consumption = Consumption*trade_scale

    #dispatch_peerG = 0.0*zeros(genset)
    for g in gensetP
        dispatch_peerG[g] = Generation[generators[g].location]
    end

    ###### Calculate Network Charge #####
    SMP = generators[root].cost[end-1]
    status, DLMP, pgextra, NodeInfo, LineInfo, DLMPInfo = calculateDLMP(dispatch_peerG, buses, generators, lines, SMP, gensetP, gensetU);

    ###### Update Network Charge #####
    deltaNw = 1*trade_scale # $1/MWh = $1e-3/kWh
    if status == Symbol(:Optimal)
        for w in 1:length(trades)
            if w in trades_selected
                trades[w].costNw = (DLMP[agents[trades[w].Ab].location] - DLMP[agents[trades[w].As].location])/2
            end
        end
    else
        for w in trades_selected
            trades[w].costNw += param_delta*deltaNw/trade_scale
        end
    end
    push!(dispatch_stack, dispatch)
    push!(DLMP_stack, DLMP)
    if GenCon_old == GenCon && status == Symbol(:Optimal)
        break
    else
        dispatch_old[:] = dispatch[:]
        GenCon_old[:] = GenCon[:]
    end
end
simulation_time = toc()

trades_dis = zeros(length(buses),length(buses));
NWcharge_dis = zeros(length(buses),length(buses));
Echarge_dis = zeros(length(buses),length(buses));
temp = Array{Array{Int64,1}}(length(buses),1)
for i in 1:length(buses)
    temp[i] = []
end

for w in trades_selected
    trades_dis[agents[trades[w].Ab].location,agents[trades[w].As].location] -= 1*trade_scale
    trades_dis[agents[trades[w].As].location,agents[trades[w].Ab].location] += 1*trade_scale
    Echarge_dis[agents[trades[w].Ab].location,agents[trades[w].As].location] -= trades[w].Peb
    Echarge_dis[agents[trades[w].As].location,agents[trades[w].Ab].location] += trades[w].Pes
    NWcharge_dis[agents[trades[w].Ab].location,agents[trades[w].As].location] -= round(trades[w].costNw,4)
    NWcharge_dis[agents[trades[w].As].location,agents[trades[w].Ab].location] += round(trades[w].costNw,4)
    push!(temp[agents[trades[w].Ab].location], w)
    push!(temp[agents[trades[w].As].location], w)
end

dispatch_peerG = 0.0*zeros(genset)
for g in gensetP
    dispatch_peerG[g] = Generation[generators[g].location]
end

###### Calculate Network Charge #####
SMP = generators[root].cost[end-1]
status, DLMP, pg, NodeInfo, LineInfo, DLMPInfo = calculateDLMP(dispatch_peerG, buses, generators, lines, SMP, gensetP, gensetU);

#CSV.write("PeerCentric_DLMP.csv", DLMPInfo);
#CSV.write("PeerCentric_NodeInfo.csv", NodeInfo);
#CSV.write("PeerCentric_LineInfo.csv", LineInfo);

price_mat = zeros(length(buses),length(buses))

#############################################
for w in trades_selected
    push!(trades_selected_mat[agents[trades[w].As].location,agents[trades[w].Ab].location],trades[w].tradeID)
end

price_mat = zeros(length(buses),length(buses));
pricemin = zeros(length(buses),length(buses));
pricemax = zeros(length(buses),length(buses));
for ii in 1:length(buses), jj in 1:length(buses)
    temp_mat = []
    for ww in trades_selected_mat[ii,jj]
        push!(temp_mat,trades[ww].Pes)
        pricemin[ii,jj] = minimum(temp_mat)
        pricemin[jj,ii] = -minimum(temp_mat)
        pricemax[ii,jj] = maximum(temp_mat)
        pricemax[jj,ii] = -maximum(temp_mat)
        if pricemin[ii,jj] == pricemax[ii,jj]
            price_mat[ii,jj] = pricemin[ii,jj]
            price_mat[jj,ii] = -pricemin[ii,jj]
        else
            error("pricemin != pricemax")
        end
    end
end

Nprice_mat = zeros(length(buses),length(buses));
for w in trades_selected
    Nprice_mat[agents[trades[w].Ab].location,agents[trades[w].As].location] = 0.5*round(trades[w].costNw,5)
    Nprice_mat[agents[trades[w].As].location,agents[trades[w].Ab].location] = 0.5*round(trades[w].costNw,5)
end

Nprice_mat_min = zeros(length(buses),length(buses));
Nprice_mat_max = zeros(length(buses),length(buses));
for w = 1:length(trades)
    Nprice_mat_min[agents[trades[w].Ab].location,agents[trades[w].As].location] = min(round(trades[w].costNw,5),Nprice_mat_min[agents[trades[w].Ab].location,agents[trades[w].As].location])
    Nprice_mat_min[agents[trades[w].As].location,agents[trades[w].Ab].location] = min(round(trades[w].costNw,5),Nprice_mat_min[agents[trades[w].As].location,agents[trades[w].Ab].location])
    Nprice_mat_max[agents[trades[w].Ab].location,agents[trades[w].As].location] = max(round(trades[w].costNw,5),Nprice_mat_max[agents[trades[w].Ab].location,agents[trades[w].As].location])
    Nprice_mat_max[agents[trades[w].As].location,agents[trades[w].Ab].location] = max(round(trades[w].costNw,5),Nprice_mat_max[agents[trades[w].As].location,agents[trades[w].Ab].location])
end

Eprice_mat = zeros(length(buses),length(buses));
for w in trades_selected
    Eprice_mat[agents[trades[w].Ab].location,agents[trades[w].As].location] = round(trades[w].Pes,5)
    Eprice_mat[agents[trades[w].As].location,agents[trades[w].Ab].location] = round(trades[w].Peb,5)
end

Eprice_mat_min = zeros(length(buses),length(buses));
Eprice_mat_max = zeros(length(buses),length(buses));
for w = 1:length(trades)
    Eprice_mat_min[agents[trades[w].Ab].location,agents[trades[w].As].location] = min(round(trades[w].Pes,5),Eprice_mat_min[agents[trades[w].Ab].location,agents[trades[w].As].location])
    Eprice_mat_min[agents[trades[w].As].location,agents[trades[w].Ab].location] = min(round(trades[w].Peb,5),Eprice_mat_min[agents[trades[w].As].location,agents[trades[w].Ab].location])
    Eprice_mat_max[agents[trades[w].Ab].location,agents[trades[w].As].location] = max(round(trades[w].Pes,5),Eprice_mat_max[agents[trades[w].Ab].location,agents[trades[w].As].location])
    Eprice_mat_max[agents[trades[w].As].location,agents[trades[w].Ab].location] = max(round(trades[w].Peb,5),Eprice_mat_max[agents[trades[w].As].location,agents[trades[w].Ab].location])
end

Peer_revenue = sum(Echarge_dis,2);
Consumer_cost = 0
Generator_revenue = 0
for bb in B_d
    Consumer_cost -= Peer_revenue[bb]
end
for bb in B_g
    Generator_revenue += Peer_revenue[bb]
end

line_marginrate = zeros(length(lines),1)
for ll in 1:length(lines)
    line_marginrate[ll] = round.((lines[ll].u  - sqrt((LineInfo[:fp][ll])^2 + (LineInfo[:fq][ll])^2))/ lines[ll].u,3)
end
line_loading = 1-line_marginrate
bus_voltage = sqrt.(NodeInfo[:v])

v_marginrate = round.(0.1-abs.(1-sqrt.(NodeInfo[:v])),3)/0.1

Loading = DataFrame(line_loading=vec(line_loading))
Vmargin = DataFrame(v_marginrate=vec(v_marginrate))
#####
loading = zeros(length(lines))
for l in 1:length(lines)
   loading[l] = LineInfo[:flow][l] / lines[l].u
end

flow = LineInfo[:flow]
voltage = sqrt.(NodeInfo[:v])

#######
pnm = trades_dis
demand_P2P = zeros(length(buses))
for m in B_d
    demand_P2P[m] = sum(pnm[n,m] for n in B_g)
end

EP_P = zeros(length(busset))
for w in trades_selected
   EP_P[agents[trades[w].Ab].location] += trades[w].Pes
end
EP_U = zeros(length(busset))
for m in B_d
   EP_U[m] = (buses[m].Pd-demand_P2P[m])*DLMP[m]
end

NUC = zeros(length(busset))
for w in trades_selected
   NUC[agents[trades[w].Ab].location] += 0.5*trades[w].costNw
   NUC[agents[trades[w].As].location] += 0.5*trades[w].costNw
end
UR = EP_U + NUC
CP = EP_U + EP_P + NUC

GC = zeros(length(genset))
for g in genset
   GC[g] = generators[g].cost[end-1] * round(pg[g],5)
end
sum_GC_P = sum(GC[g] for g in gensetP)
sum_GC_U = sum(GC[g] for g in gensetU)
sum_DGR = sum(EP_P)
sum_DGP = sum_DGR - sum_GC_P

UP = sum(UR) - sum_GC_U
CashFlow_mat = DataFrame(EP_P=EP_P,EP_U=EP_U,NUC=NUC,UR=UR,CP=CP)
CashFlow_sum = DataFrame(sum_CP=sum(CP), sum_EP_P=sum(EP_P), sum_EP_U=sum(EP_U), sum_NUC=sum(NUC), sum_DGR=sum_DGR, sum_DGP=sum_DGP, UP=UP)

println("partlevel = ""$partlevel")
println("status = ""$status")
println("max loading = ""$(maximum(loading)*100)""%")
println("max voltage = ""$(round(maximum(sqrt.(NodeInfo[:v])),3))")
println("min voltage = ""$(round(minimum(sqrt.(NodeInfo[:v])),3))")
println("max DLMP = ""$(round(maximum(DLMP),3))")
println("min DLMP = ""$(round(minimum(DLMP),3))")
println("simulation time = ""$simulation_time")
method = "PC"
voltagefile = "$method$(length(buses))""_""$(Int(100*partlevel))""_voltage"".csv"
loadingfile = "$method$(length(buses))""_""$(Int(100*partlevel))""_loading"".csv"
DLMPfile = "$method$(length(buses))""_""$(Int(100*partlevel))""_DLMP"".csv"
cashflowmatfile = "$method$(length(buses))""_""$(Int(100*partlevel))""_cfmat"".csv"
cashflowsumfile = "$method$(length(buses))""_""$(Int(100*partlevel))""_cfsum"".csv"
# CSVFiles.save(voltagefile, DataFrame(voltage=voltage))
# CSVFiles.save(loadingfile, DataFrame(loading=100*loading))
# CSVFiles.save(DLMPfile, DataFrame(DLMP=DLMP))
# CSVFiles.save(cashflowmatfile, CashFlow_mat)
# CSVFiles.save(cashflowsumfile, CashFlow_sum)

#######
resultfile = "$method$(length(buses))""_""$(Int(round(100*partlevel)))""_""$(Int(1000*trade_scale))""kW"".jld"
#@save "$resultfile"

#######################################
trades_selected_n = Array{Array{Int64,1}}(length(buses),1)
trades_n = Array{Array{Int64,1}}(length(buses),1)
for i in 1:length(buses)
    trades_selected_n[i] = []
    trades_n[i] = []
end
for w in trades_selected
    push!(trades_selected_n[trades[w].Ab],w)
    push!(trades_selected_n[trades[w].As],w)
end
for w in 1:length(trades)
    push!(trades_n[trades[w].Ab],w)
    push!(trades_n[trades[w].As],w)
end
Erho = zeros(length(buses),1)
Ecn = zeros(length(buses),1)
Tvol = zeros(length(buses),1)
NUC_total = zeros(length(buses),1)
for i in 1:length(buses)
    if ~isempty(trades_selected_n[i])
        Tvol[i] = length(trades_selected_n[i])
        Ecn[i] = trade_scale*sum(0.5*(DLMP[trades[w].Ab]-DLMP[trades[w].As]) for w in trades_selected_n[i]) / Tvol[i]
        NUC_total[i] = trade_scale*sum(0.5*(DLMP[trades[w].Ab]-DLMP[trades[w].As]) for w in trades_selected_n[i])
        Erho[i] = sum(trades[w].Pes for w in trades_selected_n[i]) / Tvol[i]
    end
end
println("NUC total = ",sum(NUC_total))
println("Tvol = ", sum(Tvol))
#######################################


text:
module OCIM0
#=
This module serves to load the OCIM0 matrix and grid.
These are loaded from the public, persistant URL in FigShare.
The Julia JLD2 format version of the OCIM0 was created with the code
in the GitHub repository https://github.com/briochemc/OceanCirculations.
=#

using SparseArrays          # For sparse matrix
using DataDeps              # For storage location of data
using Downloads
using JLD2                  # For saving circulation as JLD2 format
using CodecZlib             # for JLD2 compression JLD2
using Unitful               # for units
using Reexport
@reexport using OceanGrids            # To store the grid

function fallback_download(remotepath, localdir)
    @assert(isdir(localdir))
    filename = basename(remotepath)  # only works for URLs with filename as last part of name
    localpath = joinpath(localdir, filename)
    Downloads.download(remotepath, localpath)
    return localpath
end

# OCIM0 URLs
function url(; version="")
    url_start = "https://files.figshare.com/"
    return "$(url_start)28336086/OCIM0.1.jld2"
end

# OCIM0 Hashes
sha(; version="") = "4aaaccead0422534ff1adbd208962a8dce7f6cd1af710dcb60faff236863b266"

# Create registry entry for OCIM0 in JLD2 format
function register_OCIM0(; version="")
    register(
        DataDep(
            "AIBECS-OCIM0.1",
            """
            References:
            - $(citation())
            """,
            url(),
            sha(),
            fetch_method = fallback_download
        )
    )
    return nothing
end

"""
    load

Returns the grid and the transport matrix (in that order).

### Usage

```
julia> grd, T = OCIM0.load()
```

See *DeVries and Primeau* (2011) and *Primeau et al.* (2013) for more details.
"""
function load(; version="")
    register_OCIM0()
    jld2_file = @datadep_str string("AIBECS-OCIM0.1/", "OCIM0.1.jld2")
    @info """You are about to use the OCIM0.1 model.
          If you use it for research, please cite:

          $(citation())

          You can find the corresponding BibTeX entries in the CITATION.bib file
          at the root of the AIBECS.jl package repository.
          (Look for the "DeVries_Primeau_2011" and "Primeau_etal_2013" keys.)
          """
    jldopen(jld2_file) do file
        file["grid"], ustrip.(file["T"])
    end
end

citation() = """
             - Primeau, F. W., Holzer, M., and DeVries, T. (2013), Southern Ocean nutrient trapping and the efficiency of the biological pump, J. Geophys. Res. Oceans, 118, 2547–2564, doi:10.1002/jgrc.20181.
             - DeVries, T. and F. Primeau, 2011: Dynamically and Observationally Constrained Estimates of Water-Mass Distributions and Ages in the Global Ocean. J. Phys. Oceanogr., 41, 2381–2401, doi:10.1175/JPO-D-10-05011.1
             """

versions() = [""]

end # end module

export OCIM0


text:
using MultivariateStats
#using PDMats
using Distances

type emulator_gp_type{T<:Real}
  sigmasq_cor::T
  rho::T
  kernel::Function
  input_train::Array{T,2}
  input_means::Array{T,1}
  input_covar::PDMat
  output_train::Array{T,2}
  output_means::Array{T,1}
  Sigma_train::Array{PDMat,1}
end

type emulator_gp_savable_type{T<:Real}
  sigmasq_cor::T
  rho::T
  kernel::String
  input_train::Array{T,2}
  input_means::Array{T,1}
  input_covar::PDMat
  output_train::Array{T,2}
  output_means::Array{T,1}
  Sigma_train::Array{Array{T,2},1} # JLD doesn't like saving
end
function emulator_gp_savable_type{T<:Real}(emu::emulator_gp_type{T})
  emulator_gp_savable_type(emu.sigmasq_cor,emu.rho,string(emu.kernel),emu.input_train,emu.input_means,emu.input_covar,emu.output_train,emu.output_means,Matrix{T}[emu.Sigma_train[i].mat for i in 1:length(emu.Sigma_train)])
end
function emulator_gp_type{T<:Real}(emu::emulator_gp_savable_type{T})
  emulator_gp_type(emu.sigmasq_cor,emu.rho,eval(Symbol(emu.kernel)),emu.input_train,emu.input_means,emu.input_covar,emu.output_train,emu.output_means,PDMat[PDMat(emu.Sigma_train[i]) for i in 1:length(emu.Sigma_train)])
end

function kernel_exp(d::Float64; rho::Float64 = 1.0)
   exp(-d/rho)
end

function kernel_exp(d::Array{Float64}; rho::Float64 = 1.0)
   exp(-d/rho)
end

function make_kernel_data{T<:Real}(x::AbstractArray{T,2}, kernel::Function;
        #param_means::AbstractArray{T,1} = vec(mean(x,2)), param_covar::AbstractArray{T,2} = Base.covm(x',param_means'),
        param_means::AbstractArray{T,1} = vec(mean(x,2)), param_covar::PDMat = PDMat(Base.covm(x',param_means')),
			sigmasq_obs::Array{T,1} = zeros(T,size(x,2)), sigmasq_cor::T = one(T), rho::T = one(T))
  param_dist_sq = pairwise(SqMahalanobis(inv(param_covar).mat),x.-param_means)
  #K = param_dist_sq
  K  = diagm(sigmasq_obs)
  K += sigmasq_cor*kernel(param_dist_sq,rho=rho^2)
  return K
end

function make_kernel_data{T<:Real}(x::AbstractArray{T,2}, xstar::AbstractArray{T,2}, kernel::Function;
        #param_means::AbstractArray{T,1} = vec(mean(x,2)), param_covar::AbstractArray{T,2} = Base.covm(x',param_means'),
        param_means::AbstractArray{T,1} = vec(mean(x,2)), param_covar::PDMat = PDMat(Base.covm(x',param_means')),
			  sigmasq_cor::T = one(T), rho::T = one(T))
  param_dist_sq = pairwise(SqMahalanobis(inv(param_covar).mat),x.-param_means,xstar.-param_means)
  K = sigmasq_cor*kernel(param_dist_sq,rho=rho^2)
  return K
end

function train_gp{T<:Real}(x::AbstractArray{T,2}, y::AbstractArray{T,2}, sigmasq_y::AbstractArray{T,2}; kernel::Function = kernel_exp,
               sigmasq_cor::T = one(T), rho::T = one(T) )
  xm = vec(mean(x,2))
  xc = PDMat(Base.covm(x',xm')+diagm(1e-6*ones(size(x,1))))
  num_outputs = size(y,1)
  Sigma_tr = Array(PDMat,num_outputs)
  for i in 1:num_outputs   # TODO WARNING Update how sigmasq_obs is calculated
    Sigma_tr[i] = PDMat(make_kernel_data(x,kernel,param_means=xm,param_covar=xc,sigmasq_obs=sigmasq_y[i,:],sigmasq_cor=sigmasq_cor,rho=rho))
  end
  ym = vec(mean(y,2))
  return emulator_gp_type(sigmasq_cor,rho,kernel,x,xm,xc,y,ym,Sigma_tr)
end

function marginal_gp(emulator::emulator_gp_type, index::Integer)
  -0.5*( invquad(emulator.Sigma_train[index], emulator.output_train[index,:].-emulator.output_means[index]) + logdet(emulator.Sigma_train[index]) + size(emulator.input_train,2)*log(2pi) )
end

function marginal_gp(emulator::emulator_gp_type)
  logprob = 0.
  for i in 1:size(emulator.output_train,1)
    logprob += marginal_gp(emulator,i)
  end
  return logprob
end

function optimize_gp_cor{T<:Real}(x::AbstractArray{T,2}, y::AbstractArray{T,2}, sigmasq_y::AbstractArray{T,2}; kernel::Function = kernel_exp,
               sigmasq_cor::T = one(T), rho::T = one(T) )
  function calc_marginal(param::Array{T,1})
    emu = train_gp(x,y,sigmasq_y,kernel=kernel,sigmasq_cor=param[1],rho=param[2])
    marginal_gp(emu)
  end
  corlist = [0.001,0.003,0.01,0.03,0.1,0.3,1.0,3.0,10.]
  rholist = [0.001,0.003,0.01,0.03,0.1,0.3,1.0,3.0,10.]
  cor_best = 0.0
  rho_best = 0.0
  marginal_best = -Inf
  for i in 1:length(corlist)
  for j in 1:length(rholist)
      marg = calc_marginal([corlist[i],rholist[j]])
      println("# sigma_c=",corlist[i], " rho=",rholist[j]," logp=",marg)
      if marg > marginal_best
        cor_best = corlist[i]
        rho_best = rholist[j]
        marginal_best = marg
      end
    end
 end
 #return (cor_best, rho_best)
 train_gp(x,y,sigmasq_y,kernel=kernel,sigmasq_cor=cor_best,rho=rho_best)
end

function predict_gp{T<:Real}(emulator::emulator_gp_type, x_star::AbstractArray{T,2}; sigmasq_obs::Array{T,2} = zeros(T,(size(emulator.output_train,1),size(x_star,2))) )
  Sigma_c = make_kernel_data(emulator.input_train,x_star,emulator.kernel,param_means=emulator.input_means,param_covar=emulator.input_covar,sigmasq_cor=emulator.sigmasq_cor,rho=emulator.rho)
  Sigma_pred = make_kernel_data(x_star,emulator.kernel,param_means=emulator.input_means,param_covar=emulator.input_covar,sigmasq_cor=emulator.sigmasq_cor,rho=emulator.rho)
  num_outputs = size(emulator.output_train,1)
  predict_mean = Array(T,(num_outputs,size(x_star,2)))
  predict_covar = Array(T,(num_outputs,size(x_star,2),size(x_star,2)))
  for i in 1:num_outputs
    predict_mean[i,:] = Sigma_c' * (emulator.Sigma_train[i] \ (emulator.output_train[i,:].-emulator.output_means[i])) .+ emulator.output_means[i]
    # TODO WARNING Update how sigmasq_obs is calculated
    predict_covar[i,:,:] = Sigma_pred + diagm(sigmasq_obs[i,:]) - Sigma_c' * (emulator.Sigma_train[i] \ Sigma_c)
  end
  return (predict_mean, predict_covar)
end

function test_emulator()
  num_train = 100
  num_test = 101
  num_inputs = 3
  num_outputs = 2
  x = rand(num_inputs,num_train)
  xs = rand(num_inputs,num_test)
  A = zeros(num_outputs,num_inputs)
  for i in 1:min(num_outputs,num_inputs)
    A[i,i] = 1.0
  end
  A += 0.1randn(num_outputs,num_inputs)
  b = 3.0*rand(num_outputs)
  sigma = 0.1
  y  = b.+A*x .+sigma*rand(num_outputs,num_train)
  ys = b.+A*xs.+sigma*rand(num_outputs,num_test)
  emu = ABC.train_gp(x,y,sigma^2*reshape(ones(length(y)),size(y)))

  mar = ABC.marginal_gp(emu)
  (em_y, em_sigma) = ABC.predict_gp(emu,xs)

  resid = em_y-ys
  resid_std = similar(resid)
  for i in 1:size(ys,1)
    resid_std[i,:] = (em_y[i,:]-ys[i,:])./diag(em_sigma[i,:,:])
  end
  resid_std

  index_x = 1
  index_y = 1
  plot(vec(x[index_x,:]),vec(y[index_y,:]),"c.")
  plot(vec(xs[index_x,:]),vec(ys[index_y,:]),"b.")
  plt[:errorbar](vec(xs[index_x,:]),vec(em_y[index_y,:]),yerr=diag(em_sigma[index_y,:,:]),fmt=".")

  # plot residuals
  for index_x in 1:size(xs,1)
     for index_y in 1:size(ys,1)
       #plot(vec(xs[index_x,:]),vec(em_y[index_y,:]).-vec(ys[index_y,:]),".")
       plt[:errorbar](vec(xs[index_x,:]),vec(em_y[index_y,:]).-vec(ys[index_y,:]),yerr=diag(em_sigma[index_y,:,:]),fmt=".")
    end
  end

end

#=
(x_train,y_train) = ABC.make_training_data(pop_out.accept_log,pop_out.reject_log,10,100)
emu = ABC.optimize_gp_cor(x_train,y_train,y_train.+1.0)
x_pred = x_train
(y_pred,sigma_y) = ABC.predict_gp(emu,x_pred)
plot(x_train[1,:],y_train[3,:],"b.")
plot(x_pred[2,:],y_pred[3,:],"r.")

=#

function test_abc_emulator(accept_log::ABC.abc_log_type, reject_log::ABC.abc_log_type, generation::Integer, num_use::Integer = 10)
  (x_train,y_train) = make_training_data(accept_log, reject_log, generation, num_use)
  emu = train_gp(x_train,y_train,y_train)
  mar = marginal_gp(emu)
  println("# logprob = ",mar)
  x_pred = x_train
  y_pred = predict_gp(emu,x_pred)
  return (x_train,y_train,x_pred,y_pred)
end


text:
using Gnuplot, Random
Random.seed!(145)
x, y, yerr = 1:2:20, 5*rand(10), 0.4*abs.(randn(10))
@gp x y yerr "w boxerrorbars t 'test' lc 'black'"
@gp :- xrange= (0,22) yrange = (0,6)
save(term="pngcairo size 600,400", output="./code/plt_errs1.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
using TeraSort
yarncm = startcluster(false)
try
    @everywhere using TeraSort
    terasort()
finally
    stopcluster(yarncm)
end


text:
ENV["PLOTS_DEFAULT_BACKEND"] = "Plotly"
data = readdlm("../data/rvs_0001.txt");

want_to_plot_pgram = false
if want_to_plot_pgram
  if !isdefined(:Plots) using Plots end
  if !isdefined(:LombScargle) using LombScargle end
end

num_planets = 0
include("../src/rv_model.jl")    # provides logtarget

# Now with Keplerian model & GP noise
import RvModelKeplerianGPNoise
RvModelKeplerianGPNoise.set_times(data[:,1]);
RvModelKeplerianGPNoise.set_obs(data[:,2]);
RvModelKeplerianGPNoise.set_sigma_obs(data[:,3]);
num_offsets = 1
param_init_kepler_gp = zeros(num_planets*RvModelKeplerianGPNoise.num_param_per_planet+num_offsets+RvModelKeplerianGPNoise.num_noise_param())
for p in 1:num_planets
  (P,K,M0) = RvModelCircularGPNoise.extract_PKM0(param_init_circ_gp, plid=p)
  e = 0.01
  w = 0.0
  RvModelKeplerianGPNoise.set_PKewM0(param_init_kepler_gp, P, K, e, w, M0; plid=p)
end
RvModelKeplerianGPNoise.set_rvoffset(param_init_kepler_gp, 0.0)
RvModelKeplerianGPNoise.set_noise_param(param_init_kepler_gp, [0.0] )

if !isdefined(:QuadGK) using QuadGK end
logtarget_ref = RvModelKeplerianGPNoise.logtarget(param_init_kepler_gp)

param_kepler_gp = param_init_kepler_gp
function integrand_rvoffset(x::Real)
  RvModelKeplerianGPNoise.set_rvoffset(param_kepler_gp, x)
  exp(RvModelKeplerianGPNoise.logtarget(param_kepler_gp)-logtarget_ref)
end
function integrand_jitter(x::Real)
  RvModelKeplerianGPNoise.set_noise_param(param_kepler_gp, [x])
  exp(RvModelKeplerianGPNoise.logtarget(param_kepler_gp)-logtarget_ref)
 end

function marginalized_over_rvoffset(x::Real)
  RvModelKeplerianGPNoise.set_noise_param(param_kepler_gp, [x])
  f(y) = integrand_rvoffset(y)
  QuadGK.quadgk(f,-1000.0,1000.0)
end

g(z) = marginalized_over_rvoffset(z::Real)[1]

logevidence = log(QuadGK.quadgk(g,0.0,99.0)[1])
result = (logtarget_ref+logevidence)/log(10.0)


text:
include("getset.jl")

# Get all of the effects on all vehicles for the current state.
function get_effects(t, X, V, D, U, scenario)

    # Get all of the effects. Start with the vehicle body. Planets can consume
    # vehicle effects. Components can consume vehicle and planet effects.

    E = Dict{String, Dict{String, Tuple}}() # Instead of starting fresh, use the existing.
    c = 0 # TODO: Use something like component.id for this?
    for outer c = 1:length(scenario.vehicles)
        E[scenario.vehicles[c].name] = Dict{String, Tuple}()
        E[scenario.vehicles[c].name][scenario.vehicles[c].body.name] =
            scenario.vehicles[c].body.effects(
                t, scenario.vehicles[c].body.constants, X[c], nothing, V[c])
    end

    # Next, we'll get the effect of the environment on each body. While looping
    # over the bodies, we don't want the effects from the first body to be
    # visible to the second, so we'll need to make a copy here. We only need a
    # copy and not a deep copy, because the old effects won't be changed; we
    # simply need a dictionary that doesn't know about the new ones.
    effects_bus_copy = copy(E)

    # Get the planet's effects, given the body's effects.
    for vehicle in scenario.vehicles
        if scenario.environment.effects != nothing
            E[vehicle.name][scenario.environment.name] = scenario.environment.effects(t,
                                     scenario.environment.constants,
                                     scenario.environment.state, # TODO: Read from X once environments are allowed to have continuous-time states.
                                     draw(scenario.environment.rand, :effects), # TODO: This draw doesn't belong here. Same below.
                                     # TODO: Environments could have implicit variables.
                                     effects_bus_copy[vehicle.name],
                                     effects_bus_copy)
        end
    end

    # Get the actuators' effects, given the bodies' effects and planet's effects.
    effects_bus_copy = copy(E)
    for vehicle in scenario.vehicles

        # Get all of the component and computer board effects.
        for component in vehicle.components
            c += 1
            if component.effects != nothing
                E[vehicle.name][component.name] = component.effects(t, component.constants,
                                                                    X[c], D[c], V[c],
                                                                    U[vehicle.name][component.name],
                                                                    effects_bus_copy[vehicle.name], effects_bus_copy)
            end
        end
        for computer in vehicle.computers
            c += 1
            if computer.board.effects != nothing
                E[vehicle.name][computer.name] = computer.board.effects(t, computer.board.constants,
                                                                        X[c], D[c], V[c],
                                                                        U[vehicle.name][computer.name][computer.board.name],
                                                                        effects_bus_copy[vehicle.name], effects_bus_copy)
            end
        end

    end

    return E

end

# Calculate the Jacobian of the constraint function wrt the state using central finite differences.
function dgdX(t, X, V, D, U, scenario, α=1e-6) # Use a ! in the name? It modifies X, but should return it to its original.
    c  = zeros(scenario.dims.nv) # Constraints
    J  = zeros(scenario.dims.nv, scenario.dims.nxd) # Jacobian
    for k = 1:scenario.dims.nx
        z = scenario.dims.x_to_xd[k] # Next unused column of J
        if z != 0 # If this k corresponds to a continuous-time state...
            xo = getk(X, k) # Store the original value.
            setk!(X, k, xo + α) # Perturb it.
            stackem!(c, dae(t, X, V, D, U, scenario), 0) # Get constraints.
            J[:,z] = c[:] # Copy from the vector to the right place in the Jacobian. Do this with views instead?
            setk!(X, k, xo - α) # And now do the negative.
            stackem!(c, dae(t, X, V, D, U, scenario), 0)
            J[:,z] = (J[:,z] - c) / (2α)
            setk!(X, k, xo) # Return x to what it was.
        end
    end
    return J
end

# Calculate the Jacobian of the derivatives function wrt the implicit variables using central finite differences.
function dfdV(t, X, V, D, U, scenario, α=1e-6) # Use a ! in the name? It modifies V, but should return it to its original.
    xd = zeros(scenario.dims.nxd)
    J  = zeros(scenario.dims.nxd, scenario.dims.nv)
    for k = 1:scenario.dims.nv
        vo = getk(V, k) # Store the original value
        setk!(V, k, vo + α) # Perturb in the positive direction.
        stackem!(xd, ode(t, X, V, D, U, scenario), 0) # Replace with a version that stacks directly into J[:,k].
        J[:,k] = xd[:] # Copy from the vector to the right place in the Jacobian. Do this with views instead?
        setk!(V, k, vo - α) # And now the negative.
        stackem!(xd, ode(t, X, V, D, U, scenario), 0)
        J[:,k] = (J[:,k] - xd) / (2α)
        setk!(V, k, vo) # Return V to what it was.
    end
    return J
end

# Levenberg-Marquardt root solving for the implicity parameters, V, that satisfy the constraint functions.
function solve_for_V(t, X, V, D, U, scenario, Gx, λ=1e-6, tol=1e-12, k_max=1000)
    Xd = Vector{Any}() # This is just here for scope.
    xd = zeros(scenario.dims.nxd) # Preallocated space for "stacked" derivative
    k  = 0
    while true

        Xd = ode(t, X, V, D, U, scenario) # Get the big set derivatives.
        stackem!(xd, Xd, 0) # Create a stack of Float64s for it.
        cd = Gx * xd # Calculate the time-derivative of the constraints.
        if all(abs.(cd) .< tol); break; end # See if the constraint function is satisfied.
        if k == k_max; error("Could not satisfy the constraints functions after k iterations on the implicit variables."); end
        k += 1 # Ok, we're committing to this.

        # Calculate the Jacobian of the time-derivative of the constraints wrt the
        # implicit variables and use a Levenberg-Marquardt search for the implicit
        # parameters driving the time-derivative of the constraints to zero.
        J   =  Gx * dfdV(t, X, V, D, U, scenario)
        JtJ =  transpose(J) * J
        A   =  JtJ + λ * diagm(0 => diag(JtJ))
        σ   =  -(pinv(A) * (transpose(J) * cd)) # Use pinv in case of mischief.

        # Update each element of the implicit variables.
        for c = 1:length(σ); setk!(V, c, getk(V, c) + σ[c]); end # Could make a single updatek!(V, c, foo) function.

    end
    (V, Xd)
end

# Get the derivatives of the entire system when there are algebraic constraints.
function minor(t, X, V, D, U, scenario)

    # If there are algebraic constraints, solve for the implicit forces.
    if scenario.dims.nv > 0

        # Let g be the constraint function.
        #
        # If 0 = g(t, x), then 0 = d/dt g, and by the chain rule 0 = d/dx g * dx/dt.
        # We have a function for dx/dt given the implicit variables, v: dx/dt = f(t, x, v).
        # We can figure out the Jacobian, Gx = d/dx g, using finite differencing.
        # Then, we can solve for v s.t. Gx * f(t, x, v) = 0.
        #
        # Constraints and continuous equations will always be solved together. Should the contraints
        # be the second output of the `derivatives` function? No, that would make everyone have to
        # output a tuple. Plus, it's a toss up between contraints functions having duplicated work
        # from the derivatives function and the derivatives function being commonly "heavier" than
        # the constraints function (and being often unnecessary).

        # First, solve the constraints function for whatever part of V is observable.
        # V = solve_for_V_index_0(t, X, V, D, U, scenario) # Use SVD and if all S == 0, then there are no index 0 constraints.

        # Now differentiate the constraints function wrt the state.
        Gx = dgdX(t, X, V, D, U, scenario, 1e-9)

        # Now serach for implicit variables that put the state time derivative in the null
        # space of the constraint Jacobian. This solves for "index 1" constraints.
        V, Xd = solve_for_V(t, X, V, D, U, scenario, Gx)

    # Otherwise, just pass through whatever nothingness we're using for them.
    else
        Xd = ode(t, X, V, D, U, scenario)
    end

    (Xd, V)
end

# Get the vector of constraints.
function dae(t, X, V, D, U, scenario)

    # Get all of the effects corresponding to this state.
    E = get_effects(t, X, V, D, U, scenario)

    # Create the set of constraints that we will write to below.
    constraints = Vector{Float64}()

    # Get the constraints for each physical thing. We need to do this in the same
    # order as the states and draws arrays were built.
    c = 0
    for vehicle in scenario.vehicles
        c += 1
        if vehicle.body.constraints != nothing
            append!(constraints, vehicle.body.constraints(t, vehicle.body.constants,
                                                          X[c], D[c], V[c],
                                                          E[vehicle.name], E))
        end
    end

    for vehicle in scenario.vehicles
        for component in vehicle.components
            c += 1
            if component.constraints != nothing
                asdf = component.constraints(t, component.constants,
                    X[c], D[c], V[c],
                    U[vehicle.name][component.name],
                    E[vehicle.name], E)
                append!(constraints, asdf)
            end
        end
        for computer in vehicle.computers
            c += 1
            if computer.board.constraints != nothing
                append!(constraints, computer.board.constraints(t, computer.board.constants,
                                                                X[c], D[c], V[c],
                                                                U[vehicle.name][computer.name][computer.board.name],
                                                                E[vehicle.name], E))
            end
        end

    end

    return constraints

end

# Get the vector of derivatives.
function ode(t, X, V, D, U, scenario)

    # Get all of the effects corresponding to this state.
    E = get_effects(t, X, V, D, U, scenario)

    # Create the set of derivatives that we will write to below.
    derivs = Vector{Any}(undef, length(X)) # Would use eltype(X) but can't if we want to use nothing.
    fill!(derivs, nothing)

    # Get the derivatives for each physical thing. We need to do this in the same
    # order as the states and draws arrays were built. The derivative and the state
    # should be of the same type.
    c = 0
    for vehicle in scenario.vehicles
        c += 1
        if vehicle.body.derivatives != nothing
            derivs[c] = vehicle.body.derivatives(t, vehicle.body.constants,
                                                 X[c], D[c], V[c],
                                                 E[vehicle.name], E)
        end
    end

    for vehicle in scenario.vehicles
        for component in vehicle.components
            c += 1
            if component.derivatives != nothing
                derivs[c] = component.derivatives(t, component.constants,
                                                  X[c], D[c], V[c],
                                                  U[vehicle.name][component.name],
                                                  E[vehicle.name], E)
            end
        end
        for computer in vehicle.computers
            c += 1
            if computer.board.derivatives != nothing
                derivs[c] = computer.board.derivatives(t, computer.board.constants,
                                                       X[c], D[c], V[c],
                                                       U[vehicle.name][computer.name][computer.board.name],
                                                       E[vehicle.name], E)
            end
        end

    end

    return derivs

end


text:
Endpoint("/examples") do request::HTTP.Request
    readstring(joinpath(dirname(@__FILE__),"examples.html"))
end

Endpoint("/examples/pages") do request::HTTP.Request
    readstring(joinpath(dirname(@__FILE__),"pages.html"))
end

include("plotly.jl")
include("requests.jl")
# include("mwe.jl")

function examples()

@async Pages.start()
sleep(2.0)
Pages.launch("http://localhost:$(Pages.port)/examples")

end


text:
module TestImpurity
using Test
using StrayCopses
using StrayCopses: information_gain, information_impurity, gini_index, gini_impurity

@testset "Basic information gain" begin
    @test information_gain([3, 3], 6) == 1.0
    @test information_gain([2, 6], 8) ≈ 0.8112781244591328
    information_before = information_gain([5, 9], 14)
    @test information_impurity(information_before, [3, 3], [2, 6], 6, 8, 14) ≈ 0.04812703040826927
end

@testset "Basic gini index" begin
    @test gini_index([4, 2], 6) ≈ 0.555555555555555555555
    @test gini_index([2, 2], 4) ≈ 0.5
    @test gini_index([2, 0], 2) ≈ 1
    gini_before = gini_index([4, 2], 6)
    @test gini_impurity(gini_before, [2, 2], [2, 0], 4, 2, 6) ≈ 0.111111111111111111111
end

@testset "Edge case gini index" begin
    gini_before = gini_index([2, 2], 4)
    @test gini_impurity(gini_before, [0, 0], [2, 2], 0, 4, 4) ≈ 0.
    @test gini_index([0, 0], 0) ≈ 0.
end

end # module


text:
using Documenter, JGE

makedocs(;
    modules = [JGE],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        collapselevel = 1,
        sidebar_sitename = false
    ),
    pages = [
        "Home" => "index.md",
        "JGE" => "jge.md",
        "Developer" => [
            "dev.md",
            "Core" => "dev/core.md",
            "Audio" => "dev/audio.md",
            "Graphics" => "dev/graphics.md",
            "Networks" => "dev/networks.md",
            "Profiling" => "dev/profiling.md",
            "Debugging" => "dev/debugging.md",
            "Build" => "dev/build.md",
        ]
    ],
    repo = "https://github.com/serenity4/JGE.jl/blob/{commit}{path}#L{line}",
    sitename = "JGE.jl",
    authors = "Cédric Belmant",
)

deploydocs(;
    repo = "github.com/serenity4/JGE.jl",
)


text:
"""
    downward_recursion(x, fn,  n)

Takes as an input value of Boys function (`fn`) for given `n` value
and returns value of Boys function for one lower value.

Fₙ(x) -> Fₙ₋₁(x)
"""
function downward_recursion(x, fn,  n)
    upper = 2x * fn + exp(-x)
    lower = 2(n-1) + 1
    return upper / lower
end


text:
module SpringCollab2020BubblePushField

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/bubblePushField" BubblePushField(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight,
	strength::Number=1.0, upwardStrength::Number=1.0, direction::String="Right", water::Bool=true, flag::String="bubble_push_field", activationMode::String="Always")

const placements = Ahorn.PlacementDict(
	"Bubble Column (Spring Collab 2020)" => Ahorn.EntityPlacement(
		BubblePushField,
		"rectangle"
	)
)

Ahorn.editingOptions(entity::BubblePushField) = Dict{String,Any}(
	"direction" => ["Up", "Down", "Left", "Right"],
	"activationMode" => ["Always", "OnlyWhenFlagActive", "OnlyWhenFlagInactive"]
)

Ahorn.minimumSize(entity::BubblePushField) = 8, 8
Ahorn.resizable(entity::BubblePushField) = true, true

Ahorn.selection(entity::BubblePushField) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::BubblePushField, room::Maple.Room)
	x = Int(get(entity.data, "x", 0))
	y = Int(get(entity.data, "y", 0))

	width = Int(get(entity.data, "width", 32))
	height = Int(get(entity.data, "height", 32))

	Ahorn.drawRectangle(ctx, 0, 0, width, height, (0.7, 0.28, 0.0, 0.34), (1.0, 1.0, 1.0, 0.5))
end

end


text:
# TODO this file contains many duplications with the inlining analysis code, factor them out

import Core.Compiler:
    MethodInstance, InferenceResult, Signature, ConstResult,
    MethodResultPure, MethodMatchInfo, UnionSplitInfo, ConstCallInfo, InvokeCallInfo,
    call_sig, argtypes_to_type, is_builtin, is_return_type, istopfunction, validate_sparams,
    specialize_method, invoke_rewrite

const Linfo = Union{MethodInstance,InferenceResult}
struct CallInfo
    linfos::Vector{Linfo}
    nothrow::Bool
end

function resolve_call(ir::IRCode, stmt::Expr, @nospecialize(info))
    sig = call_sig(ir, stmt)
    if sig === nothing
        return missing
    end
    # TODO handle _apply_iterate
    if is_builtin(sig) && sig.f !== invoke
        return false
    end
    # handling corresponding to late_inline_special_case!
    (; f, argtypes) = sig
    if length(argtypes) == 3 && istopfunction(f, :!==)
        return true
    elseif length(argtypes) == 3 && istopfunction(f, :(>:))
        return true
    elseif f === TypeVar && 2 ≤ length(argtypes) ≤ 4 && (argtypes[2] ⊑ Symbol)
        return true
    elseif f === UnionAll && length(argtypes) == 3 && (argtypes[2] ⊑ TypeVar)
        return true
    elseif is_return_type(f)
        return true
    end
    if info isa MethodResultPure
        return true
    elseif info === false
        return missing
    end
    # TODO handle OpaqueClosureCallInfo
    if sig.f === invoke
        isa(info, InvokeCallInfo) || return missing
        return analyze_invoke_call(sig, info)
    elseif isa(info, ConstCallInfo)
        return analyze_const_call(sig, info)
    elseif isa(info, MethodMatchInfo)
        infos = MethodMatchInfo[info]
    elseif isa(info, UnionSplitInfo)
        infos = info.matches
    else # isa(info, ReturnTypeCallInfo), etc.
        return missing
    end
    return analyze_call(sig, infos)
end

function analyze_invoke_call(sig::Signature, info::InvokeCallInfo)
    match = info.match
    if !match.fully_covers
        # TODO: We could union split out the signature check and continue on
        return missing
    end
    result = info.result
    if isa(result, InferenceResult)
        return CallInfo(Linfo[result], true)
    else
        argtypes = invoke_rewrite(sig.argtypes)
        mi = analyze_match(match, length(argtypes))
        mi === nothing && return missing
        return CallInfo(Linfo[mi], true)
    end
end

function analyze_const_call(sig::Signature, cinfo::ConstCallInfo)
    linfos = Linfo[]
    (; call, results) = cinfo
    infos = isa(call, MethodMatchInfo) ? MethodMatchInfo[call] : call.matches
    local nothrow = true # required to account for potential escape via MethodError
    local j = 0
    for i in 1:length(infos)
        meth = infos[i].results
        nothrow &= !meth.ambig
        nmatch = Core.Compiler.length(meth)
        if nmatch == 0 # No applicable methods
            # mark this call may potentially throw, and the try next union split
            nothrow = false
            continue
        end
        for i = 1:nmatch
            j += 1
            result = results[j]
            match = Core.Compiler.getindex(meth, i)
            if result === nothing
                mi = analyze_match(match, length(sig.argtypes))
                mi === nothing && return missing
                push!(linfos, mi)
            elseif isa(result, ConstResult)
                # TODO we may want to feedback information that this call always throws if !isdefined(result, :result)
                push!(linfos, result.mi)
            else
                push!(linfos, result)
            end
            nothrow &= match.fully_covers
        end
    end
    return CallInfo(linfos, nothrow)
end

function analyze_call(sig::Signature, infos::Vector{MethodMatchInfo})
    linfos = Linfo[]
    local nothrow = true # required to account for potential escape via MethodError
    for i in 1:length(infos)
        meth = infos[i].results
        nothrow &= !meth.ambig
        nmatch = Core.Compiler.length(meth)
        if nmatch == 0 # No applicable methods
            # mark this call may potentially throw, and the try next union split
            nothrow = false
            continue
        end
        for i = 1:nmatch
            match = Core.Compiler.getindex(meth, i)
            mi = analyze_match(match, length(sig.argtypes))
            mi === nothing && return missing
            push!(linfos, mi)
            nothrow &= match.fully_covers
        end
    end
    return CallInfo(linfos, nothrow)
end

function analyze_match(match::MethodMatch, npassedargs::Int)
    method = match.method
    na = Int(method.nargs)
    if na != npassedargs && !(na > 0 && method.isva)
        # we have a method match only because an earlier
        # inference step shortened our call args list, even
        # though we have too many arguments to actually
        # call this function
        return nothing
    end

    # Bail out if any static parameters are left as TypeVar
    # COMBAK is this needed for escape analysis?
    validate_sparams(match.sparams) || return nothing

    # See if there exists a specialization for this method signature
    mi = specialize_method(match; preexisting=true) # Union{Nothing, MethodInstance}
    return mi
end

max_stars_repo_name: JuliaLabs/julia

max_issues_repo_name: JuliaLabs/julia

max_forks_repo_name: JuliaLabs/julia


text:
struct SGD_prox_iterable{R<:Real,C<:RealOrComplex{R},Tx<:AbstractArray{C},Tf,Tg}
    F::Array{Tf}            # smooth term
    g::Tg                   # nonsmooth term
    x0::Tx                  # initial point
    N::Int                  # of data points in the finite sum problem
    L::Maybe{Union{Array{R},R}}  # Lipschitz moduli of nabla f_i
    μ::Maybe{Union{Array{R},R}}  # convexity moduli of the gradients
    γ::Maybe{R}             # stepsize
    plus::Bool              # plus version (diminishing stepsize)
end

mutable struct SGD_prox_state{R<:Real,Tx}
    γ::R                    # stepsize
    z::Tx
    cind::Int               # current interation index
    idxr::Int               # current index
    # some extra placeholders
    ∇f_temp::Tx             # placeholder for gradients
    temp::Tx
end

function SGD_prox_state(γ::R, z::Tx, cind) where {R,Tx}
    return SGD_prox_state{R,Tx}(γ, z, cind, Int(0), copy(z), copy(z))
end

function Base.iterate(iter::SGD_prox_iterable{R}) where {R}
    N = iter.N
    ind = collect(1:N)
    # updating the stepsize
    if iter.γ === nothing
        if iter.L === nothing
            @warn "smoothness or convexity parameter absent"
            return nothing
        else
            L_M = maximum(iter.L)
            γ = 1/ (2*L_M)
        end
    else
        γ = iter.γ # provided γ
    end
    if iter.plus
        println("diminishing stepsize version")
    end
    # initializing the vectors
    ∇f_temp = zero(iter.x0)
    z = zero(iter.x0)
    z .= iter.x0
    cind = 0
    for i = 1:N
        ∇f, ~ = gradient(iter.F[i], z)
        if iter.plus
            cind += 1
            γ = 0.1/(1+cind/N)
        end
        ∇f ./= N
        ∇f .*= - γ
        ∇f .+= z
        CIAOAlgorithms.prox!(z, iter.g, ∇f, γ)
    end

    state = SGD_prox_state(γ, z, cind)
    return state, state
end

function Base.iterate(iter::SGD_prox_iterable{R}, state::SGD_prox_state{R}) where {R}
    # The inner cycle

    for i=1:iter.N # just for speed in implementation
        state.cind += 1
        if iter.plus
            state.γ = 0.1/(1+state.cind/iter.N)
        end

        gradient!(state.∇f_temp, iter.F[i], state.z)

        state.∇f_temp .*= - state.γ
        state.∇f_temp ./= iter.N
        state.∇f_temp .+= state.z

        CIAOAlgorithms.prox!(state.z, iter.g, state.∇f_temp, state.γ)
    end

    return state, state
end

solution(state::SGD_prox_state) = state.z


text:
# ---
# title: Shear
# cover: shear.gif
# description: shear the input image horizontally or vertically
# ---

# [`ShearX`](@ref)/[`ShearY`](@ref) can be used to shear the input image horizontally/vertically.
# The input to the constructor can be a scalar or a vector. In the case of a vector, the transformation
# will be a stochastic process.

using Augmentor
using ImageShow, ImageCore
using Random
Random.seed!(0)
img_in = testpattern(RGB, ratio=0.5)

mosaicview(
    ## deterministic transformation
    augment(img_in, ShearX(20)),
    augment(img_in, ShearY(20)),

    ## random transformation
    augment(img_in, ShearX(-20:20)),
    augment(img_in, ShearY(-20:20));

    fillvalue=colorant"white", nrow=2, npad=10
)

# Note that the output image size will be changed after transformation, [`CropNative`](@ref) can be
# particalually useful to preserve the image size.

mosaicview(
    augment(img_in, ShearX(10)),
    augment(img_in, ShearX(10) |> CropNative(axes(img_in)));
    fillvalue=colorant"white", nrow=1, npad=10
)

# ## References

#md # ```@docs
#md # ShearX
#md # ShearY
#md # ```

## save covers #src
using ImageMagick #src
include(joinpath("..", "assets", "utilities.jl")) #src
cover = make_gif(testpattern(RGB, ratio=0.5), Either(ShearX(-10:10), ShearY(-10:10)), 5) #src
ImageMagick.save("shear.gif", cover; fps=1) #src


text:
"""
Abstract type for storing Neuroimaing data.

All other neuroimaging types inherit from this type.
All neuroimaing types support the following functions:

* `samplingrate()`
* `channelnames()`
* `remove_channel!()`
* `keep_channel!()`
* `trim_channel()`
* `highpass_filter()`
* `lowpass_filter()`
* `data()`
* `plot()`

# Examples
```julia
data = # load your neuroimaging data
samplingrate(data)  # Returns the sampling rate
channelnames(data)  # Returns the channel names
```
"""
abstract type NeuroimagingMeasurement end

#######################################
#
# Filtering
#
# We provide maximum flexibility to users by providing filt and filtfilt
# to be performed on any NeuroimagingMeasurement.
# On top of this flexible default, methods are provided for each subtype
# of EEG measurements in their associated directories.
#
#######################################

function filt(
    s::NeuroimagingMeasurement,
    fcoef::Union{FilterCoefficients,AbstractVector{F}},
) where {F<:AbstractFloat}
    s2 = deepcopy(s)
    s2.data = DSP.filt(fcoef, data(s2))
    return s2
end

function filtfilt(
    s::NeuroimagingMeasurement,
    fcoef::Union{FilterCoefficients,AbstractVector{F}},
) where {F<:AbstractFloat}
    s2 = deepcopy(s)
    s2.data = DSP.filtfilt(fcoef, data(s2))
    return s2
end

function filter(
    s::NeuroimagingMeasurement,
    fcoef::Union{FilterCoefficients,AbstractVector{F}},
    phase::AbstractString,
) where {F<:AbstractFloat}
    if phase == "zero-double"

        s = filtfilt(s, fcoef)

    elseif phase == "zero"

        # # append filterdelay as zeros
        tau = filterdelay(fcoef)
        if ndims(s.data) > 1 # todo: this feels hacky
            # we want to add zeros at the end to get a longer signal for the filter, so that we can then correct the filter delay
            s_dim2 = size(s.data)[2]
        else
            s_dim2 = 1
        end
        s.data = vcat(s.data, zeros(tau * 2, s_dim2))

        # # filter
        s = filt(s, fcoef)

        # # fix grpdelay
        s.data = s.data[tau+1:end-tau, :]

    else
        throw(
            ArgumentError(
                "Unknown filter phase argument $(phase). Must be one of zero-double or zero or ...",
            ),
        )
    end
    return s
end

#######################################
#
# Plotting
#
# Use recipes for plotting.
# Here a base recipe is defined.
# Specific types can overide this.
#
#######################################

"""
    plot(n::NeuroimagingMeasurement)
    plot(n::NeuroimagingMeasurement, c::AbstractString)
    plot(n::NeuroimagingMeasurement, c::AbstractVector{AbstractString})

Plot the time series of the neuroimaing measurement.
If specified, only a selection of channels will be plotted.

# Examples
```julia
measurement = # load your neuroimaging data
plot(measurement)
plot(measurement, "TP7")
plot(measurement, ["TP7", "Cz"])
```
"""
@recipe function plot(s::NeuroimagingMeasurement)

    time_s = times(Float64, s)

    if size(data(s), 2) == 1

        labs = pop!(plotattributes, :label, channelnames(s)[1])

        RecipesBase.@series begin
            xguide := "Time (s)"
            yguide := "Amplitude (uV)"
            label := labs

            time_s, data(s)
        end

    else

        signals = data(deepcopy(s))
        variances = var(signals, dims = 1)
        mean_variance = Statistics.mean(variances)
        for c = 1:size(signals, 2)
            signals[:, c] = signals[:, c] .- Statistics.mean(signals[:, c])
            signals[:, c] = signals[:, c] ./ (mean_variance ./ 4) .+ (c - 1)
        end

        RecipesBase.@series begin
            seriestype := :path
            xguide := "Time (s)"
            yguide := "Amplitude (uV)"
            label := false
            yformatter := y -> ""
            yticks := (0:length(channelnames(s))-1, channelnames(s))

            time_s, signals
        end
    end
end

@recipe function plot(s::NeuroimagingMeasurement, c::S) where {S<:AbstractString}
    keep_channel!(deepcopy(s), [c])
end

@recipe function plot(
    s::NeuroimagingMeasurement,
    c::AbstractVector{S},
) where {S<:AbstractString}
    keep_channel!(deepcopy(s), c)
end

#######################################
#
# Internal functions
#
#######################################

"""
Internal function to find indices for channel names
"""
function _channel_indices(
    s::NeuroimagingMeasurement,
    channels::AbstractVector{S};
    warn_on_missing = true,
) where {S<:AbstractString}
    c_idx = Int[something(findfirst(isequal(c1), channelnames(s)), 0) for c1 in channels]
    if warn_on_missing
        if any(c_idx .== 0)
            throw(KeyError("Requested channel does not exist in $(channelnames(s))"))
        end
    end
    return c_idx
end


text:
function AngularSeparationDetector(arg0::PVCoordinatesProvider, arg1::PVCoordinatesProvider, arg2::jdouble)
    return AngularSeparationDetector((PVCoordinatesProvider, PVCoordinatesProvider, jdouble), arg0, arg1, arg2)
end

function g(obj::AngularSeparationDetector, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function get_beacon(obj::AngularSeparationDetector)
    return jcall(obj, "getBeacon", PVCoordinatesProvider, ())
end

function get_observer(obj::AngularSeparationDetector)
    return jcall(obj, "getObserver", PVCoordinatesProvider, ())
end

function get_proximity_angle(obj::AngularSeparationDetector)
    return jcall(obj, "getProximityAngle", jdouble, ())
end


text:
function extendSeries(X::Array{Float64}, method::String, length::String, n::Int)

  (N, ) = size(X)

  # determine final length 'n' of series after extension
  if (length == "arbitrary")
    if (n <= N)
     error("Invalid argument: 'n' must be greater than length of series when length='arbitrary'.")
    end
  elseif (length == "powerof2")
    k = N/(2^n)
    if (round(k) == k)
      error("Invalid argument: length of series should not be multiple of 2^j when length='powerof2'.")
    else
      n = ceil(k)*2^n
    end
  elseif (length == "double")
    error("The integer parameter is arbitrary, you probably made some error.")
  end

  if isa(X, Vector)
    # extend the series to length 'n'
    (method == "periodic") ? X = repeat(X, outer=[convert(Int,ceil(n/N))])[1:n] : ()
    (method == "reflection") ? X = repeat([X,X[end:-1:1]], outer=[convert(Int,ceil(n/N))])[1:n] : ()
    (method == "zeros") ? X = [X, zeros(N-n)] : ()
    (method == "mean") ? X = [X, mean(X).*ones(N-n)] : ()
    # if(method == "reflection.inverse") X = apply(X, 2, function(x,n,N){rep(c(x,2*x[N]-rev(x)),length=n)}, n=n, N=N)
  else
    # extend the series to length 'n'
    (method == "periodic") ? X = mapslices(x -> repeat(x, outer=[convert(Int,ceil(n/N))])[1:n], X, 1) : ()
    (method == "reflection") ? X = mapslices(x -> repeat(x, outer=[convert(Int,ceil(n/N))])[1:n], [X,X[end:-1:1]], 2) : ()
    (method == "zeros") ? X = mapslices(x -> [x, zeros(N-n)], X, 2) : ()
    (method == "mean") ? X = mapslices(x -> [x, mean(x).*ones(N-n)], X, 2) : ()
    # if(method == "reflection.inverse") X = apply(X, 2, function(x,n,N){rep(c(x,2*x[N]-rev(x)),length=n)}, n=n, N=N)
  end

  return X
end

function extendSeries(X::Array{Float64}, method::String, length::String)
  (N, ) = size(X)

  if length == "double"
    n = 2*N
  else
    error("The selected method does not correspond, you probably forgot to include another parameter.")
  end

  if isa(X, Vector)
    # extend the series to length 'n'
    (method == "periodic") ? X = repeat(X, outer=[convert(Int,ceil(n/N))])[1:n] : ()
    (method == "reflection") ? X = repeat([X,X[end:-1:1]], outer=[convert(Int,ceil(n/N))])[1:n] : ()
    (method == "zeros") ? X = [X, zeros(N-n)] : ()
    (method == "mean") ? X = [X, mean(X).*ones(N-n)] : ()
    # if(method == "reflection.inverse") X = apply(X, 2, function(x,n,N){rep(c(x,2*x[N]-rev(x)),length=n)}, n=n, N=N)
  else
    # extend the series to length 'n'
    (method == "periodic") ? X = mapslices(x -> repeat(x, outer=[convert(Int,ceil(n/N))])[1:n], X, 1) : ()
    (method == "reflection") ? X = mapslices(x -> repeat(x, outer=[convert(Int,ceil(n/N))])[1:n], [X,X[end:-1:1]], 2) : ()
    (method == "zeros") ? X = mapslices(x -> [x, zeros(N-n)], X, 2) : ()
    (method == "mean") ? X = mapslices(x -> [x, mean(x).*ones(N-n)], X, 2) : ()
    # if(method == "reflection.inverse") X = apply(X, 2, function(x,n,N){rep(c(x,2*x[N]-rev(x)),length=n)}, n=n, N=N)
  end

  return X
end


text:
function layoutable(::Type{Slider}, fig_or_scene; bbox = nothing, kwargs...)

    topscene = get_topscene(fig_or_scene)

    default_attrs = default_attributes(Slider, topscene).attributes
    theme_attrs = subtheme(topscene, :Slider)
    attrs = merge!(merge!(Attributes(kwargs), theme_attrs), default_attrs)

    decorations = Dict{Symbol, Any}()

    @extract attrs (
        halign, valign, horizontal, linewidth, snap,
        startvalue, value, color_active, color_active_dimmed, color_inactive
    )

    sliderrange = attrs.range

    protrusions = Node(GridLayoutBase.RectSides{Float32}(0, 0, 0, 0))
    layoutobservables = LayoutObservables{Slider}(attrs.width, attrs.height, attrs.tellwidth, attrs.tellheight,
        halign, valign, attrs.alignmode; suggestedbbox = bbox, protrusions = protrusions)

    onany(linewidth, horizontal) do lw, horizontal
        if horizontal
            layoutobservables.autosize[] = (nothing, Float32(lw))
        else
            layoutobservables.autosize[] = (Float32(lw), nothing)
        end
    end

    sliderbox = lift(identity, layoutobservables.computedbbox)

    endpoints = lift(sliderbox, horizontal) do bb, horizontal

        h = height(bb)
        w = width(bb)

        if horizontal
            y = bottom(bb) + h / 2
            [Point2f0(left(bb) + h/2, y),
             Point2f0(right(bb) - h/2, y)]
        else
            x = left(bb) + w / 2
            [Point2f0(x, bottom(bb) + w/2),
             Point2f0(x, top(bb) - w/2)]
        end
    end

    # this is the index of the selected value in the slider's range
    # selected_index = Node(1)
    # add the selected index to the attributes so it can be manipulated later
    attrs.selected_index = 1
    selected_index = attrs.selected_index

    # the fraction on the slider corresponding to the selected_index
    # this is only used after dragging
    sliderfraction = lift(selected_index, sliderrange) do i, r
        (i - 1) / (length(r) - 1)
    end

    dragging = Node(false)

    # what the slider actually displays currently (also during dragging when
    # the slider position is in an "invalid" position given the slider's range)
    displayed_sliderfraction = Node(0.0)

    on(sliderfraction) do frac
        # only update displayed fraction through sliderfraction if not dragging
        # dragging overrides the value so there is clear mouse interaction
        if !dragging[]
            displayed_sliderfraction[] = frac
        end
    end

    # when the range is changed, switch to closest value
    on(sliderrange) do rng
        selected_index[] = closest_index(rng, value[])
    end

    on(selected_index) do i
        value[] = sliderrange[][i]
    end

    # initialize slider value with closest from range
    selected_index[] = closest_index(sliderrange[], startvalue[])

    middlepoint = lift(endpoints, displayed_sliderfraction) do ep, sf
        Point2f0(ep[1] .+ sf .* (ep[2] .- ep[1]))
    end

    linepoints = lift(endpoints, middlepoint) do eps, middle
        [eps[1], middle, middle, eps[2]]
    end

    linecolors = lift(color_active_dimmed, color_inactive) do ca, ci
        [ca, ci]
    end

    endbuttons = scatter!(topscene, endpoints, color = linecolors,
        markersize = linewidth, strokewidth = 0, raw = true, inspectable = false)
    decorations[:endbuttons] = endbuttons

    linesegs = linesegments!(topscene, linepoints, color = linecolors,
        linewidth = linewidth, raw = true, inspectable = false)
    decorations[:linesegments] = linesegs

    button_magnification = Node(1.0)
    buttonsize = @lift($linewidth * $button_magnification)
    button = scatter!(topscene, middlepoint, color = color_active, strokewidth = 0,
        markersize = buttonsize, raw = true, inspectable = false)
    decorations[:button] = button

    mouseevents = addmouseevents!(topscene, layoutobservables.computedbbox)

    onmouseleftdrag(mouseevents) do event
        dragging[] = true
        dif = event.px - event.prev_px
        fraction = clamp(if horizontal[]
            (event.px[1] - endpoints[][1][1]) / (endpoints[][2][1] - endpoints[][1][1])
        else
            (event.px[2] - endpoints[][1][2]) / (endpoints[][2][2] - endpoints[][1][2])
        end, 0, 1)

        newindex = closest_fractionindex(sliderrange[], fraction)
        if snap[]
            fraction = (newindex - 1) / (length(sliderrange[]) - 1)
        end
        displayed_sliderfraction[] = fraction

        if selected_index[] != newindex
            selected_index[] = newindex
        end

        return true
    end

    onmouseleftdragstop(mouseevents) do event
        dragging[] = false
        # adjust slider to closest legal value
        sliderfraction[] = sliderfraction[]
        linecolors[] = [color_active_dimmed[], color_inactive[]]
        return true
    end

    onmouseleftdown(mouseevents) do event
        pos = event.px
        dim = horizontal[] ? 1 : 2
        frac = (pos[dim] - endpoints[][1][dim]) / (endpoints[][2][dim] - endpoints[][1][dim])
        selected_index[] = closest_fractionindex(sliderrange[], frac)
        # linecolors[] = [color_active[], color_inactive[]]
        return true
    end

    onmouseleftdoubleclick(mouseevents) do event
        selected_index[] = closest_index(sliderrange[], startvalue[])
        return true
    end

    onmouseenter(mouseevents) do event
        button_magnification[] = 1.25
        return false
    end

    onmouseout(mouseevents) do event
        button_magnification[] = 1.0
        linecolors[] = [color_active_dimmed[], color_inactive[]]
        return false
    end

    # trigger autosize through linewidth for first layout
    linewidth[] = linewidth[]

    Slider(fig_or_scene, layoutobservables, attrs, decorations)
end

function valueindex(sliderrange, value)
    for (i, val) in enumerate(sliderrange)
        if val == value
            return i
        end
    end
    nothing
end

function closest_fractionindex(sliderrange, fraction)
    n = length(sliderrange)
    onestepfrac = 1 / (n - 1)
    i = round(Int, fraction / onestepfrac) + 1
    min(max(i, 1), n)
end

function closest_index(sliderrange, value)
    for (i, val) in enumerate(sliderrange)
        if val == value
            return i
        end
    end
    # if the value wasn't found this way try inexact
    closest_index_inexact(sliderrange, value)
end

function closest_index_inexact(sliderrange, value)
    distance = Inf
    selected_i = 0
    for (i, val) in enumerate(sliderrange)
        newdist = abs(val - value)
        if newdist < distance
            distance = newdist
            selected_i = i
        end
    end
    selected_i
end

"""
Set the `slider` to the value in the slider's range that is closest to `value` and return this value.
"""
function set_close_to!(slider::Slider, value)
    closest = closest_index(slider.range[], value)
    slider.selected_index = closest
    slider.range[][closest]
end