text:
# weave("linear_sys_tmp.jl", doctype="github", out_path="build")
using Plots
default(grid=false)
cd(@__DIR__)
using Parameters, JacProp, OrdinaryDiffEq, LTVModels, LTVModelsBase, Random, LinearAlgebra, Statistics, DSP
using Flux: params, jacobian
using Flux.Optimise: Param, optimiser, expdecay

function simaverage(x, n)
    N,p = size(x)
    res = map(1:p) do i
        mean(reshape(x[:,i], N÷n,:), dims=2)
    end
    hcat(res...)
end

@with_kw struct LinearSys <: AbstractSystem
    A
    B
    N     = 1000
    nx    = size(A,1)
    nu    = size(B,2)
    h     = 0.02
    σ0    = 0
    sind  = 1:nx
    uind  = nx+1:(nx+nu)
    s1ind = (nx+nu+1):(nx+nu+nx)
end

function LinearSys(seed; nx = 10, nu = nx, h=0.02, kwargs...)
    Random.seed!(seed)
    A = randn(nx,nx)
    A = A-A'        # skew-symmetric = pure imaginary eigenvalues
    A = A - h*I     # Make 'slightly' stable
    A = exp(h*A)   # discrete time
    B = h*randn(nx,nu)
    LinearSys(;A=A, B=B, nx=nx, nu=nu, h=h, kwargs...)
end

function generate_data(sys::LinearSys, seed, validation=false)
    Parameters.@unpack A,B,N, nx, nu, h, σ0 = sys
    Random.seed!(seed)
    u      = filt(ones(5),[5], 10randn(N+2,nu))'
    t      = h:h:N*h+h
    x0     = randn(nx)
    x      = zeros(nx,N+1)
    x[:,1] = x0
    for i = 1:N-1
        x[:,i+1] = A*x[:,i] + B*u[:,i]
    end

    validation || (x .+= σ0 * randn(size(x)))
    u = u[:,1:N]
    @assert all(isfinite, u)
    x,u
end

function true_jacobian(sys::LinearSys, x, u)
    [sys.A sys.B]
end

function callbacker(epoch, loss,d,trace,model,mt)
    i = length(trace) + epoch - 1
    function ()
        l = sum(d->Flux.data(loss(d...)),d)
        increment!(trace,epoch,l)
        # i % 500 == 0 && println(@sprintf("Loss: %.4f", l))
    end
end

num_params = 20
wdecay     = 0
stepsize   = 0.02
const sys  = LinearSys(1, nx = 10, N=200, h=0.02, σ0 = 0.01)
true_jacobian(x,u) = true_jacobian(sys,x,u)
nu         = sys.nu
nx         = sys.nx
kalmanopts = [:P => 10, :R2 => 1000I, :σdivider => 20]

#' Generate validation data
function valdata()
    vx,vu,vy = Vector{Float64}[],Vector{Float64}[],Vector{Float64}[]
    for i = 20:60
        x,u = generate_data(sys,i, true)
        for j in 10:5:(sys.N-1)
            push!(vx, x[:,j])
            push!(vy, x[:,j+1])
            push!(vu, u[:,j])
        end
    end
    hcat(vx...),hcat(vu...),hcat(vy...)
end
vx,vu,vy = valdata()
vt = Trajectory(vx,vu,vy)

#' Generate training trajectories
trajs = [Trajectory(generate_data(sys, i)...) for i = 1:3]

## Monte-Carlo evaluation
num_montecarlo = 4
it = num_montecarlo ÷ 2
res = map(1:num_montecarlo) do it
    trajs = [Trajectory(generate_data(sys, i)...) for i = it .+ (1:2)]

    #' ## Without jacprop
    Random.seed!(1)
    # wdecay = exp10.(range(-6, stop=1, length=num_montecarlo))[it]
    wdecay = 1e-4
    models     = [DiffSystem(nx,nu,num_params, a) for a in default_activations]
    opts       = [[ADAM.(params.(models), stepsize, ϵ=1e-2); [expdecay(Param(p), wdecay) for p in params(models[i]) if p isa AbstractMatrix]] for i = 1:length(models)]
    trainer  = ModelTrainer(;models=models, opts=opts, losses=JacProp.loss.(models), cb=callbacker, σdivider = wdecay)
    for i = 1:2
        trainer(trajs[i], epochs=0, jacprop=0, useprior=false)
    end
    # trainer(epochs=2000)
    # serializesave("trainer", trainer)
    # trainer = deserialize("trainer")
    #' ## With jacprop

    # Random.seed!(1)
    # models     = [DiffSystem(nx,nu,num_params, a) for a in default_activations]
    # opts       = ADAM.(params.(models), stepsize, decay=0.0005)
    # trainerjn  = ModelTrainer(;models=models, opts=opts, losses=JacProp.loss.(models), cb=callbacker, kalmanopts...)
    # for i = 1:2
    #     trainerjn(trajs[i], epochs=0, jacprop=1, useprior=false)
    # end
    # trainerjn(epochs=1000, jacprop=1)
    # serialize("trainerjn", trainerjn)
    # trainerjn = deserialize("trainerjn")

    #' ## With AD-jacprop
    Random.seed!(1)
    # λ = exp10.(range(-4, stop=3, length=num_montecarlo))[it]
    λ = 0.15
    model     = JacProp.ADDiffSystem(nx,nu,num_params,tanh) # TODO: tanh has no effect
    opt       = LTVModels.ADAMOptimizer.(model.w, α = stepsize)
    trainerad = ADModelTrainer(;model=model, opt=opt, λ=λ, testdata = vt)
    for i = 1:2
        trainerad(trajs[i], epochs=0)
    end
    # trainerad(epochs=2000)
    # serializesave("trainerad", trainerad)
    # trainerad = deserialize("trainerad")

    trainer, trainerad
end

# ress = getindex.(res,1)
# resad = getindex.(res,2)
resdiff = getindex.(res,1)
resaddiff = getindex.(res,2)
# λs = getfield.(ress, :σdivider)
# λad = getfield.(resad, :λ)

# Threads.@threads for trainer=ress       train!(trainer, epochs=2500) end
# @info("Done ress")
# Threads.@threads for trainer=resad      train!(trainer, epochs=2500) end
# @info("Done resad")
Threads.@threads for trainer=resdiff    train!(trainer, epochs=300) end
@info("Done resdiff")
Threads.@threads for trainer=resaddiff  train!(trainer, epochs=300) end
@info("Done resaddiff")

# serialize("results_linear_eigvals_final", (resdiff,resaddiff))
# (resdiff,resaddiff) = deserialize("results_linear_eigvals_final")

##

using JacProp: eval_pred, eval_sim, eval_jac, eval_jac2
λd = getfield.(resdiff, :σdivider)
λadd = getfield.(resaddiff, :λ)
xv = [λd λadd]
labelvec = ["f" "g"]
infostring = @sprintf("Num hidden: %d, sigma: %2.2f, Montecarlo: %d", num_params, sys.σ0, num_montecarlo)
simvals = [Trajectory(generate_data(sys, i, true)...) for i = 4:6]
pred  = [eval_pred.(resdiff, (vt,)) eval_pred.(resaddiff, (vt,))]
sim   = simaverage(vcat([[eval_sim.(resdiff, (simval,)) eval_sim.(resaddiff, (simval,))] for simval in simvals]...), 3)
# sim = [x > 10 ? NaN : x for x in sim]
jac = [eval_jac2.(resdiff, (vt,), true_jacobian,3) eval_jac2.(resaddiff, (vt,), true_jacobian,3)]

## Plot error vs λ regularization parameter
plot(xv,pred, layout=(3,1), subplot=1, ylabel="Prediction RMS", lab=["Standard" "Jacobian propagation"], background_color_legend=false, xscale=:log10,yscale=:log10, size=(1200,1100),legend=false)
plot!(xv,min.(sim,10), subplot=2, ylabel="Simulation RMS", lab=["Standard" "Jacobian propagation"], background_color_legend=false, xscale=:log10,yscale=:log10,legend=false)
plot!(xv,jac, subplot=3, ylabel="Jacobian Error", lab=["Standard" "Jacobian propagation"], background_color_legend=false, xscale=:log10,yscale=:log10,legend=false)

plot()
traceplot!.(resdiff, m=(:blue, 1), markerstrokealpha=0)
traceplot!.(resaddiff, m=(:red, 1), markerstrokealpha=0)
plot!(getfield.(resaddiff, :tracev), m=(:green, 1), markerstrokealpha=0)
gui()

##
using StatPlots
xticks = (1:2, ["Standard" "Jacprop"])
xvals = [1 2].*ones(num_montecarlo)
common = (marker_z=xv, legend=false, xticks=xticks)
vio1 = violin(xvals,pred; title="Prediction RMS",ylims=(0,maximum(pred)), common...)
# vio2 = violin(xvals,sim; title="Simulation RMS",  common...)
vio3 = violin(xvals,(jac); title="Jacobian Error",ylims=(0,maximum(jac)), common...)
plot(vio1,vio3, colorbar=false, layout=(1,2))
# savefig3("/local/home/fredrikb/phdthesis/blackbox/figs/boxplot_linear.tex")
##

# common = (marker_z=xv, legend=false, xticks=false)
# vio1 = violin([1],pred[:,1]; ylabel="Prediction RMS",  side=:left, ylims=(0,maximum(pred)), common...)
# violin!([1],pred[:,2]; ylabel="Prediction RMS",  side=:right, ylims=(0,maximum(pred)), common...)
# vio3 = violin([1],jac[:,1]; ylabel="Jacobian Error", side=:left, ylims=(0,maximum(jac)), common...)
# violin!([1],jac[:,2]; ylabel="Jacobian Error", side=:right, ylims=(0,maximum(jac)), common...)
# plot(vio1,vio3, colorbar=false, layout=(1,2))

## Visualize result
using Plots.PlotMeasures
F = font(18, "times")
fontopts = [(:xticks, 0:0.5:1), (:yticks, -0.3:0.3:0.3), (:xlims, (0,1.4)), (:ylims, (-0.5, 0.5)), (:grid, false), (:markeralpha, 0.2), (:ds, 2), (:markerstrokealpha, 0), (:titlefont, F), (:tickfont, F), (:xtickfont, F), (:ytickfont, F), (:guidefont, F), (:xguidefont, F), (:yguidefont, F)]

# mutregplot(trainer, vt, true_jacobian, title="Witout jacprop", subplot=1, layout=(2,2), reuse=false, useprior=false, showltv=false, legend=false, xaxis=(1:3), xaxis=(1:3))
# mutregplot!(trainerjn, vt, true_jacobian, title="With jacprop", subplot=2, link=:y, useprior=false, showltv=false, legend=false, xaxis=(1:3))
# traceplot!(trainer, subplot=3, title="Training error", xlabel="Epoch", legend=false)
# traceplot!(trainerjn, subplot=4, title="Training error", xlabel="Epoch", legend=false)

#' The top row shows the error (Frobenius norm) in the Jacbians for several points sampled randomly in the state space. The bottow row shows the training errors. The training errors are lower without jacprop, but he greater error in the Jacobians for the validation data indicates overfitting, which is prevented by jacprop.

eigvalplot(resdiff[5].models, vt, true_jacobian; title="Weight decay", layout=(1,2), subplot=1, size=(800,290), fontopts...)
# eigvalplot!(trainerjn.models, vt, true_jacobian;  title="Jacprop", subplot=2, fontopts...)
eigvalplot!(resaddiff[7].model, vt, true_jacobian;  title="Jacprop", subplot=2, fontopts...)
##
# gui()
# savefig("/local/home/fredrikb/phdthesis/blackbox/figs/jacpropeig2.pdf")
# pdftopng -r 300 -aa yes jacpropeig2.pdf jacpropeig2.png

error()
#' # Weight decay
#' ## Weight decay off
Random.seed!(1)
models     = [System(nx,nu,num_params, a) for a in default_activations]
opts       = ADAM.(params.(models), stepsize, decay=0.0005)
trainers  = ModelTrainer(;models=models, opts=opts, losses=JacProp.loss.(models), cb=callbacker, kalmanopts...)
for i = 1:3
    trainers(trajs[i], epochs=0, jacprop=0, useprior=false)
end
trainers(epochs=1000, jacprop=1)

Random.seed!(1)
models     = [DiffSystem(nx,nu,num_params, a) for a in default_activations]
opts       = ADAM.(params.(models), stepsize, decay=0.0005)
trainerds  = ModelTrainer(;models=models, opts=opts, losses=JacProp.loss.(models), cb=callbacker, kalmanopts...)
for i = 1:3
    trainerds(trajs[i], epochs=0, jacprop=0, useprior=false)
end
trainerds(epochs=1000, jacprop=1)

#' ## Weight decay on
wdecay = 0.1
Random.seed!(1)
models     = [System(nx,nu,num_params, a) for a in default_activations]
opts       = [[ADAM(params(models[i]), stepsize, decay=0.0005); [expdecay(Param(p), wdecay) for p in params(models[i]) if p isa AbstractMatrix]] for i = 1:length(models)]
trainerswd  = ModelTrainer(;models=models, opts=opts, losses=JacProp.loss.(models), cb=callbacker, kalmanopts...)
for i = 1:3
    trainerswd(trajs[i], epochs=0, jacprop=0, useprior=false)
end
trainerswd(epochs=1000, jacprop=1)

Random.seed!(1)
models     = [DiffSystem(nx,nu,num_params, a) for a in default_activations]
opts       = [[ADAM(params(models[i]), stepsize, decay=0.0005); [expdecay(Param(p), wdecay) for p in params(models[i]) if p isa AbstractMatrix]] for i = 1:length(models)]
trainerdswd  = ModelTrainer(;models=models, opts=opts, losses=JacProp.loss.(models), cb=callbacker, kalmanopts...)
for i = 1:3
    trainerdswd(trajs[i], epochs=0, jacprop=0, useprior=false)
end
trainerdswd(epochs=1000, jacprop=1)

#' Visualize result
#+ fig_ext="png"

eigvalplot(trainers.models, vt, true_jacobian; title="\$f, \\lambda=0\$", layout=4, subplot=1, size=(920,700), fontopts...)
eigvalplot!(trainerds.models, vt, true_jacobian;  title="\$g, \\lambda=0\$", subplot=2, fontopts...)
eigvalplot!(trainerswd.models, vt, true_jacobian;  title="\$f, \\lambda=$wdecay\$", subplot=3, fontopts...)
eigvalplot!(trainerdswd.models, vt, true_jacobian;  title="\$g, \\lambda=$wdecay\$", subplot=4, fontopts...)
eigvalplot(trainers.models, vt, true_jacobian;  title="", axis=false, fontopts...)
# plot!(link=:both)
gui()
# savefig("/local/home/fredrikb/papers/nn_prior/figs/all_zoom.pdf")
#' # Different activation functions

# ui = display_modeltrainer(trainerdswd, size=(800,600))

for acti in 1:4
    plot(layout=4, link=:both)
    actstr = string(JacProp.default_activations[acti])
    actstr = actstr[1:2] == "NN" ? actstr[7:end] : actstr
    for tri in 1:4
        tr = [trainers, trainerds, trainerswd, trainerdswd][tri]
        eigvalplot!([tr.models[acti]], vt, true_jacobian, title="$actstr, $(tri%2==0 ? "g" : "f"), wdecay=$(tri>2 ? wdecay : 0)", subplot=tri, markeralpha=0.05, ds=3, markerstrokealpha=0)
    end
    gui()
end


text:
using MathProgBase
using DataFrames

import Mimi.metainfo
import Mimi: CertainScalarParameter, CertainArrayParameter
import Base.*, Base.-, Base.+

export LinearProgrammingHall, LinearProgrammingShaft, LinearProgrammingRoom, LinearProgrammingHouse
export hallsingle, hall_relabel
export shaftsingle
export roomdiagonal, roomsingle, roomintersect, room_relabel, room_relabel_parameter
export varsum, fromindex
export setobjective!, setconstraint!, setconstraintoffset!, clearconstraint!, setlower!, setupper!
export gethouse, constraining, houseoptimize, summarizeparameters, findinfeasiblepair, varlengths

# A hallway is a vector of variables
type LinearProgrammingHall
    component::Symbol
    name::Symbol
    f::Vector{Float64}
end

function hallsingle(model::Model, component::Symbol, name::Symbol, gen::Function)
    LinearProgrammingHall(component, name, vectorsingle(getdims(model, component, name), gen))
end

"""Connect a derivative to another component: change the variable component and name to another component."""
function hall_relabel(hall::LinearProgrammingHall, from::Symbol, tocomponent::Symbol, toname::Symbol)
    @assert hall.name == from "Name mismatch in hall_relabel: $(hall.name) <> $from"

    LinearProgrammingHall(tocomponent, toname, hall.f)
end

function -(hall::LinearProgrammingHall)
    LinearProgrammingHall(hall.component, hall.name, -hall.f)
end

function +(hall1::LinearProgrammingHall, hall2::LinearProgrammingHall; skipnamecheck=false)
    if !skipnamecheck
        @assert hall1.name == hall2.name "Hall + Hall name mismatch: $(hall1.name) <> $(hall2.name); use hall_relabel?"
    end
    LinearProgrammingHall(hall1.component, hall1.name, hall1.f + hall2.f)
end

function -(hall1::LinearProgrammingHall, hall2::LinearProgrammingHall; skipnamecheck=false)
    if !skipnamecheck
        @assert hall1.name == hall2.name "Hall - Hall name mismatch: $(hall1.name) <> $(hall2.name); use hall_relabel?"
    end
    LinearProgrammingHall(hall1.component, hall1.name, hall1.f - hall2.f)
end

# A shaft is a transpose of a hallway, used for parameters
type LinearProgrammingShaft
    component::Symbol
    name::Symbol
    x::Vector{Float64}
end

function shaftsingle(model::Model, component::Symbol, name::Symbol, gen::Function)
    LinearProgrammingShaft(component, name, vectorsingle(getdims(model, component, name), gen))
end

function -(shaft::LinearProgrammingShaft)
    LinearProgrammingShaft(shaft.component, shaft.name, -shaft.x)
end

function +(shaft1::LinearProgrammingShaft, shaft2::LinearProgrammingShaft; skipnamecheck=false)
    if !skipnamecheck
        @assert shaft1.name == shaft2.name "Shaft + Shaft name mismatch: $(shaft1.name) <> $(shaft2.name); use shaft_relabel?"
    end
    LinearProgrammingShaft(shaft1.component, shaft1.name, shaft1.x + shaft2.x)
end

type LinearProgrammingRoom
    varcomponent::Symbol
    variable::Symbol
    paramcomponent::Symbol
    parameter::Symbol
    A::SparseMatrixCSC{Float64, Int64}
end

"""
Fill in just the diagonal
"""
function roomdiagonal(model::Model, component::Symbol, variable::Symbol, parameter::Symbol, gen::Function)
    dimsvar = getdims(model, component, variable)
    dimspar = getdims(model, component, parameter)
    @assert dimsvar == dimspar "Variable and parameter in roomdiagonal do not have the same dimensions: $dimsvar <> $dimspar"
    LinearProgrammingRoom(component, variable, component, parameter, matrixdiagonal(dimsvar, gen))
end

"""
Fill in every element
"""
function roomsingle(model::Model, component::Symbol, variable::Symbol, parameter::Symbol, gen::Function)
    dimsvar = getdims(model, component, variable)
    dimspar = getdims(model, component, parameter)
    LinearProgrammingRoom(component, variable, component, parameter, matrixsingle(dimsvar, dimspar, gen))
end

"""
Call gen for each shared index, passing an array to be filled for unshared indexes.
This version assumes both variables come from the same component.
"""
function roomintersect(model::Model, component::Symbol, variable1::Symbol, variable2::Symbol, gen::Function)
    dims1 = getdims(model, component, variable1)
    dims2 = getdims(model, component, variable2)
    dimnames1 = getdimnames(model, component, variable1)
    dimnames2 = getdimnames(model, component, variable2)
    LinearProgrammingRoom(component, variable1, component, variable2, matrixintersect(dims1, dims2, dimnames1, dimnames2, gen))
end

"""
Call gen for each shared index, passing an array to be filled for unshared indexes.
This version allows the variables to come from different component.
"""
function roomintersect(model::Model, component1::Symbol, variable1::Symbol, component2::Symbol, variable2::Symbol, gen::Function)
    dims1 = getdims(model, component1, variable1)
    dims2 = getdims(model, component2, variable2)
    dimnames1 = getdimnames(model, component1, variable1)
    dimnames2 = getdimnames(model, component2, variable2)
    LinearProgrammingRoom(component1, variable1, component2, variable2, matrixintersect(dims1, dims2, dimnames1, dimnames2, gen))
end

"""Connect a gradient to another component: change the variable component and name to another component."""
function room_relabel(room::LinearProgrammingRoom, from::Symbol, tocomponent::Symbol, toname::Symbol)
    @assert room.variable == from "Variable name mismatch in room_relabel: $(room.variable) <> $from"

    LinearProgrammingRoom(tocomponent, toname, room.paramcomponent, room.parameter, room.A)
end

"""Connect a gradient to another component: change the variable component and name to another component."""
function room_relabel_parameter(room::LinearProgrammingRoom, from::Symbol, tocomponent::Symbol, toname::Symbol)
    @assert room.parameter == from "Parameter name mismatch in room_relabel_parameter: $(room.variable) <> $from"

    LinearProgrammingRoom(room.varcomponent, room.variable, tocomponent, toname, room.A)
end

function *(hall::LinearProgrammingHall, room::LinearProgrammingRoom; skipnamecheck=false)
    if !skipnamecheck
        @assert hall.name == room.variable "Hall * Room name mismatch: $(hall.name) <> $(room.variable); use room_relabel?"
    end
    LinearProgrammingHall(room.paramcomponent, room.parameter, vec(transpose(hall.f) * room.A))
end

# Returns a hall, since it then contains variable/constraint information
function *(room::LinearProgrammingRoom, shaft::LinearProgrammingShaft; skipnamecheck=false)
    if !skipnamecheck
        @assert shaft.name == room.parameter "Room * Shaft name mismatch: $(room.parameter) <> $(hall.name); use room_relabel?"
    end
    LinearProgrammingHall(room.varcomponent, room.variable, room.A * shaft.x)
end

function *(room1::LinearProgrammingRoom, room2::LinearProgrammingRoom; skipnamecheck=false)
    if !skipnamecheck
        @assert room1.parameter == room2.variable "Room * Room name mismatch: $(room1.parameter) <> $(room2.variable); use room_relabel?"
    end
    LinearProgrammingRoom(room1.varcomponent, room1.variable, room2.paramcomponent, room2.parameter, room1.A * room2.A)
end

function +(room1::LinearProgrammingRoom, room2::LinearProgrammingRoom; skipnamecheck=false)
    if !skipnamecheck
        @assert room1.parameter == room2.parameter "Room + Room parameter name mismatch: $(room1.parameter) <> $(room2.parameter); use room_relabel?"
        @assert room1.variable == room2.variable "Room + Room variable name mismatch: $(room1.variable) <> $(room2.variable); use room_relabel?"
    end
    LinearProgrammingRoom(room1.varcomponent, room1.variable, room1.paramcomponent, room1.parameter, room1.A + room2.A)
end

function -(room::LinearProgrammingRoom)
    LinearProgrammingRoom(room.varcomponent, room.variable, room.paramcomponent, room.parameter, -room.A)
end

function varsum(room::LinearProgrammingRoom)
    LinearProgrammingHall(room.paramcomponent, room.parameter, vec(sum(room.A, 1)))
end

type LinearProgrammingHouse
    model::Model
    paramcomps::Vector{Symbol}
    parameters::Vector{Symbol}
    constcomps::Vector{Symbol}
    constraints::Vector{Symbol}
    constdictionary::Dict{Symbol, Symbol}
    lowers::Vector{Float64}
    uppers::Vector{Float64}
    f::Vector{Float64}
    A::SparseMatrixCSC{Float64, Int64}
    b::Vector{Float64}
end

function LinearProgrammingHouse(model::Model, paramcomps::Vector{Symbol}, parameters::Vector{Symbol}, constcomps::Vector{Symbol}, constraints::Vector{Symbol}, constdictionary::Dict{Symbol, Symbol}=Dict{Symbol, Symbol}())
    paramlen = sum(varlengths(model, paramcomps, parameters))
    variablelen = sum(varlengths(model, constcomps, constraints, constdictionary))
    A = spzeros(variablelen, paramlen)
    f = zeros(paramlen)
    b = zeros(variablelen)
    LinearProgrammingHouse(model, paramcomps, parameters, constcomps, constraints, constdictionary, zeros(length(f)), Inf * ones(length(f)), f, A, b)
end

function setobjective!(house::LinearProgrammingHouse, hall::LinearProgrammingHall)
    @assert hall.name in house.parameters "$(hall.name) not a known parameter"
    kk = findfirst((house.paramcomps .== hall.component) & (house.parameters .== hall.name))
    paramspans = varlengths(house.model, house.paramcomps, house.parameters)
    @assert length(hall.f) == paramspans[kk] "Length of parameter $(hall.name) unexpected: $(length(hall.f)) <> $(paramspans[kk])"
    before = sum(paramspans[1:kk-1])
    for ii in 1:paramspans[kk]
        #@assert house.f[before+ii] == 0 "Overwrite existing gradient in setobjective!"
        house.f[before+ii] = hall.f[ii]
    end
end

function setconstraint!(house::LinearProgrammingHouse, room::LinearProgrammingRoom)
    @assert room.variable in house.constraints "$(room.variable) not a known variable"
    @assert room.parameter in house.parameters "$(room.parameter) not a known parameter"

    kk = findfirst((house.paramcomps .== room.paramcomponent) & (house.parameters .== room.parameter))
    ll = findfirst((house.constcomps .== room.varcomponent) & (house.constraints .== room.variable))
    paramspans = varlengths(house.model, house.paramcomps, house.parameters)
    constspans = varlengths(house.model, house.constcomps, house.constraints, house.constdictionary)
    @assert size(room.A, 1) == constspans[ll] "Length of variable $(room.variable) unexpected: $(size(room.A, 1)) <> $(constspans[ll])"
    @assert size(room.A, 2) == paramspans[kk] "Length of parameter $(room.parameter) unexpected: $(size(room.A, 2)) <> $(paramspans[kk])"

    parambefore = sum(paramspans[1:kk-1])
    constbefore = sum(constspans[1:ll-1])

    house.A[constbefore+1:constbefore+constspans[ll], parambefore+1:parambefore+paramspans[kk]] = room.A
end

function setconstraintoffset!(house::LinearProgrammingHouse, hall::LinearProgrammingHall)
    @assert hall.name in house.constraints "$(hall.name) not a known variable"
    kk = findfirst((house.constcomps .== hall.component) & (house.constraints .== hall.name))
    constspans = varlengths(house.model, house.constcomps, house.constraints, house.constdictionary)
    @assert length(hall.f) == constspans[kk] "Length of parameter $(hall.name) unexpected: $(length(hall.f)) <> $(constspans[kk])"
    before = sum(constspans[1:kk-1])
    for ii in 1:constspans[kk]
        #@assert house.b[before+ii] == 0 "Overwrite existing gradient in setobjective!"
        house.b[before+ii] = hall.f[ii]
    end
end

function clearconstraint!(house::LinearProgrammingHouse, component::Symbol, variable::Symbol)
    @assert variable in house.constraints "$(variable) not a known variable"

    ll = findfirst((house.constcomps .== component) & (house.constraints .== variable))
    constspans = varlengths(house.model, house.constcomps, house.constraints, house.constdictionary)

    constbefore = sum(constspans[1:ll-1])

    house.A[constbefore+1:constbefore+constspans[ll], :] = 0
    house.b[constbefore+1:constbefore+constspans[ll]] = 0
end

function setlower!(house::LinearProgrammingHouse, hall::LinearProgrammingHall)
    @assert hall.name in house.parameters "$(hall.name) not a known parameter"
    kk = findfirst((house.paramcomps .== hall.component) & (house.parameters .== hall.name))
    paramspans = varlengths(house.model, house.paramcomps, house.parameters)
    @assert length(hall.f) == paramspans[kk] "Length of parameter $(hall.name) unexpected: $(length(hall.f)) <> $(paramspans[kk])"
    before = sum(paramspans[1:kk-1])
    for ii in 1:paramspans[kk]
        #@assert house.b[before+ii] == 0 "Overwrite existing gradient in setobjective!"
        house.lowers[before+ii] = hall.f[ii]
    end
end

function setupper!(house::LinearProgrammingHouse, hall::LinearProgrammingHall)
    @assert hall.name in house.parameters "$(hall.name) not a known parameter"
    kk = findfirst((house.paramcomps .== hall.component) & (house.parameters .== hall.name))
    paramspans = varlengths(house.model, house.paramcomps, house.parameters)
    @assert length(hall.f) == paramspans[kk] "Length of parameter $(hall.name) unexpected: $(length(hall.f)) <> $(paramspans[kk])"
    before = sum(paramspans[1:kk-1])
    for ii in 1:paramspans[kk]
        #@assert house.b[before+ii] == 0 "Overwrite existing gradient in setobjective!"
        house.uppers[before+ii] = hall.f[ii]
    end
end

function gethouse(house::LinearProgrammingHouse, rr::Int64, cc::Int64)
    # Determine the row and column names
    varii = findlast(cumsum(varlengths(house.model, house.constcomps, house.constraints, house.constdictionary)) .< rr) + 1
    parii = findlast(cumsum(varlengths(house.model, house.paramcomps, house.parameters)) .< cc) + 1
    rrrelative = rr - sum(varlengths(house.model, house.constcomps, house.constraints, house.constdictionary)[1:varii-1])
    ccrelative = cc - sum(varlengths(house.model, house.paramcomps, house.parameters)[1:parii-1])
    vardims = getdims(house.model, house.constcomps[varii], get(house.constdictionary, house.constraints[varii], house.constraints[varii]))
    pardims = getdims(house.model, house.paramcomps[parii], house.parameters[parii])

    println("$(house.constcomps[varii]).$(house.constraints[varii])$(toindex(rrrelative, vardims)), $(house.paramcomps[parii]).$(house.parameters[parii])$(toindex(ccrelative, pardims)) = $(house.A[rr, cc])")
end

function constraining(house::LinearProgrammingHouse, solution::Vector{Float64})
    # Determine which constraint (if any) is stopping an increase or decrease of each
    df = DataFrame(solution=solution)
    df[:component] = :na
    df[:parameter] = :na
    df[:abovefail] = ""
    df[:belowfail] = ""

    # Produce names for all constraints
    varlens = varlengths(house.model, house.constcomps, house.constraints, house.constdictionary)
    names = ["" for ii in 1:sum(varlens)]
    for kk in 1:length(house.constcomps)
        ii0 = sum(varlens[1:kk-1])
        for ii in 1:varlens[kk]
            names[ii0 + ii] = "$(house.constraints[kk]).$ii"
        end
    end

    varlens = varlengths(house.model, house.paramcomps, house.parameters)
    baseconsts = house.A * solution

    println("Ignore:")
    println(join(names[find(baseconsts .> house.b)], ", "))
    ignore = baseconsts .> house.b

    for kk in 1:length(house.paramcomps)
        ii0 = sum(varlens[1:kk-1])
        for ii in 1:varlens[kk]
            df[ii0 + ii, :component] = house.paramcomps[kk]
            df[ii0 + ii, :parameter] = house.parameters[kk]

            newconst = baseconsts + house.A[:, ii0 + ii] * 1e-6
            df[ii0 + ii, :abovefail] = join(names[find((newconst .> house.b) & !ignore)], ", ")

            newconst = baseconsts - house.A[:, ii0 + ii] * 1e-6
            df[ii0 + ii, :belowfail] = join(names[find((newconst .> house.b) & !ignore)], ", ")
        end
    end

    df
end

houseoptimize(house::LinearProgrammingHouse) = linprog(-house.f, house.A, '<', house.b, house.lowers, house.uppers)
houseoptimize(house::LinearProgrammingHouse, solver) = linprog(-house.f, house.A, '<', house.b, house.lowers, house.uppers, solver)

function findinfeasiblepair(house::LinearProgrammingHouse, solver)
    top = findinfeasiblepairhelper(house, solver, "top", 2, length(house.b))
    bottom = findinfeasiblepairhelper(house, solver, "bottom", 1, length(house.b) - 1)

    [top, bottom]
end

"""
Look for the row that makes the problem feasible
"""
function findinfeasiblepairhelper(house::LinearProgrammingHouse, solver, direction, searchtop, searchbottom)
    searchpoint = round(Int, (searchtop + searchbottom) / 2 + .01) # Always round up
    println([direction, searchtop, searchbottom])
    isfeasible = checkfeasibility(house, solver, direction, searchpoint)

    if !isfeasible
        # Look for a shorter span
        if searchpoint == searchbottom
            if direction == "bottom" || checkfeasibility(house, solver, direction, searchtop)
                searchbottom
            else
                searchtop
            end
        elseif direction == "top"
            findinfeasiblepairhelper(house, solver, direction, searchtop, searchpoint)
        else
            findinfeasiblepairhelper(house, solver, direction, searchpoint, searchbottom)
        end
    else
        # Look for a longer span
        if searchpoint == searchbottom
            searchtop
        elseif (direction == "top")
            findinfeasiblepairhelper(house, solver, direction, searchpoint, searchbottom)
        else
            findinfeasiblepairhelper(house, solver, direction, searchtop, searchpoint)
        end
    end
end

function checkfeasibility(house::LinearProgrammingHouse, solver, direction, searchpoint)
    if (direction == "top")
        sol = linprog(-house.f, house.A[1:searchpoint, :], '<', house.b[1:searchpoint], house.lowers, house.uppers, solver)
    else
        sol = linprog(-house.f, house.A[searchpoint:end, :], '<', house.b[searchpoint:end], house.lowers, house.uppers, solver)
    end

    sol.status != :Infeasible
end

function summarizeparameters(house::LinearProgrammingHouse, solution::Vector{Float64})
    # Look at parameter values
    varlens = varlengths(house.model, house.paramcomps, house.parameters)
    for ii in 1:length(house.parameters)
        println(house.parameters[ii])
        index1 = sum(varlens[1:ii-1]) + 1
        index2 = sum(varlens[1:ii])

        values = solution[index1:index2]

        if (sum(values .!= 0) == 0)
            println("All zero.")
        else
            println(values[1:min(100, index2 - index1 + 1)])
            println("Sum: $(sum(values))")
        end
    end
end

function getconstraintsolution(house, sol, constraint)
    constvalues = house.A * sol.sol

    varlens = varlengths(m, house.constcomps, house.constraints, house.constdictionary)

    ii = find(house.constraints .== constraint)[1]

    index1 = sum(varlens[1:ii-1]) + 1
    index2 = sum(varlens[1:ii])

    constvalues[index1:index2]
end

## Helpers

function rangeof(m::Model, name, components, names, vardictionary::Dict{Symbol, Symbol}=Dict{Symbol, Symbol}())
    varlens = varlengths(m, components, names, vardictionary)
    kk = findfirst(name .== names)
    sum(varlens[1:kk-1])+1:sum(varlens[1:kk])
end

"Translate an offset value (+1) to an index vector."
function toindex(ii::Int64, dims::Vector{Int64})
    indexes = Vector{Int64}(length(dims))
    offset = ii - 1
    for dd in 1:length(dims)
        indexes[dd] = offset % dims[dd] + 1
        offset = floor(Int64, offset / dims[dd])
    end

    return indexes
end

"Translate an index vector to an offset (+1)."
function fromindex(index::Vector{Int64}, dims::Vector{Int64})
    offset = index[end]
    for ii in length(dims)-1:-1:1
        offset = (offset - 1) * dims[ii] + index[ii]
    end

    return offset
end

"Return a vector of the indices defining the parameter or variable."
function getdims(model::Model, component::Symbol, name::Symbol)
    if name in keys(model.parameters)
        if isa(model.parameters[name], CertainScalarParameter)
            Int64[1]
        elseif isa(model.parameters[name], CertainArrayParameter)
            Int64[size(model.parameters[name].values)...]
        end
    else
        meta = metainfo.getallcomps()
        convert(Vector{Int64}, map(dim -> model.indices_counts[dim], meta[(:Main, component)].variables[name].dimensions))
    end
end

"""
Return the symbols representing each of the dimensions for this variable or parameter.
"""
function getdimnames(model::Model, component::Symbol, name::Symbol)
    meta = metainfo.getallcomps()
    if name in keys(meta[(:Main, component)].parameters)
        convert(Vector{Symbol}, meta[(:Main, component)].parameters[name].dimensions)
    else
        convert(Vector{Symbol}, meta[(:Main, component)].variables[name].dimensions)
    end
end

"Return the total span occupied by each variable or parameter."
function varlengths(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, vardictionary::Dict{Symbol, Symbol}=Dict{Symbol, Symbol}())
    Int64[prod(getdims(model, components[ii], get(vardictionary, names[ii], names[ii]))) for ii in 1:length(components)]
end

## Matrix methods

"Construct a matrix with the given dimensions, calling gen for each element."
function vectorsingle(dims::Vector{Int64}, gen)
    f = Vector{Float64}(prod(dims))
    for ii in 1:length(f)
        f[ii] = gen(toindex(ii, dims)...)
    end

    f
end

"Call the generate function for all indices along the diagonal."
function matrixdiagonal(dims::Vector{Int64}, gen)
    dimlen = prod(dims)
    A = spzeros(dimlen, dimlen)
    for ii in 1:dimlen
        A[ii, ii] = gen(toindex(ii, dims)...)
    end

    A
end

"Call the generate function for every element."
function matrixsingle(vardims::Vector{Int64}, pardims::Vector{Int64}, gen)
    vardimlen = prod(vardims)
    pardimlen = prod(pardims)
    A = spzeros(vardimlen, pardimlen)
    for ii in 1:vardimlen
        for jj in 1:pardimlen
            A[ii, jj] = gen(toindex(ii, vardims)..., toindex(jj, pardims)...)
        end
    end

    A
end

"""
Call the generate function with all combinations of the shared
indices; shared dimensions must come in the same order and at the end
of the dimensions lists.
"""
function matrixintersect(rowdims::Vector{Int64}, coldims::Vector{Int64}, rowdimnames::Vector{Symbol}, coldimnames::Vector{Symbol}, gen::Function)
    A = spzeros(prod(rowdims), prod(coldims))

    # Determine shared dimensions: counting from 0 = end
    numshared = 0
    while numshared < min(length(rowdims), length(coldims))
        if rowdims[end - numshared] != coldims[end - numshared]
            break
        end
        if rowdimnames[end - numshared] != coldimnames[end - numshared]
            break
        end
        numshared += 1
    end

    sharedims = rowdims[end - numshared + 1:end]
    sharedimslen = prod(sharedims)
    for ii in 1:sharedimslen
        index = toindex(ii, sharedims)
        topii = fromindex([ones(Int64, length(rowdims) - numshared); index], rowdims)
        bottomii = fromindex([rowdims[1:length(rowdims) - numshared]; index], rowdims)
        leftii = fromindex([ones(Int64, length(coldims) - numshared); index], coldims)
        rightii = fromindex([coldims[1:length(coldims) - numshared]; index], coldims)
        subA = sub(A, topii:bottomii, leftii:rightii)
        gen(subA, index...)
    end

    A
end

ext: jl


text:
"""
    union(X::T, Y::T) where {T<:BipartiteNetwork}

Union of two bipartite networks -- interactions *and* species which are present
in either networks are also present in the final network.
"""
function Base.union(X::T, Y::T) where {T<:BipartiteNetwork}
    new_t = union(species(X; dims=1), species(Y; dims=1))
    new_b = union(species(X; dims=2), species(Y; dims=2))
    int = interactions(X) ∪ interactions(Y)
    I = zeros(Int64, length(int))
    J = zeros(Int64, length(int))
    for i in eachindex(int)
        I[i] = findfirst(isequal(int[i].from), new_t)
        J[i] = findfirst(isequal(int[i].to), new_b)
    end
    new_a = sparse(I, J, true, length(new_t), length(new_b))
    return BipartiteNetwork(new_a, new_t, new_b)
end

"""
    intersect(X::T, Y::T) where {T<:BipartiteNetwork}

Intersect between two bipartite networks. The resulting network has the *species
and interactions* common to both networks. This can result in species being
disconnected, if they are found in both networks but with no operations in
common.
"""
function Base.intersect(X::T, Y::T) where {T<:BipartiteNetwork}
    new_t = intersect(species(X; dims=1), species(Y; dims=1))
    new_b = intersect(species(X; dims=2), species(Y; dims=2))
    int = interactions(X) ∩ interactions(Y)
    I = zeros(Int64, length(int))
    J = zeros(Int64, length(int))
    for i in eachindex(int)
        I[i] = findfirst(isequal(int[i].from), new_t)
        J[i] = findfirst(isequal(int[i].to), new_b)
    end
    new_a = sparse(I, J, true, length(new_t), length(new_b))
    return BipartiteNetwork(new_a, new_t, new_b)
end

"""
    union(X::T, Y::T) where {T<:UnipartiteNetwork}

Union of two unipartite networks -- interactions *and* species which are present
in either networks are also present in the final network.
"""
function Base.union(X::T, Y::T) where {T<:UnipartiteNetwork}
    new_s = union(species(X), species(Y))
    int = interactions(X) ∪ interactions(Y)
    I = zeros(Int64, length(int))
    J = zeros(Int64, length(int))
    for i in eachindex(int)
        I[i] = findfirst(isequal(int[i].from), new_s)
        J[i] = findfirst(isequal(int[i].to), new_s)
    end
    new_a = sparse(I, J, true, length(new_s), length(new_s))
    return UnipartiteNetwork(new_a, new_s)
end

"""
    intersect(X::T, Y::T) where {T<:UnipartiteNetwork}

Intersect between two unipartite networks. The resulting network has the
*species and interactions* common to both networks. This can result in species
being disconnected, if they are found in both networks but with no operations in
common.
"""
function Base.intersect(X::T, Y::T) where {T<:UnipartiteNetwork}
    new_s = intersect(species(X), species(Y))
    int = interactions(X) ∩ interactions(Y)
    I = zeros(Int64, length(int))
    J = zeros(Int64, length(int))
    for i in eachindex(int)
        I[i] = findfirst(isequal(int[i].from), new_s)
        J[i] = findfirst(isequal(int[i].to), new_s)
    end
    new_a = sparse(I, J, true, length(new_s), length(new_s))
    return UnipartiteNetwork(new_a, new_s)
end

"""
    setdiff(X::T, Y::T) where {T<:BipartiteNetwork}

Difference between two bipartite networks. This operation is *sensitive* to the
order of arguments, as the resulting network will have the species present in
the first network (and their interactions) only.
"""
function Base.setdiff(X::T, Y::T) where {T<:BipartiteNetwork}
    new_t = setdiff(species(X; dims=1), species(Y; dims=1))
    new_b = setdiff(species(X; dims=2), species(Y; dims=2))
    return X[new_t, new_b]
end

"""
    setdiff(X::T, Y::T) where {T<:UnipartiteNetwork}

Difference between two unipartite networks. This operation is *sensitive* to the
order of arguments, as the resulting network will have the species present in
the first network (and their interactions) only.
"""
function Base.setdiff(X::T, Y::T) where {T<:UnipartiteNetwork}
    new_s = setdiff(species(X; dims=1), species(Y; dims=1))
    return X[new_s]
end


text:
function compute()
N=4000000
ans=0
x=1
y=2
while x<=N
	if(x%2==0)
		 ans+=x
	end
	 x,y=y,x+y

end
ans
end
print(compute())


text:
@testset "Matrix" begin

    @testset "Matrix: Determinant" begin
        M1 = [1 0; 2 2]
        M2 = rand(3,3)
        M3 = rand(4,4)

        @test determinant(M1) == det(M1)
        @test round(determinant(M2),digits = 4) == round(det(M2),digits = 4)
        @test round(determinant(M3),digits = 4) == round(det(M3),digits = 4)
    end

    @testset "Matrix: LU Decompose" begin
        mat = [
            2  -1 -2;
            -4   6  3;
            -4  -2  8
        ]

        L = [
            1  0 0;
            -2  1 0;
            -2 -1 1
        ]
        U = [
            2 -1 -2;
            0  4 -1;
            0  0  3
        ]

	    @test lu_decompose(mat) == (L,U)
    end

	@testset "Matrix: rotation matrix" begin
		theta = pi/6
		a = [1;0]
		b = [0;1]
		@test rotation_matrix(theta)*a == [cos(theta);sin(theta)]
		@test rotation_matrix(theta)*b == [-sin(theta);cos(theta)]
	end
end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia


text:
export savitsky_golay, central_difference

#Polynomial smoothing with the Savitsky Golay filters
#
# Sources
# ---------
# Theory: http://www.ece.rutgers.edu/~orfanidi/intro2sp/orfanidis-i2sp.pdf
# Python Example: http://wiki.scipy.org/Cookbook/SavitzkyGolay
# Modified from https://github.com/BBN-Q/Qlab.jl/blob/master/src/SavitskyGolay.jl

#This can't handle NaNs or skipped data
function savitsky_golay(x::Vector, windowSize::Integer, polyOrder::Integer; deriv::Integer=0)

#Some error checking
    @assert isodd(windowSize) "Window size must be an odd integer."
    @assert polyOrder < windowSize "Polynomial order must me less than window size."

    halfWindow = round(Int64,(windowSize-1)/2)

    #Setup the S matrix of basis vectors.
    S = zeros(windowSize, polyOrder+1)
    for ct = 0:polyOrder
        S[:,ct+1] = collect(-halfWindow:halfWindow).^(ct)
    end

    #Compute the filter coefficients for all orders
    #From the scipy code it seems pinv(S) and taking rows should be enough
    G = S*pinv(S'*S)

    #Slice out the derivative order we want
    filterCoeffs = G[:,deriv+1] * factorial(deriv);

    #Pad the signal with the endpoints and convolve with filter
    paddedX = [x[1]*ones(halfWindow); x; x[end]*ones(halfWindow)]
    y = conv(filterCoeffs[end:-1:1], paddedX)

    #Return the valid midsection
    return y[2*halfWindow+1:end-2*halfWindow]

end

#=
Outlier removal
=#

function outlier_removal_min(cov1,min_ex,ind_range=(1,length(cov1)))

    cov2=cov1[ind_range[1]:ind_range[2]]

    outlier_inds=cov2 .< percentile(cov2,min_ex)

    itp_out = interpolate((find(.!outlier_inds),), cov2[.!outlier_inds], Gridded(Linear()))

    for i in find(outlier_inds)[2:(end-1)]
        cov2[i] = itp_out[i]
    end

    cov1[ind_range[1]:ind_range[2]]=cov2
    nothing
end

function outlier_removal_max(cov1,max_ex,ind_range=(1,length(cov1)))

    cov2=cov1[ind_range[1]:ind_range[2]]

    outlier_inds=cov2 .> percentile(cov2,max_ex)

    itp_out = interpolate((find(.!outlier_inds),), cov2[.!outlier_inds], Gridded(Linear()))

    for i in find(outlier_inds)[2:(end-1)]
        cov2[i] = itp_out[i]
    end

    cov1[ind_range[1]:ind_range[2]]=cov2
    nothing
end
function outlier_removal_twosided(cov1,min_ex,max_ex,ind_range=(1,length(cov1)))

    cov2=cov1[ind_range[1]:ind_range[2]]

    outlier_inds1=cov2 .< percentile(cov2,min_ex)
    outlier_inds2=cov2 .> percentile(cov2,max_ex)

    outlier_inds = outlier_inds1 .| outlier_inds2

    itp_out = interpolate((find(.!outlier_inds),), cov2[.!outlier_inds], Gridded(Linear()))

    for i in find(outlier_inds)[2:(end-1)]
        cov2[i] = itp_out[i]
    end

    cov1[ind_range[1]:ind_range[2]]=cov2
    nothing
end

#Central Differences
#=
One implementation only calculates if points are tracked on either side
=#

function central_difference(x::Array{T,1}) where T
    y=zeros(Float64,length(x))

    for i=2:length(y)-1
        y[i]=(x[i+1] - x[i-1]) / 2
    end
    y
end

function central_difference(x::Array{T,1},tracked) where T
    y=zeros(Float64,length(x))
    diff_tracked=deepcopy(tracked)

    for i=2:length(y)-1
        if (tracked[i-1]&tracked[i+1])
            y[i]=(x[i+1] - x[i-1]) / 2
        else
            diff_tracked[i] = false
        end
    end

    tracked=diff_tracked

    y
end

function count_consecutive_bits(x,con_thres)
    interp_inds = Array{Int64,1}()

    i = 2
    while (i < (length(x)-con_thres))
        if !(x[i]) #false detected
            j = 1
            while (j<con_thres)
                if (x[i+j])

                    for k=1:j
                        push!(interp_inds,i+k)
                    end
                    break
                end
                j = j + 1
            end
            i = i+j
        else
            i = i+1
        end
    end

    interp_inds
end

function remove_two_sided_outlier(x::Array{T,1},t,min_p,max_p) where T

    min_x = percentile(x[t],min_p)
    max_x = percentile(x[t],max_p)

    for i=1:length(t)

        if t[i]
            if (x[i]<min_x)|(x[i]>max_x)
                t[i] = false
            end
        end
    end
    nothing
end

function interpolate_forward(x::Array{T,1},t,interp_ids::Array{Int64,1}) where T

    interp_x = interpolate((findall(t),),x[t],Gridded(Linear()))

    for i in interp_ids
        try
            x[i] = interp_x(i)
        catch
        end
    end
end

function reset_tracked(t,interp_ids::Array{Int64,1})
   t[interp_ids] .= true
end


text:
function reverse_direction(d::Integer)
    if d == NORTH
        return SOUTH
    elseif d == SOUTH
        return NORTH
    elseif d == WEST
        return EAST
    elseif d == EAST
        return WEST
    else
        error("Reverse direction for $d not defined")
    end
end

function neighbor_cell(i::Integer, j::Integer, out_direction::Integer)
    if out_direction == NORTH
        return (i - 1, j)
    elseif out_direction == EAST
        return (i, j + 1)
    elseif out_direction == SOUTH
        return (i + 1, j)
    elseif out_direction == WEST
        return (i, j - 1)
    else
        error("Out-direction $out_direction is not valid")
    end
end

function vertices_on_cell(g::FlatlandGraph, i::Integer, j::Integer)
    return [
        get_vertex(g, (i, j, direction, REAL)) for
        direction in CARDINAL_POINTS if haskey(g, (i, j, direction, REAL))
    ]
end

function mirror_vertex(g::FlatlandGraph, v::Integer)
    (i, j, direction, kind) = get_label(g, v)
    if kind == REAL
        i2, j2 = neighbor_cell(i, j, reverse_direction(direction))
        return get_vertex(g, (i2, j2, reverse_direction(direction), REAL))
    end
end

function transition_exists(
    transition_map::String, in_direction::Integer, out_direction::Integer
)
    return transition_map[4(in_direction - 1) + out_direction] == '1'
end

function direction_exists(transition_map::String, in_direction::Integer)
    return any(
        transition_exists(transition_map, in_direction, out_direction) for
        out_direction in CARDINAL_POINTS
    )
end


text:
# Defines the UncertainBound type to handle uncertain untegration start and end points

# Define which quantiles are stored in UncertainBound objects
const BOUND_QUANTILES = 0.1:0.2:0.9 |> collect
const IDX_BOUND_MEDIAN = 3 # the median is at index 3

"""
    UncertainBound{T, N}

Used to define uncertain integration bounds.

Holds `N` samples of start and end points of type `T` for numeric integration.

# Fields

- `left::Particles{T, N}`: samples of starting point
- `right::Particles{T, N}`: samples of end point

# Constructors

    UncertainBound(left::S, right::T, N::Int=10_000) where {S <: ContinuousUnivariateDistribution, T <: ContinuousUnivariateDistribution}

Create an `UncertainBound` from two distributions `left` and `right` defining uncertain start and end points of the integration window.

**Example**

Start point falls in the range [1, 2] with uniform probability, end point falls in the range [5, 6] with uniform probability:

```jldoctest UBexample
julia> using Distributions

julia> using Random: seed!; seed!(1);

julia> UncertainBound(Uniform(1, 2), Uniform(5, 6), 20_000)
UncertainBound{Float64, 20000}(start = 1.5 ± 0.29, end = 5.5 ± 0.29)
```

---

    UncertainBound(pos::T, width::ContinuousUnivariateDistribution, uc::UncertainCurve{T, N}) where {T, N}

Create an `UncertainBound` from a position `pos`, a distribution for the integration window width `width`, and an uncertain curve object `uc`.
Useful for integrating symmetric peaks.
The constructor will find start and end points for the integration window such that, for each sample in `uc`, the integration window will
be symmetric around the peak that falls in the interval `pos` ± (sample of `width`)/2, i.e. the integration window will "follow" the peak
position from draw to draw (note that, due to noise, the peak may change its position from draw to draw).

**Example**

Uncertainty of width is modeled with a scaled and shifted Beta(2, 2) distribution:

```jldoctest UBexample
julia> uc = begin # create uncertain curve with one symmetric peak
           x = 0:0.1:10;
           y = @. exp(-(x - 5)^2)
           add_noise(Curve(x, y), MvGaussianNoiseModel(0.1, 0.03, 0.5))
       end;

julia> ub = UncertainBound(5., scale_shift_beta(2, 2, 3.5, 4.0), uc)
UncertainBound{Float64, 10000}(start = 3.12 ± 0.063, end = 6.87 ± 0.064)
```

---

    UncertainBound(pos::Vector{T}, width::ContinuousUnivariateDistribution, uc::UncertainCurve{T, N}) where {T, N}

Create several `UncertainBound` objects with identical integration window widths in each draw.
Useful for integration of several symmetric peaks for which the same width may be assumed.

**Example**

Uncertainty of width is modeled with a scaled and shifted Beta(2, 2) distribution:

```jldoctest UBexample
julia> uc = begin # create uncertain curve with one symmetric peak
           x = 0:0.1:10;
           y = @. exp(-(x - 3)^2/0.15) + exp(-(x - 7)^2/0.15)
           add_noise(Curve(x, y), MvGaussianNoiseModel(0.1, 0.03, 0.5))
       end;

julia> ubs = UncertainBound([3., 7.], scale_shift_beta(2, 2, 1.3, 1.5), uc)
2-element Vector{UncertainBound{Float64, 10000}}:
 UncertainBound{Float64, 10000}(start = 2.3 ± 0.022, end = 3.7 ± 0.022)
 UncertainBound{Float64, 10000}(start = 6.3 ± 0.022, end = 7.7 ± 0.022)
```
"""
struct UncertainBound{T, N}
    left::Particles{T, N}
    right::Particles{T, N}
    # quantiles for local baseline subtraction and plotting
    _left_quantiles
    _right_quantiles
end

function Base.show(io::IO, c::UncertainBound{T, N}) where {T, N}
    print(io, "UncertainBound{$T, $N}(start = $(c.left), end = $(c.right))")
end

# Constructors

# Base constructor that injects quantiles
function UncertainBound(left::Particles{T, N}, right::Particles{T, N}) where {T, N}
    qs = BOUND_QUANTILES
    lqs, rqs = [[quantile(lr, q) for q in qs] for lr in (left, right)]
    return UncertainBound{T, N}(left, right, lqs, rqs)
end

# Create a left/right bound
function UncertainBound(left::S, right::T, N::Int=10_000) where {S <: ContinuousUnivariateDistribution, T <: ContinuousUnivariateDistribution}
    left  = Particles(N, left)
    right = Particles(N, right)
    return UncertainBound(left, right)
end

# Create several, correlated width bounds
function UncertainBound(
    pos::Vector{T},
    width::ContinuousUnivariateDistribution,
    uc::UncertainCurve{T, N}
) where {T, N}

    M = length(pos)
    left = Array{T}(undef, M, N)
    right = Array{T}(undef, M, N)
    width = Particles(N, width)

    for i ∈ 1:M
        pᵢ = pos[i]
        for j ∈ 1:N
            cⱼ = get_draw(j, uc) # sample j from curve
            wⱼ = get_draw(j, width) # sample j from width
            wⱼ < zero(T) && replace("""Negative integration window width encountered ($wⱼ),
                                       make sure the support region of your width distribution
                                       does not allow negative values!""", "\n" => " ") |> ArgumentError |> throw
            left[i, j], right[i, j] = left_right_from_peak(uc.x, cⱼ.y, pᵢ, wⱼ)
        end
    end

    return [UncertainBound(Particles(left[i, :]), Particles(right[i, :])) for i in 1:M]
end

# Create single width bound
function UncertainBound(
    pos::T,
    width::ContinuousUnivariateDistribution,
    uc::UncertainCurve{T, N}
) where {T, N}
    bnd = UncertainBound([pos], width, uc)
    return bnd[1]
end

"""
    get_draw(n, bnd::UncertainBound)

Retrieve the `n`th sample of the samples stored in `UncertainBound` `bnd`.
"""
get_draw(n, bnd::UncertainBound) = [get_draw(n, bnd.left), get_draw(n, bnd.right)]

"""
    mean(uc::UncertainBound)

Retrieve the mean of the `UncertainBound` start and end point.
"""
Statistics.mean(bnd::UncertainBound) = [Statistics.mean(bnd.left), Statistics.mean(bnd.right)]

"""
    scale_shift_beta(α, β, a, b)

Create a scaled and shifted `Beta(α, β)` distribution.
Samples fall in the interval [`a`, `b`].
"""
function scale_shift_beta(α, β, a, b)
    return LocationScale(a, b - a, Beta(α, β))
end


text:
# File containing functions to update particle positions based on different forces

include("../../src/part/UpdateParticles.jl")  # loads functions to be tested

# needs tests :)

max_stars_repo_name: pseastham/eFEMpart

max_issues_repo_name: pseastham/eFEMpart


text:
using Random

abstract type AbstractEnv end
# もしくは `using ReinforcementLearningBase` するか (要検討)
abstract type Node end

mutable struct InputNode{T<:Number} <: Node
    id::Int
    state::T
    targets::Array{Int, 1}
    weights::Array{T, 1}
end

mutable struct MiddleNode{T<:Number} <: Node
    id::Int
    state::T
    targets::Array{Int, 1}
    weights::Array{T, 1}
    bias::T
    activation_function::Function
end

mutable struct OutputNode{T<:Number} <: Node
    id::Int
    state::T
    bias::T
    activation_function::Function
end

mutable struct Model{T<:Number}
    nodes::Array{Node, 1}
    input_size::Int
    output_size::Int
    nodes_num::Int
    weights_num::Int
    score::T
end

# モデルを作る
function Model(input_size::Int, output_size::Int; T=Float64, rate=0.5, rng=MersenneTwister())
    nodes = Array{Node, 1}()
    nodes_num = input_size + output_size
    weights_num = 0

    # InputNodeの設定
    for i in 1:input_size
        # 出力へつなげたり繋げなかったりする
        targets = Array{Int, 1}()
        for j in (input_size+1):nodes_num
            if rand(rng) < rate
                push!(targets, j)
                weights_num += 1
            end
        end
        node = InputNode(i, T(0), targets, 2rand(rng, T, length(targets)).-1)
        push!(nodes, node)
    end

    # OutputNodeの設定
    for i in (input_size+1):nodes_num
        node = OutputNode(i, T(0), 2rand(T)-1, x -> x)
        push!(nodes, node)
    end

    return Model{T}(nodes, input_size, output_size, nodes_num, weights_num, T(0))
end

Model{T}(input_size, output_size; kwargs...) where {T} = Model(input_size, output_size; T=T, kwargs...)

# InputNode以外の全ノードのstateをbiasにする
function reset!(model::Model)
    for i in (model.input_size+1):length(model.nodes)
        model.nodes[i].state = model.nodes[i].bias
    end
    nothing
end

# モデルに入力を入れて計算する
function (model::Model{T})(input::Array{T, 1}) where T <: Number
    reset!(model)
    middlenode_ids = Array{Int, 1}()

    # InputNodeからの出力を計算
    for i in 1:model.input_size
        model.nodes[i].state = input[i]
        for j in 1:length(model.nodes[i].targets)
            target = model.nodes[i].targets[j]
            model.nodes[target].state += model.nodes[i].weights[j] * model.nodes[i].state
            if typeof(model.nodes[j]) <: MiddleNode
                push!(target, middlenode_ids)
            end
        end
    end

    # MiddleNodeからの出力を計算
    while length(middlenode_ids) > 0
        id = popfirst!(middlenode_ids)
        for j in 1:length(model.nodes[id].targets)
            target = model.nodes[i].targets[j]
            model.nodes[target].state += model.nodes[id].activation_function(model.nodes[id].weights[j] * model.nodes[id].state)
        end
        if typeof(model.nodes[j]) == MiddleNode
            push!(j, middlenode_ids)
        end
    end

    return [model.nodes[i].activation_function(model.nodes[i].state) for i in (model.input_size+1):(model.input_size+model.output_size)]
end

mutable struct Species
    representative_id::Int
    nodes_num::Int
    weights_num::Int
    model_ids::Array{Int, 1}
end

mutable struct Population{T<:Number}
    models::Array{Model{T}, 1}
    species_set::Array{Species, 1}
    species_num::Int
    generation::Int
end

function node_distance(model1::Model, model2::Model; compatibility_weight_coefficient=0.5, compatibility_disjoint_coefficient=0.5)
    disjoint_nodes = abs(model1.nodes_num - model2.nodes_num) * compatibility_disjoint_coefficient

    bias_distance = 0.0
    for node1 in model1.nodes
        if typeof(node1) <: InputNode
            continue
        end
        for node2 in model2.nodes
            if typeof(node2) <: InputNode
                continue
            end
            if node1.id == node2.id
                bias_distance += abs(node1.bias - node2.bias)
                break
            end
        end
    end
    bias_distance *= compatibility_weight_coefficient

    return (disjoint_nodes+bias_distance)/max(model1.nodes_num, model2.nodes_num)
end

function connection_distance(model1::Model, model2::Model; compatibility_weight_coefficient=0.5, compatibility_disjoint_coefficient=0.5)
    disjoint_connections = abs(model1.weights_num - model2.weights_num) * compatibility_disjoint_coefficient

    # nodeの数の6乗ぐらいの計算量のクソアルゴリズムなので改善したい
    weight_distance = 0.0
    for node1 in model1.nodes
        if typeof(node1) <: OutputNode
            continue
        end
        for node2 in model2.nodes
            if typeof(node2) <: OutputNode
                continue
            end
            if node1.id == node2.id

                for i in 1:length(node1.targets)
                    target1 = node1.targets[i]
                    for j in 1:length(node2.targets)
                        target2 = node2.targets[j]
                        if target1 == target2
                            weight_distance += abs(node1.weights[i] - node2.weights[j])
                            break
                        end
                    end
                end
                break

            end
        end
    end
    weight_distance *= compatibility_weight_coefficient

    return (disjoint_connections + weight_distance) / max(model1.weights_num, model2.weights_num)
end

# 種を分ける (最初の世代の場合)
function Population(models::Array{Model{T}, 1}; compatibility_threshold=3.0, compatibility_disjoint_coefficient=0.5, compatibility_weight_coefficient=0.5) where T
    species_set_toReturn = [Species(1, models[1].nodes_num, models[1].weights_num, [1])]
    species_set_toExamine = copy(species_set_toReturn)
    unexamined = [i for i in 2:length(models)]

    while length(unexamined) > 0
        current_species = pop!(species_set_toExamine)
        for i in copy(unexamined)
            distance = node_distance(models[current_species.representative_id], models[i],
                            compatibility_disjoint_coefficient=compatibility_disjoint_coefficient,
                            compatibility_weight_coefficient=compatibility_weight_coefficient)
                        + connection_distance(models[current_species.representative_id], models[i],
                            compatibility_disjoint_coefficient=compatibility_disjoint_coefficient,
                            compatibility_weight_coefficient=compatibility_weight_coefficient)

            if (distance > compatibility_threshold && length(species_set_toExamine) == 0)
                new_species = Species(i, models[i].nodes_num, models[i].weights_num, [i])
                push!(species_set_toReturn, new_species)
                push!(species_set_toExamine, new_species)
                deleteat!(unexamined, findfirst(x->x==i, unexamined))
            elseif (distance <= compatibility_threshold)
                push!(species_set_toReturn[end].model_ids, i)
                deleteat!(unexamined, findfirst(x->x==i, unexamined))
            end
        end
    end

    return Population{T}(models, species_set_toReturn, length(species_set_toReturn), 1)
end

# 種を分ける（最初の世代じゃない場合）
function Population(models::Array{Model{T}, 1}, previous_population::Population;
    compatibility_threshold=3.0, compatibility_disjoint_coefficient=1.0, compatibility_weight_coefficient=0.5) where T

    species_set_toReturn = Array{Species, 1}()
    previous_species_set_toExamine = copy(previous_population.species_set)
    unexamined = [i for i in 1:length(models)]

    # 前の世代の種に含まれるかどうか調査する
    for current_species in previous_species_set_toExamine
        isnew = true

        for i in unexamined
            distance = node_distance(previous_population.models[current_species.representative_id], models[i],
                            compatibility_disjoint_coefficient=compatibility_disjoint_coefficient,
                            compatibility_weight_coefficient=compatibility_weight_coefficient)
                        + connection_distance(previous_population.models[current_species.representative_id], models[i],
                            compatibility_disjoint_coefficient=compatibility_disjoint_coefficient,
                            compatibility_weight_coefficient=compatibility_weight_coefficient)

            if (distance <= compatibility_threshold && isnew)
                new_species = Species(i, models[i].nodes_num, models[i].weights_num, [i])
                push!(species_set_toReturn, new_species)
                deleteat!(unexamined, findfirst(x->x==i, unexamined))
                isnew = false
            elseif (distance <= compatibility_threshold)
                push!(species_set_toReturn[end].model_ids, i)
                deleteat!(unexamined, findfirst(x->x==i, unexamined))
            end
        end
    end

    # 新種を分ける
    species_set_toExamine = Array{Species, 1}()
    while length(unexamined) > 0
        new_species = Species(unexamined[1], models[unexamined[1]].nodes_num, models[unexamined[1]].weights_num, [unexamined[1]])
        popfirst!(unexamined)

        for i in unexamined
            distance = node_distance(models[new_species.representative_id], models[i],
                            compatibility_disjoint_coefficient=compatibility_disjoint_coefficient,
                            compatibility_weight_coefficient=compatibility_weight_coefficient)
                        + connection_distance(models[new_species.representative_id], models[i],
                            compatibility_disjoint_coefficient=compatibility_disjoint_coefficient,
                            compatibility_weight_coefficient=compatibility_weight_coefficient)

            if (distance <= compatibility_threshold)
                push!(new_species.model_ids, i)
                deleteat!(unexamined, findfirst(x->x==i, unexamined))
            end
        end
        push!(species_set_toReturn, new_species)
    end

    return Population{T}(models, species_set_toReturn, length(species_set_toReturn), previous_population.generation+1)
end

# 成績を評価
function evaluate!(model::Model, env::AbstractEnv)
    while !get_terminal(env)
        env(model(get_state(env)))
        model.score += get_reward(env)
    end
    nothing
end

function evaluate!(popu::Population, env::AbstractEnv)
    bestmodel_id = popu.model_ids[1]

    for model in popu.models
        evaluate!(model, env)
        if popu.models[bestmodel_id].score < model.score
            bestmodel_id = model.id
        end
    end

    return bestmodel_id
end

function crossover(parent1::Model, parent2::Model)
    child = parent1.nodes_num > parent2.nodes_num ? parent1 : parent2
    small_parent = parent1.nodes_num < parent2.nodes_num ? parent1 : parent2

    # 二つの親の InputNode に共通している weight をランダムに選択する
    for i in 1:child.input_size
        weight_ind = findall(x -> (small_parent.nodes[i].targets[x] ∈ child.nodes[i].targets), 1:length(small_parent.nodes[i].targets))
        for j in weight_ind
            if rand() < 0.5
                child.nodes[i].weights[j] = parent1.nodes[i].weights[j]
            else
                child.nodes[i].weights[j] = parent2.nodes[i].weights[j]
            end
        end
    end

    # 二つの親の MiddleNode に共通している weight をランダムに選択する
    for i in (child.input_size+child.output_size+1):child.nodes_num
        weight_ind = findall(x -> (small_parent.nodes[i].targets[x] ∈ child.nodes[i].targets), 1:length(small_parent.nodes[i].targets))
        for j in weight_ind
            if rand() < 0.5
                child.nodes[i].weights[j] = parent1.nodes[i].weights[j]
            else
                child.nodes[i].weights[j] = parent2.nodes[i].weights[j]
            end
        end
    end

    # 小さい parent の InputNode にしかない weight を追加する
    for i in 1:child.input_size
        adding_target = findall(x -> !(x ∈ child.nodes[i].targets), small_parent.nodes[i].targets)
        weight_ind = findall(x -> !(small_parent.nodes[i].targets[x] ∈ child.nodes[i].targets), 1:length(small_parent.nodes[i].targets))
        append!(child.nodes[i].targets, adding_target)
        append!(child.nodes[i].weight, small_parent.nodes[i].weight[weight_ind])
    end

    # 小さい parent の MiddleNode にしかない weight を追加する
    for i in (child.input_size+child.output_size+1):child.nodes_num
        adding_target = findall(x -> !(x ∈ child.nodes[i].targets), small_parent.nodes[i].targets)
        weight_ind = findall(x -> !(small_parent.nodes[i].targets[x] ∈ child.nodes[i].targets), 1:length(small_parent.nodes[i].targets))
        append!(child.nodes[i].targets, adding_target)
        append!(child.nodes[i].weight, small_parent.nodes[i].weight[weight_ind])
    end

    return child
end


text:
module JuliaNantes2019

end


text:
"""
    UniformMesh{T}
    UniformMesh(start::T, stop::T, length::Int) where {T}

1D uniform mesh data.

# Arguments
- `start::T` : beginning of the mesh
- `stop::T` : end of the mesh
- `length::Int` : number of cells of the mesh

# Implementation
- `step::T`  : size step
- `points::Vector{T}`: Array with node positions
- `width::T` : Distance between left and right edges.

"""
struct UniformMesh{T}
    step::T
    points::Vector{T}
    width::T

    function UniformMesh(start::T, stop::T, length::Int) where {T}
        pdeb = range(start, stop = stop, length = length + 1)
        points = pdeb[1:end-1]
        step_loc = step(pdeb)
        width = stop - start
        new{T}(step_loc, points, width)
    end
end
"""
    Base.step(mesh::UniformMesh)

Get the step of the mesh

# Argument
- `mesh::UniformMesh` : the mesh

# Return
- `step` : the step of the mesh that is the difference between two contiguous points
"""
Base.step(mesh::UniformMesh) = mesh.step

start(mesh::UniformMesh) = mesh.points[1]
stop(mesh::UniformMesh) = mesh.points[end] + mesh.step
"""
    Base.length(mesh::UniformMesh)

Get the length of the mesh

# Argument
- `mesh::UniformMesh` : the mesh

# Return
- `length` : the length of the mesh that is the number of points and cells
"""
Base.length(mesh::UniformMesh) = length(mesh.points)
"""
    points(mesh::UniformMesh)

Get the points of the mesh

# Argument
- `mesh::UniformMesh` : the mesh

# Return
- `points` : the points of the mesh that is the vector of all points of the mesh except the last
"""
points(mesh::UniformMesh) = mesh.points
"""
    width(mesh::UniformMesh)

Get the width of the mesh

# Argument
- `mesh::UniformMesh` : the mesh

# Return
- `width` : the width that is step*length or distance between left and right edges.
"""
width(mesh::UniformMesh) = mesh.width
"""
    vec_k_fft(mesh::UniformMesh{T}) where{T}

Get the fft coefficients of the mesh

# Argument
- `mesh::UniformMesh{T}` : the mesh

# Return
- fft coefficients
"""
function vec_k_fft(mesh::UniformMesh{T}) where {T}
    midx = div(length(mesh), 2)
    k = 2T(pi) / (width(mesh))
    res = k * vcat(0:midx-1, -midx:-1)
    return res
end

# @inline function newval(valref, val, borne, width)
#     diff = val - valref
#     if abs(diff) > borne
#         s = sign(diff)
#         divdiff, moddiff = divrem(diff, width)
#         if abs(moddiff) > borne
#             val -= (divdiff + s) * width
#         else
#             val -= divdiff * width
#         end
#     end
#     return val
# end

@inline function newval(valref, val, borne, width)
    diff = val - valref
    return abs(diff) >= borne ? valref + mod(diff +borne, width) - borne : val
end
# @inline function newval(valref::OpTuple{N,T}, val::OpTuple{N,T}, borne::OpTuple{N,T}, width::OpTuple{N,T}) where{N,T}
#     return OpTuple(ntuple(x -> newval(valref[x],val[x],borne[x],width[x]),N))
# end
function traitmodend!(lg::T2, f::Array{T2,N}) where {N, T2<:Union{OpTuple{N,<:Number},Number}} # T2 must OpTuple or Number
    f .= mod.(f, lg)
end
function traitmodbegin!(lg::T2, f::Array{T2,N}) where {N, T2<:Union{OpTuple{N,<:Number},Number}} # T2 must OpTuple or Number
    sz = size(f)
    borne = lg / 2
    flagmod = false
    list = [CartesianIndex(ntuple(x->1,N))]
    flagfait = zeros(Bool, sz)
    incr = map(x -> CartesianIndex(ntuple(y -> y == x, N)), 1:N)
    while !isempty(list)
        indice = popfirst!(list)
#        @show indice
        valref = f[indice]
        for i = 1:N
            newindice = indice + incr[i]
            if newindice[i] <= sz[i]
                if !flagfait[newindice]
                    val = f[newindice]
                    nval = newval(valref, val, borne, lg)
                    if nval != val
                        f[newindice] = nval
                        flagmod = true
                    end
                    flagfait[newindice] = true
                    push!(list, newindice)
                end
            end
        end
    end
    return flagmod
end

traitmodbegin!(mesh::UniformMesh{T}, f::Array{T,N}) where {T,N} = traitmodbegin!(width(mesh),f)
# traitmodbegin!(mesh::NTuple{N,UniformMesh{T}}, f::Array{OpTuple{T,N},N}) where {T,N} = traitmodbegin!(OpTuple(width.(mesh)),f)

# function traitmodbegin!(mesh::Union{NTuple{N,UniformMesh{T}},UniformMesh{T}}, f::Array{T,N}) where {T,N}
#     sz = size(f)
#     wdth = width.(mesh)
#     borne = wdth ./ 2
#     flagmod = false
#     list = [ones(Int, N)]
#     flagfait = zeros(Bool, sz)
#     incr = map(x -> map(y -> y == x, 1:N), 1:N)
#     while !isempty(list)
#         indice = popfirst!(list)
# #        @show indice
#         valref = f[indice...]
#         for i = 1:N
#             newindice = indice + incr[i]
#             if newindice[i] <= sz[i]
#                 if !flagfait[newindice...]
#                     val = f[newindice...]
#                     nval = newval(valref, val, borne, wdth)
#                     if nval != val
#                         f[newindice...] = nval
#                         flagmod = true
#                     end
#                     flagfait[newindice...] = true
#                     push!(list, newindice)
#                 end
#             end
#         end
#     end
#     return flagmod
# end
function traitmodend!(mesh::UniformMesh{T}, res::Array{T,N}) where {T,N}
    strt = start(mesh)
    res .-= strt
    traitmodend!(width(mesh), res)
    res .+= strt
end
# function traitmodend!(t_mesh::NTuple{N, UniformMesh{T}}, res::Array{OpTuple{N,T},N}) where {T,N}
#     strt = OpTuple(start.(t_mesh))
#     res .-= strt
#     traitmodend!(OpTuple(width.(t_mesh)), res)
#     res .+= strt
# end

stdtomesh(mesh::UniformMesh{T}, v) where{T}=mesh.step*v .+ mesh.points[1]
meshtostd(mesh::UniformMesh{T}, v) where{T}=(v .- mesh.points[1]) / mesh.step

# function gettuple_x(t_mesh::NTuple{N,UniformMesh{T}}) where {T,N}
#     sz = length.(t_mesh)
#     ret = ntuple(x -> zeros(T, sz), N)
#     for ind in CartesianIndices(sz)
#         for i = 1:N
#             ret[i][ind] = t_mesh[i].points[ind.I[i]]
#         end
#     end
#     return ret
# end


text:
using DifferentialEquations
using MATLAB
using ProgressMeter
using BSON: @save, @load

n_iter = Int(1e5)
barrier = 4
H_start = 0.125
H_end = 0.4999
HH = range(H_start, stop = H_end, length = n_iter)
exit_time = zeros(n_iter)
max_time = 1e5
tspan = (0.0, max_time)
length_max = 100000

ODE1(z, w) = conj(im * w .* (1 ./ (z .^ 2 - w .^ 2) + 1 ./ (1 + w .^ 2)))
ODE2(z, w) = conj(im * z .* (1 ./ (w .^ 2 - z .^ 2) + 1 ./ (1 + z .^ 2)))
condition_escape(u, t, integrator) =
    u[5] > length_max || maximum([abs(u[1]), abs(u[2]), abs(u[3]), abs(u[4])]) > barrier
function condition_hits_X(u, t, integrator) # Event when event_f(u,t) == 0
    u[4]
end

function affect_update_iterator_X!(integrator)
    integrator.u[5] = integrator.u[5] + 1
end
affect_stop!(integrator) = terminate!(integrator)

callback_escape = DiscreteCallback(condition_escape, affect_stop!)
callback_hits_PSS_X =
    ContinuousCallback(condition_hits_X, affect_update_iterator_X!, rootfind = false)
cb = CallbackSet(callback_hits_PSS_X, callback_escape)
function Eq_of_M(du, u, p, t)
    du[1] = real(ODE1(u[1] + im * u[2], u[3] + im * u[4]))
    du[2] = imag(ODE1(u[1] + im * u[2], u[3] + im * u[4]))
    du[3] = real(ODE2(u[1] + im * u[2], u[3] + im * u[4]))
    du[4] = imag(ODE2(u[1] + im * u[2], u[3] + im * u[4]))
    du[5] = 0
    return
end
# refernce
# z=X+iP W=Q+iY . Plane Q=0, P=0
# z=u[1]+i u[2] ; w=u[3]+i u[4]
function Yfind_Aref(Q, P, H_Aref)
    H = (2 * H_Aref)^2
    try
        Y = sqrt(
            (1 + (-1) * H * ((-1) + P^2)^2)^(-1) * (
                P^2 +
                (-1) * Q^2 +
                H * ((-1) + P^2)^2 * ((-1) + Q^2) +
                sqrt(
                    H * ((-1) + P^2)^4 +
                    (-4) *
                    ((-1) + H * ((-1) + P^2)^2) *
                    ((-1) * P^2 + H * ((-1) + P^2)^2) *
                    Q^2,
                )
            ),
        )
    catch
        Y = zeros(0)
    end
end
n_iter_Q = 2;#50
n_iter_P = 2
N = n_iter_P * n_iter_Q
Q_end = 1e-5
P_end = 1e-6

ArrP = range(0, stop = P_end, length = n_iter_P)
ArrQ = range(0, stop = Q_end, length = n_iter_Q)
mesh = [(Q, P) for P in ArrP, Q in ArrQ]
mesh_list = reshape(mesh, 1, :)
hits = zeros(n_iter)

# @showprogress for j = 1:n_iter
#     H = HH[j]
#     for k = 1:N
#         Q,P=mesh_list[k]
#         # println(mesh_list[k])
#         Y_0 = Yfind_Aref(Q, P, H)
#         u0 = [0;P;Q;Y_0;0]
#         prob = ODEProblem(Eq_of_M, u0, tspan)
#         t,A=solve(prob, Tsit5(),reltol=1e-6,abstol=1e-8,maxiters=1e15,callback=cb)
#         # t, A = solve(prob, Tsit5(), maxiters = 1e15, callback = cb)
#         # t,A=solve(prob, Tsit5(),reltol=1e-6,abstol=1e-9,maxiters=1e15)
#         if maximum([abs(A[1,end]),abs(A[2,end]),abs(A[3,end]),abs(A[4,end])]) > barrier
#             hits[j]=hits[j]+1
#         end
#     end
# end
# mat"figure();set(gcf, 'Position',  [0, 0, 1500, 1500]); hold on;"
# normalized_hits=hits/N
# mat"plot($HH,$normalized_hits,'b.','MarkerSize',10)"

# mat"axis([ $H_start, $H_end, 0,1])"

@showprogress for j in 1:n_iter
    H = HH[j]
    if H>.195
    Y_end = Yfind_Aref(Q_end, P_end, H)

    u0 = [0; P_end; Q_end; Y_end; 0]
    prob = ODEProblem(Eq_of_M, u0, tspan)
    sol= solve(prob, Tsit5(), reltol = 1e-7, abstol = 1e-10, maxiters = 1e15, callback = cb)
        # t, A = solve(prob, Tsit5(), maxiters = 1e15, callback = cb)
        # t,A=solve(prob, Tsit5(),reltol=1e-6,abstol=1e-9,maxiters=1e15)

        if sol.u[end][5]>=length_max-1
            exit_time[j] = max_time
        else
            exit_time[j] = sol.t[end]
        end
    else
        exit_time[j] = max_time
    end
end
mat"figure();set(gcf, 'Position',  [0, 0, 1500, 1500]); hold on;"
logT = log10.(exit_time)
# normalized_hits=hits/N
mat"plot($HH,$logT,'k.','MarkerSize',3)"
min=minimum(logT)
max=maximum(logT)
mat"axis([ .125, .4999, $min,$max])"


text:
using Core.Compiler: CodeInfo
import Core.Compiler: widenconst

struct CompClosure; opaque; end # TODO: Is this a YAKC?
(::CompClosure)(x) = error("Hello")

# TODO: Are these lattice elements or info components?
struct AbstractCompClosure
    order::Int
    seq::Int
    primal_info::Any
    prev_seq_infos::Vector
end

function Base.show(io::IO, cc::AbstractCompClosure)
    print(io, "∂⃖", superscript(cc.order), subscript(cc.seq))
    if isa(cc.primal_info, MethodMatchInfo) &&
       length(cc.primal_info.results.matches) == 1
        match = cc.primal_info.results.matches[1]
        print(io, match.method.name)
    else
        printstyled(io, "!!!", typeof(cc.primal_info), color=:red)
    end
end

widenconst(::AbstractCompClosure) = CompClosure

struct RRuleInfo
    rrule_rt
    info
end

struct CompClosInfo
    clos::AbstractCompClosure
    infos::Vector
end

struct PrimClosInfo
    next::Any
end

struct PrimClosure
    name::Symbol
    order::Int
    seq::Int
    dual::Any
    info_below::Any
    info_carried::Any
end
widenconst(::PrimClosure) = CompClosure

function Base.show(io::IO, pc::PrimClosure)
    printstyled(IOContext(io, :color=>true), "∂⃖", superscript(pc.order), subscript(pc.seq), pc.name; color=:green)
end

struct RecurseInfo
    info
end

struct ReifyInfo
    info
end

# Helpers
tuple_type_fields(rt) = isa(rt, PartialStruct) ? rt.fields : widenconst(rt).parameters