text:
# ------------------------------------------------------------------------------------------
# # Funciones
#
# Temas:
# 1. Cómo declara una función
# 2. Duck-typing en Julia
# 3. Funciones mutantes vs. no-mutantes
# 4. Broadcasting
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Cómo declarar una función
# Julia nos permite definir una función de varias maneras. La primera requiere de las
# palabras reservadas `function` y `end`
# ------------------------------------------------------------------------------------------

function sayhi(name)
    println("Hi $name, it's great to see you!")
end

function f(x)
    x^2
end

# ------------------------------------------------------------------------------------------
# Y las podemos llamar así
# ------------------------------------------------------------------------------------------

sayhi("C-3PO")

f(42)

# ------------------------------------------------------------------------------------------
# Alternativamente, las podemos declarar en una sóla línea
# ------------------------------------------------------------------------------------------

sayhi2(name) = println("Hi $name, it's great to see you!")

f2(x) = x^2

sayhi2("R2D2")

f2(42)

# ------------------------------------------------------------------------------------------
# Finalmente, pudimos declararlas como funciones "anónimas"
# ------------------------------------------------------------------------------------------

sayhi3 = name -> println("Hi $name, it's great to see you!")

f3 = x -> x^2

sayhi3("Chewbacca")

f3(42)

# ------------------------------------------------------------------------------------------
# ## Duck-typing en Julia
# *"If it quacks like a duck, it's a duck."*
# *"Si suena como pato, es un pato"*
#

# En Julia, las funciones operaran con el cualquier valor que haga sentido.

# Por ejemplo, `sayhi` funciona con el nombre de este personaje de tele, escrito como
# entero...
# ------------------------------------------------------------------------------------------

sayhi(55595472)

# ------------------------------------------------------------------------------------------
# Y `f` va a funcionar en una matriz.
# ------------------------------------------------------------------------------------------

A = rand(3, 3)
A

f(A)

# ------------------------------------------------------------------------------------------
# `f` funcionará con "hi" porque `*` para inputs de cadenas como concatenación.
# ------------------------------------------------------------------------------------------

f("hi")

# ------------------------------------------------------------------------------------------
# Por el otro lado, `f` no funcionará sobre un vector. A diferencia de `A^2`, la cual es una
# operación bien definida, el signifiado de `v^2` para un vector, `v`, es ambigua.
# ------------------------------------------------------------------------------------------

v = rand(3)

f(v)

# ------------------------------------------------------------------------------------------
# ## Funciones mutantes vs no-mutantes
#
# Por convención, funciones seguidas por un `!` alteran, o bien mutan, sus contenidos y las
# que carecen de un `!` no lo hacen.
#
# Por ejemplo, `sort` y `sort!`.
#
# ------------------------------------------------------------------------------------------

v = [3, 5, 2]

sort(v)

v

# ------------------------------------------------------------------------------------------
# `sort(v)` regresa el arreglo ordenado de `v`, pero `v` no cambia.

#
# Por otro lado, si corremos `sort!(v)`, el contenido del arreglo es ordenado dentro de `v`.
# ------------------------------------------------------------------------------------------

sort!(v)

v

# ------------------------------------------------------------------------------------------
# ## Broadcasting
#
# Si ponemos `.` entre el nombre de la funcion y su lista de argumento,

# le estamos diciendo a la función que se "difunda"/haga broadcasting sobre los elementos
# del input.

#
# Primero veamos la diferencia entre `f()` y `f.()`.

#
# Primero definimos una nueva matriz `A` que hará la diferencia fácil de observar
# ------------------------------------------------------------------------------------------

A = [i + 3*j for j in 0:2, i in 1:3]

f(A)

# ------------------------------------------------------------------------------------------
# Cómo se vio antes, para una matriz, `A`,
# ```
# f(A) = A^2 = A * A
# ```
#
# `f.(A)` por el otro lado va a regresar un objeto que contiene el cuadrado de `A[i, j]` en
# su entrada correspondiente.
# ------------------------------------------------------------------------------------------

B = f.(A)

A[2, 2]

A[2, 2]^2

B[2, 2]

# ------------------------------------------------------------------------------------------
# Esto significa que para `v`, `f.(v)` está definido, pero no para `f(v)` :
# ------------------------------------------------------------------------------------------

v = [1, 2, 3]

f.(v)

# ------------------------------------------------------------------------------------------
# ### Ejercicios
#
# 6.1 En vez de broadcastear `f` sobre `v`,  pudimos haber hecho `v .^ 2`.
#
# Sin declarar una nueva funcion, agrega 1 a cada elemento de una matriz de `3x3` llena de
# `0`'s.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
#
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 6.3 Una cifra de César recorre cada letra un número determinado de plazas más adelante en
# el abecedario. Un corrimiento, o shift, de 1 manda "A" a "B". Escribe una función llamada
# `cesar` que toma una cadena como input y un corrimiento y regresa una cadena desencriptada
# tal que obtengas
#
# ```julia
# cesar("abc", 1)
# "bcd"
#
# cesar("hello", 4)
# "lipps"
# ```
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
#
# Plot the equilibrium solutions for normalized CR3BP dynamics.
#

"""
Plot specified lagrange points in the rotating
reference frame of CR3BP system μ.
"""
function lagrangeplot(μ::Real, L=1:5; kwargs...)

    defaults = (; title  = "Nondimensional Lagrange Points",
                  xlabel = "X (DU)", ylabel="Y (DU)",
                  labels = ["Body 1" "Body 2" [string("L",i) for i ∈ L]...],
                  legend = :topleft, dpi=150)
    options  = merge(defaults, kwargs)

    lagrange_points = lagrange(μ)

    fig = scatter([-μ], [0]; markersize=10, markercolor=:lightblue, label="Body 1")
    scatter!(fig, [1-μ], [0]; markersize=6, markercolor=:gray, label="Body 2")

    colors = (:red, :orange, :tan, :cyan, :indigo)
    for point ∈ zip(lagrange_points, 1:length(lagrange_points))
        p, i = point
        scatter!(fig, [p[1]], [p[2]];
                markershape=:x, markercolor=colors[i], label=string("L", i))
    end

    scatter!(fig; options...)
    for i ∈ 1:min(length(fig.series_list), length(options.labels))
        fig.series_list[i].plotattributes[:label] = options.labels[i]
    end

    fig

end

"""
Plot specified lagrange points in the rotating
reference frame of CR3BP system μ.
"""
function lagrangeplot(sys::CircularRestrictedThreeBodySystem, L=1:5; normalize = true, exclude_z = true, kwargs...)

    defaults = (; title  = sys.name != "" ? "$(sys.name) Lagrange Points" : "Lagrange Plots",
                  xlabel = "X ($(string(normalized_length_unit(sys))))", ylabel="Y ($(string(normalized_length_unit(sys))))",
                  labels = [string("L",i) for i ∈ L],
                  legend = :topleft, dpi=150, aspect_ratio = 1)
    options  = merge(defaults, kwargs)

    μ = normalized_mass_parameter(sys)
    lagrange_points = lagrange(μ)

    colors = (:red, :orange, :tan, :cyan, :indigo)
    fig = scatter(; options...)
    for point ∈ zip(lagrange_points, 1:length(lagrange_points))
        p, i = point
        if !normalize
            p .= ustrip.(lengthunit(sys), p .* normalized_length_unit(sys))
        end
        if exclude_z
            scatter!(fig, [p[1]], [p[2]]; markershape=:x, markercolor=colors[i], markerstrokewidth=2, label=string("L", i))
        else
            scatter!(fig, [p[1]], [p[2], [p[3]]]; markershape=:x, markercolor=colors[i], markerstrokewidth=2, label=string("L", i))
        end
    end

    for i ∈ 1:min(length(fig.series_list), length(options.labels))
        fig.series_list[i].plotattributes[:label] = options.labels[i]
    end

    return fig

end

"""
Plot specified lagrange points in the rotating
reference frame of CR3BP system μ to the last figure.
"""
function lagrangeplot!(fig, sys::CircularRestrictedThreeBodySystem, L=1:5; normalize = true, exclude_z = true, kwargs...)

    defaults = (; )
    options  = merge(defaults, kwargs)

    μ = normalized_mass_parameter(sys)
    lagrange_points = lagrange(μ)

    colors = (:red, :orange, :tan, :cyan, :indigo)
    scatter!(fig; options...)
    for point ∈ zip(lagrange_points, 1:length(lagrange_points))
        p, i = point
        if !normalize
            p .= ustrip.(lengthunit(sys), p .* normalized_length_unit(sys))
        end
        if exclude_z
            scatter!(fig, [p[1]], [p[2]]; markershape=:x, markercolor=colors[i], markerstrokewidth=2, label=string("L", i))
        else
            scatter!(fig, [p[1]], [p[2], [p[3]]]; markershape=:x, markercolor=colors[i], markerstrokewidth=2, label=string("L", i))
        end
    end

    return fig

end

ext: jl


text: Dict{String, Any}("int64-max" => Dict{String, Any}("value" => "9223372036854775807", "type" => "integer"), "int64-max-neg" => Dict{String, Any}("value" => "-9223372036854775808", "type" => "integer"))


text:
push!(LOAD_PATH,"../src")

using Test
using TypesParticles
using TypesBasis
using BuildBasis
using InputParser
using ElecElecRepulsion

"""
		run_buildelecelecrepul_test()

""" function run_buildelecelecrepul_test()
    @testset "Electron-Electron Repulsion" begin
       system,basisset = getcalcsetup("../examples/HatreeFock.in")
       basisfunc,numelec = buildbasisfunc(system,basisset);
       elecelecrepcalc = buildelecelecrepulsion(numelec,basisfunc,basisset);
       elecelecreptest = gettestdata();
       tolerance = 1.0e-4;
       for i=1:4
           for j=1:4
               for k=1:4
                   for l=1:4
                       @test sign(elecelecrepcalc[i,j,k,l]) == sign(elecelecreptest[i,j,k,l]);
                       diff = abs(elecelecrepcalc[i,j,k,l]-elecelecreptest[i,j,k,l]);
                       @test diff <= tolerance;
	           end
               end #k=1:4
           end
       end #i=1:4

    end#testset

end #run_elecelecrepul_test()

"""
    Test data

"""
function gettestdata()
   testdata = zeros(4,4,4,4);

   testdata[:,:,1,1]= [1.07657   0.57847   0.77610   0.53268;
   0.57847   0.58740   0.48970   0.55254;
   0.77610   0.48970   0.89630   0.51160;
   0.53268   0.55254   0.51160   0.55730];

   testdata[:,:,2,1] =[0.57847   0.32942   0.42735   0.30366;
   0.32942   0.36130   0.28576   0.34043;
   0.42735   0.28576   0.51160   0.30076;
   0.30366   0.34043   0.30076   0.34503];

   testdata[:,:,3,1] = [0.77610   0.42735   0.61361   0.40139;
   0.42735   0.44585   0.40139   0.43066;
   0.61361   0.40139   0.77610   0.42735;
   0.40139   0.43066   0.42735   0.44585];

   testdata[:,:,4,1] = [0.53268   0.30366   0.40139   0.28144;
   0.30366   0.33360   0.27019   0.31697;
   0.40139   0.27019   0.48970   0.28576;
   0.28144   0.31697   0.28576   0.32383];

   testdata[:,:,1,2] = [0.57847   0.32942   0.42735   0.30366;
   0.32942   0.36130   0.28576   0.34043;
   0.42735   0.28576   0.51160   0.30076;
   0.30366   0.34043   0.30076   0.34503];

   testdata[:,:,2,2] = [0.58740   0.36130   0.44585   0.33360;
   0.36130   0.45315   0.32383   0.42858;
   0.44585   0.32383   0.55730   0.34503;
   0.33360   0.42858   0.34503   0.43911];

   testdata[:,:,3,2] = [0.48970   0.28576   0.40139   0.27019;
   0.28576   0.32383   0.28144   0.31697;
   0.40139   0.28144   0.53268   0.30366;
   0.27019   0.31697   0.30366   0.33360];

   testdata[:,:,4,2] = [0.55254   0.34043   0.43066   0.31697;
   0.34043   0.42858   0.31697   0.41183;
   0.43066   0.31697   0.55254   0.34043;
   0.31697   0.41183   0.34043   0.42858];

   testdata[:,:,1,3] = [0.77610   0.42735   0.61361   0.40139;
   0.42735   0.44585   0.40139   0.43066;
   0.61361   0.40139   0.77610   0.42735;
   0.40139   0.43066   0.42735   0.44585];

   testdata[:,:,2,3] = [0.48970   0.28576   0.40139   0.27019;
   0.28576   0.32383   0.28144   0.31697;
   0.40139   0.28144   0.53268   0.30366;
   0.27019   0.31697   0.30366   0.33360];

   testdata[:,:,3,3] = [0.89630   0.51160   0.77610   0.48970;
   0.51160   0.55730   0.53268   0.55254;
   0.77610   0.53268   1.07657   0.57847;
   0.48970   0.55254   0.57847   0.58740];

   testdata[:,:,4,3] = [0.51160   0.30076   0.42735   0.28576;
   0.30076   0.34503   0.30366   0.34043;
   0.42735   0.30366   0.57847   0.32942;
   0.28576   0.34043   0.32942   0.36130];

   testdata[:,:,1,4] = [0.53268   0.30366   0.40139   0.28144;
   0.30366   0.33360   0.27019   0.31697;
   0.40139   0.27019   0.48970   0.28576;
   0.28144   0.31697   0.28576   0.32383];

   testdata[:,:,2,4] = [0.55254   0.34043   0.43066   0.31697;
   0.34043   0.42858   0.31697   0.41183;
   0.43066   0.31697   0.55254   0.34043;
   0.31697   0.41183   0.34043   0.42858];

   testdata[:,:,3,4] = [0.51160   0.30076   0.42735   0.28576;
   0.30076   0.34503   0.30366   0.34043;
   0.42735   0.30366   0.57847   0.32942;
   0.28576   0.34043   0.32942   0.36130];

   testdata[:,:,4,4] = [0.55730   0.34503   0.44585   0.32383;
   0.34503   0.43911   0.33360   0.42858;
   0.44585   0.33360   0.58740   0.36130;
   0.32383   0.42858   0.36130   0.45315];

    return testdata
end #gettestdata


text:
"""
    weights(model)

Load the pre-trained weights for `model` using the stored artifacts.
"""
function weights(model)
  try
    path = joinpath(@artifact_str(model), "$model.bson")
    return BSON.load(path, @__MODULE__)[:weights]
  catch e
    throw(ArgumentError("No pre-trained weights available for $model."))
  end
end

"""
    loadpretrain!(model, name)

Load the pre-trained weight artifacts matching `<name>.bson` into `model`.
"""
loadpretrain!(model, name) = Flux.loadparams!(model, weights(name))


text:
using GLNS

datafile = joinpath( @__DIR__, "..", "examples", "test.gtsp" )
@time GLNS.solver( datafile )


text:
@testset "permutation" begin

N = [3,4]
M = prod(N)
Z = zeros(M,M)

@testset "create_πcmp" begin
    permute = [3,1,2]
    scale = [-1,1,-1]
    Kf = length(permute)
    KfM = Kf * M

    # order_cmpfirst = false
    P = create_πcmp(N, permute, scale=scale, order_cmpfirst=false)
    P2 = zeros(KfM,KfM)
    Id = Matrix{eltype(scale)}(I, M, M)

    for nw = 1:Kf
        iblk = permute[nw]
        jblk = nw
        I = (iblk-1)*M+1:iblk*M
        J = (jblk-1)*M+1:jblk*M
        P2[I,J] = scale[nw] * Id
    end
    @test P2 == P

    # order_cmpfirst = true
    P_cmpfirst = create_πcmp(N, permute, scale=scale)
    p = spzeros(eltype(scale), Kf, Kf)
    for nw = 1:Kf
        p[permute[nw], nw] = scale[nw]
    end
    P2_cmpfirst = blockdiag(ntuple(k->p, M)...)
    @test P2_cmpfirst == P_cmpfirst
end  # @testset "create_πcmp"

end  # @testset "permutation"


text:
#=
The area of a circle is defined as πr^2. Estimate π to 3 decimal places using a Monte Carlo method.

Hint: The basic equation of a circle is x2 + y2 = r2
=#
using Test
include("Solutions/problem14_estimate_pi.jl")

@test abs(estimate_pi() - pi) <= 0.001
@test abs(estimate_pi(number_of_decimal = 5) - pi) <= 0.00001

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Cone
# FIXME bugged, so not included!

"""
    Cone{radius_ratio} <: SimpleShape

A [`SimpleShape`](@ref) representing a truncated cone of unitary height, unitary base radius, and `radius_ratio` upper circle diameter.

The `radius_ratio` is part of the type signature so that a unitary shape can be defined for the type. The origin point is located at the center of the circular base.

# Members

- `transformation::Transformation`: the `Transformation` associated with the cylinder.
- `material::Material`: the [`Material`](@ref) of the cylinder.
"""
struct Cone{radius_ratio} <: SimpleShape
    transformation::Transformation
    material::Material
	function Cone(transformation::Transformation, material::Material, radius_ratio::Float32)
		abs(radius_ratio) <= 1f0 || throw(ArgumentError("the `radius_ratio` must have an absolute value less than 1: got $radius_ratio"))
		new{radius_ratio}(transformation, material)
	end
end

function Cone(; transformation::Transformation = Transformation(), material::Material = Material(), radius_ratio::Float32 = 0f0)
	Cone(transformation, material, radius_ratio)
end

@doc """
    Cone(transformation::Transformation, material::Material)

Constructor for a [`Cone`](@ref) instance.
""" Cone(::Transformation, ::Material)

@doc """
    Cone(transformation::Transformation = Transformation(),
           material::Material = Material())

Constructor for a [`Cone`](@ref) instance.
""" Cone(; ::Transformation, ::Material)

function get_t(::Type{Cone{RR}}, ray::Ray) where {RR}
	res = filter(t -> ray.tmin < t < ray.tmax, get_all_ts(Cone{RR}, ray))
	isempty(res) && return Inf32
	minimum(res)
end

function get_all_ts(::Type{Cone{RR}}, ray::Ray) where {RR}
    res = Vector{Float32}()
	â  = Normal{true}(0, 0, 1f0)
	dx, dy, dz = s⃗ = ray.dir
	ox, oy, oz = e⃗ = convert(Vec, ray.origin)
	R₀ = 1f0
	R₁ = RR - 1f0
	Z₀ = e⃗ ⋅ â
	Z₁ = s⃗ ⋅ â

	P⃗  = s⃗ - Z₁ * â
	Q⃗  = e⃗ - Z₀ * â

	B₂ = norm²(P⃗) - R₁^2 * Z₁^2
	B₁ = Z₁ * R₁ * (R₁ * Z₀ + R₀) - P⃗ ⋅ Q⃗
	B₀ = norm²(Q⃗) - (R₁ * Z₀ + R₀)^2

	if B₂ ≈ 0f0 && !(B₁ ≈ 0f0)
		t = B₀/2B₁
		0 <= Z₁ * t + Z₀ <= 1 && push!(res, t)
	elseif (Δ = B₁^2 - B₀ * B₂) >= 0f0
		t1 = (B₁ + sqrt(Δ))/B₂
		t2 = (B₁ - sqrt(Δ))/B₂
		0 <= Z₁ * t1 + Z₀ <= 1 && push!(res, t1)
		0 <= Z₁ * t2 + Z₀ <= 1 && push!(res, t2)
	end

    # check if caps are hit
    tz1, tz2 = (( 1f0 - oz) / dz, - oz / dz)

	if !iseven(length(res))
		(ox + tz1 * dx)^2 + (oy + tz1 * dy)^2 <= RR^2 ? push!(res, tz1) : push!(res, tz2)
		@assert (check = (ox + tz2 * dx)^2 + (oy + tz2 * dy)^2) <=  1f0 "it was actually $check !<= 1f0"
	elseif (ox + tz2 * dx)^2 + (oy + tz2 * dy)^2 <=  1f0
		append!(res, [tz1, tz2])
	end
	return res
end

function get_normal(::Type{Cone{RR}}, point::Point, ray::Ray) where {RR}
    z = ray.origin[3]
	abs_z = abs(z)
    # if it comes from the upper cap
	# or if it comes from the base
    (abs_z ≈ 1f0 || abs_z ≈ 0f0) && return -sign(z) * sign(ray.dir.z) * NORMAL_Z

    # if it comes from the side
	z = 1/(1f0 - RR) * sign(1/(1 - RR) - point[3])
    x, y = point.v[1:2]
	normal = Normal(x, y, z) |> normalize
    (normal ⋅ ray.dir < 0) ? normal : -normal
end

function get_uv(::Type{<:Cone}, point::Point)
    x, y, z = point.v
    z ≈ 0f0 && return Vec2D(x, (1 + y/2) * 0.25f0)
    return Vec2D(x, (3 + y/2) * 0.25f0)
end


text:
"""
    find_tokens(s, templates)

Go through a text left to right, one (valid) char at the time and keep track
of sequences of chars that match specific tokens. The list of tokens found is
returned.

**Arguments**

* `s`: the initial text
* `templates`: dictionary of possible tokens
"""
function find_tokens(
            s::SS,
            templates::Dict{Char, Vector{Pair{TokenFinder, Symbol}}}
            )::Vector{Token}

    tokens = Token[]
    isempty(s) && return tokens

    # start with a LINERETURN to allow proper treatment of special lines
    # such as blockquote lines or items that could be right at the start
    # see "process_line_return!"
    push!(tokens, Token(:SOS, subs(parent_string(s), from(s))))
    head_idx = firstindex(s)
    end_idx  = lastindex(s)

    @inbounds while head_idx <= end_idx
        head_char = s[head_idx]
        if haskey(templates, head_char)
            # Look at each possible finder sequentially
            for (tf, case) in templates[head_char]
                # ------------------------------------
                # exact match of a given fixed pattern
                # --> we form a candidate substring with a fixed number of characters
                # and try to see if it matches a fixed rule. Possibly the substring
                # contains an extra character for rules where we must match only when
                # the next character is or isn't something.
                if (tf.steps >= 0)
                    tail_idx = nextind(s, head_idx, tf.steps)
                    at_eos   = false
                    if tail_idx == nextind(s, end_idx)
                        tail_idx = end_idx
                        at_eos = true
                    end
                    (tail_idx > end_idx) && continue

                    # if there is space, consider the substring and verify whether it matches
                    candidate       = subs(s, head_idx, tail_idx)
                    matches, offset = fixed_lookahead(tf, candidate, at_eos)

                    # if it matches, form the token and break the for loop: no need to check
                    # other cases.
                    if matches
                        head_idx = prevind(s, tail_idx, offset)
                        token    = Token(case, chop(candidate, tail=offset))
                        push!(tokens, token)
                        break
                    end
                # -----------------------------------------
                # rule-based match: greedy catch until fail
                # --> we gradually form a candidate substring of increasing length until
                # the next character doesn't meet the condition.
                else
                    nchars    = 1
                    tail_idx  = head_idx
                    probe_idx = nextind(s, head_idx)
                    probe_idx > end_idx && continue
                    probe_char::Char = s[probe_idx]

                    # while the condition holds, consume get next char
                    while greedy_lookahead(tf, nchars, probe_char)
                        tail_idx   = probe_idx
                        probe_idx  = nextind(s, probe_idx)
                        (probe_idx > end_idx) && break
                        probe_char = s[probe_idx]
                        nchars    += 1
                    end

                    # if we took in at least a char, validate then form the token
                    if tail_idx > head_idx
                        candidate = subs(s, head_idx, tail_idx)
                        # check if the backward validator is happy otherwise skip
                        check(tf, candidate) || continue
                        # if it's happy move head and push the token
                        head_idx = tail_idx
                        token = Token(case, candidate)
                        push!(tokens, token)
                        break
                    end
                end
            end
        end
        head_idx = nextind(s, head_idx)
    end

    # finally push the end token on the stack observe that it can overlap a token
    # that would be at the end of the string.
    eos = Token(:EOS, subs(s, end_idx))
    push!(tokens, eos)

    # discard header tokens that are not at the start of a line or
    # only preceded by whitespaces
    process_header_tokens!(tokens)
    # validate or drop emphasis tokens
    process_emphasis_tokens!(tokens)
    # discard autolink_close tokens which are preceded by a space
    process_autolink_close_tokens!(tokens)
    return tokens
end

@inline find_tokens(s::String, templates) = find_tokens(subs(s), templates)

"""
    process_header_tokens!(tokens)

Discard header tokens that are not at the start of a line or only preceded by
whitespaces.
"""
function process_header_tokens!(tokens::Vector{Token})
    remove = Int[]
    @inbounds for (i, t) in enumerate(tokens)
        if t.name in MD_HEADERS
            ss = until_previous_line_return(t)
            isempty(strip(ss)) || push!(remove, i)
        end
    end
    deleteat!(tokens, remove)
end

"""
    process_emphasis_tokens!(tokens)

Process emphasis token candidates and either take them or discard them if
they don't look correct.

* `sTs` with token `T` is `s` is a space
* `xTs` with token `T` is a valid CLOSE if `x` is a character and `s` a space
* `sTx` with token `T` is a valid OPEN if `x` is a character and `s` a space
* `xTy` with token `T` is a valid MIXED if `x`, `y` are characters
"""
function process_emphasis_tokens!(tokens::Vector{Token})
    isempty(tokens) && return
    remove = Int[]
    ps = parent_string(first(tokens))
    N  = lastindex(ps)
    @inbounds for (i, t) in enumerate(tokens)
        if t.name in (:EM, :STRONG, :EM_STRONG)
            prev_char = previous_chars(t)
            next_char = next_chars(t)
            # if the token is surrounded by spaces, discard it
            bad = !isempty(prev_char) && first(prev_char) in SPACE_CHAR &&
                  !isempty(next_char) && first(next_char) in SPACE_CHAR
            # sTs
            if bad
                push!(remove, i)
            # Tx or sTx
            elseif isempty(prev_char) || first(prev_char) in SPACE_CHAR
                n = Symbol(string(t.name) * "_OPEN")
                tokens[i] = Token(n, t.ss)
            # xT or xTs
            elseif isempty(next_char) || first(next_char) in SPACE_CHAR
                n = Symbol(string(t.name) * "_CLOSE")
                tokens[i] = Token(n, t.ss)
            else # xTy
                n = Symbol(string(t.name) * "_MX")
                tokens[i] = Token(n, t.ss)
            end
        end
    end
    deleteat!(tokens, remove)
    return
end

"""
    process_autolink_close_tokens!(tokens)

Discard :AUTOLINK_CLOSE that are preceded by a space.
"""
function process_autolink_close_tokens!(tokens::Vector{Token})
    isempty(tokens) && return
    remove = Int[]
    @inbounds for (i, t) in enumerate(tokens)
        t.name == :AUTOLINK_CLOSE || continue
        c = previous_chars(t)
        (isempty(c) || first(c) ∈ SPACE_CHAR) && push!(remove, i)
    end
    deleteat!(tokens, remove)
    return
end


text:
# main-sc for movies
gmtset(FONT_TAG="14p,Helvetica-Bold")
grdimage("USEast_Coast.nc", proj="GMOVIE_COL0/MOVIE_COL1/160/210/55/0/36/34/MOVIE_WIDTH+", region=:g, I="int_US.nc", C="globe_US.cpt", X=0, Y=0)
plot("flight_path.txt", lw=1)


text:
function change!(v::GraphEdge, style::Symbol, to::Any, frames=:automatic)
end

function change_prop!(v::GraphEdge, property, to::Any; frames=:automatic)
end

function settheme(v::GraphEdge, theme; frames=:automatic)
end

function Javis.setopacity(v::GraphEdge, opacity; frames=:automatic)
end

function setscale(v::GraphEdge, scale; frames=:automatic)
end

function highlightedge(v::GraphEdge; frames=:automatic)
end


text:
"""
    `calculated_attributes!(trait::Type{<: AbstractPlot}, plot)`
trait version of calculated_attributes
"""
calculated_attributes!(trait, plot) = nothing

"""
    `calculated_attributes!(plot::AbstractPlot)`
Fill in values that can only be calculated when we have all other attributes filled
"""
calculated_attributes!(plot::T) where T = calculated_attributes!(T, plot)

"""
    image(x, y, image)
    image(image)

Plots an image on range `x, y` (defaults to dimensions).

"""
@recipe(Image, x, y, image) do scene
    Attributes(;
        default_theme(scene)...,
        colormap = [:black, :white],
        colorrange = automatic,
        interpolate = true,
        fxaa = false,
        lowclip = nothing,
        highclip = nothing,
        inspectable = theme(scene, :inspectable)
    )
end

"""
    heatmap(x, y, values)
    heatmap(values)

Plots a heatmap as an image on `x, y` (defaults to interpretation as dimensions).

"""
@recipe(Heatmap, x, y, values) do scene
    Attributes(;
        default_theme(scene)...,
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        linewidth = 0.0,
        interpolate = false,
        levels = 1,
        fxaa = true,
        lowclip = nothing,
        highclip = nothing,
        inspectable = theme(scene, :inspectable)
    )
end

"""
    volume(volume_data)
    volume(x, y, z, volume_data)

Plots a volume, with optional physical dimensions `x, y, z`.
Available algorithms are:
* `:iso` => IsoValue
* `:absorption` => Absorption
* `:mip` => MaximumIntensityProjection
* `:absorptionrgba` => AbsorptionRGBA
* `:additive` => AdditiveRGBA
* `:indexedabsorption` => IndexedAbsorptionRGBA
"""
@recipe(Volume, x, y, z, volume) do scene
    Attributes(;
        default_theme(scene)...,
        algorithm = :mip,
        isovalue = 0.5,
        isorange = 0.05,
        color = nothing,
        colormap = theme(scene, :colormap),
        colorrange = (0, 1),
        fxaa = true,
        inspectable = theme(scene, :inspectable)
    )
end

"""
    surface(x, y, z)

Plots a surface, where `(x, y)`  define a grid whose heights are the entries in `z`.
`x` and `y` may be `Vectors` which define a regular grid, **or** `Matrices` which define an irregular grid.

"""
@recipe(Surface, x, y, z) do scene
    Attributes(;
        default_theme(scene)...,
        color = nothing,
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        shading = true,
        fxaa = true,
        lowclip = nothing,
        highclip = nothing,
        invert_normals = false,
        inspectable = theme(scene, :inspectable)
    )
end

"""
    lines(positions)
    lines(x, y)
    lines(x, y, z)

Creates a connected line plot for each element in `(x, y, z)`, `(x, y)` or `positions`.

!!! tip
    You can separate segments by inserting `NaN`s.
"""
@recipe(Lines, positions) do scene
    Attributes(;
        default_theme(scene)...,
        linewidth = theme(scene, :linewidth),
        color = theme(scene, :linecolor),
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        linestyle = nothing,
        fxaa = false,
        cycle = [:color],
        inspectable = theme(scene, :inspectable)
    )
end

"""
    linesegments(positions)
    linesegments(x, y)
    linesegments(x, y, z)

Plots a line for each pair of points in `(x, y, z)`, `(x, y)`, or `positions`.

"""
@recipe(LineSegments, positions) do scene
    default_theme(scene, Lines)
end

# alternatively, mesh3d? Or having only mesh instead of poly + mesh and figure out 2d/3d via dispatch
"""
    mesh(x, y, z)
    mesh(mesh_object)
    mesh(x, y, z, faces)
    mesh(xyz, faces)

Plots a 3D or 2D mesh. Supported `mesh_object`s include `Mesh` types from [GeometryBasics.jl](https://github.com/JuliaGeometry/GeometryBasics.jl).

"""
@recipe(Mesh, mesh) do scene
    Attributes(;
        default_theme(scene)...,
        color = :black,
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        interpolate = false,
        shading = true,
        fxaa = true,
        inspectable = theme(scene, :inspectable),
        cycle = [:color => :patchcolor],
    )
end

"""
    scatter(positions)
    scatter(x, y)
    scatter(x, y, z)

Plots a marker for each element in `(x, y, z)`, `(x, y)`, or `positions`.

"""
@recipe(Scatter, positions) do scene
    Attributes(;
        default_theme(scene)...,
        color = theme(scene, :markercolor),
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        marker = theme(scene, :marker),
        markersize = theme(scene, :markersize),

        strokecolor = theme(scene, :markerstrokecolor),
        strokewidth = theme(scene, :markerstrokewidth),
        glowcolor = (:black, 0.0),
        glowwidth = 0.0,

        rotations = Billboard(),
        marker_offset = automatic,
        transform_marker = false, # Applies the plots transformation to marker
        distancefield = nothing,
        uv_offset_width = (0.0, 0.0, 0.0, 0.0),
        markerspace = Pixel,
        fxaa = false,
        cycle = [:color],
        inspectable = theme(scene, :inspectable)
    )
end

"""
    meshscatter(positions)
    meshscatter(x, y)
    meshscatter(x, y, z)

Plots a mesh for each element in `(x, y, z)`, `(x, y)`, or `positions` (similar to `scatter`).
`markersize` is a scaling applied to the primitive passed as `marker`.

"""
@recipe(MeshScatter, positions) do scene
    Attributes(;
        default_theme(scene)...,
        color = :black,
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        marker = :Sphere,
        markersize = 0.1,
        rotations = 0.0,
        # markerspace = relative,
        shading = true,
        fxaa = true,
        inspectable = theme(scene, :inspectable),
        cycle = [:color],
    )
end

"""
    text(string)

Plots a text.

"""
@recipe(Text, text) do scene
    Attributes(;
        default_theme(scene)...,
        color = theme(scene, :textcolor),
        font = theme(scene, :font),
        strokecolor = (:black, 0.0),
        strokewidth = 0,
        align = (:left, :bottom),
        rotation = 0.0,
        textsize = 20,
        position = (0.0, 0.0),
        justification = automatic,
        lineheight = 1.0,
        space = :screen, # or :data
        offset = (0.0, 0.0),
        inspectable = theme(scene, :inspectable)
    )
end


text:
Random.seed!(0)

#sig = [1 2; 2 2]
#sig = [1 1; 2 3]
#sig = [2 3]
#sig = [2 2]
sig = [3 2; 2 3]

S = random_S0Graph(sig)
T = complement(S)

w = random_bounded(S.n)

@time opt1 = dsw(S, w, eps=solver_eps)[1]
@time opt2 = dsw_via_complement(complement(S), w, eps=solver_eps)[1]
@test opt1 ≈ opt2  atol=tol

max_stars_repo_name: dstahlke/NoncommutativeGraphs

max_issues_repo_name: dstahlke/NoncommutativeGraphs

max_forks_repo_name: dstahlke/NoncommutativeGraphs


text:
# Define the dynamics model of the model.
struct BicycleGame{N,M,P,LF,LR,SVu,SVx,SVz,T} <: AbstractGameModel
    n::Int  # Number of states
    m::Int  # Number of controls
	p::Int  # Number of players
	ni::Vector{Int}  # Number of states for each player
	mi::Vector{Int}  # Number of controls for each player
	pu::SVu # Indices of the each player's controls
	px::SVx # Indices of the each player's x and y positions
	pz::SVz # Indices of the each player's states
	lf::T # Length from center of mass to front axle
	lr::T # Length from center of mass to rear axle
end

function BicycleGame(;p::Int=2, lf::T=0.05, lr::T=0.05) where {T}
	# p = number of players
	n = 4p
	m = 2p
	pu = [SVector{2,Int}([i + (j-1)*p for j=1:2]) for i=1:p]
	px = [SVector{2,Int}([i + (j-1)*p for j=1:2]) for i=1:p]
	pz = [SVector{4,Int}([i + (j-1)*p for j=1:4]) for i=1:p]
	TYPE = typeof.((pu,px,pz))
	ni = 4*ones(Int,p)
	mi = 2*ones(Int,p)
	return BicycleGame{n,m,p,lf,lr,TYPE...,T}(n,m,p,ni,mi,pu,px,pz,lf,lr)
end

@generated function RobotDynamics.dynamics(model::BicycleGame{N,M,P,lf,lr}, x, u) where {N,M,P,lf,lr}
	# https://archit-rstg.medium.com/two-to-four-bicycle-model-for-car-898063e87074
	# X = [x,y,v,ψ]
	# U = [a,δ]
	# β = atan(lr*tan(δ), lr+lf)
	# ̇X = [v*cos(β+ψ), v*cos(β+ψ), a, v*sin(β)/lr]

	L = :(lr+lf)
	xd  = [:(x[2P+$i]*cos(atan($lr*tan(u[P+$i]), $L) + x[3P+$i])) for i=1:P]
	yd  = [:(x[2P+$i]*sin(atan($lr*tan(u[P+$i]), $L) + x[3P+$i])) for i=1:P]
	vd  = [:(u[$i]) for i=1:P]
	ψd  = [:(x[2P+$i]*sin(atan($lr*tan(u[P+$i]), $L))/$lr) for i=1:P]
	return :(SVector{$N}($(xd...), $(yd...), $(vd...), $(ψd...)))
end

dim(model::BicycleGame) = 2


text:
# Device intrinsics, corresponding with CUDA extensions to the C language

# TODO: "CUDA C programming guide" > "C language extensions" lists mathematical functions,
#       without mentioning libdevice. Is this implied, by NVCC always using libdevice,
#       or are there some natively-supported math functions as well?

for i in ["memory_shared", "indexing", "synchronization",
          "warp_vote", "warp_shuffle", "output"]
    include(joinpath("intrinsics", "$i.jl"))
end


text:
#=
Modia main module.

* Developers: Hilding Elmqvist, Mogram AB, Toivo Henningsson, Lund and Martin Otter, DLR
* Copyright (c) 2016-2018: Hilding Elmqvist, Toivo Henningsson, Martin Otter
* License: MIT (expat)
=#

import ModiaMath

doc"""
Modia - Dynamic Modeling and Simulation in Julia

To define a model:

```julia
  @model FirstOrder begin
     x = Variable(start=1)   # start means x(0)
     T = Parameter(0.5)      # Time constant
     u = 2.0                 # Same as Parameter(2.0)
  @equations begin
     T*der(x) + x = u        # der() means time derivative
     end
  end;
```

To simulate a model:
```julia
  result = simulate(FirstOrder, 2);
  @show result["x"][end];
  ModiaMath.plot(result, "x")
```

To run examples:
```julia
  include("$(Modia.ModiaDir)/examples/runexamples.jl")
```

To run tests:
```julia
  include("$(Modia.ModiaDir)/test/runtests.jl")
```
For more information, see (https://github.com/ModiaSim/Modia.jl/blob/master/README.md)
"""
module Modia

const Version = "0.2.0-beta.1"
const Date = "2018-08-03"

#println(" \n\nWelcome to Modia - Dynamic MODeling and Simulation in julIA")
print(" \n\nWelcome to ")
print("Mod")
print_with_color(:red, "ia", bold=true)
print(" - ")
print_with_color(:light_black, "Dynamic ")
print("Mod")
print_with_color(:light_black, "eling and Simulation with Jul")
print_with_color(:red, "ia", bold=true)
println()

println("Version $Version ($Date)")
println("Type \"?Modia\" for help.\n\n")

const ModiaDir = Pkg.dir("Modia")
export ModiaDir

export @model, simulateModel, simulate, checkSimulation, simulateMultiModeModel
export Variable, Float, Boolean, Integ, Str, Parameter, Var, Par, undefined
export Variability, constant, parameter, discrete, continuous
export Property, general, symmetric, orthogonal, rotationGroup3D
export skew, skewCoords
export allInstances
export @component, addComponent!

include("language/ModiaLogging.jl")
include("language/Instantiation.jl")
include("language/Execution.jl")

include("symbolic/DAEquations/Synchronous.jl") # Before models/Electric, etc

include("models/models.jl")  # Before symbolic because MultiBody is used in BasicStructuralTransform

include("symbolic/symbolic.jl")

using .Instantiation

include("models/ModiaBase.jl")

end


text:
module TestExamples2DUnstructuredQuad

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d")

@testset "Unstructured Curve Mesh for Euler" begin
  @testset "elixir_euler_unstructured_quad.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_unstructured_quad.jl"),
      l2   = [7.844127595829963e-5, 0.00011715233141491677, 0.00011715235371677694, 0.0002587793057642201],
      linf = [0.007935754054894772, 0.011476103071539345, 0.01147611262431103, 0.025080478753547464])
  end

  @testset "elixir_euler_unstructured_quad_periodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_unstructured_quad_periodic.jl"),
      l2   = [0.00010720127811104554, 0.00012453133183340216, 0.00012453133183106435, 0.0002862523679720807],
      linf = [0.005838648336111252, 0.008981476526964016, 0.008981476526997323, 0.019987154066631874])
  end
end

end # module


text:
# This file contains some aliases and generic methods to use neural networks.

export logψ, ∇logψ, ∇logψ!, logψ_and_∇logψ, logψ_and_∇logψ!, grad_cache

"""
    grad_cache([T=out_type(net)], net::N) -> tuple

Creates a tuple holding all derivatives of the network N.
If the type is not specified, it defaults to the output type of the network.
"""
grad_cache(net::NeuralNetwork) = grad_cache(out_type(net), net)
grad_cache(T::Type{<:Number}, net::NeuralNetwork) = begin
    is_analytic(net) && return RealDerivative(T, net)
    return WirtingerDerivative(T, net)
end

# Define various aliases used a bit everywhere in networks and when testing
"""
    logψ(net, σ) -> Number

Evaluates the neural network `net` for the given input `σ`.
This will return the value ⟨σ|ψ⟩ if the network represents a pure state (a ket)
or ⟨σᵣ|ρ̂|σᵪ⟩, where σ=(σᵣ,σᵪ), if the network encodes a mixed state.

This is equivalent to `net(σ)`.

The numerical type of the output is `out_type(net)`.

If `net isa CachedNet` then the computation will be performed efficiently
with minimal allocations.
"""
function logψ end

@specialize_vararg 5 @inline logψ(net::NeuralNetwork, σ...) = net(σ...)
# Standard probability
@inline log_prob_ψ(net, σ)          = 2.0*real(net(σ))

@specialize_vararg 5 @inline ∇logψ(args...)   = logψ_and_∇logψ(args...)[2]
@specialize_vararg 5 @inline ∇logψ!(args...)  = logψ_and_∇logψ!(args...)[2]

"""
    logψ_and_∇logψ(net, σ)

Evaluates the neural network `net` for the given input `σ` and returns both the
value and it's gradient w.r.t. to the nework's parameters.

!!  The gradient will be allocated. Don't use this in a hot loop and prefer the
    in-place version `logψ_and_∇logψ!`. The gradient is always allocated by calling
    `grad_cache(net)`.

!!  If `net isa CachedNet` and the gradient has been hand-coded, the hand-coded
    version will be used, otherwise `Zygote` will be used to generate the
    gradient with AD. This case is much slower and incurs a long precompilation
    time.

See also: `logψ`, `logψ_and_∇logψ!`
"""
function logψ_and_∇logψ(net::NeuralNetwork, σ) #::Vararg{N,V})where {N,V}
    der = grad_cache(net)
    y, der = logψ_and_∇logψ!(der, net, σ)
    return y, der
end

# Autodiff version of the gradient: if you are hitting this, expect
# bad performance. Maybe you should hand-code the gradient of your model?
function logψ_and_∇logψ!(der, net::NeuralNetwork, σ)

    # Old implementation. Probably must be resuscitated one day, but not compatible
    # with our implementation of gradient caches and trainable.
    #=
    # If it is a KetNet we should not use σ...
    if isa(net, KetNeuralNetwork)
        y, back = pullback(net -> net(σ), net)
    else
        # Zygote's autodiff: generate the pullback
        y, back = pullback(net -> net(σ...), net)
    end

    # This computes the gradient, which is the conjugate of the derivative
    _der = back(Int8(1))[1]
    =#

    pars = params(net)
    fun = isa(net, KetNeuralNetwork) ? ()->net(σ) : ()->net(σ...)
    y, back = pullback(fun, pars)
    _der = back(Int8(1))

    # Convert the AD derivative to our contiguous type.
    apply_recurse!(fields(der), _der, net, conj!)
    return y, der
end

# Common Operations on networks that should always be defined
# If you hit this section it probably means that you forgot to define one
# of those functions on the type of your neural network.
"""
    input_type(net) -> Type

Returns the numerical `eltype` of the input of the network.
By default it returns `real(eltype(trainable_first(net)))`, but you should
ensure this result is correct (or specialize input_type for your network) so
that the computation is type stable
"""
input_type(net::NeuralNetwork) = real(eltype(trainable_first(net)))

"""
    out_type(net) -> Type

Returns the numerical `eltype` of the output of the network.
"""
out_type(net::NeuralNetwork) = error("Not Implemented")

"""
    out_similar(net) -> zero(T)
    out_similar(net, dims...) -> zeros(T, dims...)

Returns the result of similar for the output of the network.
This is either a scalar or a row vector.
"""
@inline out_similar(net::NeuralNetwork) = zero(out_type(net))
@inline out_similar(net::NeuralNetwork, dims::Vararg{T,N}) where {T<:Integer,N} = zeros(out_type(net), 1, dims...)

"""
    is_analytic(net) -> Bool

Returns true if the network is an R-analytic function (almost everywhere).

Please note that this is a weaker condition than holomorphicity, which is
equivalent to C-analytic.

Ref:
    [1] K. Kreuz-Delgado The Complex Gradient Operator and the CR-Calculus,
        arXiv:0906.4835
"""
is_analytic(net::NeuralNetwork) = false

"""
    num_params(net) -> Net

Returns the total number of parameters of the neural network.
"""
num_params(net::NeuralNetwork) = trainable_length(net)

# TODO does this even make sense?!
# the idea was that a shallow-copy of the weights of the net is not
# even a copy....
Base.copy(net::NeuralNetwork) = net

## For Matrix networks evaluated squared
abstract type MatrixNeuralNetwork <: NeuralNetwork end
abstract type PureNeuralNetwork <: NeuralNetwork end
abstract type KetNeuralNetwork <: NeuralNetwork end

@inline log_prob_ψ(net::MatrixNeuralNetwork,
                        σ::ADoubleStateOrBatchOrVec) = 2.0*real(net(σ))
# For Matrix network evaluated on diagonal
@inline log_prob_ψ(net::MatrixNeuralNetwork,
                        σ::AStateOrBatchOrVec) = real(net(σ))


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
function add_export_constraints(m, p; _n="")

    ##Constraint (8e): Production export and curtailment no greater than production
    @constraint(m, [t in p.techs.elec, ts in p.time_steps_with_grid],
        p.production_factor[t,ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t,ts]
        >= sum(m[Symbol("dvProductionToGrid"*_n)][t, u, ts] for u in p.export_bins_by_tech[t]) +
           m[Symbol("dvCurtail"*_n)][t, ts]
    )

    binNEM = 0
    binWHL = 0
    NEM_benefit = 0
    EXC_benefit = 0
    WHL_benefit = 0
    NEM_techs = String[t for t in p.techs.elec if :NEM in p.export_bins_by_tech[t]]
    WHL_techs = String[t for t in p.techs.elec if :WHL in p.export_bins_by_tech[t]]

    if !isempty(NEM_techs)
        # Constraint (9c): Net metering only -- can't sell more than you purchase
        # hours_per_timestep is cancelled on both sides, but used for unit consistency (convert power to energy)
        @constraint(m,
            p.hours_per_timestep * sum( m[Symbol("dvProductionToGrid"*_n)][t, :NEM, ts]
            for t in NEM_techs, ts in p.time_steps)
            <= p.hours_per_timestep * sum( m[Symbol("dvGridPurchase"*_n)][ts, tier]
                for ts in p.time_steps, tier in 1:p.s.electric_tariff.n_energy_tiers)
        )

        if p.s.electric_utility.net_metering_limit_kw == p.s.electric_utility.interconnection_limit_kw && isempty(WHL_techs)
            # no need for binNEM nor binWHL
            binNEM = 1
            @constraint(m,
                sum(m[Symbol("dvSize"*_n)][t] for t in NEM_techs) <= p.s.electric_utility.interconnection_limit_kw
            )
            NEM_benefit = @expression(m, p.pwf_e * p.hours_per_timestep *
                sum( sum(p.s.electric_tariff.export_rates[:NEM][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :NEM, ts]
                    for t in p.techs_by_exportbin[:NEM]) for ts in p.time_steps)
            )
            if :EXC in p.s.electric_tariff.export_bins
                EXC_benefit = @expression(m, p.pwf_e * p.hours_per_timestep *
                    sum( sum(p.s.electric_tariff.export_rates[:EXC][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :EXC, ts]
                        for t in p.techs_by_exportbin[:EXC]) for ts in p.time_steps)
                )
            end
        else
            if !(isempty(_n))
                @error """Binaries decisions for net metering capacity limit is not implemented for multinode models to keep
                            them linear. Please set the net metering limit to zero or equal to the interconnection limit."""
            end

            binNEM = @variable(m, binary = true)
            @warn "Adding binary variable for net metering choice. Some solvers are slow with binaries."

            # Good to bound the benefit
            max_bene = sum([ld*rate for (ld,rate) in zip(p.s.electric_load.loads_kw, p.s.electric_tariff.export_rates[:NEM])])*10
            NEM_benefit = @variable(m, lower_bound = max_bene)

            # If choosing to take advantage of NEM, must have total capacity less than net_metering_limit_kw
            @constraint(m,
                binNEM => {sum(m[Symbol("dvSize"*_n)][t] for t in NEM_techs) <= p.s.electric_utility.net_metering_limit_kw}
            )
            @constraint(m,
                !binNEM => {sum(m[Symbol("dvSize"*_n)][t] for t in NEM_techs) <= p.s.electric_utility.interconnection_limit_kw}
            )

            # binary choice for NEM benefit
            @constraint(m,
                binNEM => {NEM_benefit >= p.pwf_e * p.hours_per_timestep *
                    sum( sum(p.s.electric_tariff.export_rates[:NEM][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :NEM, ts]
                        for t in p.techs_by_exportbin[:NEM]) for ts in p.time_steps)
                }
            )
            @constraint(m, !binNEM => {NEM_benefit >= 0})

            EXC_benefit = 0
            if :EXC in p.s.electric_tariff.export_bins
                EXC_benefit = @variable(m, lower_bound = max_bene)
                @constraint(m,
                    binNEM => {EXC_benefit >= p.pwf_e * p.hours_per_timestep *
                        sum( sum(p.s.electric_tariff.export_rates[:EXC][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :EXC, ts]
                            for t in p.techs_by_exportbin[:EXC]) for ts in p.time_steps)
                    }
                )
                @constraint(m, !binNEM => {EXC_benefit >= 0})
            end
        end
    end

    if !isempty(WHL_techs)

        if typeof(binNEM) <: Real  # no need for wholesale binary
            binWHL = 1
            WHL_benefit = @expression(m, p.pwf_e * p.hours_per_timestep *
                sum( sum(p.s.electric_tariff.export_rates[:WHL][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :WHL, ts]
                        for t in p.techs_by_exportbin[:WHL]) for ts in p.time_steps)
            )
        else
            binWHL = @variable(m, binary = true)
            @warn "Adding binary variable for wholesale export choice. Some solvers are slow with binaries."
            max_bene = sum([ld*rate for (ld,rate) in zip(p.s.electric_load.loads_kw, p.s.electric_tariff.export_rates[:WHL])])*10
            WHL_benefit = @variable(m, lower_bound = max_bene)

            @constraint(m, binNEM + binWHL == 1)  # can either NEM or WHL export, not both

            @constraint(m,
                binWHL => {WHL_benefit >= p.pwf_e * p.hours_per_timestep *
                    sum( sum(p.s.electric_tariff.export_rates[:WHL][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :WHL, ts]
                            for t in p.techs_by_exportbin[:WHL]) for ts in p.time_steps)
                }
            )
            @constraint(m, !binWHL => {WHL_benefit >= 0})
        end
    end

    # register the benefits in the model
    m[Symbol("NEM_benefit"*_n)] = NEM_benefit
    m[Symbol("EXC_benefit"*_n)] = EXC_benefit
    m[Symbol("WHL_benefit"*_n)] = WHL_benefit
    nothing
end

"""
    add_monthly_peak_constraint(m, p; _n="")

Only used if ElectricTariff has monthly demand rates.
Sets dvPeakDemandMonth to greater than dvGridPurchase across each month.
If the monthly demand rate is tiered than also adds binMonthlyDemandTier and constraints.
"""
function add_monthly_peak_constraint(m, p; _n="")

	## Constraint (11d): Monthly peak demand is >= demand at each hour in the month
    if (!isempty(p.techs.chp)) && !(p.s.chp.reduces_demand_charges)
        @constraint(m, [mth in p.months, ts in p.s.electric_tariff.time_steps_monthly[mth]],
            sum(m[Symbol("dvPeakDemandMonth"*_n)][mth, t] for t in 1:p.s.electric_tariff.n_monthly_demand_tiers)
            >= sum(m[Symbol("dvGridPurchase"*_n)][ts, tier] for tier in 1:p.s.electric_tariff.n_energy_tiers) +
            sum(p.production_factor[t, ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t, ts] for t in p.techs.chp) -
            sum(sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.s.storage.types) for t in p.techs.chp) -
            sum(sum(m[Symbol("dvProductionToGrid")][t,u,ts] for u in p.export_bins_by_tech[t]) for t in p.techs.chp)

        )
    else
        @constraint(m, [mth in p.months, ts in p.s.electric_tariff.time_steps_monthly[mth]],
            sum(m[Symbol("dvPeakDemandMonth"*_n)][mth, t] for t in 1:p.s.electric_tariff.n_monthly_demand_tiers)
                >= sum(m[Symbol("dvGridPurchase"*_n)][ts, tier] for tier in 1:p.s.electric_tariff.n_energy_tiers)
        )
    end

    if p.s.electric_tariff.n_monthly_demand_tiers > 1  # only need binaries if more than one tier
        @warn "Adding binary variables to model monthly demand tiers."
        ntiers = p.s.electric_tariff.n_monthly_demand_tiers
        dv = "binMonthlyDemandTier" * _n
        m[Symbol(dv)] = @variable(m, [p.months, 1:ntiers], binary = true, base_name = dv)
        b = m[Symbol(dv)]
        # Upper bound on peak electrical power demand by month, tier; if tier is selected (0 o.w.)
        @constraint(m, [mth in p.months, tier in 1:ntiers],
            m[Symbol("dvPeakDemandMonth"*_n)][mth, tier] <= p.s.electric_tariff.monthly_demand_tier_limits[tier] *
                b[mth, tier]
        )

        # Monthly peak electrical power demand tier ordering
        @constraint(m, [mth in p.months, tier in 2:ntiers], b[mth, tier] <= b[mth, tier-1])

        # One monthly peak electrical power demand tier must be full before next one is active
        @constraint(m, [mth in p.months, tier in 2:ntiers],
        b[mth, tier] * p.s.electric_tariff.monthly_demand_tier_limits[tier-1] <=
            m[Symbol("dvPeakDemandMonth"*_n)][mth, tier-1]
        )
        # TODO implement NewMaxDemandMonthsInTier, which adds mth index to monthly_demand_tier_limits
    end
end

function add_tou_peak_constraint(m, p; _n="")
    ## Constraint (12d): Ratchet peak demand is >= demand at each hour in the ratchet`
    @constraint(m, [r in p.ratchets, ts in p.s.electric_tariff.tou_demand_ratchet_timesteps[r]],
        sum(m[Symbol("dvPeakDemandTOU"*_n)][r, tier] for tier in 1:p.s.electric_tariff.n_tou_demand_tiers) >=
        sum(m[Symbol("dvGridPurchase"*_n)][ts, tier] for tier in 1:p.s.electric_tariff.n_energy_tiers)
    )

    if p.s.electric_tariff.n_tou_demand_tiers > 1
        @warn "Adding binary variables to model TOU demand tiers."
        ntiers = p.s.electric_tariff.n_tou_demand_tiers
        dv = "binTOUDemandTier" * _n
        m[Symbol(dv)] = @variable(m, [p.ratchets, 1:ntiers], binary = true, base_name = dv)
        b = m[Symbol(dv)]

        # Upper bound on peak electrical power demand by tier, by ratchet, if tier is selected (0 o.w.)
        @constraint(m, [r in p.ratchets, tier in 1:ntiers],
            m[Symbol("dvPeakDemandTOU"*_n)][r, tier] <= p.s.electric_tariff.tou_demand_tier_limits[tier] * b[r, tier]
        )

        # Ratchet peak electrical power ratchet tier ordering
        @constraint(m, [r in p.ratchets, tier in 2:ntiers],
            b[r, tier] <= b[r, tier-1]
        )

        # One ratchet peak electrical power demand tier must be full before next one is active
        @constraint(m, [r in p.ratchets, tier in 2:ntiers],
            b[r, tier] * p.s.electric_tariff.tou_demand_tier_limits[tier-1]
            <= m[Symbol("dvPeakDemandTOU"*_n)][r, tier-1]
        )

        # Ratchet peak demand is >= demand at each hour in the ratchet
        @constraint(m, [r in p.ratchets, ts in p.s.electric_tariff.tou_demand_ratchet_timesteps[r]],
            sum(m[Symbol("dvPeakDemandTOU"*_n)][r, tier] for tier in 1:ntiers)
            >= sum( m[Symbol("dvGridPurchase"*_n)][ts, tier] for tier in 1:p.s.electric_tariff.n_energy_tiers )
        )
    end
    # TODO implement NewMaxDemandInTier
end

function add_mincharge_constraint(m, p; _n="")
    @constraint(m,
        m[Symbol("MinChargeAdder"*_n)] >= m[Symbol("TotalMinCharge"*_n)] - ( m[Symbol("TotalEnergyChargesUtil"*_n)] +
        m[Symbol("TotalDemandCharges"*_n)] + m[Symbol("TotalExportBenefit"*_n)] + m[Symbol("TotalFixedCharges"*_n)] )
    )
end

function add_simultaneous_export_import_constraint(m, p; _n="")
    @constraint(m, NoGridPurchasesBinary[ts in p.time_steps],
        m[Symbol("binNoGridPurchases"*_n)][ts] => {
          sum(m[Symbol("dvGridPurchase"*_n)][ts, tier] for tier in 1:p.s.electric_tariff.n_energy_tiers) +
          sum(m[Symbol("dvGridToStorage"*_n)][b, ts] for b in p.s.storage.types) <= 0
        }
    )
    @constraint(m, ExportOnlyAfterSiteLoadMetCon[ts in p.time_steps],
        !m[Symbol("binNoGridPurchases"*_n)][ts] => {
            sum(m[Symbol("dvProductionToGrid"*_n)][t,u,ts] for t in p.techs.elec, u in p.export_bins_by_tech[t]) <= 0
        }
    )
end

"""
    add_energy_tier_constraints(m, p; _n="")

Only necessary if n_energy_tiers > 1
"""
function add_energy_tier_constraints(m, p; _n="")
    @warn "Adding binary variables to model energy cost tiers."
    ntiers = p.s.electric_tariff.n_energy_tiers
    dv = "binEnergyTier" * _n
    m[Symbol(dv)] = @variable(m, [p.months, 1:ntiers], binary = true, base_name = dv)
    b = m[Symbol(dv)]
    ##Constraint (10a): Usage limits by pricing tier, by month
    @constraint(m, [mth in p.months, tier in 1:p.s.electric_tariff.n_energy_tiers],
        p.hours_per_timestep * sum( m[Symbol("dvGridPurchase"*_n)][ts, tier] for ts in p.s.electric_tariff.time_steps_monthly[mth] )
        <= b[mth, tier] * p.s.electric_tariff.energy_tier_limits[tier]
    )
    ##Constraint (10b): Ordering of pricing tiers
    @constraint(m, [mth in p.months, tier in 2:p.s.electric_tariff.n_energy_tiers],
        b[mth, tier] - b[mth, tier-1] <= 0
    )
    ## Constraint (10c): One tier must be full before any usage in next tier
    @constraint(m, [mth in p.months, tier in 2:p.s.electric_tariff.n_energy_tiers],
        b[mth, tier] * p.s.electric_tariff.energy_tier_limits[tier-1] -
        sum( m[Symbol("dvGridPurchase"*_n)][ts, tier-1] for ts in p.s.electric_tariff.time_steps_monthly[mth])
        <= 0
    )
    # TODO implement NewMaxUsageInTier
end

"""
    add_demand_lookback_constraints(m, p; _n="")

Only necessary if ElectricTariff.demand_lookback_percent > 0
"""
function add_demand_lookback_constraints(m, p; _n="")
    dv = "dvPeakDemandLookback" * _n
    m[Symbol(dv)] = @variable(m, [p.months], base_name = dv, lower_bound = 0)

	if p.s.electric_tariff.demand_lookback_range != 0  # then the dvPeakDemandLookback varies by month

		##Constraint (12e): dvPeakDemandLookback is the highest peak demand in DemandLookbackMonths
		for mth in p.months
			if mth > p.s.electric_tariff.demand_lookback_range
				@constraint(m, [lm in 1:p.s.electric_tariff.demand_lookback_range, ts in p.s.electric_tariff.time_steps_monthly[mth - lm]],
					m[Symbol(dv)][mth] ≥ sum( m[Symbol("dvGridPurchase"*_n)][ts, tier]
                                              for tier in 1:p.s.electric_tariff.n_energy_tiers )
				)
			else  # need to handle rollover months
				for lm in 1:p.s.electric_tariff.demand_lookback_range
					lkbkmonth = mth - lm
					if lkbkmonth ≤ 0
						lkbkmonth += 12
					end
					@constraint(m, [ts in p.s.electric_tariff.time_steps_monthly[lkbkmonth]],
						m[Symbol(dv)][mth] ≥ sum( m[Symbol("dvGridPurchase"*_n)][ts, tier]
                                                  for tier in 1:p.s.electric_tariff.n_energy_tiers )
					)
				end
			end
		end

		##Constraint (12f): Ratchet peak demand charge is bounded below by lookback
		@constraint(m, [mth in p.months],
			sum( m[Symbol("dvPeakDemandMonth"*_n)][mth, tier] for tier in 1:p.s.electric_tariff.n_monthly_demand_tiers ) >=
			p.s.electric_tariff.demand_lookback_percent * m[Symbol(dv)][mth]
		)

	else  # dvPeakDemandLookback does not vary by month

		##Constraint (12e): dvPeakDemandLookback is the highest peak demand in DemandLookbackMonths
		@constraint(m, [lm in p.s.electric_tariff.demand_lookback_months],
			m[Symbol(dv)][1] >= sum(m[Symbol("dvPeakDemandMonth"*_n)][lm, tier] for tier in 1:p.s.electric_tariff.n_monthly_demand_tiers)
		)

		##Constraint (12f): Ratchet peak demand charge is bounded below by lookback
		@constraint(m, [mth in p.months],
			sum( m[Symbol("dvPeakDemandMonth"*_n)][mth, tier] for tier in 1:p.s.electric_tariff.n_monthly_demand_tiers ) >=
			p.s.electric_tariff.demand_lookback_percent * m[Symbol(dv)][1]
		)
	end
end

function add_coincident_peak_charge_constraints(m, p; _n="")
	## Constraint (14a): in each coincident peak period, charged CP demand is the max of demand in all CP timesteps
    dv = "dvPeakDemandCP" * _n
    m[Symbol(dv)] = @variable(m, [p.s.electric_tariff.coincpeak_periods], lower_bound = 0, base_name = dv)
	@constraint(m,
        [prd in p.s.electric_tariff.coincpeak_periods,
         ts in p.s.electric_tariff.coincident_peak_load_active_timesteps[prd]],
		m[Symbol("dvPeakDemandCP"*_n)][prd] >= sum(m[Symbol("dvGridPurchase"*_n)][ts, tier]
                                                                      for tier in 1:p.s.electric_tariff.n_energy_tiers)
	)
end

function add_elec_utility_expressions(m, p; _n="")

    if !isempty(p.s.electric_tariff.export_bins) && !isempty(p.techs.all)
        # NOTE: levelization_factor is baked into dvProductionToGrid
        m[Symbol("TotalExportBenefit"*_n)] = m[Symbol("NEM_benefit"*_n)] + m[Symbol("WHL_benefit"*_n)] +
                                             m[Symbol("EXC_benefit"*_n)]
    else
        m[Symbol("TotalExportBenefit"*_n)] = 0
    end

    m[Symbol("TotalEnergyChargesUtil"*_n)] = @expression(m, p.pwf_e * p.hours_per_timestep *
        sum( p.s.electric_tariff.energy_rates[ts, tier] * m[Symbol("dvGridPurchase"*_n)][ts, tier]
            for ts in p.time_steps, tier in 1:p.s.electric_tariff.n_energy_tiers)
    )

    if !isempty(p.s.electric_tariff.tou_demand_rates)
        m[Symbol("DemandTOUCharges"*_n)] = @expression(m,
            p.pwf_e * sum( p.s.electric_tariff.tou_demand_rates[r] * m[Symbol("dvPeakDemandTOU"*_n)][r, tier]
            for r in p.ratchets, tier in p.s.electric_tariff.n_tou_demand_tiers)
        )
    else
        m[Symbol("DemandTOUCharges"*_n)] = 0
    end

    if !isempty(p.s.electric_tariff.monthly_demand_rates)
        m[Symbol("DemandFlatCharges"*_n)] = @expression(m, p.pwf_e *
            sum( p.s.electric_tariff.monthly_demand_rates[mth, t] * m[Symbol("dvPeakDemandMonth"*_n)][mth, t]
                for mth in p.months, t in 1:p.s.electric_tariff.n_monthly_demand_tiers)
        )
    else
        m[Symbol("DemandFlatCharges"*_n)] = 0
    end

    m[Symbol("TotalDemandCharges"*_n)] = m[Symbol("DemandTOUCharges"*_n)] + m[Symbol("DemandFlatCharges"*_n)]

    m[Symbol("TotalFixedCharges"*_n)] = p.pwf_e * p.s.electric_tariff.fixed_monthly_charge * 12

    if p.s.electric_tariff.annual_min_charge > 12 * p.s.electric_tariff.min_monthly_charge
        m[Symbol("TotalMinCharge"*_n)] = p.s.electric_tariff.annual_min_charge
    else
        m[Symbol("TotalMinCharge"*_n)] = 12 * p.s.electric_tariff.min_monthly_charge
    end

	if m[Symbol("TotalMinCharge"*_n)] >= 1e-2
		add_mincharge_constraint(m, p)
	else
		@constraint(m, m[Symbol("MinChargeAdder"*_n)] == 0)
	end

    if !isempty(p.s.electric_tariff.coincpeak_periods)
        m[Symbol("TotalCPCharges"*_n)] = @expression(m, p.pwf_e *
            sum( p.s.electric_tariff.coincident_peak_load_charge_per_kw[prd] * m[Symbol("dvPeakDemandCP"*_n)][prd]
                for prd in p.s.electric_tariff.coincpeak_periods )
        )
    else
        m[Symbol("TotalCPCharges"*_n)] = 0
    end

    m[Symbol("TotalElecBill"*_n)] = (
        m[Symbol("TotalEnergyChargesUtil"*_n)]
        + m[Symbol("TotalDemandCharges"*_n)]
        + m[Symbol("TotalExportBenefit"*_n)]
        + m[Symbol("TotalFixedCharges"*_n)]
        + m[Symbol("TotalCPCharges"*_n)]
        + 0.999 * m[Symbol("MinChargeAdder"*_n)]
    )
    #= Note: 0.999 * MinChargeAdder in Objective b/c when
        TotalMinCharge > (TotalEnergyCharges + TotalDemandCharges + TotalExportBenefit + TotalFixedCharges)
		it is arbitrary where the min charge ends up (eg. could be in TotalDemandCharges or MinChargeAdder).
		0.001 * MinChargeAdder is added back into LCC when writing to results.
    =#
    nothing
end


text:
import CSV
using Decimals

@enum INGDirection direction_in=1 direction_out=2

struct INGTransaction
    date::Date
    countername::String
    account::String
    counteraccount::String
    code::String
    direction::INGDirection
    amount::Decimal
    type::String
    description::String
end

function load_ing_transactions_from_file(filePath::String)::Array{INGTransaction}
    csvrows = CSV.File(
      filePath,
      skipto=2,
      quotechar='"',
      escapechar='"',
      delim=",",
      header=["date", "countername", "account", "counteraccount",
              "code", "direction", "amount", "type",
              "description"],
      types=Dict("date" => String, "countername" => String, "account" => String,
                 "counteraccount" => String, "direction" => String, "amount" => String,
                 "type" => String, "description" => String))
    transactions = map(row -> INGTransaction(
        Date(parse(Int, row.date[1:4]), parse(Int, row.date[5:6]), parse(Int, row.date[7:8])),
        row.countername,
        row.account,
        coalesce(row.counteraccount, ""),
        row.code,
        if (row.direction == "Debit") direction_out else direction_in end,
        parse(Decimal, replace(row.amount, "," => ".")) * if (row.direction == "Debit") parse(Decimal, "-1.0") else parse(Decimal,"1.0") end,
        row.type,
        row.description
    ), csvrows)

    return if first(transactions).date > last(transactions).date
        reverse(transactions)
    else
        transactions
    end
end

function categorise(tr::INGTransaction, config::INGConfig)::String
    for matchrule in config.categorymatchrules
        if (typeof(match(Regex(matchrule.matchers.countername), tr.countername)) != Nothing
            && typeof(match(Regex(matchrule.matchers.counteraccount), tr.counteraccount)) != Nothing
            && typeof(match(Regex(matchrule.matchers.description), tr.description)) != Nothing)
            return matchrule.category
        end
    end
    error("Couldn't match transaction " * string(tr))
end

function convert_ing_transactions(transactions::Array{INGTransaction}, config::INGConfig)::Array{Transaction}
    balance::Dict{Account, Float64} = Dict(
        account_main => 0.0,
        account_net => 0.0
    )
    for accountname in config.accountnames
        balance[Account(accountname)] = 0.0
    end
    new_transactions::Array{Transaction} = []
    for tr in transactions
        category = categorise(tr, config)
        change::Dict{Account, Float64} = Dict(
            account_main => Float64(tr.amount)
        )
        if category in config.accountnames
            change[Account(category)] = Float64(-tr.amount)
        end
        if category in config.accountnames && !(category in config.netignoreaccounts)
            change[account_net] = 0.0
        else
            change[account_net] = change[account_main]
        end

        new_balance = Dict()
        for account in keys(balance)
            new_balance[account] = balance[account]
        end
        for account in keys(change)
            new_balance[account] = get(new_balance, account, 0) + change[account]
        end
        push!(new_transactions, Transaction(
            tr.date,
            category,
            tr.countername,
            tr.description,
            change,
            new_balance
           ))
        balance = new_balance
    end
    return new_transactions
end

max_stars_repo_name: CrazedProgrammer/BalancePlot

max_issues_repo_name: CrazedProgrammer/BalancePlot

max_forks_repo_name: CrazedProgrammer/BalancePlot


text:
#!/usr/bin/env julia

using EwaldSum, Test, SafeTestsets

@time begin
@time @safetestset "Madelung" begin include("test_vmad.jl") end
end


text:
function search(tkn, query; limit = 10, offset = 0, type = "all")
    @assert type ∈  ["all", "artists", "tracks", "playlists", "curators", "albums", "stations", "cities"]
    ret = getparsed(tkn, ["search"], ["limit=$limit", "offset=$offset", "type=$type"])
    if haskey(ret, "obj")
        return ret["obj"]
    else
        return ret
    end
end

function albumrequest(token::Token, variable, id, maxtries = 5;
        type=nothing, platform = nothing, status = nothing,
        parameters::Vector{String})
    allowedvars = ["charts", "get-ids", "playlists", "stats", "tracks", "tunefind", "metadata"]
    variable = lowercase(variable)
    variable ∈ allowedvars || @warn "Variable is unknown. Trying to parse..."
    path = ["album"]
    id = string(id)
    if variable == "charts"
        required =  [id, type, variable]
        any(x -> x === nothing, required) && @error string("Missing required argument: type")
        append!(path, required)
    elseif variable === "get-ids"
        required = [type, id, variable]
        any(x -> x === nothing, required) && @error string("Missing required argument: type")
        append!(path, required)
    elseif variable === "metadata"
        push!(path, id)
    elseif variable == "playlists"
        required = [id, platform, status, variable]
        any(x -> x === nothing, required) && @error string("Missing required argument(s): platform, status")
        append!(path, required)
    elseif variable == "stats"
        required = [id, platform, variable]
        any(x -> x === nothing, required) && @error string("Missing required argument: platform")
        append!(path, required)
    elseif variable ∈ ["tracks", "tunefind"]
        append!(path,  [id, variable])
    else
        vars_all = [id, type, platform, status, variable]
        vars_provided = vars_all[vars_all .!== nothing]
        append!(path, vars_provided)
        println(string("Parsed URL: ", buildrequesturl(path)))
    end
    url_req = buildrequesturl(path)
    req = Request(token, url_req, maxtries)
    setparameters!(req, parameters)
    return req
end

function curatorrequest(token::Token, variable, maxtries = 5;
            platform = nothing, id = nothing, parameters::Vector{String})
            # TODO input checks
            @assert platform !== nothing
            path = ["curator"]
            if variable == "metadata"
                @assert id !== nothing
                append!(path, [platform, id])
            elseif variable == "lists"
                append!(path, [platform, variable])
            elseif variable == "playlists"
                @assert id !== nothing
                append!(path, [platform, id, variable])
            else
                vars_all = [platform, id, variable]
                vars_provided = vars_all[vars_all .!== nothing]
                append!(path, vars_provided)
                println(string("Parsed URL: ", buildrequesturl(path)))
            end
            url_req = buildrequesturl(path)
            req = Request(token, url_req, maxtries)
            setparameters!(req, parameters)
            return req
end

function playlistrequest(token::Token, variable, maxtries = 5;
            id = nothing, platform = nothing, type = nothing, span = nothing,
            parameters::Vector{String} = [""])
            allowedvars = ["metadata", "playlist-evolution", "journey-progression", "lists", "snapshot", "stats", "tracks"]
            variable = lowercase(variable)
            variable ∈ allowedvars || @warn "Variable is unknown. Trying to parse..."
            path = ["playlist"]
            if variable == "metadata"
                @assert platform !== nothing && id !== nothing
                id = string(id)
                append!(path, [platform, id])
            elseif variable == "playlist-evolution"
                @assert type !== nothing && id !== nothing
                id = string(id)
                append!(path, ["by", type, id, variable])
            elseif variable == "journey-progression"
                @assert platform !== nothing && id !== nothing && type !== nothing
                id = string(id)
                append!(path, [platform, id, variable, type])
            elseif variable == "lists"
                @assert platform !== nothing
                append!(path, [platform, variable])
            elseif variable == "snapshot"
                @assert platform !== nothing && id !== nothing
                id = string(id)
                append!(path, [platform, id, variable])
            elseif variable == "stats"
                @assert platform !== nothing && id !== nothing
                id = string(id)
                append!(path, [platform, id, variable])
            elseif variable == "tracks"
                @assert platform !== nothing && id !== nothing && span !== nothing
                id = string(id)
                append!(path, [platform, id, span, variable])
            else
                id !== nothing && (id = string(id))
                vars_all = [platform, id, span, variable, type]
                vars_provided = vars_all[vars_all .!== nothing]
                append!(path, vars_provided)
                println(string("Parsed URL: ", buildrequesturl(path)))
            end
            url_req = buildrequesturl(path)
            req = Request(token, url_req, maxtries)
            setparameters!(req, parameters)
            return req
end


text:
module 03_example_par_study_three_columns

function independent_LSM_ratios()
  reutrn ImageToEIS.par_study(
          Dict(
              "matrix_template" => three_column_domain_matrix,
              #
              "trials_count" => 1,
              #
              #
              "LSM_ratio1" => 0.5,
              "LSM_ratio2" => 0.5,
              "LSM_ratio3" => 0.5,
              #
              "porosity1" => 0.5,
              "porosity2" => 0.5,
              "porosity3" => 0.5,
              #
              "positions_of_contacts" => (10, 20),
              "height_of_contacts" => 5,
              #
              "column_width" => 6,
              "height" => 30
          ),
          save_to_file = ""
  )
end

function LSM_ratios_as_tuple()
  return ImageToEIS.par_study(
                 Dict(
                     "matrix_template" => ImageToEIS.three_column_domain_LSM_ratios,
                     #
                     "trials_count" => 20,
                     #
                     #
                     "LSM_ratios" => [
                                       (0.5, 0.5, 0.5),
                                       (0.0, 0.5, 0.5),
                                       (0.0, 0.0, 0.5)
                                       ],
                     #
                     "porosity1" => 0.5,
                     "porosity2" => 0.5,
                     "porosity3" => 0.5,
                     #
                     "positions_of_contacts" => (10, 30),
                     "height_of_contacts" => 5,
                     #
                     "column_width" => 5,
                     "height" => 60
                 ),
                 save_to_file = "jojo"
       )
end

end # of module

max_stars_repo_name: Masicko/ImageToEIS

max_issues_repo_name: Masicko/ImageToEIS

max_forks_repo_name: Masicko/ImageToEIS


text:
"""
    lu_decompose(mat)
Decomposes a `n x n` non singular matrix into a lower triangular matrix (L) and an upper triangular matrix (U)
"""
function lu_decompose(mat)
    n = first(size(mat))
    L = zeros(n, n)
    U = zeros(n, n)

    for i in 1:n
        for j in i:n
            s = 0
            for k in 1:i
                s += L[i, k] * U[k, j]
            end
            U[i, j] = mat[i, j] - s
        end

        for k in i:n
            if i == k
                L[i, i] = 1
            else
                s = 0
                for j in 1:i
                    s += L[k, j] * U[j, i]
                end
                L[k, i] = (mat[k, i] - s) / U[i, i]
            end
        end
    end

    return L, U
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
using LinearAlgebra
using StaticArrays
using Distributions
using PyPlot
const Vec2f = SVector{2, Float64}
const Vec4f = SVector{4, Float64}
const VecFrame = SVector{18, Float64}

include("read_scenario.jl")

@inline function distance(p::Vec2f, q::Vec2f)
    return sqrt((p[1]-q[1])^2+(p[2]-q[2])^2)
end

function predict_vehicle(scenario::Vector{VecFrame}, center::Vec2f, vehicle_id::Int64, frame_id::Int64, start_frame)
    frame = get_frame(scenario, frame_id, start_frame)
    n = length(frame)
    next_lane_id = 0
    for i = 1:n
        if convert(Int64, frame[i][1]) == vehicle_id
            if distance(Vec2f(frame[i][5], frame[i][6]), center) < 150
                next_lane_id = frame[i][14]
                break
            end
        end
    end
    return next_lane_id
end

function cal_time(p::Vec4f, q::Vec4f)
    v = (sqrt(p[3]^2 + p[4]^2) + sqrt(q[3]^2 + q[4]^2)) / 2
    d = distance(Vec2f(p[1], p[2]), Vec2f(q[1], q[2]))
    return d/v
end

function cal_time(waypoints::Vector{Vec4f})
    n = length(waypoints)
    cum_time = [0.0]
    t = 0.0
    for i = 1:(n-1)
        t += cal_time(waypoints[i], waypoints[i+1])
        push!(cum_time, t)
    end
    return cum_time
end

function get_next_position(cum_time::Vector{Float64}, waypoints::Vector{Vec4f}, next_time_stamp::Float64)
    n = length(waypoints)
    t = 0
    for i = 1:(n-1)
        if cum_time[i] < next_time_stamp < cum_time[i+1]
            t = i + 1
            break
        end
    end
    return waypoints[1:t]
end

function get_next_position(cum_time::Vector{Float64}, waypoints::Vector{Vec4f}, next_time_stamp::Float64)
    n = length(waypoints)
    t = 0
    for i = 1:(n-1)
        if cum_time[i] < next_time_stamp < cum_time[i+1]
            t = i
            break
        end
    end
    println("the final time is ", cum_time[t])
    return waypoints[1:t], waypoints[t:end]
end

function get_next_position(cum_time::Vector{Float64}, waypoints::Vector{Vec4f}, next_time_stamp::Float64, s_goal::Vec4f)
    n = length(waypoints)
    t = 0
    reach = false
    println(length(waypoints))
    println(length(cum_time))
    for i = 1:(n-1)
        if cum_time[i] < next_time_stamp
            if distance(Vec2f(waypoints[i][1], waypoints[i][2]), Vec2f(s_goal[1], s_goal[2])) < 10
                t = i
                reach = true
                println("the final time is ", cum_time[t])
                break
            end
        end
        if cum_time[i] < next_time_stamp < cum_time[i+1]
            t = i
            break
        end
    end
    print(t)
    return waypoints[1:t], waypoints[t:end], reach
end

function get_traj(cum_time::Vector{Float64}, waypoints::Vector{Vec4f}, cur_time_stamp::Float64, next_time_stamp::Float64)
    n = length(waypoints)
    t = 0
    next_t = 0
    if cur_time_stamp - 0 < 1e-5
        t = 1
    end
    for i = 1:(n-1)
        if cum_time[i] < cur_time_stamp < cum_time[i+1]
            t = i
        end
        if cum_time[i] < next_time_stamp < cum_time[i+1]
            next_t = i
        end
    end
    return waypoints[t:next_t]
end

function check_traj(cum_time::Vector{Float64}, waypoints, time, frame_id, s_goal, scenario, start_frame)
    lane12 = (lane1_c + lane2_c)/2
    lane23 = (lane2_c + lane3_c)/2
    lane34 = (lane3_c + lane4_c)/2
    lane45 = (lane4_c + lane5_c)/2
    lane5edge = lane5_c + (lane5_c - lane45)
    traj = get_traj(cum_time, waypoints, time-0.1, time)
    # println(traj)
    s_init = traj[1]
    frame1 = get_frame(scenario, frame_id, start_frame)

    center = Vec2f(s_init[1], s_init[2])
    if 0 < s_init[1] < lane12
        lane_id = 1
    elseif lane12 < s_init[1] < lane23
        lane_id = 2
    elseif lane23 < s_init[1] < lane34
        lane_id = 3
    elseif lane34 < s_init[1] < lane45
        lane_id = 4
    elseif lane45 < s_init[1] < lane5edge
        lane_id = 5
    end
    neighbors1, front_neighbors1 = find_neighbors(frame1, center, lane_id, 150.0)
    obstacle_set = set_obstacles(neighbors1)
    x_min, x_max, v_min, v_max = get_boundary(neighbors1, center, lane_id, s_goal)
    W = World(x_min, x_max, v_min, v_max, obstacle_set)
    if isValid(W, traj)
        return true
    else
        return false
    end
end

function softmax(l1::Float64, l2::Float64, l3::Float64)
    sum = exp(l1) + exp(l2) + exp(l3)
    return convert(Float64, exp(l1)/sum), convert(Float64, exp(l2)/sum), convert(Float64, exp(l3)/sum)
end

function softmax(l1::Float64, l2::Float64)
    sum = exp(l1) + exp(l2)
    return convert(Float64, exp(l1)/sum), convert(Float64, exp(l2)/sum)
end


text:
# Subspace matrix describing the coupling between partial waves k1 and
# k2.
type SubOperator
    m::AbstractMatrix
    k1::Integer
    k2::Integer
end
type Operator
    o::Vector{SubOperator}
end

export SubOperator, Operator


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct MeterInfo <: SwaggerModel
    MeterId::Any # spec type: Union{ Nothing, String } # spec name: MeterId
    MeterName::Any # spec type: Union{ Nothing, String } # spec name: MeterName
    MeterCategory::Any # spec type: Union{ Nothing, String } # spec name: MeterCategory
    MeterSubCategory::Any # spec type: Union{ Nothing, String } # spec name: MeterSubCategory
    Unit::Any # spec type: Union{ Nothing, String } # spec name: Unit
    MeterTags::Any # spec type: Union{ Nothing, Vector{String} } # spec name: MeterTags
    MeterRegion::Any # spec type: Union{ Nothing, String } # spec name: MeterRegion
    MeterRates::Any # spec type: Union{ Nothing, Dict{String, Float32} } # spec name: MeterRates
    EffectiveDate::Any # spec type: Union{ Nothing, DateTime } # spec name: EffectiveDate
    IncludedQuantity::Any # spec type: Union{ Nothing, Float32 } # spec name: IncludedQuantity

    function MeterInfo(;MeterId=nothing, MeterName=nothing, MeterCategory=nothing, MeterSubCategory=nothing, Unit=nothing, MeterTags=nothing, MeterRegion=nothing, MeterRates=nothing, EffectiveDate=nothing, IncludedQuantity=nothing)
        o = new()
        validate_property(MeterInfo, Symbol("MeterId"), MeterId)
        setfield!(o, Symbol("MeterId"), MeterId)
        validate_property(MeterInfo, Symbol("MeterName"), MeterName)
        setfield!(o, Symbol("MeterName"), MeterName)
        validate_property(MeterInfo, Symbol("MeterCategory"), MeterCategory)
        setfield!(o, Symbol("MeterCategory"), MeterCategory)
        validate_property(MeterInfo, Symbol("MeterSubCategory"), MeterSubCategory)
        setfield!(o, Symbol("MeterSubCategory"), MeterSubCategory)
        validate_property(MeterInfo, Symbol("Unit"), Unit)
        setfield!(o, Symbol("Unit"), Unit)
        validate_property(MeterInfo, Symbol("MeterTags"), MeterTags)
        setfield!(o, Symbol("MeterTags"), MeterTags)
        validate_property(MeterInfo, Symbol("MeterRegion"), MeterRegion)
        setfield!(o, Symbol("MeterRegion"), MeterRegion)
        validate_property(MeterInfo, Symbol("MeterRates"), MeterRates)
        setfield!(o, Symbol("MeterRates"), MeterRates)
        validate_property(MeterInfo, Symbol("EffectiveDate"), EffectiveDate)
        setfield!(o, Symbol("EffectiveDate"), EffectiveDate)
        validate_property(MeterInfo, Symbol("IncludedQuantity"), IncludedQuantity)
        setfield!(o, Symbol("IncludedQuantity"), IncludedQuantity)
        o
    end
end # type MeterInfo

const _property_map_MeterInfo = Dict{Symbol,Symbol}(Symbol("MeterId")=>Symbol("MeterId"), Symbol("MeterName")=>Symbol("MeterName"), Symbol("MeterCategory")=>Symbol("MeterCategory"), Symbol("MeterSubCategory")=>Symbol("MeterSubCategory"), Symbol("Unit")=>Symbol("Unit"), Symbol("MeterTags")=>Symbol("MeterTags"), Symbol("MeterRegion")=>Symbol("MeterRegion"), Symbol("MeterRates")=>Symbol("MeterRates"), Symbol("EffectiveDate")=>Symbol("EffectiveDate"), Symbol("IncludedQuantity")=>Symbol("IncludedQuantity"))
const _property_types_MeterInfo = Dict{Symbol,String}(Symbol("MeterId")=>"String", Symbol("MeterName")=>"String", Symbol("MeterCategory")=>"String", Symbol("MeterSubCategory")=>"String", Symbol("Unit")=>"String", Symbol("MeterTags")=>"Vector{String}", Symbol("MeterRegion")=>"String", Symbol("MeterRates")=>"Dict{String, Float32}", Symbol("EffectiveDate")=>"DateTime", Symbol("IncludedQuantity")=>"Float32")
Base.propertynames(::Type{ MeterInfo }) = collect(keys(_property_map_MeterInfo))
Swagger.property_type(::Type{ MeterInfo }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_MeterInfo[name]))}
Swagger.field_name(::Type{ MeterInfo }, property_name::Symbol) =  _property_map_MeterInfo[property_name]

function check_required(o::MeterInfo)
    true
end

function validate_property(::Type{ MeterInfo }, name::Symbol, val)
end


text:
using Test
using CellMLToolkit
using DifferentialEquations, Plots

path = @__DIR__
ml = CellModel(path * "/../models/beeler_reuter_1977.cellml.xml")

# @test length(ml.eqs) == 8
# @test ml.iv.op.name == :time

# eqs, vs = CellMLToolkit.flat_equations(ml)
# @test length(vs) == 8

# @test find_V(ml).op.name == :V

prob = ODEProblem(ml, (0,10000.0))
sol1 = solve(prob, Euler(), dt=0.01, saveat=1.0)
sol2 = solve(prob, TRBDF2(), dtmax=0.5, saveat=1.0)
V1 = map(x -> x[1], sol1.u)
V2 = map(x -> x[1], sol2.u)
err = sum(abs.(V1 .- V2)) / length(V1)
@test err < 0.1

# prob = ODEProblem(ml, (0,10000.0); jac=true)
# sol3 = solve(prob, TRBDF2(), dtmax=0.5, saveat=1.0)
# V3 = map(x -> x[1], sol2.u)
# err = sum(abs.(V1 .- V3)) / length(V1)
# @test err < 0.1

p = list_params(ml)
update_list!(p, "IstimPeriod", 280.0)
prob = ODEProblem(ml, (0,10000.0); jac=false, p=p)
sol4 = solve(prob, TRBDF2(), dtmax=0.5, saveat=1.0)
V4 = map(x -> x[1], sol2.u)
err = sum(abs.(V1[1:250] .- V4[1:250])) / 250
@test err < 0.1

##############################################################################

function test_plot()
    ml = CellModel(path * "/../models/beeler_reuter_1977.cellml.xml")
    tspan = (0, 5000.0)
    prob = ODEProblem(ml, tspan)
    sol = solve(prob, TRBDF2(), dtmax=0.5)
    x = Array(sol)
    plot(x[1,:])
end


text:
#=
    TODO: how to address good ideal?

=#

"""
    An interface to a field of rational functions and a main package structure

    A convenient way to construct a field is to provide a set of generators

    TODO
"""
mutable struct RationalFunctionField
    generating_set       #
    groebner_set         # lazy field, not to be addressed directly
end

function RationalFunctionField(genset)
    RationalFunctionField(genset, [])
end

"""
    Returns generators of the given rational function field

    If the relative groebner ideal was already computed,
    groebner coefficients are returned
"""
function AbstractAlgebra.gens(FF::RationalFunctionField)
    if isempty(FF.groebner_set)
        return FF.generating_set
    else
        return obtain_generators(FF.groebner_set)
    end
end

"""
    Assesses the membership of the given `elem` in the provided field `FF`
    with the given probability of a correct answer `p`

    Correctness of this algorithm relies on *Theorem 3.3* from
        *Structural identifiability via input-output projections, Pogudin et alia*
"""
function contains_randomized(
            FF::RationalFunctionField, elem;
            p=0.99)

    # Theorem 3.3 from
    # Structural identifiability via input-output projections, Pogudin et alia

    # choose appropriate to the given 'p' random upper bound
    n = length(gens(parent(numerator(elem))))
    d = maximum([
            maximum([degrees(numerator(f))..., degrees(denominator(f))...])
            for f in FF.generating_set
    ])
    M = round(Int, 6d^(n+3) / (1 - p))

    @info "Random values bound M = $M"

    # idealize the generators and the elem in the same rings
    I = idealize(FF.generating_set)
    i = idealize(elem, I)

    # evaluate at a random point
    G = GroebnerEvaluator(I)
    g = GroebnerEvaluator(i)

    p = generate_point(G, M=M)

    gb = evaluate(G, p)
    elem = evaluate(g, p)

    @debug "Evaluated eLement is $elem"
    @debug "Evaluated groebner basis is $gb"

    Is = Singular.Ideal(I.evalring, gb)

    # check ideal membership
    return iszero(Singular.reduce(elem..., GroebnerBasis.f4(Is)))
end

"""
    Assesses the membership of the given `elem` in the provided field `FF`

    If the groebner ideal corresponding to FF generators is not computed yet,
    this might take some time
"""
function contains_using_groebner(FF::RationalFunctionField, elem)

    if isempty(FF.groebner_set)
        compute_groebner!(FF)
    end

    # TODO: think of something more elegant
    Is, yoverx, basepolyring, yoverxs, basepolyrings = groebner_ideal_to_singular(FF.groebner_set)

    f = idealize_element(elem, basepolyring, basepolyrings, yoverx, yoverxs)

    fs = Ideal(yoverxs, f)
    Is = Ideal(yoverxs, Is)
    Singular.contains(Is, fs)
end

function compute_groebner!(FF::RationalFunctionField)
    FF.groebner_set = groebner_ideal(gens(FF))
end

function simple_preprocess(generators)
    tmp = filter((c) -> !isconstant(c), generators)
    sort!(tmp, by=(c) -> sum(degrees(c)), rev=true)
    ans = [ first(tmp) ]
    for elem in tmp
        field = RationalFunctionField(ans)
        compute_groebner!(field)
        if !contains(field, elem, p=1)
            push!(ans, elem)
        end
    end
    ans = [ x * inv(Nemo.lc(numerator(x))) for x in ans ]
    ans
end

"""
    Simplifies generators of the given RationalFunctionField
    by computing the Groebner Basis of the corresponding good ideal
    and applying some heuristics to its coefficients after
"""
function simplify!(
            FF::RationalFunctionField;
            backend_algorithm=simple_preprocess)
    FF.generating_set = backend_algorithm(copy(gens(FF)))
end

"""
    Assesses the field membership of `elem` in `FF` by checking
    if `elem` lies in the "good" ideal generated by the given field generators

    `p` - probability of a correct answer (both positive and negative)

    If p is set to 1, deterministic algorithm is used

"""
function Base.contains(FF::RationalFunctionField, elem; p=0.99)
    if isone(p)
        contains_using_groebner(FF, elem)
    else
        contains_randomized(FF, elem, p=p)
    end
end


text:
using PSCOPF

using Test
using Dates
using DataStructures

@testset verbose=true "test_starts" begin

    TS = [DateTime("2015-01-01T11:00:00"), DateTime("2015-01-01T11:15:00"), DateTime("2015-01-01T11:30:00"), DateTime("2015-01-01T11:45:00")]

    initial_state = SortedDict{String, PSCOPF.GeneratorState}(
            "gen1" => PSCOPF.ON,
            "gen2" => PSCOPF.OFF,
            "gen3" => PSCOPF.OFF,
        )

    commitments = SortedDict{Tuple{String, Dates.DateTime}, PSCOPF.GeneratorState}(
        ("gen1", TS[1]) => PSCOPF.ON,
        ("gen1", TS[2]) => PSCOPF.ON,
        ("gen1", TS[3]) => PSCOPF.OFF,
        ("gen1", TS[4]) => PSCOPF.ON,#start

        ("gen2", TS[1]) => PSCOPF.ON,#start
        ("gen2", TS[2]) => PSCOPF.ON,
        ("gen2", TS[3]) => PSCOPF.ON,
        ("gen2", TS[4]) => PSCOPF.ON,

        ("gen3", TS[1]) => PSCOPF.OFF,
        ("gen3", TS[2]) => PSCOPF.ON,#start
        ("gen3", TS[3]) => PSCOPF.OFF,
        ("gen3", TS[4]) => PSCOPF.ON,#start
    )

    expected = Set{Tuple{String,Dates.DateTime}}([
            ("gen1", TS[4]), ("gen2", TS[1]), ("gen3", TS[2]), ("gen3", TS[4])
        ])

    @testset "test_starts_testcase" begin

        @testset "test_starts_simple" begin
            result = PSCOPF.get_starts(commitments, initial_state)
            @test result == expected
        end

        @testset "test_starts_from_schedule" begin
            schedule = PSCOPF.Schedule(PSCOPF.Utilitary(), DateTime("2015-01-01T06:00:00"), SortedDict(
                    "gen1" => PSCOPF.GeneratorSchedule("gen1",
                        SortedDict(TS[1]=> PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.ON, SortedDict()),
                                    TS[2] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.ON, SortedDict()),
                                    TS[3] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.OFF, SortedDict()),
                                    TS[4] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.ON, SortedDict()),
                                    ),
                        SortedDict(),
                        ),
                    "gen2" => PSCOPF.GeneratorSchedule("gen2",
                        SortedDict(TS[1]=> PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.ON, SortedDict()),
                                    TS[2] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.ON, SortedDict()),
                                    TS[3] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.ON, SortedDict()),
                                    TS[4] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.ON, SortedDict()),
                                    ),
                        SortedDict(),
                        ),
                    "gen3" => PSCOPF.GeneratorSchedule("gen3",
                        SortedDict(TS[1]=> PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.OFF, SortedDict()),
                                    TS[2] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.ON, SortedDict()),
                                    TS[3] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.OFF, SortedDict()),
                                    TS[4] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.ON, SortedDict()),
                                    ),
                        SortedDict(),
                        )
                    )
            )

            result = PSCOPF.get_starts(schedule, initial_state)
            @test result == expected
        end

    end

    @testset "test_starts_if_different_timesteps_in_comitments" begin
        initial_state = SortedDict{String, PSCOPF.GeneratorState}(
            "gen1" => PSCOPF.ON,
            "gen2" => PSCOPF.OFF,
            "gen3" => PSCOPF.OFF,
        )

        commitments = SortedDict{Tuple{String, Dates.DateTime}, PSCOPF.GeneratorState}(
            ("gen1", TS[1]) => PSCOPF.ON,
            ("gen1", TS[2]) => PSCOPF.ON,
            ("gen1", TS[3]) => PSCOPF.OFF,
            ("gen1", TS[4]) => PSCOPF.ON,#start

            ("gen2", Hour(1) + TS[1]) => PSCOPF.ON,#start
            ("gen2", Hour(1) + TS[2]) => PSCOPF.ON,
            ("gen2", Hour(1) + TS[3]) => PSCOPF.ON,
            ("gen2", Hour(1) + TS[4]) => PSCOPF.ON,

            ("gen3", TS[1]) => PSCOPF.OFF,
            ("gen3", TS[2]) => PSCOPF.ON,#start
            ("gen3", TS[3]) => PSCOPF.OFF,
            ("gen3", TS[4]) => PSCOPF.ON,#start
        )

        expected = Set{Tuple{String,Dates.DateTime}}([
            ("gen1", TS[4]), ("gen2", Hour(1) + TS[1]), ("gen3", TS[2]), ("gen3", TS[4])
        ])

        result = PSCOPF.get_starts(commitments, initial_state)
        @test result == expected
    end

    @testset "test_starts_from_schedule_ignores_non_definitive_starts" begin
        initial_state = SortedDict{String, PSCOPF.GeneratorState}(
            "unit_1" => PSCOPF.OFF,
        )

        schedule = PSCOPF.Schedule(PSCOPF.Utilitary(), DateTime("2015-01-01T06:00:00"), SortedDict(
                "unit_1" => PSCOPF.GeneratorSchedule("unit_1",
                    SortedDict(TS[1]=> PSCOPF.UncertainValue{PSCOPF.GeneratorState}(missing, SortedDict("S1"=>PSCOPF.ON)), # non definitive start
                                TS[2] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(missing, SortedDict("S1"=>PSCOPF.ON)),
                                TS[3] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(missing, SortedDict("S1"=>PSCOPF.ON)),
                                TS[4] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(missing, SortedDict("S1"=>PSCOPF.ON)),
                                ),
                    SortedDict(),
                    )
                )
        )

        result = PSCOPF.get_starts(schedule, initial_state)
        @test isempty(result)
    end

    @testset "test_starts_from_schedule_breaks_at_first_non_definitive_value" begin
        initial_state = SortedDict{String, PSCOPF.GeneratorState}(
            "unit_1" => PSCOPF.OFF,
        )

        schedule = PSCOPF.Schedule(PSCOPF.Utilitary(), DateTime("2015-01-01T06:00:00"), SortedDict(
                "unit_1" => PSCOPF.GeneratorSchedule("unit_1",
                    SortedDict(TS[1]=> PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.ON, SortedDict()), # definitive start
                                TS[2] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(missing, SortedDict()), # Non-definitive value => later values will be ignored
                                TS[3] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.OFF, SortedDict()),
                                TS[4] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.ON, SortedDict()), # ignored start
                                ),
                    SortedDict(),
                    ),
                )
        )

        # ("unit_1", TS[4]) is not reported
        #  because it is preceded by non-definitive commitment at TS[2]
        #  even if TS[3] is definitive
        #  This should not happen in PSCOPF launches
        expected = Set{Tuple{String,Dates.DateTime}}([
            ("unit_1", TS[1])
        ])

        result = PSCOPF.get_starts(schedule, initial_state)
        @test result == expected
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf