text:
using Clang.Generators

cd(@__DIR__)

@add_def off_t

function main()
    options = load_options(joinpath(@__DIR__, "wrap.toml"))

    args = get_default_args()
    includedir = joinpath(dirname(@__DIR__), "cmt", "include")
    push!(args, "-I$includedir")

    headers = [joinpath(includedir, "cmt", "cmt.h")]

    ctx = create_context(headers, args, options)

    build!(ctx, BUILDSTAGE_NO_PRINTING)

    # tweak exprs
    for node in get_nodes(ctx.dag)
        exprs = get_exprs(node)
        for (i, expr) in enumerate(exprs)
            exprs[i] = make_opaque(expr)
        end
    end

    build!(ctx, BUILDSTAGE_PRINTING_ONLY)
end

# cmt declares all Ns/Mt objects as `typedef void`, instead of using opaque types.
# make those definitions opaque instead to avoid potential conversiosn between them.
#
# TODO: fix cmt? consider LLVM's approach; look for LLVMModuleRef
function make_opaque(x::Expr)
    Meta.isexpr(x, :const) || return x
    y = x.args[1]
    Meta.isexpr(y, :(=)) || return x
    var, val = y.args
    if (startswith(String(var), "Mt") || startswith(String(var), "Ns")) && val == :Cvoid
        return :(struct $var end)
    end
    return x
end

isinteractive() || main()


text:
################################################################################
# Simulator Data
################################################################################

simdata_dim(mechanism::Mechanism) = sum(simdata_dim.(mech.contacts))
simdata_dim(contact::ContactConstraint) = sum(simdata_dim.(contact.constraints))
simdata_dim(model::NonlinearContact) = 5 # [friction_coefficient, p, contact_radius]
simdata_dim(model::LinearContact) = 5 # [friction_coefficient, p, contact_radius]
simdata_dim(model::ImpactContact) = 4 # [p, contact_radius]

function set_simulator_data!(model::NonlinearContact, data::AbstractVector)
	model.friction_coefficient = data[1]
    model.collision.contact_radius = norm(data[SVector{3,Int}(2:4)]) #TODO: Confirm this is correct
    model.collision.contact_origin = data[SVector{3,Int}(5:7)]
    return nothing
end

function set_simulator_data!(model::LinearContact, data::AbstractVector)
	model.friction_coefficient = data[1]
    model.collision.contact_radius = norm(data[SVector{3,Int}(2:4)]) #TODO: Confirm this is correct
    model.collision.contact_origin = data[SVector{3,Int}(5:7)]
    return nothing
end

function set_simulator_data!(model::ImpactContact, data::AbstractVector)
    model.collision.contact_radius = norm(data[SVector{3,Int}(2:4)]) #TODO: Confirm this is correct
    model.collision.contact_origin = data[SVector{3,Int}(5:7)]
    return nothing
end

function set_simulator_data!(mechanism::Mechanism, data::AbstractVector)
    c = 0
    for contact in mechanism.contacts
		for element in contact.constraints
			N = simdata_dim(element)
	        set_simulator_data!(element, data[c .+ (1:N)]); c += N
		end
    end
    return nothing
end

function get_simulator_data(model::NonlinearContact)
	return [model.friction_coefficient; model.collision.contact_radius; model.collision.contact_origin]
end

function get_simulator_data(modelset::LinearContact)
	return [model.friction_coefficient; model.collision.contact_radius; model.collision.contact_origin]
end

function get_simulator_data(model::ImpactContact)
	return [model.collision.contact_radius; model.collision.contact_origin]
end

function get_simulator_data(mechanism::Mechanism{T}) where T
    data = zeros(T,simdata_dim(mechanism))
	e = 0
    for contact in mechanism.contacts
		for element in contact.constraints
			N = simdata_dim(element)
        	data[e .+ (1:N)] = get_simulator_data(model); e += N
		end
    end
    return data
end

################################################################################
# Simulator Data Jacobian
################################################################################

function simdata_jacobian(mechanism::Mechanism{T,Nn,Ne,Nb,Ni}) where {T,Nn,Ne,Nb,Ni}
	nsol = solution_dimension(mechanism)
    ndata = simdata_dim(mechanism)
    data = get_simulator_data(mechanism)

	A = zeros(nsol, ndata)
	njoints = joint_dimension(mechanism)
	n = sum(length.(mech.joints.values))
	rb = njoints # row contact
	ri = njoints + 6Nb # row body
	c = 0 # column
	for contact in mechanism.contacts
		for element in contact.constraints
			N = length(contact)
			N½ = Int(N/2)
			Ns = simdata_dim(element)
			∇contact, ∇body = ∂g∂simdata(mechanism, contact) # assumes only one element per contact
			A[rb + 3 .+ (1:3), c .+ (1:Ns)] -= ∇body # only works for one body for now
			A[ri + N½ .+ (1:N½), c .+ (1:Ns)] -= ∇contact
			ri += N
			c += Ns
		end
	end
	return A
end

function ∂g∂simdata(mechanism, contact::ContactConstraint{T,N,Nc,Cs}) where {T,N,Nc,Cs<:NonlinearContact{T,N}}
    model = contact.model
	p = model.collision.contact_origin
	offset = model.collision.contact_normal' * model.collision.contact_radius
    body = get_body(mechanism, contact.parent_id)
    x2, v25, q2, ϕ25 = current_configuration_velocity(body.state)
    x3, q3 = next_configuration(body.state, mechanism.timestep)
	s = contact.impulses_dual[2]
	γ = contact.impulses[2]

	# Contribution to Injointonstraint
	∇friction_coefficient = SA[0,γ[1],0,0]
	∇off = [-model.collision.contact_normal; szeros(T,1,3); -model.contact_tangent * skew(vector_rotate(ϕ25, q3))]
	∇p = [model.collision.contact_normal * rotation_matrix(q3); szeros(T,1,3); model.contact_tangent * skew(vector_rotate(ϕ25, q3)) * rotation_matrix(q3)]
	∇contact = [∇friction_coefficient ∇off ∇p]

	# Contribution to Body dynamics
	∇friction_coefficient = szeros(T,3)
	X = force_mapping(model, x3, q3)
	# this what we differentiate: Qᵀγ = - skew(p - vector_rotate(offset, inv(q3))) * VRmat(q3) * LᵀVᵀmat(q3) * X' * γ
	∇off = - ∂skew∂p(VRmat(q3) * LᵀVᵀmat(q3) * X' * γ) * -rotation_matrix(inv(q3))
	∇p = - ∂skew∂p(VRmat(q3) * LᵀVᵀmat(q3) * X' * γ)
	∇body = [∇friction_coefficient ∇off ∇p]
	return ∇contact, ∇body
end

################################################################################
# Generate & Save Dataset
################################################################################
function datafilename(model::Symbol; kwargs...)
    eval(Symbol(:datafilename, model))(; kwargs...)
end

function datafilenamesphere(; N::Int=10, friction_coefficient=0.1, radius=0.5)
    "sphere_dim_N_$(N)_friction_coefficient_$(friction_coefficient)_radius_$(radius).jld2"
end

function datafilenameblock2d(; N::Int=10, friction_coefficient=0.1, radius=0.05, side=0.50)
    "block2d_dim_N_$(N)_friction_coefficient_$(friction_coefficient)_radius_$(radius)_side_$(side).jld2"
end

function datafilenamebox(; N::Int=10, friction_coefficient=0.1, radius=0., side=0.50)
    "box_dim_N_$(N)_friction_coefficient_$(friction_coefficient)_radius_$(radius)_side_$(side).jld2"
end

function datafilenamehardwarebox(; N::Int=10, s=1)
    "hardwarebox_dim_N_$(N)_sample_S_$(S).jld2"
end

function initial_state(model::Symbol; kwargs...)
    eval(Symbol(:initial_state_, model))(; kwargs...)
end

function initial_state_sphere(;
	xlims=[[0,0,0], [1,1,0.2]],
	vlims=[-1ones(3), 1ones(3)],
	ωlims=[-5ones(3), 5ones(3)],)
	x = xlims[1] + rand(3) .* (xlims[2] - xlims[1])
	v = vlims[1] + rand(3) .* (vlims[2] - vlims[1])
	ω = ωlims[1] + rand(3) .* (ωlims[2] - ωlims[1])
	return Dict(:x => x, :v => v , :ω => ω)
end

function initial_state_block2d(;
		xlims=[[0,0.2], [1,0.4]],
        vlims=[-ones(2), ones(2)],
		θlims=[-π, π],
        ωlims=[-10, 10],)
	x = xlims[1] + rand(2) .* (xlims[2] - xlims[1])
	v = vlims[1] + rand(2) .* (vlims[2] - vlims[1])
	θ = θlims[1] + rand() * (θlims[2] - θlims[1])
	ω = ωlims[1] + rand() * (ωlims[2] - ωlims[1])
	return Dict(:x => x, :v => v , :θ => θ, :ω => ω)
end

function initial_state_box(;
		xlims=[[0,0,0.2], [1,1,0.4]],
        vlims=[-2ones(3), [2,2,-1.]],
        ωlims=[-6ones(3), 6ones(3)],)
	x = xlims[1] + rand(3) .* (xlims[2] - xlims[1])
	v = vlims[1] + rand(3) .* (vlims[2] - vlims[1])
	ω = ωlims[1] + rand(3) .* (ωlims[2] - ωlims[1])
	return Dict(:x => x, :v => v , :q => normalize(rand(Quaternion{Float64})), :ω => ω)
end

function build_pairs(mechanism::Mechanism, trajs::AbstractVector)
    pairs = []
    for traj in trajs
        push!(pairs, build_pairs(mechanism, traj)...)
    end
    return pairs
end

function build_pairs(mechanism::Mechanism, traj::Storage{T,N}) where {T,N}
    pairs = []
    z = get_maximal_state(traj)
    for t = 1:N-1
        z1 = z[t]
        z2 = z[t+1]
        pair = [z1, z2]
        push!(pairs, pair)
    end
    return pairs
end

function generate_dataset(model::Symbol;
		N::Int=10, H=2.0, timestep=0.05, g=-9.81,
		opts=SolverOptions(btol=1e-6, rtol=1e-6),
		init_kwargs=Dict(), # xlims, vlims, θlims, ωlims...
		mech_kwargs=Dict(), # friction_coefficient, radius, side...
		sleep_ratio = 0.0,
		show_contact = true,
		)
    mechanism = get_mechanism(model, timestep=timestep, gravity=gravity; mech_kwargs...);
    trajs = []
    for i = 1:N
		state = initial_state(model; init_kwargs...)
        initialize!(mechanism, model; state...)
        storage = simulate!(mechanism, H, record=true, opts=opts)
        push!(trajs, storage)
        visualize(mechanism, storage, vis=vis, show_contact=show_contact)
		sleep(H*sleep_ratio)
    end

	data = get_simulator_data(mechanism)
    params = Dict(:N => N, :H => H, :timestep => timestep, :g => g, :data => data)
    pairs = build_pairs(mechanism, trajs)
    jldsave(joinpath(@__DIR__, "data", "dataset", datafilename(model; N = N, mech_kwargs...));
        params=params, trajs=trajs, pairs=pairs)
    return nothing
end

################################################################################
# Load Dataset
################################################################################
function open_dataset(model::Symbol; kwargs...)
    dataset = jldopen(joinpath(@__DIR__, "data", "dataset", datafilename(model; kwargs...)))
    params = dataset["params"]
    trajs = dataset["trajs"]
    pairs = dataset["pairs"]
    JLD2.close(dataset)
    return params, trajs, pairs
end

################################################################################
# Optimization Loss: Evaluation & Gradient
################################################################################
function getSimulatorMaxGradients(mechanism::Mechanism{T,Nn,Ne,Nb,Ni}) where {T,Nn,Ne,Nb,Ni}
	timestep= mechanism.timestep
	nu = input_dimension(mechanism)
	nsd = simdata_dim(mechanism)
	njoints = joint_dimension(mechanism)
	datamat = simdata_jacobian(mechanism)
	solmat = full_matrix(mechanism.system)

	data_jacobian = - solmat \ datamat
	∇data_vϕ = data_jacobian[njoints .+ (1:6Nb),:]
	data_jacobian̄ = zeros(13Nb,nsd)
	for (i, body) in enumerate(mechanism.bodies)
		# Fill in gradients of v25, ϕ25
		data_jacobian̄[13*(i-1) .+ [4:6; 11:13],:] += ∇data_vϕ[6*(i-1) .+ (1:6),:]

		# Fill in gradients of x3, q3
		x2 = body.state.x2
		q2 = body.state.q2
		v25 = body.state.vsol[2]
		ϕ25 = body.state.ϕsol[2]
		data_jacobian̄[13*(i-1) .+ (1:3),:] += linear_integrator_jacobian_velocity(x2, v25, timestep) * ∇data_vϕ[6*(i-1) .+ (1:3),:]
		data_jacobian̄[13*(i-1) .+ (7:10),:] += rotational_integrator_jacobian_velocity(q2, ϕ25, timestep) * ∇data_vϕ[6*(i-1) .+ (4:6),:]
	end
	return data_jacobian̄
end

function getSimulatorMaxGradients!(mechanism::Mechanism{T,Nn,Ne,Nb,Ni}, z::AbstractVector{T}, u::AbstractVector{T};
		opts=SolverOptions()) where {T,Nn,Ne,Nb,Ni}
	step!(mechanism, z, u, opts=opts)
	data_jacobian̄ = getSimulatorMaxGradients(mechanism)
	return data_jacobian̄
end

function clean_loss(model::Symbol, pairs, data; timestep=0.05, g=-9.81,
		opts=SolverOptions(btol=1e-6, rtol=1e-6), n_sample = 20, rot = 0.0)
	mechanism = get_mechanism(model, timestep=timestep, gravity=gravity)
	nsd = simdata_dim(mechanism)
	set_simulator_data!(mechanism, data)

	l = 0.0
	∇ = zeros(nsd)
    H = zeros(nsd, nsd)
	n = length(pairs)

	offset = Int(floor(rot))
	Random.seed!(0)
	mask = shuffle!(Vector(1:n))
	mask = [mask; mask]
	mask = mask[offset%n .+ (1:n_sample)]
    for i in mask
        li, ∇i, Hi = loss(mechanism, pairs[i], opts=opts)
        l += li
        ∇ += ∇i
		H += Hi
    end
    return l / n_sample, ∇ ./ n_sample, H ./ n_sample
end

function loss(mechanism::Mechanism, pair; opts=SolverOptions(btol=1e-6, rtol=1e-6))
	nu = input_dimension(mechanism)
	nz = maximal_dimension(mechanism)
	u = zeros(nu)
    z1 = pair[1]
    z2true = pair[2]
    z2 = step!(mechanism, z1, u, opts=opts)
    ∇z2 = getSimulatorMaxGradients(mechanism)
	# @show scn.(z2 - z2true)
	Q = Diagonal([ones(3); 1e-1ones(3); ones(4); 1e-1ones(3)])
    l = 0.5 * (z2 - z2true)'* Q *(z2 - z2true)
    ∇ = - ∇z2' * Q * (z2 - z2true)
	H = ∇z2'* Q * ∇z2 # quasi newton approx oft he Hessian
    return l, ∇, H
end

ext: jl


text:
"""
    vernacular(t::NCBITaxon)

This function will return `nothing` if no vernacular name is known, and an array
of names if found. It searches the "common name" and "genbank common name"
category of the NCBI taxonomy name table.
"""
function vernacular(t::NCBITaxon)
    x = NCBITaxonomy.names_table[findall(NCBITaxonomy.names_table.tax_id .== t.id),:]
    p = findall(!isnothing, indexin(x.class, [NCBITaxonomy.class_common_name, NCBITaxonomy.class_genbank_common_name]))
    return length(p) == 0 ? nothing : x.name[p]
end

"""
    synonyms(t::NCBITaxon)

This function will return `nothing` if no synonyms exist, and an array of names
if they do. It returns all of the
"""
function synonyms(t::NCBITaxon)
    x = NCBITaxonomy.names_table[findall(NCBITaxonomy.names_table.tax_id .== t.id),:]
    p = findall(isequal(NCBITaxonomy.class_synonym), x.class)
    return length(p) == 0 ? nothing : x.name[p]
end

"""
    authority(t::NCBITaxon)

This function will return `nothing` if no authority exist, and a string with the
authority if found.
"""
function authority(t::NCBITaxon)
    x = NCBITaxonomy.names_table[findall(NCBITaxonomy.names_table.tax_id .== t.id),:]
    p = findall(isequal(NCBITaxonomy.class_authority), x.class)
    return length(p) == 0 ? nothing : first(x.name[p])
end


text:
####
#### Scatter2D
####

"""
   set_label!(obj, v)

Internal function to set the label associated with an object `obj` to
string `v`.
"""
set_label!(o::Union{Hist2D, Fill2D, Scatter3D}, v::String) =
    (o.label = v;)

# label of a drawing (cf legend)
function set_labels!(o::Union{Scatter2D, Bar2D}, v::Vector{String})
    length(v) == o.nobj || throw(
        DimensionMismatch("labels // dimensions don't match")
    )
    o.labels = v
    return
end
# unusual calls...
set_labels!(o::Scatter2D, s::String) = set_labels!(o, fill(s, length(o.linestyles)))
set_labels!(o::Bar2D,     s::String) = set_labels!(o, fill(s, length(o.barstyles)))

####
#### Hist2D / Bar
####

set_bins!(o::Hist2D, v::Int) = (o.bins = v;)

function set_scaling!(o::Hist2D, v::String)
   o.scaling = get(HIST2D_SCALING, v) do
      throw(
            OptionValueError("scaling", v)
      )
   end
   return
end

set_horiz!(o::Union{Hist2D,Bar2D,Boxplot}, v::Bool) = (o.horiz = v)

####
#### Fill2D
####

set_xmin!(o::Fill2D, v::F64) = (o.xmin = v;)
set_xmax!(o::Fill2D, v::F64) = (o.xmax = v;)

####
#### Bar2D
####

set_stacked!(o::Bar2D, v::Bool) = (o.stacked = v;)

####
#### Heatmap
####

set_cmap!(h::Heatmap, c::Vector{<:Color}) = (h.cmap = c;)


text:
# Physical units.

# mass:        g / mol
# length:      nm
# time:        ps
# momentum:    g nm / ps mol
# energy:      kJ / mol
# temperature: K
# charge:      e

const deg_to_rad = pi / 180.0
const to_kilo = 1e-3
const to_nano = 1e9
const to_pico = 1e12

# NIST Special Publication 811 (2008), Appendix B.8.
const A_to_nm = 0.1
const kcal_to_kJ = 4.184

# 2019 SI definitions from BIPM.
const e_to_C = 1.602_176_634e-19

# 2018 CODATA recommended values from NIST.
const N_A = 6.022_140_76e23 # 1 / mol
const planck = 6.626_070_15e-34 * to_kilo * to_pico * N_A # kJ ps / mol
const hbar = planck / 2pi
const kB = 1.380_649e-23 * to_kilo * N_A # kJ / mol K

"""
    beta(T::Float64)

Compute the reciprocal temperature for `T`.
"""
beta(T::Float64) = 1.0 / (kB * T) # mol / kJ


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libmicrohttpd_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libmicrohttpd")
JLLWrappers.@generate_main_file("libmicrohttpd", UUID("9663f319-c281-55ef-904e-227b7281b1a1"))
end  # module libmicrohttpd_jll


text:
#######################################################################
# Coverage.jl
# Take Julia test coverage results and bundle them up in JSONs
# https://github.com/JuliaCI/Coverage.jl
#######################################################################

using Coverage, Compat, Compat.Test
using Suppressor

@testset "Coverage" begin

@testset "iscovfile" begin
    # test our filename matching. These aren't exported functions but it's probably
    # a good idea to have explicit tests for them, as they're used to match files
    # that get deleted
    @test Coverage.iscovfile("test.jl.cov")
    @test Coverage.iscovfile("test.jl.2934.cov")
    @test Coverage.iscovfile("/home/somebody/test.jl.2934.cov")
    @test !Coverage.iscovfile("test.ji.2934.cov")
    @test !Coverage.iscovfile("test.jl.2934.cove")
    @test !Coverage.iscovfile("test.jicov")
    @test !Coverage.iscovfile("test.c.cov")
    @test Coverage.iscovfile("test.jl.cov", "test.jl")
    @test !Coverage.iscovfile("test.jl.cov", "other.jl")
    @test Coverage.iscovfile("test.jl.8392.cov", "test.jl")
    @test Coverage.iscovfile("/somedir/test.jl.8392.cov", "/somedir/test.jl")
    @test !Coverage.iscovfile("/otherdir/test.jl.cov", "/somedir/test.jl")
end

@testset "isfuncexpr" begin
    @test !Coverage.isfuncexpr("2")
end

@testset "Processing coverage" begin
    cd(dirname(@__DIR__)) do
        datadir = joinpath("test", "data")
        # Process a saved set of coverage data...
        r = process_file(joinpath(datadir,"Coverage.jl"))
        # ... and memory data
        analyze_malloc(datadir)
        lcov = IOBuffer()
        # we only have a single file, but we want to test on the Vector of file results
        LCOV.write(lcov, FileCoverage[r])
        expected = read(joinpath(datadir, "expected.info"), String)
        if Compat.Sys.iswindows()
            expected = replace(expected, "SF:test/data/Coverage.jl\n" => "SF:test\\data\\Coverage.jl\n")
        end
        @test String(take!(lcov)) == expected

        # Test a file from scratch
        srcname = joinpath("test", "data","testparser.jl")
        covname = srcname*".cov"
        # clean out any previous coverage files. Don't use clean_folder because we
        # need to preserve the pre-baked coverage file Coverage.jl.cov
        clean_file(srcname)
        cmdstr = "include(\"$(escape_string(srcname))\"); using Compat, Compat.Test; @test f2(2) == 4"
        run(`$(Compat.Sys.BINDIR)/julia --code-coverage=user -e $cmdstr`)
        r = process_file(srcname, datadir)

        # Parsing seems to have changed slightly in Julia (or JuliaParser?) between v0.6 and v0.7.
        # Line 10 is the end of a @doc string (`""" ->`), and Line 11 is an expression (`f6(x) = 6x`)
        # In v0.6, the zero count goes to line 10, and in v0.7, the zero count goes (more correctly?)
        # to line 11
        # NOTE: The commit that made this change in Base was backported to 0.6.1, which necessitates
        # another VERSION check.

        if (VERSION.major == 0 && VERSION.minor == 7 && VERSION < v"0.7.0-DEV.468") || VERSION < v"0.6.1-pre.93"
            target = Union{Int64,Nothing}[nothing, 1, nothing, 0, nothing, 0, nothing, nothing, nothing, 0, nothing, nothing, nothing, nothing, nothing, nothing, 0, nothing, nothing, 0, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing]
        elseif (VERSION.major == 0 && VERSION.minor == 6)
            target = Union{Int64,Void}[nothing, 1, nothing, 0, nothing, 0, nothing, nothing, nothing, nothing, 0, nothing, nothing, nothing, nothing, nothing, 0, nothing, nothing, 0, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing]
        elseif VERSION < v"1.0.0"
            target = Union{Int64,Nothing}[nothing, 2, nothing, 0, nothing, 0, nothing, nothing, nothing, nothing, 0, nothing, nothing, nothing, nothing, nothing, 0, nothing, nothing, 0, nothing, nothing, 3, nothing, nothing, nothing, nothing, nothing, nothing]
        else
            target = Union{Int64,Nothing}[nothing, 2, nothing, 0, nothing, 0, nothing, nothing, nothing, nothing, 0, nothing, nothing, nothing, nothing, nothing, 0, nothing, nothing, 0, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing]
        end
        @test r.coverage == target

        covtarget = (sum(x->x != nothing && x > 0, target), sum(x->x != nothing, target))
        @test get_summary(r) == covtarget
        @test get_summary(process_folder(datadir)) != covtarget

        # Handle an empty coverage vector
        emptycov = FileCoverage("", "", [])
        @test get_summary(emptycov) == (0, 0)

        #json_data = Codecov.build_json_data(Codecov.process_folder("data"))
        #@test typeof(json_data["coverage"]["data/Coverage.jl"]) == Array{Union{Int64,Nothing},1}
        if VERSION < v"0.7.0"
            close(open("fakefile",true,true,true,false,false))
        else
            close(open("fakefile",read=true,write=true,create=true,truncate=false,append=false))
        end
        @test isempty(Coverage.process_cov("fakefile",datadir))
        rm("fakefile")
    end
end

@testset "codecovio.jl" begin
    """
    extracts the api URL from stdout in a codecov.io submit call
    very helpful for testing codecovio.jl
    """
    function extract_codecov_url(fun)
        data = @capture_out fun()
        lines = split(data, "\n")

        url = "None"
        get_next = false
        for line in lines
            if get_next
                url = line
                get_next = false
            end
            if occursin("Codecov.io API URL", line)
                get_next = true
            end
        end

        # println("url: $(url)")
        @assert url != "None" "unable to find codecov api url in stdout, check for changes in codecovio.jl"
        return url
    end

    # empty file coverage for testing
    fcs = FileCoverage[]

    # set up base system ENV vars for testing
    withenv(
        "CODECOV_URL" => nothing,
        "CODECOV_TOKEN" => nothing,
        "TRAVIS" => nothing,
        "TRAVIS_BRANCH" => nothing,
        "TRAVIS_COMMIT" => nothing,
        "TRAVIS_PULL_REQUEST" => nothing,
        "TRAVIS_JOB_ID" => nothing,
        "TRAVIS_REPO_SLUG" => nothing,
        "TRAVIS_JOB_NUMBER" => nothing,
        "APPVEYOR" => nothing,
        "APPVEYOR_PULL_REQUEST_NUMBER" => nothing,
        "APPVEYOR_ACCOUNT_NAME" => nothing,
        "APPVEYOR_PROJECT_SLUG" => nothing,
        "APPVEYOR_BUILD_VERSION" => nothing,
        "APPVEYOR_REPO_BRANCH" => nothing,
        "APPVEYOR_REPO_COMMIT" => nothing,
        "APPVEYOR_REPO_NAME" => nothing,
        "APPVEYOR_JOB_ID" => nothing,
        ) do

        # test local submission process

        # default values
        codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit_local(fcs; dry_run = true) )
        @test occursin("codecov.io", codecov_url)
        @test occursin("commit", codecov_url)
        @test occursin("branch", codecov_url)
        @test !occursin("service", codecov_url)

        # env var url override
        withenv( "CODECOV_URL" => "https://enterprise-codecov-1.com" ) do

            codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit_local(fcs; dry_run = true) )
            @test occursin("enterprise-codecov-1.com", codecov_url)
            @test occursin("commit", codecov_url)
            @test occursin("branch", codecov_url)
            @test !occursin("service", codecov_url)

            # function argument url override
            codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit_local(fcs; dry_run = true, codecov_url = "https://enterprise-codecov-2.com") )
            @test occursin("enterprise-codecov-2.com", codecov_url)
            @test occursin("commit", codecov_url)
            @test occursin("branch", codecov_url)
            @test !occursin("service", codecov_url)

            # env var token
            withenv( "CODECOV_TOKEN" => "token_name_1" ) do

                codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit_local(fcs; dry_run = true) )
                @test occursin("enterprise-codecov-1.com", codecov_url)
                @test occursin("token=token_name_1", codecov_url)
                @test !occursin("service", codecov_url)

                # function argument token url override
                codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit_local(fcs; dry_run = true, token="token_name_2") )
                @test occursin("enterprise-codecov-1.com", codecov_url)
                @test occursin("token=token_name_2", codecov_url)
                @test !occursin("service", codecov_url)
            end
        end

        # test faulty non-CI submission

        @test_throws ErrorException extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true) )

        # test travis-ci submission process

        # set up travis env
        withenv(
            "TRAVIS" => "true",
            "TRAVIS_BRANCH" => "t_branch",
            "TRAVIS_COMMIT" => "t_commit",
            "TRAVIS_PULL_REQUEST" => "t_pr",
            "TRAVIS_JOB_ID" => "t_job_id",
            "TRAVIS_REPO_SLUG" => "t_slug",
            "TRAVIS_JOB_NUMBER" => "t_job_num",
            ) do

            # default values
            codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true) )
            @test occursin("codecov.io", codecov_url)
            @test occursin("service=travis-org", codecov_url)
            @test occursin("branch=t_branch", codecov_url)
            @test occursin("commit=t_commit", codecov_url)
            @test occursin("pull_request=t_pr", codecov_url)
            @test occursin("job=t_job_id", codecov_url)
            @test occursin("slug=t_slug", codecov_url)
            @test occursin("build=t_job_num", codecov_url)

            # env var url override
            withenv( "CODECOV_URL" => "https://enterprise-codecov-1.com" ) do

                codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true) )
                @test occursin("enterprise-codecov-1.com", codecov_url)
                @test occursin("service=travis-org", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("pull_request=t_pr", codecov_url)
                @test occursin("job=t_job_id", codecov_url)
                @test occursin("slug=t_slug", codecov_url)
                @test occursin("build=t_job_num", codecov_url)

                # function argument url override
                codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true, codecov_url = "https://enterprise-codecov-2.com") )
                @test occursin("enterprise-codecov-2.com", codecov_url)
                @test occursin("service=travis-org", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("pull_request=t_pr", codecov_url)
                @test occursin("job=t_job_id", codecov_url)
                @test occursin("slug=t_slug", codecov_url)
                @test occursin("build=t_job_num", codecov_url)

                # env var token
                withenv( "CODECOV_TOKEN" => "token_name_1" ) do

                    codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true) )
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("token=token_name_1", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("pull_request=t_pr", codecov_url)
                    @test occursin("job=t_job_id", codecov_url)
                    @test occursin("slug=t_slug", codecov_url)
                    @test occursin("build=t_job_num", codecov_url)

                    # function argument token url override
                    codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true, token="token_name_2") )
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("token=token_name_2", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("pull_request=t_pr", codecov_url)
                    @test occursin("job=t_job_id", codecov_url)
                    @test occursin("slug=t_slug", codecov_url)
                    @test occursin("build=t_job_num", codecov_url)
                end
            end
        end

        # test appveyor submission process

        # set up appveyor env
        withenv(
            "APPVEYOR" => "true",
            "APPVEYOR_PULL_REQUEST_NUMBER" => "t_pr",
            "APPVEYOR_ACCOUNT_NAME" => "t_account",
            "APPVEYOR_PROJECT_SLUG" => "t_slug",
            "APPVEYOR_BUILD_VERSION" => "t_version",
            "APPVEYOR_REPO_BRANCH" => "t_branch",
            "APPVEYOR_REPO_COMMIT" => "t_commit",
            "APPVEYOR_REPO_NAME" => "t_repo",
            "APPVEYOR_JOB_ID" => "t_job_num",
            ) do

                # default values
                codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true) )
                @test occursin("codecov.io", codecov_url)
                @test occursin("service=appveyor", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("pull_request=t_pr", codecov_url)
                @test occursin("job=t_account%2Ft_slug%2Ft_version", codecov_url)
                @test occursin("build=t_job_num", codecov_url)

                # env var url override
                withenv( "CODECOV_URL" => "https://enterprise-codecov-1.com" ) do

                    codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true) )
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("service=appveyor", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("pull_request=t_pr", codecov_url)
                    @test occursin("job=t_account%2Ft_slug%2Ft_version", codecov_url)
                    @test occursin("build=t_job_num", codecov_url)

                    # function argument url override
                    codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true, codecov_url = "https://enterprise-codecov-2.com") )
                    @test occursin("enterprise-codecov-2.com", codecov_url)
                    @test occursin("service=appveyor", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("pull_request=t_pr", codecov_url)
                    @test occursin("job=t_account%2Ft_slug%2Ft_version", codecov_url)
                    @test occursin("build=t_job_num", codecov_url)

                    # env var token
                    withenv( "CODECOV_TOKEN" => "token_name_1" ) do

                        codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true) )
                        @test occursin("enterprise-codecov-1.com", codecov_url)
                        @test occursin("token=token_name_1", codecov_url)
                        @test occursin("branch=t_branch", codecov_url)
                        @test occursin("commit=t_commit", codecov_url)
                        @test occursin("pull_request=t_pr", codecov_url)
                        @test occursin("job=t_account%2Ft_slug%2Ft_version", codecov_url)
                        @test occursin("build=t_job_num", codecov_url)

                        # function argument token url override
                        codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true, token="token_name_2") )
                        @test occursin("enterprise-codecov-1.com", codecov_url)
                        @test occursin("token=token_name_2", codecov_url)
                        @test occursin("branch=t_branch", codecov_url)
                        @test occursin("commit=t_commit", codecov_url)
                        @test occursin("pull_request=t_pr", codecov_url)
                        @test occursin("job=t_account%2Ft_slug%2Ft_version", codecov_url)
                        @test occursin("build=t_job_num", codecov_url)
                    end
                end
            end

        # test circle ci submission process

        # set up circle ci env
        withenv(
            "CIRCLECI" => "true",
            "CIRCLE_PR_NUMBER" => "t_pr",
            "CIRCLE_PROJECT_USERNAME" => "t_proj",
            "CIRCLE_BRANCH" => "t_branch",
            "CIRCLE_SHA1" => "t_commit",
            "CIRCLE_PROJECT_REPONAME" => "t_repo",
            "CIRCLE_BUILD_URL" => "t_url",
            "CIRCLE_BUILD_NUM" => "t_num",
            ) do

            # default values
            codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true) )
            @test occursin("codecov.io", codecov_url)
            @test occursin("service=circleci", codecov_url)
            @test occursin("branch=t_branch", codecov_url)
            @test occursin("commit=t_commit", codecov_url)
            @test occursin("pull_request=t_pr", codecov_url)
            @test occursin("build_url=t_url", codecov_url)
            @test occursin("build=t_num", codecov_url)

            # env var url override
            withenv( "CODECOV_URL" => "https://enterprise-codecov-1.com" ) do

                codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true) )
                @test occursin("enterprise-codecov-1.com", codecov_url)
                @test occursin("service=circleci", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("pull_request=t_pr", codecov_url)
                @test occursin("build_url=t_url", codecov_url)
                @test occursin("build=t_num", codecov_url)

                # function argument url override
                codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true, codecov_url = "https://enterprise-codecov-2.com") )
                @test occursin("enterprise-codecov-2.com", codecov_url)
                @test occursin("service=circleci", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("pull_request=t_pr", codecov_url)
                @test occursin("build_url=t_url", codecov_url)
                @test occursin("build=t_num", codecov_url)

                # env var token
                withenv( "CODECOV_TOKEN" => "token_name_1" ) do

                    codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true) )
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("token=token_name_1", codecov_url)
                    @test occursin("service=circleci", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("pull_request=t_pr", codecov_url)
                    @test occursin("build_url=t_url", codecov_url)
                    @test occursin("build=t_num", codecov_url)

                    # function argument token url override
                    codecov_url = extract_codecov_url( () -> Coverage.Codecov.submit(fcs; dry_run = true, token="token_name_2") )
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("service=circleci", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("pull_request=t_pr", codecov_url)
                    @test occursin("build_url=t_url", codecov_url)
                    @test occursin("build=t_num", codecov_url)
                end
            end
        end
    end
end

end


text:
module  DifferentiableFunctionsTests

using Test
using ProximalBase
using Random
using LinearAlgebra

Random.seed!(1)

@testset "L2Loss" begin

  x = randn(5)
  y = randn(5)
  f1 = L2Loss()
  f2 = L2Loss(y)

  @test value(f1, x) ≈ sum(abs2,x)/2.
  @test value(f2, x) ≈ sum(abs2,x-y)/2.

  hat_x = similar(x)
  @test value_and_gradient!(f1, hat_x, x) ≈ sum(abs2,x)/2.
  @test hat_x ≈ x
  hat_x = similar(x)
  @test value_and_gradient!(f2, hat_x, x) ≈ sum(abs2,x-y)/2.
  @test hat_x ≈ x-y

end

@testset "Quadratic Function" begin

  A = randn(10, 10)
  A = A + A'
  b = randn(10)
  c = 1.

  x = randn(10)

  q1 = QuadraticFunction(A)
  q2 = QuadraticFunction(A, b)
  q3 = QuadraticFunction(A, b, c)

  @test value(q1, x) ≈ dot(x, A*x)/2.
  @test value(q2, x) ≈ dot(x, A*x)/2. + dot(x, b)
  @test value(q3, x) ≈ dot(x, A*x)/2. + dot(x, b) + c

  hat_x = similar(x)
  @test value_and_gradient!(q1, hat_x, x) ≈ dot(x, A*x)/2.
  @test hat_x ≈ A*x

  @test value_and_gradient!(q2, hat_x, x) ≈ dot(x, A*x)/2. + dot(x, b)
  @test hat_x ≈ A*x + b

  @test value_and_gradient!(q3, hat_x, x) ≈ dot(x, A*x)/2. + dot(x, b) + c
  @test hat_x ≈ A*x + b
end

@testset "Least Squares Loss" begin

  y = randn(10)
  X = randn(10, 5)
  b = rand(5)
  f = LeastSquaresLoss(y, X)

  @test value(f, b) ≈ norm(y-X*b)^2 / 2.
  grad_out = zeros(5)
  @test value_and_gradient!(f, grad_out, b) ≈ norm(y-X*b)^2 / 2.
  @test grad_out ≈ -X'*(y - X*b)

end

end # module


text:
using IcanProj
using NetCDF
using DataFrames
using PyPlot
using Statistics

# Settings

cfg = 32

variable = "rnet"

file_fine = "/data04/jmg/fsm_simulations/netcdf/fsmres_forest/results_$(cfg)/$(variable)_1km.nc"

file_coarse = "/data04/jmg/fsm_simulations/netcdf/fsmres_forest/results_$(cfg)/$(variable)_50km.nc"

#file_fine = "/data04/jmg/fsm_simulations/netcdf/forcings_st/$(variable)_1km.nc"

#file_coarse = "/data04/jmg/fsm_simulations/netcdf/forcings_st/$(variable)_50km.nc"

# Load data

df_links = link_results(file_fine, file_coarse)

data_coarse, data_aggregated, ngrids  = unify_results(file_fine, file_coarse, df_links, variable)

# Compute metrics

rmse = sqrt.(mean((data_coarse .- data_aggregated).^2 , dims = 1))

meanref = mean(data_aggregated, dims = 1)

meancmp = mean(data_coarse, dims = 1)

nrmse = rmse ./ meanref

bias = (meancmp .- meanref)

nse = 1 .- var(data_coarse .- data_aggregated, dims = 1) ./ var(data_aggregated .- mean(data_aggregated, dims = 1), dims = 1)

# Project to map

rmse_map = project_results(rmse[:], df_links)

meanref_map = project_results(meanref[:], df_links)

meancmp_map = project_results(meancmp[:], df_links)

nrmse_map = project_results(nrmse[:], df_links)

bias_map = project_results(bias[:], df_links)

nse_map = project_results(nse[:], df_links)

# Plot results

figure()
imshow(rmse_map)
cb = colorbar()
cb[:set_label]("RMSE (-)")
title(variable)

figure()
imshow(bias_map)
cb = colorbar()
cb[:set_label]("BIAS (-)")
title(variable)


text:
@testset "SymMatrix" begin
    Q = SymMatrix([1, 2, 3], 2)
    symmetric_setindex!(Q, 4, 1, 1)
    @test Q.Q == [4, 2, 3]
    symmetric_setindex!(Q, 5, 1, 2)
    @test Q.Q == [4, 5, 3]
    symmetric_setindex!(Q, 6, 2, 2)
    @test Q.Q == [4, 5, 6]
end

@testset "MomentMatrix" begin
    Mod.@polyvar x y
    @test_throws ArgumentError moment_matrix(measure([1], [x]), [y])
    μ = measure(1:3, [x^2, x*y, y^2])
    X = [x, y]
    for ν in (moment_matrix(μ, X), moment_matrix((i, j) -> i + j - 1, X))
        @test ν.Q[1:4] == [1, 2, 2, 3]
        @test ν.Q[1, 1] == 1
        @test ν.Q[1, 2] == 2
        @test ν.Q[2, 1] == 2
        @test ν.Q[2, 2] == 3
        @test variables(ν)[1] == x
        @test variables(ν)[2] == y
        @test nvariables(ν) == 2
    end
    @test_throws ArgumentError moment_matrix(measure([1], [x]), [y])
end

# [HL05] Henrion, D. & Lasserre, J-B.
# Detecting Global Optimality and Extracting Solutions of GloptiPoly 2005

@testset "[HL05] Section 2.3" begin
    Mod.@polyvar x y
    η = AtomicMeasure([x, y], WeightedDiracMeasure.([[1, 2], [2, 2], [2, 3]], [0.4132, 0.3391, 0.2477]))
    μ = measure(η, [x^4, x^3*y, x^2*y^2, x*y^3, y^4, x^3, x^2*y, x*y^2, y^3, x^2, x*y, y^2, x, y, 1])
    ν = moment_matrix(μ, [1, x, y, x^2, x*y, y^2])
    for lrc in (SVDChol(), ShiftChol(1e-14))
        atoms = extractatoms(ν, 1e-4, lrc)
        @test atoms !== nothing
        @test atoms ≈ η
    end
end

function testelements(X, Y, atol)
    @test length(X) == length(Y)
    for y in Y
        @test any(x -> isapprox(x, y, atol=atol), X)
    end
end
@testset "[HL05] Section 3.3.1" begin
    Mod.@polyvar x y z
    U = [1 0 0 0 0 0;
         0 1 0 0 0 0;
         0 0 1 0 0 0;
         2 0 0 0 0 0; # z^2 = 2z
         0 0 0 1 0 0;
         0 2 0 0 0 0; # y^2 = 2y
         0 0 0 0 1 0;
         0 0 0 0 0 1;
         0 0 2 0 0 0] # x^2 = 2x
    # β will be [z, y, x, y*z, x*z, x*y]
    x = monovec([z, y, x, z^2, y*z, y^2, x*z, x*y, x^2])
    # x*β contains x*y*z, x^2*z, x^2*y which are not present so it show fail
    V = MultivariateMoments.build_system(U', x, sqrt(eps(Float64)))
    @test iszerodimensional(V)
    testelements(V, [[2.0, 2.0, 2.0], [2.0, 2.0, 0.0], [2.0, 0.0, 2.0], [0.0, 2.0, 2.0], [2.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0], [0.0, 0.0, 0.0]], 1e-11)
end

@testset "[HL05] Section 4" begin
    Mod.@polyvar x y
    s3 = sqrt(1/3)
    η = AtomicMeasure([x, y], WeightedDiracMeasure.([[-s3, -s3], [-s3, s3], [s3, -s3], [s3, s3]], [0.25, 0.25, 0.25, 0.25]))
    μ = measure([1/9,     0,     1/9,     0, 1/9,   0,     0,     0,   0, 1/3,   0, 1/3, 0, 0, 1],
                [x^4, x^3*y, x^2*y^2, x*y^3, y^4, x^3, x^2*y, x*y^2, y^3, x^2, x*y, y^2, x, y, 1])
    ν = moment_matrix(μ, [1, x, y, x^2, x*y, y^2])
    for lrc in (SVDChol(), ShiftChol(1e-16))
        atoms = extractatoms(ν, 1e-16, lrc)
        @test atoms !== nothing
        @test atoms ≈ η
    end
end

@testset "[LJP17] Example ?" begin
    Mod.@polyvar x[1:2]
    η = AtomicMeasure(x, [WeightedDiracMeasure([1., 0.], 2.)])
    μ = measure([2., 0.0, 1e-6],
                monomials(x, 2))
    ν = moment_matrix(μ, monomials(x, 1))
    # All singular values will be at least 1e-6 > 1e-12 it won't eliminate any row
    @test extractatoms(ν, 1e-12, ShiftChol(1e-6)) === nothing
    # The following tests that the method does not error if ranktol eliminates everything
    # In particular, this tests that the function equation(i) do not call sum when r equal to 0
    # this that throws an ArgumentError as details in src/extract.jl
    @test extractatoms(ν, 1.0) === nothing
    atoms = extractatoms(ν, 1e-5)
    @test atoms !== nothing
    @test atoms ≈ η
end

@testset "[LJP17] Example ?" begin
    Mod.@polyvar x[1:2]
    η = AtomicMeasure(x, [WeightedDiracMeasure([1., 0.], 2.53267)])
    ν = moment_matrix([2.53267 -0.0 -5.36283e-19; -0.0 -5.36283e-19 -0.0; -5.36283e-19 -0.0 7.44133e-6], monomials(x, 2))
    atoms = extractatoms(ν, 1e-5)
    @test atoms !== nothing
    @test atoms ≈ η
end

@testset "[LJP17] Example ?" begin
    Mod.@polyvar x[1:2]
    # η1 is slightly closer to ν than η2
    η1 = AtomicMeasure(x, [WeightedDiracMeasure([1., 0.2278331868065880], 14.81069999810012)])
    η2 = AtomicMeasure(x, [WeightedDiracMeasure([1., 0.2278466167919198], 14.81065219789812)])
    η3 = AtomicMeasure(x, [WeightedDiracMeasure([1., 0.2278545192047005], 14.81062406873638)])
    η4 = AtomicMeasure(x, [WeightedDiracMeasure([1., 0.2278512018654204], 14.81063587727590)])
    η5 = AtomicMeasure(x, [WeightedDiracMeasure([1., 11.15240711893039],   6.038060135782876e-19),
                           WeightedDiracMeasure([1., 0.991590438508714],   8.48980575400066e-6),
                           WeightedDiracMeasure([1., 0.5482585022039119],  0.0011737977570447154),
                           WeightedDiracMeasure([1., -0.2808541188399676], 9.72284626287903e-5),
                           WeightedDiracMeasure([1., 0.22799858558381308], 14.781765717198045),
                           WeightedDiracMeasure([1., 0.12343866254667915], 0.02765476684449834)])
    ν = moment_matrix([14.8107 3.37426 0.769196 0.175447 0.0400656 0.00917416 0.00211654; 3.37426 0.769196 0.175447 0.0400656 0.00917416 0.00211654 0.000498924; 0.769196 0.175447 0.0400656 0.00917416 0.00211654 0.000498924 0.000125543; 0.175447 0.0400656 0.00917416 0.00211654 0.000498924 0.000125543 3.76638e-5; 0.0400656 0.00917416 0.00211654 0.000498924 0.000125543 3.76638e-5 1.62883e-5; 0.00917416 0.00211654 0.000498924 0.000125543 3.76638e-5 1.62883e-5 1.07817e-5; 0.00211654 0.000498924 0.000125543 3.76638e-5 1.62883e-5 1.07817e-5 1.10658e-5], monomials(x, 6))
    for (ranktol, η) in ((1e-3, η1), (1e-4, η1), (1e-5, η2), (1e-6, η3), (1e-7, η4), (1e-8, η5))
        # With 1e-3 and 1e-4, the rank is detected to be 1
        # With 1e-5, the rank is detected to be 2
        # With 1e-6, the rank is detected to be 3
        # With 1e-7, the rank is detected to be 5
        # With 1e-8, the rank is detected to be 6
        atoms = extractatoms(ν, ranktol)
        @test atoms !== nothing
        @test atoms ≈ η
    end
    # With 1e-9, the rank is detected to be 7
    atoms = extractatoms(ν, 1e-9)
    @test atoms === nothing
end

@testset "[LJP17] applied to [JCG14] Example 6.1" begin
    Mod.@polyvar x[1:4]
    η = AtomicMeasure(x, [WeightedDiracMeasure([1.0, 4.78736282579504, 1.24375738760842, -1.4231836829978], 0.039112791390926646)])
    ν = moment_matrix([0.0397951 0.187094 0.0489553 -0.0551816; 0.187094 0.896353 0.232962 -0.265564; 0.0489553 0.232962 0.0614682 -0.0676226; -0.0551816 -0.265564 -0.0676226 0.0837186], monomials(x, 1))
    atoms = extractatoms(ν, 6e-3)
    @test atoms !== nothing
    @test atoms ≈ η
    atoms = extractatoms(ν, 8e-4)
    @test atoms === nothing
end

@testset "Large norm" begin
    # If the norm of `M` is given to `rref!` instead of `√||M||`, `extractatoms` will error.
    Mod.@polyvar x[1:2]
    ν = moment_matrix([2749.376669556701 -800.152792847183; -800.152792847183 586.8034549325414], monomials(x, 1))
    @test nothing === extractatoms(ν, 1e-2)
end


text:
######################################################################
# Computing with polynomials
######################################################################

"""
Decide whether a given monic polynomial is a square, and if so compute a square root.

This test can fail in positive characteristic.
"""

function issquare(P::Nemo.PolyElem)
	d = degree(P)
	if isodd(d)
		return (false, P)
	else
		dP = Nemo.derivative(P)
		sqrt = Nemo.gcd(dP, P)
		if sqrt^2 == P
			return (true, sqrt)
		else
			return (false, sqrt)
		end
	end
end

######################################################################
# Computing with power series
######################################################################

"""
Get the polynomial associated with a power series.
"""
function convert(R::PolyRing{T}, F::AbsSeriesElem{T}) where T
	d = precision(F)
	pol = zero(R)
	for i = 0: d-1
		setcoeff!(pol, i, coeff(F, i))
	end
	return pol
end

"""
Get the power series associated with a polynomial.
"""
function convert(S::SeriesRing{T}, P::PolyElem{T}) where T
	d = min(degree(P), Nemo.max_precision(S)-1)
	series = zero(S)
	for i = 0:d
		setcoeff!(series, i, coeff(P, i))
	end
	return series
end

"""
Transform the coefficients of an absolute power series.
"""
function map(F::AbsSeriesElem, f::Function)
	R = parent(F)
	G = zero(R)
	d = precision(F)
	set_prec!(G, d)
	for i = 0: d-1
		setcoeff!(G, i, f(i, coeff(F, i)))
	end
	return G
end

"""
Compute the derivative of an absolute capped power series.
"""
function derivative(F::AbsSeriesElem)
	return shift_right(map(F, (i, c)-> i * c), 1)
end

"""
Compute the integral of an absolute capped power series, with constant term zero.

This does not change the maximum precision of the parent power series ring, however.
"""
function integral(F::AbsSeriesElem{T}) where T<:FieldElem
	return shift_left(map(F, (i, c)-> c // (i+1)), 1)
end

"""
Compose two absolute capped power series.
"""
function compose(F::AbsSeriesElem{T}, G::AbsSeriesElem{T}) where T
	R = parent(F)
	K = base_ring(R)
	A, _ = PolynomialRing(K, "x")
	Fpol = convert(A, F)
	Gpol = convert(A, G)
	comp = compose(Fpol, Gpol)
	return convert(R, comp)
end

######################################################################
# Berlekamp-Massey
######################################################################

"""
Berlekamp-Massey algorithm to compute the rational reconstruction of a given polynomial of degree 2M - 1 (or 2M - 2) mod X^2M.
"""

function berlekamp_massey(a::PolyElem{T}) where T<:FieldElem
    M = 1 + div(degree(a), 2)
    A = parent(a)
    x = gen(A)

	#Extended euclidean algorithm : we stop when the remainder has degree < M
    fjm2 = a
    fjm1 = x^(2*M)
	sjm1 = 0
    sjm2 = 1
    j = 0
    while degree(fjm1) >= M
        j += 1
        q, f = divrem(fjm2, fjm1)
        s = sjm2 - q * sjm1
		fjm1, fjm2 = f, fjm1
		sjm1, sjm2 = s, sjm1
	end

	t = sjm1
	u = fjm1
    den = inv(Nemo.lead(t)) * t  # make monic
	num = inv(Nemo.lead(u)) * u
    return num, den
end

######################################################################
# BMSS algorithm
######################################################################

"""
Compute the rational fraction of the normalized rational isogeny between elliptic curves of the form y^2 = f(x).
``E1`` and ``E2`` of *odd* degree ``deg``.

Assuming a rational normalized separable isogeny of degree ``deg`` exists between
``E1`` and
``E2``, kernelpoly returns the squarefree polynomial vanishing on the
abscissae of its kernel. unsafe_kernelpoly returns its square if the isogeny has odd degree, with some correction made on two-torsion points otherwise.

This algorithm works only if the characteristic is zero or greater than 4*deg - 1.
"""

function unsafe_kernelpoly(E1::AbstractWeierstrass{T}, E2::AbstractWeierstrass{T},
	deg::Integer) where T<:FieldElem

	K = base_ring(E1)
	#We hope no one will use this with function fields...
	if isa(K, Nemo.FinField)
		p = Nemo.characteristic(K)
	else
		p = 0
	end

	((p > 0) & (p <= 4*deg-1)) &&
		throw(ArgumentError("BMSS algorithm only works for characteristic 0 or greater than 4*deg - 1."))

    (a1, a2, a3, a4, a6) = a_invariants(E1)
    (b1, b2, b3, b4, b6) = a_invariants(E2)

    ((a1 != 0) | (a3 != 0) | (b1 != 0) | (b3 != 0)) &&
        throw(ArgumentError("Curves must have a model of the form y^2 = f(x)."))

	#We need precision 2*deg + 1
	R, x = PowerSeriesRing(K, 2*deg +1, "x", model=:capped_absolute)
    G = a6 * x^3 + a4 * x^2 + a2 * x + one(R)
    H = b6 * x^3 + b4 * x^2 + b2 * x + one(R)

    # solve the differential equation G(x) T'^2 = (T/x) H(T)
    sol = _BMSS_diffeq(G, H, 2*deg + 1)

    # We recover the rational fraction using the relation T = D(1/x) / N(1/x)
    U = shift_right(sol, 1)   #why ?
	A, _ = PolynomialRing(K, "x")
	Upol = convert(A, U)
    num, den = berlekamp_massey(Upol)
	N = reverse(den) #maybe not the right numerator
	D = reverse(num)

    # If the points of abscissa 0 are in the kernel, correct the degree of D
    gap = degree(N) - degree(D) - 1
    (gap > 0) && (D = shift_left(D, gap))

	#Making denominator monic
	lambda = inv(Nemo.lead(D))
	D = lambda * D
	N = lambda * N
    return (N, D)
end

raw"""
    Compute a power-series solution to the differential equation used by
    the BMSS algorithm. The differential equation is

    .. math::  P(T,x) = (T/x) H(T) - G(x) {T'}^2 = 0

    with initial conditions `T = O(x)`. The output is truncated to
    precision ``prec``.

    :param G: a power series with non-zero constant coefficient.
    :param H: a power series with non-zero constant coefficient.
    :param prec: an integer denoting the truncation order of the solution.

    :returns: the solution to the differential equation.
    :rtype: Power series

    :raises ZeroDivisionError: when the characteristic is smaller than
        ``2*prec-2`` and not 0.
    :raises ValueError: when ``G`` or ``H`` have constant term \neq 1.

    ALGORITHM:

    It is a Newton iteration. After some substitution we get

    .. math::

       T_0 = ax + O(x^2)\\
       T_{i+1} = T_i + T_i' \sqrt{G} \sqrt{x} \int \frac{k_i(x)}{2\sqrt{x}}
       + O(x^{2^i+1}),

    where

    .. math::

       k_i(x) = \frac{P(T_i, x)} {{T_i'}^2 G \sqrt{G}}.
"""

function _BMSS_diffeq(G::AbsSeriesElem{T}, H::AbsSeriesElem{T}, prec::Integer) where T<:FieldElem

    # The power series ring
    R = parent(G)
	one = Nemo.one(R)
	x = gen(R)
	K = base_ring(R)
	#We hope no one will use this with function fields...
	if isa(K, Nemo.FinField)
    	p = Nemo.characteristic(base_ring(R))
	else
		p = 0
	end

	#Checking the arguments are sane
    ((coeff(G,0) != 1) | (coeff(H,0) != 1)) &&
        throw(ArgumentError("Arguments must have constant coefficient one."))

    (0 < p & p <= 2*prec-3) &&
        throw(DivideError("Characteristic must be greater than 2*prec - 3 in order to compute a solution to precision 'prec'."))

    # the precision to which the solution T is known
    d = 1
    # 1/(T'^2 G) to precision d ; at the beginning we set T = x * H/G
    diffT2G = R(1) + Nemo.O(x)
    # Sqrt(G) and Sqrt(1/G) to precision d
    sqG = R(1) + Nemo.O(x)
	invsqG = R(1) + Nemo.O(x)

    sol = x + Nemo.O(x^2)

    while d < prec-1
        # update diffT, sqG and invsqG to precision d
        # (nothing changes in the first iteration)
        diffT2G = diffT2G * (2 - G * derivative(sol)^2 * diffT2G)
        sqG = inv(K(2)) * (sqG + G * invsqG * (2 - sqG * invsqG))
        invsqG = invsqG * (2 - sqG * invsqG)

        # double the current precision
        d = min(2*d, prec-1)
        set_prec!(diffT2G, d)
        set_prec!(sqG, d)
        set_prec!(invsqG, d)
        set_prec!(sol, d+1)

        f = (compose(H, sol) * shift_right(sol, 1) - G * derivative(sol)^2) * diffT2G * invsqG
        # F = 1/2 sqrt(x) integral(k/sqrt(x))
        F = shift_left( map(f, (i, c) -> c // (2*i+1)), 1)

        # update the solution
        sol += derivative(sol) * sqG * F
	end
    return sol
end

"""
Sanity checks
"""

function kernelpoly(E1::AbstractWeierstrass{T}, E2::AbstractWeierstrass{T},
	deg::Integer) where T<:FieldElem
	N, D = unsafe_kernelpoly(E1, E2, deg)

	evenpart = Nemo.gcd(divisionpolynomial(E1, 2), D)
    oddpart = Nemo.divexact(D, evenpart)
    check, sqrtoddpart = issquare(oddpart)
    if !check
        throw(ArgumentError("The two curves are not linked by a rational normalized isogeny of this degree"))
	end
    return evenpart * sqrtoddpart
end


text:
# script that wants to numerically test chatteris diaconis (misspelled with 99% prob)
# for the estimates of beta, fitness,ergm (many names..) in the DirBin1 case
using Utilities,AReg,StaticNets,JLD,MLBase,StatsBase,DynNets
using PyCall; pygui(:qt); using PyPlot

Nsample = 100
sizes = [50 100  200] #400 800 1200]
Nsizes = length(sizes)
estParSiz = Array{Array{Float64,2},1}(Nsizes)
dgpParSiz = Array{Array{Float64,1},1}(Nsizes)

for k=1:Nsizes
    ##Define the Dgp model
    N = sizes[k]
    deltaN = 5
    maxDeg = sqrt(N)
    tmpdegsIO = round.([ linspace(deltaN,maxDeg,N);Vector(linspace(deltaN,maxDeg,N))])
    dgpPar,dgpDegs, = estimate(ErgmDirBin1(tmpdegsIO))
    #remove infinities
    dgpPar[.!isfinite(dgpPar)] = -3

    dgpParSiz[k] = dgpPar

    L = maximum(abs.(dgpPar))
    C = L
    errBnd = C*sqrt(log(N)/N )
    prob = 1 - C/(N^2)
    dgpMod = ErgmDirBin1(zeros(dgpDegs),dgpPar,Int.(1:N))
    findmax(dgpDegs)
    #sample the dgp model
    sam = sampl(dgpMod,Nsample )

    est = zeros(2N,Nsample)
    degsIO_Sample = [sumSq(sam,2);sumSq(sam,1)]
    estPar = zeros(2N,Nsample)
    #estimate
    for i=1:Nsample
        println(i)
        estPar[:,i]= estimate(fooErgmDirBin1;degIO =degsIO_Sample[:,i])[1]
    end

    estPar[.!isfinite(estPar)] = -3
    estParSiz[k] = estPar
end

k = 1

diff =[ abs.(estParSiz[k] .- dgpParSiz[k]) for k=1:Nsizes]
plt[:hist]([diff[i][:] for i = 1:Nsizes] ,15,density=true )
 title(" A First Test for Asymptotic Consistency    ")
 legend(["N = $(sizes[k])  Bnd = C(L) $(round(sqrt(log(sizes[k])/sizes[k]),3))" for k=1:Nsizes])
  xlabel("θ_est - θ_dgp  ")
 grid(which = "both" , linewidth=1,alpha = 0.5)
##
Mod = ErgmDirBin1(deg,zeros(N),groupsInds)
groupsNames = unique(groupsInds)
NG = length(groupsNames)
tmpDic = countmap(groupsInds)
#number of nodes per group
NGMemb = sortrows([[i for i in keys(tmpDic)] [  tmpDic[i] for i in keys(tmpDic) ] ])
# agrregate degree of each group
degGroups = zeros(2NG) ; for i =1:N degGroups[groupsInds[i]] += deg[i]; degGroups[NG + groupsInds[i]] += deg[N + i];  end
L = sum(degGroups,1)[1]
L/(N^2-N)
groupsPar = ones(2NG)
firstOrderCond(Mod;degIO = deg,parGroupsIO = groupsPar, groupsInds = groupsInds )
expMatrix(Mod,groupsPar)

targetErrValStaticNets = 0.005
estPar,estIt,estMod = estimate(Mod)

##
sam = sampl(estMod,1000)

meanMatSam = squeeze(mean(sam,3),3)
degsSam = expValStatsFromMat(estMod,meanMatSam)
degGroupsSam = zeros(2NG) ; for i =1:N degGroupsSam[groupsInds[i]] += degsSam[i]; degGroupsSam[NG + groupsInds[i]] += degsSam[ N + i] end
relErrSam = maximum(abs.(( degGroupsSam - degGroups )./degGroups )[degGroups.!=0])

## Test of estimates efficiency on the data

## Repeat on piero's DATA
using MAT
pieroEmidEst = matread("/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/matlabFiles/data\ Piero/forecastEMID/fitness_timeseries.mat")
pieroEmidData = matread("/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/matlabFiles/data\ Piero/forecastEMID/eMID_weeklyaggregated_binaryDIRECT_postLTRO.mat")

pieroA_T = Bool.(pieroEmidData["A"])

# load my data
halfPeriod = false
## Load dataj
fold_Path =  "/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/juliaFiles/"
loadFilePartialName = "Weekly_eMid_Data_from_"
halfPeriod? periodEndStr =  "2012_03_12_to_2015_02_27.jld": periodEndStr =  "2009_06_22_to_2015_02_27.jld"
@load(fold_Path*loadFilePartialName*periodEndStr, AeMidWeekly_T,banksIDs,inactiveBanks, YeMidWeekly_T,weekInd,datesONeMid ,degsIO_T,strIO_T)

matA_T = AeMidWeekly_T# pieroA_T
degsIO_T = [sumSq(matA_T,2);sumSq(matA_T,1)]
T = size(matA_T)[3]
allFitSS =  StaticNets.estimate( StaticNets.SnapSeqNetDirBin1(degsIO_T); identPost = false,identIter= true )
dgpPar = meanSq(allFitSS,2 )
expDgpMat = StaticNets.expMatrix(fooErgmDirBin1,dgpPar)
dgpDegs = [sumSq(expDgpMat,2); sumSq(expDgpMat,1)]
N = length(matA_T[1,:,1])
Nsample = 1000
L = maximum(abs.(dgpPar))
C = L
errBnd = C*sqrt(log(N)/N )
prob = 1 - C/(N^2)
dgpMod = ErgmDirBin1(zeros(dgpDegs),dgpPar,Int.(1:N))
findmax(dgpDegs)
#sample the dgp model
sam = sampl(dgpMod,Nsample )

est = zeros(2N,Nsample)
degsIO_Sample = [sumSq(sam,2);sumSq(sam,1)]
estPar = zeros(2N,Nsample)
estDegs = zeros(2N,Nsample)
#estimate
for i=1:Nsample
    println(i)
    estPar[:,i]= estimate(fooErgmDirBin1;degIO =degsIO_Sample[:,i])[1]
    expEstMat = StaticNets.expMatrix(fooErgmDirBin1,estPar[:,i])
    estDegs[:,i] = [sumSq(expEstMat,2); sumSq(expEstMat,1)]

end

estPar[.!isfinite(estPar)] = -5

diffdat = estDegs .- dgpDegs # abs.(estPar .- dgpPar)
plt[:hist](diffdat[:]  ,100,density=true )
 title("Sample  Static Model, Run MLE and Conpare Degrees")
 grid(which = "both" , linewidth=1,alpha = 0.5)
 xlim(-3, 3)
 xlabel("deg_dgp - deg_MLE")

sampMeanDegs = meanSq(degsIO_Sample,2)
ind1= findmin(sampMeanDegs)[2]
ind2 = findmax(sampMeanDegs)[2]
diffdat = abs.(estPar .- dgpPar)./abs.(dgpPar)
plt[:hist]([diffdat[ind1,:] diffdat[ind2,:]]  ,50,density=true )
 title("DGP fit = mean sing snap estimates  ")
 legend(["Sample Mean Deg = $(round(sampMeanDegs[ind1],2)) " ; "Sample Mean Deg = $(round(sampMeanDegs[ind2],2)) "])
  xlabel("θ_est - θ_dgp  ")
 grid(which = "both" , linewidth=1,alpha = 0.5)

degsIO_Sample[1,:]
degsIO_Sample[ind2,:]

funProb(x::Array{T,1} where T<: Real ) = 1./(1+exp.(-x))
derfunProb(x::Array{T,1} where T<: Real )  =  exp.(-x)./((1+exp.(-x)).^2)
x = Vector(-5:0.001:5)
othPar = -3:1.5:3
for i = 1:length(othPar) plot(x, funProb(x.+othPar[i] )); end
 legTex = ["θ_2 =  $(othPar[i])" for i = 1:length(othPar)]
 title( "1/(1+ exp( - θ - θ_2 ) )")
 legend(legTex)
 xlabel("θ")
 grid()

#


text:
const SpecDict = Dict{String, Set{String}}()

function changestatespec(currws)
    !inRoom(currws) && return
    roomname = getRoomfromWS[currws]
    handle = getHandlefromWS[currws]

    if !haskey(SpecDict, roomname)
        SpecDict[roomname] = Set{String}()
    end

    push!(SpecDict[roomname], handle)

    msg = Dict{String, Any}()
    msg["responsetype"] = "room-change-status-spectate"
    msg["handle"] = handle
    broadcastmsg(roomname, JSON.json(msg))

    nothing
end

function changestateplay(currws)
    !inRoom(currws) && return
    roomname = getRoomfromWS[currws]
    handle = getHandlefromWS[currws]

    if haskey(SpecDict, roomname) && (handle in SpecDict[roomname])
        pop!(SpecDict[roomname], handle)

        msg = Dict{String, Any}()
        msg["responsetype"] = "room-change-status-play"
        msg["handle"] = handle
        broadcastmsg(roomname, JSON.json(msg))
    end

    nothing
end


text:
# include("../../__lib__/std/datastructs/src/datastructs_module.jl")

# module rw
# using ..datastructs

abstract type AbstractVariable end
abstract type AbstractTape end
abstract type AbstractNode end

struct Node{T} <: AbstractNode
    parder1::T
    parder2::T
    parent1::Int
    parent2::Int

    function Node{T}(parder1::T, parder2::T, parent1::Int, parent2::Int) where {T}
        return new(parder1, parder2, parent1, parent2)
    end
end

function Node(parder1::T, parder2, parent1, parent2) where {T}
    return Node{T}(promote(parder1, parder2)..., parent1, parent2)
end

# function create_node(tape::T, dzdx, dzdy, x::V) where {T <: AbstractTape, V <: AbstractVariable}

#     node = Node(dzdx, dzdy, x.index, 0)
#     push!(tape, node)

#     return node
# end

struct Tape
    stack::datastructs.Stack{Node}
    # stack::Vector{Node}

    function Tape()
        s = datastructs.Stack{Node}()
        # s = Vector{Node}(undef, 0)
        obj = new(s)
        return obj
    end
end
function Base.push!(tape::Tape, node::Node)
    return push!( tape.stack, node )
    # push!(tape.stack, node)
    return tape
end

# function last_valid_index(vec::Vector{T}) where {T}
#     return length(vec)
# end

# struct Tape <: AbstractTape
#     stack::datastructs.Stack{Node}

#     function Tape()
#         stack = datastructs.Stack{Node}()
#         return new(stack)
#     end
# end

# function Base.push!(tape::Tape, node::Node)
#     return push!( tape.stack, node )
# end
function Base.peek(tape::Tape)
    return peek( tape.stack )
end
function Base.length(tape::Tape)
    return length( tape.stack )
end
function Base.isempty(tape::Tape)
    return isempty( tape.stack )
end
function Base.size(tape::Tape)
    return size( tape.stack )
end

macro tape(tape_name)
    ex = :($(esc(tape_name)) = Tape())
    return ex
end

# must be mutable in order to work with a hash map
# or let the indices in the struct
struct Variable{T} <: AbstractVariable
    value::T
    tape::Tape
    index::Int

    function Variable{T}(value::T, tape::Tape, index::Int) where {T}
        return new(value, tape, index)
    end
end

function create_variable(tape::Tape, value::T) where {T}

    node  = Node(value, zero(T), 0, 0)
    push!(tape, node)

    index = datastructs.last_valid_index(tape.stack)
    # index = last_valid_index(tape.stack)
    var   = Variable{T}(value, tape, index)

    return var
end

function create_temporary_variable(tape::Tape, value::T) where {T}
    index = datastructs.last_valid_index(tape.stack)
    # index = last_valid_index(tape.stack)
    var   = Variable{T}(value, tape, index)
    return var
end

macro variable(tape, varname, value)
    ex = :($(esc(varname)) = create_variable($tape, $value))
    return ex
end

function create_node(tape::Tape, dzdx, dzdy, x::V1, y::V2) where {V1 <: AbstractVariable, V2 <: AbstractVariable}
    return __create_node(tape, dzdx, dzdy, x.index, y.index)
end
function create_node(tape::Tape, dzdx, dzdy, x::Int, y::Int)
    return __create_node(tape, dzdx, dzdy, x, y)
end
function create_node(tape::Tape, dzdx, dzdy, x::Int, y::V2) where {V2 <: AbstractVariable}
    return __create_node(tape, dzdx, dzdy, x, y.index)
end

function create_node(tape::Tape, dzdx, dzdy, x::V1, y::Int) where {V1 <: AbstractVariable}
    return __create_node(tape, dzdx, dzdy, x.index, y)
end

function __create_node(tape::Tape, dzdx, dzdy, x_index::Int, y_index::Int)
    node = Node(dzdx, dzdy, x_index, y_index)
    push!(tape, node)

    return node
end

# operations:
# +, -, *, /
# v1^2
# 2^v1
# v1^v2
# min, max

# sqrt(v1)
# abs
# sin, cos, tan
# sinh, cosh, tanh
# exp, log

function Base.:+(v1::T) where {T <: Variable}
    return v1
end
function Base.:-(v1::T) where {T <: Variable}

    dzdx = -1.0
    dzdy = 0.0
    create_node(v1.tape, dzdx, dzdy, v1, 0);

    tmp = -v1.value
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.:+(v1::T, v2::T) where {T <: Variable}

    dzdx = 1.0
    dzdy = 1.0
    create_node(v1.tape, dzdx, dzdy, v1, v2);

    tmp = v1.value + v2.value
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.:+(v1::T, v2::N) where {T <: Variable, N <: Number}

    dzdx = 1.0
    dzdy = 0.0
    create_node(v1.tape, dzdx, dzdy, v1, 0);

    tmp = v1.value + v2
    var = create_temporary_variable(v1.tape, tmp)

    return var
end
function Base.:+(v1::N, v2::T) where {T <: Variable, N <: Number}
    return v2 + v1
end

function Base.:-(v1::T, v2::T) where {T <: Variable}

    dzdx = 1.0
    dzdy = -1.0
    create_node(v1.tape, dzdx, dzdy, v1, v2);

    tmp = v1.value - v2.value
    var = create_temporary_variable(v1.tape, tmp)

    return var
end
function Base.:-(v1::T, v2::N) where {T <: Variable, N <: Number}

    dzdx = 1.0
    dzdy = 0.0
    create_node(v1.tape, dzdx, dzdy, v1, 0);

    tmp = v1.value - v2
    var = create_temporary_variable(v1.tape, tmp)

    return var
end
function Base.:-(v1::N, v2::T) where {T <: Variable, N <: Number}

    dzdx = 0.0
    dzdy = -1.0
    create_node(v2.tape, dzdx, dzdy, 0, v2);

    tmp = v1 - v2.value
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.:*(v1::T, v2::T) where {T <: Variable}

    dzdx = v2.value
    dzdy = v1.value
    create_node(v1.tape, dzdx, dzdy, v1, v2);

    tmp = v1.value * v2.value
    var = create_temporary_variable(v1.tape, tmp)

    return var
end
function Base.:*(v1::T, v2::N) where {T <: Variable, N <: Number}

    dzdx = convert(typeof(v1.value), v2)
    dzdy = v1.value
    create_node(v1.tape, dzdx, dzdy, v1, 0);

    tmp = v1.value * v2
    var = create_temporary_variable(v1.tape, tmp)

    return var
end
function Base.:*(v1::N, v2::T) where {T <: Variable, N <: Number}
    return v2 * v1
end

function Base.:/(v1::T, v2::T) where {T <: Variable}

    dzdx = 1 / v2.value
    dzdy = -v1.value / v2.value^2
    create_node(v1.tape, dzdx, dzdy, v1, v2);

    tmp = v1.value / v2.value
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.:/(v1::T, v2::N) where {T <: Variable, N <: Number}
    return v1 * (1/v2)
end

function Base.:/(v1::N, v2::T) where {T <: Variable, N <: Number}
    dzdx = zero(N)
    dzdy = -v1.value / v2.value^2
    create_node(v1.tape, dzdx, dzdy, v1, v2);

    tmp = v1 / v2.value
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.:^(v1::T, v2::N) where {T <: Variable, N <: Number}

    dzdx = v2 * v1.value ^ (v2-1)
    dzdy = 0.0
    create_node(v1.tape, dzdx, dzdy, v1, 0);

    tmp = v1.value ^ v2
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.:^(v1::N, v2::T) where {T <: Variable, N <: Number}

    dzdx = log(v1) * v1 ^ v2.value
    dzdy = 0.0
    create_node(v2.tape, dzdx, dzdy, 0, v2);

    tmp = v1 ^ v2.value
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.:^(v1::T, v2::T) where {T <: Variable}

    dzdx = v2.value * v1.value ^ (v2.value - 1)
    dzdy = log(v1.value) * v1.value ^ v2.value
    create_node(v1.tape, dzdx, dzdy, v1, v2);

    tmp = v1.value ^ v2.value
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.max(v1::T, v2::T) where {T <: Variable}

    dzdx = v1.value > v2.value
    dzdy = v1.value < v2.value
    create_node(v1.tape, dzdx, dzdy, v1, v2);

    tmp = max(v1.value, v2.value)
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.max(v1::T, v2::N) where {T <: Variable, N <: Number}

    dzdx = v1.value > v2
    dzdy = v1.value < v2
    create_node(v1.tape, dzdx, dzdy, v1, 0);

    tmp = max(v1.value, v2)
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.max(v1::N, v2::T) where {T <: Variable, N <: Number}

    dzdx = v2.value > v1
    dzdy = v2.value < v1
    create_node(v2.tape, dzdx, dzdy, 0, v2);

    tmp = max(v1, v2.value)
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.min(v1::T, v2::T) where {T <: Variable}

    dzdx = v1.value < v2.value
    dzdy = v1.value > v2.value
    create_node(v1.tape, dzdx, dzdy, v1, v2);

    tmp = min(v1.value, v2.value)
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.min(v1::T, v2::N) where {T <: Variable, N <: Number}

    dzdx = v1.value < v2
    dzdy = v1.value > v2
    create_node(v1.tape, dzdx, dzdy, v1, 0);

    tmp = min(v1.value, v2)
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.min(v1::N, v2::T) where {T <: Variable, N <: Number}

    dzdx = v2.value > v1
    dzdy = v2.value < v1
    create_node(v2.tape, dzdx, dzdy, 0, v2);

    tmp = max(v1, v2.value)
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.sin(v1::T) where {T <: Variable}

    dzdx = cos(v1.value)
    dzdy = 0.0
    create_node(v1.tape, dzdx, dzdy, v1, 0);

    tmp = sin(v1.value)
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.cos(v1::T) where {T <: Variable}

    dzdx = -sin(v1.value)
    dzdy = 0.0
    create_node(v1.tape, dzdx, dzdy, v1, 0);

    tmp = cos(v1.value)
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.tan(v1::T) where {T <: Variable}

    dzdx = (1/cos(v1.value))^2
    dzdy = 0.0
    create_node(v1.tape, dzdx, dzdy, v1, 0);

    tmp = tan(v1.value)
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.abs(v1::T) where {T <: Variable}

    dzdx = sign(v1.value)
    dzdy = 0.0
    create_node(v1.tape, dzdx, dzdy, v1, 0);

    tmp = abs(v1.value)
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.sqrt(v1::T) where {T <: Variable}

    dzdx = 1 / sqrt(v1.value)
    dzdy = 0.0
    create_node(v1.tape, dzdx, dzdy, v1, 0);

    tmp = sqrt(v1.value)
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.exp(v1::T) where {T <: Variable}

    dzdx = exp(v1.value)
    dzdy = 0.0
    create_node(v1.tape, dzdx, dzdy, v1, 0);

    tmp = exp(v1.value)
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.log(v1::T) where {T <: Variable}

    dzdx = 1 / (v1.value)
    dzdy = 0.0
    create_node(v1.tape, dzdx, dzdy, v1, 0);

    tmp = log(v1.value)
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.sinh(v1::T) where {T <: Variable}

    dzdx = cosh(v1.value)
    dzdy = 0.0
    create_node(v1.tape, dzdx, dzdy, v1, 0);

    tmp = sinh(v1.value)
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.cosh(v1::T) where {T <: Variable}

    dzdx = sinh(v1.value)
    dzdy = 0.0
    create_node(v1.tape, dzdx, dzdy, v1, 0);

    tmp = cosh(v1.value)
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

function Base.tanh(v1::T) where {T <: Variable}

    dzdx = (1/cosh(v1.value))^2
    dzdy = 0.0
    create_node(v1.tape, dzdx, dzdy, v1, 0);

    tmp = tanh(v1.value)
    var = create_temporary_variable(v1.tape, tmp)

    return var
end

struct Gradient{T}
    der::Vector{T}

    function Gradient{T}(n::Int) where {T}
        v = zeros(n)
        return new(v)
    end
end

function Gradient(n::Int)
    T = Float64
    return Gradient{T}(n)
end

function grad_wrt(g::Gradient, var::Variable)
    return g.der[ var.index ]
end

function create_gradient(var::Variable)

    L = length(var.tape)
    g = Gradient(L)

    # seed
    g.der[var.index] = one(Float64) # 1.0

    # traverse the tape in reverse
    for kk = L:-1:1

        node = var.tape.stack[kk]
        d    = g.der[kk]

        # update the adjoints for its parent nodes
        if node.parent1 > 0
            g.der[ node.parent1 ] += node.parder1 * d
        end
        if node.parent2 > 0
            g.der[ node.parent2 ] += node.parder2 * d
        end
    end

    return g
end

function derivative(::ReverseMode, fcn, x0::T) where {T <: Number}
    @tape t
    var = create_variable(t, x0)
    y   = fcn(var)
    gr  = create_gradient(y)
    g   = grad_wrt(gr, var)
    return g
end

function gradient(::ReverseMode, fcn, x0::Vector{T}) where {T <: Number}
    @tape t
    L       = length(x0)
    var_vec = Vector{Variable}(undef, L)

    for ii = 1:L
        var_vec[ii] = create_variable(t, x0[ii])
    end

    y  = fcn(var_vec)
    gr = create_gradient(y)
    g  = similar(x0)

    for ii = 1:L
        g[ii] = grad_wrt(gr, var_vec[ii])
    end

    return g
end

function jacobian(::ReverseMode, fcn, x0::Vector{T}) where {T <: Number}

    # f: R^n -> R^m
    # g_ij = df_i / dx_j

    @tape t
    L       = length(x0)
    var_vec = Vector{Variable}(undef, L)

    for ii = 1:L
        var_vec[ii] = create_variable(t, x0[ii])
    end

    y  = fcn(var_vec)
    m = length(y)
    jac = zeros(T, m, L) # m-by-n
    for ii = 1:m
        gr = create_gradient(y[ii])

        for jj = 1:L
            jac[ii, jj] = grad_wrt(gr, var_vec[jj])
        end

    end

    return jac

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using PSCOPF

using Test
using JuMP
using Dates
using DataStructures

@testset verbose=true "test_tso_slacks" begin

    #=
    Currently
    Pinj = min(Plim, uncertainties)
    so Pinj cannot be fixed to a level lower than uncertainties
    ie. if limit is 25, uncertainty is 20 and demand is 15
    the expression above would cause a production of min(25,20)=20MW
    So an extra 5MW (which is not allowed in the model)
                    S1   S2
    Load            15    25
    prod            20    25
    we could have (Not what we want, injection should = min(uncertainty, Plimit)) :
    Plim         =     25
    Pinj         = [15 , 25]
    B_islim      = [1  , 0]
    Pislim_x_lim = [0  , 0]
    loss_of_load    = [0  , 0]
    We have what we want :
    Plim         =     15
    Pinj         = [15 , 15]
    B_islim      = [0  , 1]
    Pislim_x_lim = [0  , 15]
    loss_of_load    = [0  , 10]

    TS: [11h]
    S: [S1]
                      bus 1
                        |
    (limitable) wind_1_1|
    Pmin=0, Pmax=100    |
    Csta=0, Cprop=1     |     load_1
    DP=>9h30            |
      S1: 20            | S1: 15
      S2: 25            | S2: 25
                        |
    =#
    @testset "tso_cant_cap_limitable_power_by_choosing_prod_level_after_dp" begin
        TS = [DateTime("2015-01-01T11:00:00")]
        ech = DateTime("2015-01-01T09:30:00")
        network = PSCOPF.Networks.Network()
        PSCOPF.Networks.add_new_bus!(network, "bus_1")
        # Limitables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 1.,
                                                Dates.Second(90*60), Dates.Second(90*60))
        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 20.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S2", 25.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 15.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S2", 25.)
        # firmness
        firmness = PSCOPF.Firmness(
                    SortedDict{String, SortedDict{Dates.DateTime, PSCOPF.DecisionFirmness} }(),
                    SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.TO_DECIDE))
                    )
        # initial generators state : No need because all pmin=0 => ON by default
        generators_init_state = SortedDict{String, PSCOPF.GeneratorState}()

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                        generators_init_state,
                                        uncertainties, nothing)

        context.market_schedule = PSCOPF.Schedule(PSCOPF.Market(), Dates.DateTime("2015-01-01T07:00:00"), SortedDict(
            "wind_1_1" => PSCOPF.GeneratorSchedule("wind_1_1",
                SortedDict(),
                SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>20.,"S2"=>25.)))
                ),
            )
        )

        tso = PSCOPF.TSOOutFO()
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_tso_schedule!(context, ech, result, firmness, tso)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK
        @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"]) < 1e-09
        @test 10. ≈ value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S2"])
        # Limitable produces to the available level
        @test 15. ≈ PSCOPF.get_prod_value(PSCOPF.get_tso_schedule(context), "wind_1_1", TS[1], "S1")
        @test 15. ≈ PSCOPF.get_prod_value(PSCOPF.get_tso_schedule(context), "wind_1_1", TS[1], "S2")
        # Limitable was capped when prod > load (ie. S1):
        @test -5. ≈ ( PSCOPF.get_prod_value(PSCOPF.get_tso_schedule(context), "wind_1_1", TS[1], "S1")
                    - PSCOPF.get_uncertainties(uncertainties[ech], "wind_1_1", TS[1], "S1") )
        @test -10. ≈ ( PSCOPF.get_prod_value(PSCOPF.get_tso_schedule(context), "wind_1_1", TS[1], "S2")
                    - PSCOPF.get_uncertainties(uncertainties[ech], "wind_1_1", TS[1], "S2") )
        # We do not pay for capped power
        #If it was due to TSO constraints, we would have paid for the generator used instead of limitables
        @test (10*1e7 + 1e-3) ≈ value(result.objective_model.penalty) # 10 loss_of_load + 1 limitation
        @test value(result.objective_model.start_cost) < 1e-09
        @test ((20. - 15.) *1. + (25. - 15.) *1. ) ≈ value(result.objective_model.prop_cost)
    end

    #=
    It is possible to cut consumption.
    This variable can be used in two cases :
    1- due to EOD constraints : we don't have enough production capacity (illustrated in S1)
    2- due to Pmin constraints : we don't have enough demand to start a unit (illustrated in S2)
    3- due to RSO constraints (c.f. tso_cuts_consumption_due_to_rso)

    TS: [11h]
    S: [S1]
                        bus 1
                         |
    (pilotable) prod_1_1 |    load_1
    Pmin=20, Pmax=100    |  S1: 150
    Csta=100k, Cprop=1   |  S2: 15
                         |  S3: 25
    =#
    @testset "tso_cuts_consumption_due_to_pmin" begin
        TS = [DateTime("2015-01-01T11:00:00")]
        ech = DateTime("2015-01-01T07:00:00")
        network = PSCOPF.Networks.Network()
        PSCOPF.Networks.add_new_bus!(network, "bus_1")
        # Pilotables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                                20., 100.,
                                                100000., 1.,
                                                Dates.Second(0), Dates.Second(0))
        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 150.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S2", 15.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S3", 25.)
        # firmness
        firmness = PSCOPF.Firmness(
                    SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),),
                    SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE))
                    )
        # initial generators state :
        generators_init_state = SortedDict("prod_1_1" => PSCOPF.OFF)

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                        generators_init_state,
                                        uncertainties, nothing)
        context.market_schedule = PSCOPF.Schedule(PSCOPF.Market(), Dates.DateTime("2015-01-01T07:00:00"), SortedDict(
            "prod_1_1" => PSCOPF.GeneratorSchedule("prod_1_1",
                SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.ON,
                                                                                SortedDict("S1"=>PSCOPF.ON, "S2"=>PSCOPF.OFF, "S3"=>PSCOPF.ON))),
                SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>100.,"S2"=>0.,"S3"=>25.)))
                ),
            )
        )

        tso = PSCOPF.TSOOutFO()
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_tso_schedule!(context, ech, result, firmness, tso)

        # indicate use of slacks for feasibility
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK
        # S1 : prod_capacity < load => cannot satisfy demand
        @test 100. ≈ PSCOPF.get_prod_value(PSCOPF.get_tso_schedule(context), "prod_1_1", TS[1], "S1")
        @test 50. ≈ value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"])
        # S2 : load < pmin => cannot start the unit for such load
        @test PSCOPF.get_prod_value(PSCOPF.get_tso_schedule(context), "prod_1_1", TS[1], "S2") < 1e-09
        @test 15. ≈ value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S2"])
        # S3 : works fine
        @test 25. ≈ PSCOPF.get_prod_value(PSCOPF.get_tso_schedule(context), "prod_1_1", TS[1], "S3")
        @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S3"]) < 1e-09
        # penalize cutting consumption
        @test 1e7 ≈ tso.configs.loss_of_load_penalty
        @test (50. * 1e7 + 15. * 1e7 + 0. ) ≈ value(result.objective_model.penalty)
        @test (0. + 0. + 0.) ≈ value(result.objective_model.start_cost) #The market paid for starting
        @test (100. + 0. + 25. ) ≈ value(result.objective_model.prop_cost) #FIXME : pay just for the difference ? i.e. delta*prop_cost instead of injected*prop_cost
    end

    #=
    TS: [11h]
    S: [S1]
                        bus 1
                         |
    (limitable) wind_1_1 |    load_1
    Pmin=0, Pmax=100     |  S1: 15
    Csta=0, Cprop=0.     |  S2: 25
    DP => 9h30           |
         S1 : 10         |
         S1 : 10         |
                         |
    (pilotable) prod_1_1 |
    Pmin=20, Pmax=100    |
    Csta=100k, Cprop=100 |
    =#
    @testset "tso_capping_limitables_due_to_pilotable_pmin" begin
        TS = [DateTime("2015-01-01T11:00:00")]
        ech = DateTime("2015-01-01T09:30:00")
        network = PSCOPF.Networks.Network()
        PSCOPF.Networks.add_new_bus!(network, "bus_1")
        # Pilotables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 1.,
                                                Dates.Second(90*60), Dates.Second(90*60))
        # Pilotables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                                20., 100.,
                                                100000., 100.,
                                                Dates.Second(0), Dates.Second(0))
        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 10.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S2", 10.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 15.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S2", 25.)
        # firmness
        firmness = PSCOPF.Firmness(
                    SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),),
                    SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.TO_DECIDE),
                               "prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE))
                    )
        # initial generators state :
        generators_init_state = SortedDict("prod_1_1" => PSCOPF.OFF)

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                        generators_init_state,
                                        uncertainties, nothing)
        context.market_schedule = PSCOPF.Schedule(PSCOPF.Market(), Dates.DateTime("2015-01-01T07:00:00"), SortedDict(
            "wind_1_1" => PSCOPF.GeneratorSchedule("wind_1_1",
                SortedDict(),
                SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>10.,"S2"=>10.)))
            ),
            "prod_1_1" => PSCOPF.GeneratorSchedule("prod_1_1",
                SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.ON,
                                                                                SortedDict("S1"=>PSCOPF.OFF, "S2"=>PSCOPF.ON))),
                SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>0.,"S2"=>20.)))
                ),
            )
        )

        tso = PSCOPF.TSOOutFO()
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_tso_schedule!(context, ech, result, firmness, tso)

        # indicates using slacks for feasibility (because of scenario S1)
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK

        # In S1 : Load=15, wind provides 10 => still missing 5 but pmin=20
        # => we want to reduce consumption by 5
        # But limit links the two scenarios and we can only produce to limit or uncertainty levels
        @test 5. ≈ PSCOPF.get_prod_value(PSCOPF.get_tso_schedule(context), "wind_1_1", TS[1], "S1")
        @test PSCOPF.get_prod_value(PSCOPF.get_tso_schedule(context), "prod_1_1", TS[1], "S1") < 1e-09
        @test 10. ≈ value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"])

        # In S2 : Load=25, wind provides 10 => still missing 15 but pmin=20
        # pilotable produces 20 => 5 extra prod (20+10 - 25) => reduce wind by 5.
        @test 5. ≈ PSCOPF.get_prod_value(PSCOPF.get_tso_schedule(context), "wind_1_1", TS[1], "S2")
        @test 20. ≈ PSCOPF.get_prod_value(PSCOPF.get_tso_schedule(context), "prod_1_1", TS[1], "S2")
        @test -5. ≈ ( PSCOPF.get_prod_value(PSCOPF.get_tso_schedule(context), "wind_1_1", TS[1], "S2")
                    - PSCOPF.get_uncertainties(uncertainties[ech], "wind_1_1", TS[1], "S2") )
        @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S2"]) < 1e-09
    end

    #=
    Capped power is the difference between the available power in a limitable
     and the decided production allowed for that limitable unit.

    TS: [11h]
    S: [S1]
                      bus 1
                        |
    (limitable) wind_1_1|
    Pmin=0, Pmax=100    |
    Csta=0, Cprop=2     |     load_1
      S1: 50            | S1: 15
                        |
    (limitable) wind_1_2|
    Pmin=0, Pmax=100    |
    Csta=0, Cprop=1     |
      S1: 10            |
                        |
    (limitable) wind_1_3|
    Pmin=0, Pmax=100    |
    Csta=0, Cprop=1     |
      S1: 0             |
    =#
    @testset "tso_capping_in_schedule" begin
        TS = [DateTime("2015-01-01T11:00:00")]
        ech = DateTime("2015-01-01T07:00:00")
        network = PSCOPF.Networks.Network()
        PSCOPF.Networks.add_new_bus!(network, "bus_1")
        # Limitables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 2.,
                                                Dates.Second(0), Dates.Second(0))
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_2", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 1.,
                                                Dates.Second(0), Dates.Second(0))
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_3", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 1.,
                                                Dates.Second(0), Dates.Second(0))
        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 50.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_2", DateTime("2015-01-01T11:00:00"), "S1", 10.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_3", DateTime("2015-01-01T11:00:00"), "S1", 0.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 15.)
        # firmness
        firmness = PSCOPF.Firmness(
                    SortedDict{String, SortedDict{Dates.DateTime, PSCOPF.DecisionFirmness} }(),
                    SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                                "wind_1_2" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                                "wind_1_3" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE)),
                    )
        # initial generators state : No need because all pmin=0 => ON by default
        generators_init_state = SortedDict{String, PSCOPF.GeneratorState}()

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                        generators_init_state,
                                        uncertainties, nothing)
        tso = PSCOPF.TSOOutFO()
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_tso_schedule!(context, ech, result, firmness, tso)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

        # Limitable injections
        @test 15. ≈ PSCOPF.get_prod_value(context.tso_schedule, "wind_1_1", TS[1], "S1")
        @test PSCOPF.get_prod_value(context.tso_schedule, "wind_1_2", TS[1], "S1") <= 1e-09
        @test PSCOPF.get_prod_value(context.tso_schedule, "wind_1_3", TS[1], "S1") <= 1e-09

        # Capped power
        @test (50. - 15) ≈ context.tso_schedule.capping["wind_1_1", TS[1], "S1"]
        @test (10. - 0)  ≈ context.tso_schedule.capping["wind_1_2", TS[1], "S1"]
        @test (0.  - 0)  ≈ context.tso_schedule.capping["wind_1_3", TS[1], "S1"]

        #Cost
        @test value(result.objective_model.start_cost) < 1e-09
        @test value(result.objective_model.prop_cost) ≈ (
              ((50. - 15.) * 2 + (10. - 0.) * 1. + (0. - 0.) * 1.)
        )
    end

    #=
    TS: [11h, 11h15]
    S: [S1,S2]
                        bus 1                   bus 2
                        |                      |
    (limitable) wind_1_1|       "1_2"          |
    Pmin=0, Pmax=200    |                      |
    Csta=0, Cprop=1     |                      |
      S1: 160           |----------------------|
                        |        5             |
                        |                      |
    load                |                      |load
      S1: 180           |                      |  S1: 20

    available prod : 160
    total load : 200
    => need to cut 40
    optimization decides on how much load to shed on each bus:
        bus1: cuts at least 15 (cause branch is limited to 5 : 20-5 = 15)
        bus2: cuts at most 40-15=25
        bus1+bus2 = 40
    =#
    @testset "tso_cutting_load_by_bus_in_schedule" begin
        TS = [DateTime("2015-01-01T11:00:00")]
        ech = DateTime("2015-01-01T07:00:00")
        network = PSCOPF.Networks.Network()
        # Buses
        PSCOPF.Networks.add_new_bus!(network, "bus_1")
        PSCOPF.Networks.add_new_bus!(network, "bus_2")
        # Branches
        PSCOPF.Networks.add_new_branch!(network, "branch_1_2", 5.);
        # PTDF
        PSCOPF.Networks.add_ptdf_elt!(network, "branch_1_2", "bus_1", 0.5)
        PSCOPF.Networks.add_ptdf_elt!(network, "branch_1_2", "bus_2", -0.5)
        # Limitables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                0., 200.,
                                                0., 1.,
                                                Dates.Second(0), Dates.Second(0))
        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 160.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 180.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_2", DateTime("2015-01-01T11:00:00"), "S1", 20.)
        # firmness
        firmness = PSCOPF.Firmness(
                    SortedDict{String, SortedDict{Dates.DateTime, PSCOPF.DecisionFirmness} }(),
                    SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE)),
                    )
        # initial generators state : No need because all pmin=0 => ON by default
        generators_init_state = SortedDict{String, PSCOPF.GeneratorState}()

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                        generators_init_state,
                                        uncertainties, nothing)
        tso = PSCOPF.TSOOutFO()
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_tso_schedule!(context, ech, result, firmness, tso)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK

        # possible production is 160 but load is 200
        @test 160. ≈ PSCOPF.get_prod_value(context.tso_schedule, "wind_1_1", TS[1], "S1")
        @test 40. ≈ (context.tso_schedule.loss_of_load_by_bus["bus_1", TS[1], "S1"]
                    + context.tso_schedule.loss_of_load_by_bus["bus_2", TS[1], "S1"])
        # cut conso by bus
        @test context.tso_schedule.loss_of_load_by_bus["bus_1", TS[1], "S1"] <= 25
        @test context.tso_schedule.loss_of_load_by_bus["bus_2", TS[1], "S1"] >= 15
    end

    #TODO : illustrate ptdf effect

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf