text:
using DataFrames
using CSV
using SparseArrays

"""
    admittance_matrices(nodes, lines, harmonics)

Build the nodal admittance matrices (admittance laplacian) for all harmonics.

Admittance scales with frequency: X_h = X_f * h. Return a dictionary of DataFrames.
"""
function admittance_matrices(net, harmonics)
    LY = Dict()
    for h in harmonics
        LY[h] = spzeros(ComplexF64, net.n, net.n)
        # non-diagonal elements
        for line in eachrow(net.lines)
            LY[h][line.fromID, line.toID] = -1/(line.R + 1im*line.X*h)
            # nodal admittance matrix is assumed to be symmetric
            LY[h][line.toID, line.fromID] = -1/(line.R + 1im*line.X*h)
        end
        # diagonal elements
        for i in 1:net.n
            # node shunt admittances only for harmonic frequencies
            if net.nodes.X_sh[i] > 0 && h != 1
                LY[h][i, i] = -sum(LY[h][i, :]) + 1/(1im*net.nodes.X_sh[i]*h)
            else
                LY[h][i, i] = -sum(LY[h][i, :])
            end
            # Adding shunt admittances for each pi-model line connected to node i
            for line in eachrow(net.lines)
                # check if line is connected to node i
                if line.fromID == i || line.toID == i
                    LY[h][i, i] = LY[h][i, i] +
                    (line.G + 1im*h*line.B)/2
                end
            end
        end
    end
    return LY
end

"""
    init_voltages(nodes, settings)

Initialize voltages for all nodes and harmonics using guesses from on `Settings`.
"""
function init_voltages(nodes, settings)
    u = Dict()
    for h in settings.harmonics
        if h == 1
            u[h] = DataFrame(
                v = ones(size(nodes, 1))*settings.v_f,
                ϕ = ones(size(nodes, 1))*settings.ϕ_f
            )
        else
            u[h] = DataFrame(
                v = ones(size(nodes, 1))*settings.v_h,
                ϕ = ones(size(nodes, 1))*settings.ϕ_h
            )
        end
    end
    return u
end

"""
    fund_state_vec(net, u)

Take the voltage DataFrame `u` and return a vector of voltage phases and magnitudes.
"""
function fund_state_vec(net, u)
    xϕ = u[1].ϕ[2:end]
    xv = u[1].v[(net.c+1):end]
    vcat(xϕ, xv)
end

"""
    fund_mismatch(net, u, LY)

Evaluate active and reactive power mismatch at fundamental frequency.
"""
function fund_mismatch(net, u, LY)
    LY_1 = LY[1]
    u_1 = u[1].v .* exp.(1im*u[1].ϕ)
    s = (net.nodes.P + 1im*net.nodes.Q)
    mismatch = u_1 .* conj(LY_1*u_1) + s
    f = vcat(real(mismatch[2:end]), imag(mismatch[(net.c+1):end]))
    err = maximum(abs.(f))
    return f, err
end

"""
    fund_jacobian(net, u, LY)

Construct fundamental Jacobian, containing powers derived wrt voltages.
"""
function fund_jacobian(net, u, LY)
    # helpful definitions for Jacobian construction
    u_1 = u[1].v .* exp.(1im*u[1].ϕ)
    i_diag = spdiagm(sparse(LY[1] * u_1))
    u_1_diag = spdiagm(u_1)
    u_1_diag_norm = spdiagm(u_1./abs.(u_1))

    # complex power mismatch derivatives
    dSdϕ = 1im*u_1_diag*conj(i_diag - LY[1]*u_1_diag)
    dSdv = u_1_diag_norm*conj(i_diag) + u_1_diag*conj(LY[1]*u_1_diag_norm)

    # divide sub-matrices into real and imag part, cut off slack
    dPdϕ = real(dSdϕ[2:end, 2:end])
    dPdv = real(dSdv[2:end, (net.c+1):end])
    dQdϕ = imag(dSdϕ[(net.c+1):end, 2:end])
    dQdv = imag(dSdv[(net.c+1):end, (net.c+1):end])

    vcat(hcat(dPdϕ, dPdv),
         hcat(dQdϕ, dQdv))
end

"""
    update_state_vec!(J, x, f)

Perform Newton-Raphson iteration and update the fundamental state vector.
"""
function update_state_vec!(J, x, f)
    x - J\f
end

"""
    update_fund_voltages!(net, u, x)

Write new voltages from `x` to `u`.
"""
function update_fund_voltages!(net, u, x)
    u[1].ϕ[2:end] = x[1:(net.n-1)]
    u[1].v[(net.c+1):end] = x[net.n:end]
    return u
end

"""
    pf(net, settings, LY, plt_convergence = false)

Perform fundamental power flow, iterate until mismatch is smaller than threshold.
"""
function pf(net, settings, LY, plt_convergence = false)
    # initialization
    u = init_voltages(net.nodes, settings)
    x_f = fund_state_vec(net, u)
    f_f, err_f = fund_mismatch(net, u, LY)

    n_iter_f = 0
    err_f_t = Dict()
    while err_f > settings.thresh_f && n_iter_f <= settings.max_iter_f
        J_f = fund_jacobian(net, u, LY)
        x_f = update_state_vec!(J_f, x_f, f_f)
        u = update_fund_voltages!(net, u, x_f)
        f_f, err_f = fund_mismatch(net, u, LY)
        err_f_t[n_iter_f] = err_f  # currently not returned
        n_iter_f += 1
    end

    println(u[1])
    if n_iter_f < settings.max_iter_f
        println("Fundamental power flow converged after ", n_iter_f,
              " iterations (err < ", settings.thresh_f, ").")
    elseif n_iter_f == settings.max_iter_f
        println("Maximum of ", n_iter_f, " iterations reached.")
    end
    return u
end

# Harmonic Power Flow functions
"""
    import_Norton_Equivalents(nodes, settings, folder_path="devices\\")

Import NE parameters for all nonlinear devices in nodes from folder devices."""
function import_Norton_Equivalents(nodes, settings, folder_path="devices\\")
    NE = Dict()
    nl_components = unique(nodes[nodes.type .== "nonlinear", "component"])
    for device in nl_components
        NE_df = CSV.read(folder_path * device * "_NE.csv", DataFrame)
        # transform to Complex type, enough to strip first paranthesis for successful parse
        vals = mapcols!(col -> parse.(ComplexF64, strip.(col, ['('])), NE_df[:, 3:end])
        NE_device = hcat(NE_df[:,1:2], vals)
        # filter columns for considered harmonics
        NE_device = NE_device[:, Between(begin, string(maximum(settings.harmonics)*settings.base_frequency))]

        # change to pu system and choose if coupled
        if settings.coupled
            if maximum(settings.harmonics) != Int(parse(Int, names(NE_df)[end])/settings.base_frequency)
                println("Warning! Inprecise results expected, as coupled method is being used but not the full spectrum taken into account.")
            end
            I_N = Array(NE_device[NE_device.Parameter .== "I_N_c", 3:end])/settings.base_current
            LY_N_full = NE_device[NE_device.Parameter .== "Y_N_c", 2:end]
            LY_N = Array(LY_N_full[LY_N_full.Frequency .<= maximum(settings.harmonics)*settings.base_frequency, 2:end])/settings.base_admittance
        else
            I_N = Array(NE_device[NE_device.Parameter .== "I_N_uc", 3:end])/settings.base_current
            LY_N = Array(NE_device[NE_device.Parameter .== "Y_N_uc", 3:end])/settings.base_admittance
        end
        NE[device] = [I_N, LY_N]
    end
    return NE
end

"""
    current_injections(nodes, nodeID, u, NE, harmonics)

Calculate the harmonic current injections at one node.

Uses coupled or uncoupled method based on dimensions of NE."""
function current_injections(nodes, nodeID, u, NE, harmonics)
    component = nodes[nodes.ID .== nodeID, "component"][1]
    I_N, LY_N = NE[component]
    # u as dict of dfs makes building this vector a bit complicated
    u_h = vcat([u[h][nodeID, "v"] .* exp.(1im*u[h][nodeID, "ϕ"]) for h in harmonics]...)
    # coupled: Y_N is a matrix, uncoupled: vector
    if size(LY_N)[1] > 1  # coupled case
        i_inj = vec(I_N) - vec(LY_N*u_h)
    else  # uncoupled case
        i_inj = vec(I_N) - spdiagm(vec(LY_N))*u_h
    end
    return i_inj
end

"""
    current_balance(net, settings, u, LY, NE)

Evaluate current balance at all nodes and frequencies, except linear nodes at fundamental frequency.
"""
function current_balance(net, settings, u, LY, NE)
    # fundamental admittance matrix for nonlinear nodes
    LY_1_nl = LY[1][net.m:end,:]
    u_1 = u[1].v .* exp.(1im*u[1].ϕ)
    # fundamental line currents at nonlinear nodes
    dI_1 = LY_1_nl * u_1
    # harmonic admittance matrices as diagonal block matrix
    LY_h = blockdiag([LY[h] for h in settings.harmonics[2:end]]...)
    u_h = vcat([u[h][:, "v"] .* exp.(1im*u[h][:, "ϕ"]) for h in settings.harmonics[2:end]]...)
    dI_h = LY_h * u_h

    # subtract the injected currents at each nonlinear node i
    for i in net.m:net.n
        i_inj = current_injections(net.nodes, net.nodes.ID[i], u, NE, settings.harmonics)
        dI_1[i-net.m+1] += i_inj[1]  # add injections at fundamental frequency...
        # ... and at all harmonic frequencies
        for p in 0:(settings.K-1)
            dI_h[p*net.n + i] += i_inj[p+2]
        end
    end
    vcat(dI_1, dI_h)
end

"""
    harmonic_mismatch(net, settings, u, LY, NE)

Combine current and power mismatch to construct harmonic mismatch vector `f`.
"""
function harmonic_mismatch(net, settings, u, LY, NE)
    # fundamental power mismatch at linear nodes except slack
    s = net.nodes.P[2:(net.m-1)] + 1im*net.nodes.Q[2:(net.m-1)]
    u_i = u[1][2:(net.m-1), "v"].*exp.(1im*u[1][2:(net.m-1), "ϕ"])
    u_j = u[1][:, "v"].*exp.(1im*u[1][:, "ϕ"])
    LY_ij = LY[1][2:(net.m-1), :]
    # power balance
    sl = u_i.*conj(LY_ij*u_j)
    ds = s + sl
    di = current_balance(net, settings, u, LY, NE)
    # harmonic mismatch vector
    f_c = vcat(ds, di)
    f = vcat(real(f_c), imag(f_c[(net.c):end]))  # crop PV nodes
    err_h = maximum(abs.(f))
    return f, err_h
end

"""
    harmonic_state_vec(net, u, harmonics)

Construct harmonic state vector `x` from fundamental and harmonic voltages.
"""
function harmonic_state_vec(net, u, harmonics)
    # fundamental voltages, crop slack and magnitude for PV nodes
    xϕ = u[1].ϕ[2:end]
    xv = u[1].v[(net.c+1):end]
    # all voltages and angles at all harmonic frequencies
    for h in harmonics[2:end]
        xϕ = vcat(xϕ, u[h].ϕ)
        xv = vcat(xv, u[h].v)
    end
    vcat(xϕ, xv)
end

"""
    build_harmonic_jacobian(net, settings, u, LY, NE)

Construct harmonic Jacobian `J` containing all derivatives of power and current mismatches wrt voltages from harmonic state vector.
"""
function build_harmonic_jacobian(net, settings, u, LY, NE)
    m = net.m
    n = net.n
    c = net.c
    K = settings.K

    # helpful definitions for Jacobian construction
    u_vec = vcat([u[h].v .* exp.(1im*u[h].ϕ) for h in settings.harmonics]...)
    v_vec = vcat([u[h].v for h in settings.harmonics]...)
    u_diag = spdiagm(u_vec)
    u_norm = u_vec./v_vec
    u_norm_diag = spdiagm(u_norm)
    LY_diag = blockdiag([LY[h] for h in settings.harmonics]...)

    # construct Jacobian sub-matrices
    dIdv = LY_diag*u_norm_diag
    dIdϕ = 1im*LY_diag*u_diag

    # indices of first nonlinear node at each harmonic
    nl_idx_start = m:n:n*(K+1)
    nl_idx_all = vcat([nl:(nl+n-m) for nl in nl_idx_start]...)
    u_nl = u_vec[nl_idx_all]
    v_nl = v_vec[nl_idx_all]
    u_nl_norm = u_nl./v_nl

    if settings.coupled
        # iterating through blocks vertically
        for h in 0:K
            # ... and horizontally
            for p in 0:K
                # iterating through nonlinear nodes
                for i in m:n
                    # within NE "[2]" points to LY_N
                    LY_N = NE[net.nodes.component[i]][2]
                    # subtract derived current injections at respective idx
                    dIdv[h*n+i, p*n+i] -= LY_N[h+1, p+1]*u_nl_norm[(i-m+1)+p*(n-m+1)]
                    dIdϕ[h*n+i, p*n+i] -= 1im*LY_N[h+1, p+1]*u_nl[(i-m+1)+p*(n-m+1)]
                end
            end
        end
    else
        # iterating through blocks diagonally (p=h)
        for h in 0:K
            for i in m:n
                # LY_N is one-dimensional in the uncoupled case
                LY_N = NE[net.nodes.component[i]][2]
                dIdv[h*n+i, h*n+i] -= LY_N[h+1]*u_nl_norm[(i-m+1)+h*(n-m+1)]
                dIdϕ[h*n+i, h*n+i] -= 1im*LY_N[h+1]*u_nl[(i-m+1)+h*(n-m+1)]
            end
        end
    end

    # crop slack and PV nodes
    dIdϕ = dIdϕ[m:end, 2:end]
    dIdv = dIdv[m:end, (c+1):end]

    # helpful definitions for dS construction
    LY_1 = LY[1]
    u_1 = u[1].v .* exp.(1im*u[1].ϕ)
    i_diag = spdiagm(LY_1*u_1)
    u_diag = spdiagm(u_1)
    u_diag_norm = spdiagm(u_1./abs.(u_1))

    # fundamental powers derived wrt to voltage phase and magnitude
    dS1dϕ1 = 1im*u_diag*(conj(i_diag - LY_1*u_diag))
    dS1dv1 = u_diag_norm*conj(i_diag) + u_diag*conj(LY_1*u_diag_norm)

    # only fundamental power derivatives are non-zero
    dSdϕ = hcat(dS1dϕ1, zeros(n, n*K))
    dSdv = hcat(dS1dv1, zeros(n, n*K))

    # crop slack and PV nodes from power derivatives
    dPdϕ = real(dSdϕ[2:(m-1), 2:end])
    dPdv = real(dSdv[2:(m-1), (c+1):end])
    dQdϕ = imag(dSdϕ[(c+1):(m-1), 2:end])
    dQdv = imag(dSdv[(c+1):(m-1), (c+1):end])

    # combine all sub-matrices and return complete Jacobian
    vcat(hcat(dPdϕ, dPdv),
         hcat(real(dIdϕ), real(dIdv)),
         hcat(dQdϕ, dQdv),
         hcat(imag(dIdϕ), imag(dIdv)))
end

"""
    update_fund_voltages!(net, u, x, settings)

Write new voltages from `x` to `u` for all harmonics.
"""
function update_harmonic_voltages!(net, u, x, settings)
    n = net.n
    c = net.c

    # slice x to separate voltage magnitude and phase
    xϕ = x[1:(n*settings.K1-1)]
    xv = x[n*settings.K1:end]

    # xϕ = xϕ .% (2*π)  # ensure phase smaller 2π, not really necessary

    # update fundamental voltages
    u[1].ϕ[2:end] = xϕ[1:(n-1)]
    u[1].v[(c+1):end] = xv[1:(n-c)]

    # crop slack and PV nodes
    xϕ_h = xϕ[n:end]
    xv_h = xv[(n-c+1):end]

    for h in settings.harmonics[2:end]
        # enumerate harmonics for easier indexing
        i = findall(settings.harmonics[2:end] .== h)[1] - 1
        u[h].ϕ = xϕ_h[(i*n+1):(i*n+n)]
        u[h].v = xv_h[(i*n+1):(i*n+n)]
    end
    return u
end

"""
    hpf(net, settings)

Solve the harmonic power flow problem for a given `PowerGrid` and `Settings`."""
function hpf(net, settings)
    LY = admittance_matrices(net, settings.harmonics)
    u = pf(net, settings, LY)
    NE = import_Norton_Equivalents(net.nodes, settings)
    f, err_h = harmonic_mismatch(net, settings, u, LY, NE)
    x = harmonic_state_vec(net, u, settings.harmonics)
    n_iter_h = 0
    err_h_t = Dict()
    while err_h > settings.thresh_h && n_iter_h < settings.max_iter_h
        J = build_harmonic_jacobian(net, settings, u, LY, NE)
        x = update_state_vec!(J, x, f)
        u = update_harmonic_voltages!(net, u, x, settings)
        f, err_h = harmonic_mismatch(net, settings, u, LY, NE)
        err_h_t[n_iter_h] = err_h
        n_iter_h += 1
    end

    # getting rid of negative voltage magnitudes:
    for h in settings.harmonics
        u[h].ϕ[u[h].v .< 0] .+= π
        u[h].ϕ .= u[h].ϕ .% (2*π)
        u[h].v[u[h].v .< 0] = -u[h].v[u[h].v .< 0]
    end

    if n_iter_h < settings.max_iter_h
        println("Harmonic power flow converged after ", n_iter_h,
                " iterations (err < ", settings.thresh_h, ").")
    elseif n_iter_h == settings.max_iter_h
        println("Maximum of ", n_iter_h, " iterations reached. Harmonic power flow did not converge.")
    end
    return u, err_h_t, n_iter_h
end


text:
#Gemm Layer

mutable struct Gemm
    weight
    bias
    alpha
    beta
end

function Gemm(;alpha, beta=1, transA=0, transB=0, w, b)
    if transA != 0; w = transpose(w); end
    if transB != 0; x = transpose()
    Gemm(w,b,alpha,beta)
end

function (g::Gemm)(x)
    g.alpha * g.weight * x + (g.beta * g.bias)
end

#graph node to Gemm Layer
function node_to_gemm(node, weightdims, g)
    alpha = node.attribute[:alpha]
    beta = node.attribute[:beta]
    transA = node.attribute[:transA]
    transB = node.attribute[:transB]
    transB = 0
    w_name = node.input[2]
    b_name = node.input[3]
    #wsize = weightdims[w_name]
    #bsize = weightdims[b_name]
    w = g.initializer[w_name]
    b = g.initializer[b_name]
    println(transA, transB)
    layer = Gemm(alpha=alpha, beta=beta, transA=transA, transB=transB, w=w, b=b)
end

ext: jl


text:
# Autogenerated wrapper script for Triangle_jll for i686-linux-gnu
export libtriangle

JLLWrappers.@generate_wrapper_header("Triangle")
JLLWrappers.@declare_library_product(libtriangle, "libtriangle.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libtriangle,
        "lib/libtriangle.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function recover_dynamics(X, t, svname, grpname, p_ = Float32[1.3, 0.9, 0.8, 1.8])

    @info "Started run $svname $(grpname)."

    # Define the network 2->5->5->5->2
    U = FastChain(
        FastDense(2, 5, rbf),
        FastDense(5, 5, rbf),
        FastDense(5, 5, rbf),
        FastDense(5, 2),
    )
    # Get the initial parameters
    p = initial_params(U)

    # Define the hybrid model
    function ude_dynamics!(du, u, p, t, p_true)
        û = U(u, p) # Network prediction
        du[1] = p_true[1] * u[1] + û[1]
        du[2] = -p_true[4] * u[2] + û[2]
    end

    # Closure with the known parameter
    nn_dynamics!(du, u, p, t) = ude_dynamics!(du, u, p, t, p_)
    # Define the problem
    prob_nn = ODEProblem(nn_dynamics!, X[:, 1], (t[1], t[end]), p)

    ## Function to train the network
    # Define a predictor
    function predict(θ, X = X[:, 1], T = t)
        Array(
            solve(
                prob_nn,
                Vern7(),
                u0 = X,
                p = θ,
                tspan = (T[1], T[end]),
                saveat = T,
                abstol = 1e-6,
                reltol = 1e-6,
                sensealg = ForwardDiffSensitivity(),
            ),
        )
    end

    # Simple L2 loss
    function loss(θ)
        X̂ = predict(θ)
        sum(abs2, X .- X̂) + Float32(1e-4)*sum(abs2, θ)/length(θ)
    end

    # Container to track the losses
    losses = Float32[]

    # Callback to show the loss during training
    callback(θ, l) = begin
        push!(losses, l)
        false
    end

    ## Training

    # First train with ADAM for better convergence -> move the parameters into a
    # favourable starting positing for BFGS
    res1 = DiffEqFlux.sciml_train(
        loss,
        p,
        ADAM(0.1f0),
        cb = callback,
        maxiters = 200,
    )
    println("$svname $(grpname) : Training loss after $(length(losses)) iterations: $(losses[end])")
    # Train with BFGS
    res2 = DiffEqFlux.sciml_train(
        loss,
        res1.minimizer,
        BFGS(initial_stepnorm = 0.01f0),
        cb = callback,
        maxiters = 10000,
    )
    println(
        "$svname $(grpname) : Final training loss after $(length(losses)) iterations: $(losses[end])",
    )

    p_trained = res2.minimizer

    X̂ = predict(p_trained)
    Ŷ = U(X̂, p_trained)

    # Create a Basis
    @variables u[1:2]
    # Generate the basis functions, multivariate polynomials up to deg 5
    # and sine
    b = [polynomial_basis(u, 5); sin.(u)]
    basis = Basis(b, u)

    # Create an optimizer for the SINDy problem
    opt = SR3(Float32(1e-2), Float32(0.1))
    # Create the thresholds which should be used in the search process
    λ = Float32.(exp10.(-7:0.1:7))
    # Target function to choose the results from; x = L0 of coefficients and L2-Error of the model
    g(x) = x[1] < 1 ? Inf : norm(x, 2)

    # Test on uode derivative data
    Ψ = SINDy(X̂, Ŷ, basis, λ,  opt, g = g, maxiter = 50000, normalize = true, denoise = true, convergence_error = Float32(1e-10)) # Succeed

    # Extract the parameter
    p̂ = parameters(Ψ)

    # Just the equations
    b = Basis((u, p, t) -> Ψ(u, ones(Float32, length(p̂)), t), u)

    # Retune for better parameters -> we could also use DiffEqFlux or other parameter estimation tools here.
    Ψf = SINDy(
        X̂,
        Ŷ,
        b,
        STRRidge(0.01f0),
        maxiter = 100,
        convergence_error = Float32(1e-18),
    ) # Succeed

    p̂ = parameters(Ψf)

    @info "$svname $(grpname) : Found equations : $(Ψf.equations.eqs)"
    @info "$svname $(grpname) : Parameter estimation : $(p̂)"

    ## Save the results

    jldopen("$(svname)_recovery_loop.jld2", "a+") do file
        mygroup = JLD2.Group(file, grpname)
        mygroup["X"] = X
        mygroup["t"] = t
        mygroup["initial_parameters"] = p
        mygroup["trained_parameters"] = p_trained
        mygroup["losses"] = losses
        mygroup["result"] = Ψf
    end

    @info "Finished run $svname $(grpname)."

    return
end


text:
using Oceananigans.Grids: Face, Face, Center
using Oceananigans.Operators: ζ₃ᶠᶠᵃ
using Oceananigans.Fields: KernelComputedField, ComputedFieldBoundaryConditions

using KernelAbstractions: @kernel

@kernel function compute_vertical_vorticity!(ζ, grid, u, v)
    i, j, k, = @index(Global, NTuple)
    @inbounds ζ[i, j, k] = ζ₃ᶠᶠᵃ(i, j, k, grid, u, v)
end

"""
    VerticalVorticityField(model)

Returns a KernelComputedField that `compute!`s vertical vorticity in a
manner consistent with the `VectorInvariant` momentum advection scheme
for curvilinear grids.

In particular, `VerticalVorticityField` uses `ζ₃ᶠᶠᵃ`, which in turn computes
the vertical vorticity by first integrating the velocity field around the borders
of the vorticity cell to find the vertical circulation, and then dividing by the area of
the vorticity cell to compute vertical vorticity.
"""
function VerticalVorticityField(model)
    return KernelComputedField(Face, Face, Center, compute_vertical_vorticity!, model,
                        boundary_conditions = ComputedFieldBoundaryConditions(model.grid, (Face, Face, Center)),
                        computed_dependencies = (model.velocities.u, model.velocities.v))
end


text:
#################################################
#
#   Type definitions and properties
#
#################################################

"""
Joint finite distribution represented by probability mass function.
"""
mutable struct JointPMF{T<:AbstractFloat}
    varlist::SVarList
    p::Array{T}

    function JointPMF{T}(vl::SVarList, p::Array{T}) where T<:AbstractFloat
        nv = length(vl)
        ndims(p) == nv ||
            throw(DimensionMismatch(
                "ndims(p) does not match the number of variables."))
        for i = 1:nv
            size(p,i) == cardinality(vl[i]) ||
                throw(DimensionMismatch(
                    "size(p,i) does not match the corresponding space size."))
        end
        new(vl, p)
    end
end

JointPMF{T<:AbstractFloat}(vl::VarList{SVar}, p::Array{T}) = JointPMF{T}(vl, p)

function JointPMF{T<:AbstractFloat}(vnames::Vector{String}, p::Array{T})
    vlst = SVar[dvar(v, size(p,i)) for (i, v) in enumerate(vnames)]
    JointPMF{T}(VarList(vlst), p)
end

function JointPMF{T<:AbstractFloat,N}(vnames::NTuple{N,String}, p::Array{T,N})
    vlst = SVar[dvar(v, size(p,i)) for (i, v) in enumerate(vnames)]
    JointPMF{T}(VarList(vlst), p)
end

nvars(jp::JointPMF) = length(jp.varlist)
vars(jp::JointPMF) = jp.varlist
probs(jp::JointPMF) = jp.p

function show{T}(io::IO, jp::JointPMF{T})
    println(io, "JointPMF{$T} with $(nvars(jp)) variables:")
    for v in jp.varlist
        println(io, "  $(v.id): $(values(v.space))")
    end
end

#################################################
#
#   Compute marginals and conditionals
#
#################################################

function marginal_i{T}(jp::JointPMF{T}, vidx::Int)
    if nvars(jp) == 1
        vidx == 1 || throw(BoundsError("Variable index out of range."))
        return copy(probs(jp))::Vector{T}
    else
        reduc_region = Utils.range_rm(nvars(jp), vidx)
        _mp = sum(probs(jp), reduc_region)
        return reshape(_mp, size(_mp, vidx))::Vector{T}
    end
end

function marginal_i{T}(jp::JointPMF{T}, vinds)
    @assert !isa(vinds, Int)
    !isempty(vinds) || throw(ArgumentError("Input indices are empty."))

    # process the case for single variable
    if length(vinds) == 1
        return marginal_i(jp, first(vinds))
    end

    # marginalize out unspecified variables
    vinds_ = collect(vinds)::Vector{Int}
    reduc_region = Utils.range_rm(nvars(jp), vinds_)
    if isempty(reduc_region)
        # no remaining variables, just keep p
        mp = copy(probs(jp))
    else
        s_vinds = issorted(vinds_) ? vinds_ : sort(vinds_)
        # reduce along unspecified dimensions
        _mp = sum(probs(jp), reduc_region)
        # squeeze out reduced dimensions
        rsiz = ntuple(i->size(_mp, s_vinds[i]), length(s_vinds))
        mp = reshape(_mp, rsiz)
    end
    # make output, permute dimensions according to specified order
    issorted(vinds)    ? mp :
    length(vinds) == 2 ? transpose(mp) :
                         permutedims(mp, Utils.isortperm(vinds_))
end

marginal(jp::JointPMF, vids) = marginal_i(jp, indexof(vids, vars(jp)))

function conditional_i(jp::JointPMF, vinds, cinds)
    # no conditions, degenerate to marginal distribution
    if isempty(cinds)
        return marginal_i(jp, vinds)
    end

    # compute the marginal of both targets and conditions
    mp = marginal_i(jp, Int[vinds..., cinds...])

    # obtain the conditional probabilities
    mp ./= sum(mp, 1:length(vinds))
end

conditional(jp::JointPMF, vs, cs) =
    conditional_i(jp, indexof(vs, vars(jp)), indexof(cs, vars(jp)))


text:
using Plots
using Random

αs = [0.0, 0.9, 0.98]
n = 100
p = plot()

function getpricepath(α, n)
    x = zeros(n + 1)
    x[1] = 0.0
    for t in 1:n
        x[t+1] = α * x[t] + rand()
    end
    x
end

pricepaths = [getpricepath(α, n) for α in αs]

for (idx, pricepath) in enumerate(pricepaths)
    plot!(p, pricepath, title="Random price paths", label="alpha = $(αs[idx])", xlabel="Timestep", ylabel="Price")
end

p

max_stars_repo_name: JuliaAI/MLFlowClient

max_issues_repo_name: JuliaAI/MLFlowClient

max_forks_repo_name: JuliaAI/MLFlowClient


text:
include("grammar.jl")
include("Node.jl")

function no_action(rule, value, first, last, children)
  return Node(rule.name, value, first, last, children, typeof(rule))
end

or_default_action(rule, value, first, last, children) = children[1]

# Terminal
type Terminal <: Rule
  name::AbstractString
  value::AbstractString
  action

  function Terminal(name::AbstractString, value::AbstractString)
    return new(name, SubString(value, 1), no_action);
  end

  function Terminal(name::AbstractString, value::Char)
    return new(name, "$value", no_action);
  end

  function Terminal(value::AbstractString)
    return new("", SubString(value, 1), no_action)
  end
end

function show(io::IO, t::Terminal)
  print(io, "$(t.value)");
end

function parseDefinition(name::AbstractString, value::AbstractString, pdata::ParserData)
  return Terminal(name, value)
end

function parseDefinition(name::AbstractString, value::Char, pdata::ParserData)
  return Terminal(name, value)
end

# References
type ReferencedRule <: Rule
  name::AbstractString
  symbol::Symbol
  action

  function ReferencedRule(name::AbstractString, symbol::Symbol)
    return new(name, symbol, no_action)
  end

  function ReferencedRule(symbol::Symbol)
    return new("", symbol, no_action)
  end
end

function parseDefinition(name::AbstractString, sym::Symbol, pdata::ParserData)
  fn = get(pdata.parsers, sym, nothing)

  if fn !== nothing
    return fn(name, pdata, [])
  end

  return ReferencedRule(name, sym)
end

# And
type AndRule <: Rule
  name::AbstractString
  values::Array{Rule}
  action

  function AndRule(name::AbstractString, values::Array{Rule})
    return new(name, values, no_action)
  end

  function AndRule(values::Array{Rule})
    return new("", values, no_action);
  end
end

get_children(rule::AndRule) = rule.values

function +(name::AbstractString, pdata::ParserData, args::Array)
  if length(args) > 1
    # Addition can contain multiple entries
    values::Array{Rule} = [parseDefinition("$(name)_$i", arg, pdata) for (i, arg) in enumerate(args)]
    return AndRule(name, values)
  else
    # it's prefix, so it maps to one or more rule
    return OneOrMoreRule(name, parseDefinition("$(name)_values", args[1], pdata))
  end
end

function show(io::IO, rule::AndRule)
  values = [string(r) for r in rule.values]
  joinedValues = join(values, " ")
  print(io, "($(rule.name),$joinedValues,$(rule.action))");
end

# TODO: check if actually being used
+(a::Rule, b::Rule) = AndRule([a, b]);
+(a::AndRule, b::AndRule) = AndRule(append!(a.values, b.values));
+(a::AndRule, b::Rule) = AndRule(push!(a.values, b));
+(a::Rule, b::AndRule) = AndRule(push!(b.values, a));

# Or
type OrRule <: Rule
  name::AbstractString
  values::Array{Rule}
  action

  function OrRule(name::AbstractString, left::OrRule, right::OrRule)
    return new(name, append!(left.values, right.values), or_default_action)
  end

  function OrRule(name::AbstractString, left::OrRule, right::Rule)
    return new(name, push!(left.values, right), or_default_action)
  end

  function OrRule(name::AbstractString, left::Rule, right::OrRule)
    return new(name, [left, right], or_default_action)
  end

  function OrRule(name::AbstractString, left::Rule, right::Rule)
    return new(name, [left, right], or_default_action)
  end
end

get_children(rule::OrRule) = rule.values

function show(io::IO, rule::OrRule)
  values = [string(r) for r in rule.values]
  joinedValues = join(values, "|")
  print(io, "($(rule.name),$joinedValues,$(rule.action))")
end

function |(name::AbstractString, pdata::ParserData, args::Array)
  left = parseDefinition("$(name)_1", args[1], pdata)
  right = parseDefinition("$(name)_2", args[2], pdata)
  return OrRule(name, left, right)
end

function parseDefinition(name::AbstractString, range::UnitRange, pdata::ParserData)
  values = [Terminal(value) for value in range];
  return OrRule(name, values);
end

# OneOrMore
type OneOrMoreRule <: Rule
  name::AbstractString
  value::Rule
  action

  function OneOrMoreRule(name::AbstractString, value::Rule)
    return new(name, value, no_action)
  end

  function OneOrMoreRule(value::Rule)
    return new("", value, no_action);
  end
end

function show(io::IO, rule::OneOrMoreRule)
  print(io, "+($(rule.value))");
end

# ZeroOrMore
type ZeroOrMoreRule <: Rule
  name::AbstractString
  value::Rule
  action

  function ZeroOrMoreRule(name::AbstractString, value::Rule)
    return new(name, value, no_action)
  end

  function ZeroOrMoreRule(value::Rule)
    return new("", value, no_action);
  end
end

function show(io::IO, rule::ZeroOrMoreRule)
  print(io, "*($(rule.value))");
end

function *(name::AbstractString, pdata::ParserData, args::Array)
  if length(args) == 1
    # it's a prefix, so it maps to zero or more rule
    return ZeroOrMoreRule(parseDefinition(name, args[1], pdata))
  end
end

# Multiple
type MultipleRule <: Rule
  name::AbstractString
  value::Rule
  minCount::Int
  maxCount::Int
  action

  function MultipleRule(name::AbstractString, value::Rule, minCount::Int, maxCount::Int)
    return new(name, value, minCount, maxCount, no_action)
  end

  function MultipleRule(value::Rule, minCount::Int, maxCount::Int)
    return new("", value, minCount, maxCount, no_action)
  end
end

function show(io::IO, rule::MultipleRule)
  print(io, "($(rule.value)){$(rule.minCount), $(rule.maxCount)}");
end

function ^(name::AbstractString, pdata::ParserData, args::Array)
  # FIXME: not sure this is correct..
  count = args[2]
  return MultipleRule(args[1], count.args[1], count.args[2])
end

# RegEx
type RegexRule <: Rule
  name::AbstractString
  value::Regex
  action

  function RegexRule(name::AbstractString, value::Regex)
    return new(name, value, no_action)
  end

  function RegexRule(value::Regex)
    return new("", Regex("^$(value.pattern)"), no_action)
  end
end

function show(io::IO, rule::RegexRule)
  print(io, "r($(rule.value.pattern))")
end

function parseDefinition(name::AbstractString, regex::Regex, pdata::ParserData)
  # TODO: Need to do this to ensure we always match at the beginning,
  # but there should be a safer way to do this
  modRegex = Regex("^$(regex.pattern)")
  return RegexRule(name, modRegex)
end

# Optional
type OptionalRule <: Rule
  name::AbstractString
  value::Rule
  action

  function OptionalRule(name::AbstractString, value::Rule)
    return new(name, value, or_default_action)
  end

  function OptionalRule(value::Rule)
    return new("", value, or_default_action)
  end
end

function ?(name::AbstractString, pdata::ParserData, args::Array)
  return OptionalRule(parseDefinition(name, args[1], pdata))
end

# Look ahead
type LookAheadRule <: Rule
    name::AbstractString
    value::Rule
    action

    function LookAheadRule(name::AbstractString, value::Rule)
        return new(name, value, no_action)
    end
end

function >(name::AbstractString, pData::ParserData, args::Array)
    if length(args) == 1
        return LookAheadRule(name, parseDefinition("$(name)_value", args[1], pData))
    end
end

# Suppress
type SuppressRule <: Rule
  name::AbstractString
  value::Rule
  action

  function SuppressRule(name::AbstractString, value::Rule)
    return new(name, value, no_action)
  end
end

function -(name::AbstractString, pdata::ParserData, args::Array)
  if length(args) == 1
    return SuppressRule(name, parseDefinition("$(name)_value", args[1], pdata))
  end
end

# List
type ListRule <: Rule
  name::AbstractString
  entry::Rule
  delim::Rule
  min::Int
  action

  function ListRule(name::AbstractString, entry::Rule, delim::Rule)
    return new(name, entry, delim, 1, no_action)
  end

  function ListRule(name::AbstractString, entry::Rule, delim::Rule, min::Int)
    return new(name, entry, delim, min, no_action)
  end
end

function list(name::AbstractString, pdata::ParserData, args::Array)
  entry = parseDefinition("$(name)_entry", args[1], pdata)
  delim = parseDefinition("$(name)_delim", args[2], pdata)

  if length(args) > 2
    return ListRule(name, entry, delim, args[3])
  end

  return ListRule(name, entry, delim)
end

# Not
type NotRule <: Rule
  name
  entry
  action

  function NotRule(name::AbstractString, entry::Rule)
    return new(name, entry, no_action)
  end
end

function !(name::AbstractString, pdata::ParserData, args::Array)
  entry = parseDefinition("$(name)_entry", args[1], pdata)
  return NotRule(name, entry)
end

function show(io::IO, rule::NotRule)
  print(io, "!($(rule.entry))");
end

function empty(name::AbstractString, pdata::ParserData, args::Array)
  return EmptyRule(name)
end

type EndOfFileRule <: Rule
  name::AbstractString
  action

  EndOfFileRule(name::AbstractString) = new(name, no_action)
end

function eof(name::AbstractString, pdata::ParserData, args::Array)
  return EndOfFileRule(name)
end

# common parsers
type IntegerRule <: Rule
  name::AbstractString
  action

  IntegerRule(name::AbstractString) = new(name, no_action)
end

type FloatRule <: Rule
  name::AbstractString
  action

  FloatRule(name::AbstractString) = new(name, no_action)
end

function integer(name::AbstractString, pdata::ParserData, args::Array)
  return IntegerRule(name)
end

function float(name::AbstractString, pdata::ParserData, args::Array)
  return FloatRule(name)
end

macro grammar(name, definitions)
  parsers = [:+, :*, :?, :|, :-, :^, :!, :>, :list, :empty, :eof, :integer, :float]
  mapped_parsers = map_symbol_to_function(parsers)
  return parseGrammar(name, definitions, ParserData(mapped_parsers))
end

# by default don't show anything
displayValue{T <: Rule}(value, ::Type{T}) = ""

# except for terminals and regex
displayValue(value, ::Type{Terminal}) = "'$value',"
displayValue(value, ::Type{RegexRule}) = "'$value',"
displayValue(value, ::Type{IntegerRule}) = "$value,"
displayValue(value, ::Type{FloatRule}) = "$value,"

function map_symbol_to_function(lst)
  m = Dict{Symbol, Function}()
  for sym in lst
    m[sym] = eval(sym)
  end

  return m
end


text:
module DictionarySlicing

using OrderedCollections

export collectall, sliced, directslicing

"""
    collectall(sets...; maxdepth=5)

Recursively collects elements from `sets`.

Keyword arguments:
maxdepth --------- determines maximum recursion level
currentdepth ----- current recursion level
"""
function collectall(sets...; maxdepth = 5, currentdepth = 1)
	ff = []
	f1 = collect([sets...])
	f2 = collect([Base.Iterators.flatten(f1)...])

    for x in f2
		if typeof(x) <: Union{AbstractArray,AbstractRange} && currentdepth < maxdepth
			xx = try
                collect([Base.Iterators.flatten(x)...])
            catch
                collect(x)
            end
			push!(ff,collectall(xx; maxdepth, currentdepth = currentdepth+1)...)
		else
			push!(ff,x)
		end
	end

	typ = typeof(ff[1])
    try
        ff = Vector{typ}(ff)
    catch

    end

	return ff
end

"""
    sliced(D::T, idxs...; keep = :first, filter = nothing) where {T<:AbstractDict}

Returns the slices from dict `D`.

Keyword args:
keep --------- either keep the :first or :last instance of a sliced
filter ------- true/false-returning function to apply to collected indices
"""
function sliced(D::T, idxs...; keep = :first, filter = nothing) where {T<:AbstractDict}
	inputs = collectall(idxs...)
	_keys = []
	_indices = []
	for i in 1:size(inputs,1)
		if typeof(inputs[i])<:AbstractString
			push!(_keys,(i,inputs[i]))
		else
			push!(_indices,(i,inputs[i]))
		end
	end

	if filter != nothing
		try
			if keep == :lastbefore || keep == "lastbefore"
				revindices = reverse(indices)
				revindices = Base.filter(filter,revindices)
				od = OrderedDict([([keys(D)...][i],[values(D)...][i]) for i in revindices])
				rv = length(od):-1:1
				return sliced(od,rv...)
			elseif keep == :firstbefore || keep == "firstbefore"
				indices = unique(indices)
				return sliced(D,indices...; filter = filter)
			else
				indices = Base.filter(filter,indices)
			end
		catch
			error("could not filter")
		end
	end
	if keep == :last || keep == "last"
		revindices = reverse(indices)
		od = OrderedDict([([keys(D)...][i],[values(D)...][i]) for i in revindices])
		rv = length(od):-1:1
		return sliced(od,rv...)
	end
	ks = [[keys(D)...][i] for i in indices]
	vs = [[values(D)...][i] for i in indices]

    return T(ks.=>vs)
end
"""
    directslicing()

Activates indexing/slicing for OrderedDict with function indexing.
`dict(2,3)` returns an OrderedDict of the 2nd and 3rd entry.
"""
function directslicing()
	@eval(function (D::OrderedDict)(idxs...; keep = :first, filter = nothing)
		indices = collectall(idxs...)

		if filter != nothing
			try
				if keep == :lastbefore || keep == "lastbefore"
					revindices = reverse(indices)
					revindices = Base.filter(filter,revindices)
					od = OrderedDict([([keys(D)...][i],[values(D)...][i]) for i in revindices])
					rv = length(od):-1:1
					return od(rv...)
				elseif keep == :firstbefore || keep == "firstbefore"
					indices = unique(indices)
					return D(indices...; filter = filter)
				else
					indices = Base.filter(filter,indices)
				end
			catch
				error("could not filter")
			end
		end
		if keep == :last || keep == "last"
			revindices = reverse(indices)
			od = OrderedDict([([keys(D)...][i],[values(D)...][i]) for i in revindices])
			rv = length(od):-1:1
			return od(rv...)
		end

	    return OrderedDict([([keys(D)...][i],[values(D)...][i]) for i in indices])
	end)
end

end


text:
module AssociationRules

using RCall
using DataFrames

export apriori, dataset, sanitize_input!

include("arules.jl")
include("utils.jl")

init_arules()

end # module


text:
function check_why_there_is_no_sell(closed_order, id_str, open_orders, client, config, bot_path)
    t = Pinguu.examine_order(closed_order, id_str, client, config, bot_path)
    fees = config["current_fees"][config["symbol"]]
    new_id    = "SELL_" * id_str
    new_quant = floor(t["use_vol"], digits=config["precision"])

    new_price = Pinguu.price2string(round((abs(t["current_profit"]) / new_quant *
        (1 + fees["maker"] +
            parse(Float64, config["take_profit"])/100)
            * t["price_correction"]), digits=config["sell_precision"]), config["sell_precision"])
    if abs(t["current_profit"]) > config["min_qty"]
        order = try client.order_limit_sell(
            symbol=config["symbol"],
            quantity=new_quant,
            price=new_price);
            catch
            nothing
        end
        open(joinpath(bot_path, "trigger.json"), "w") do f
            JSON.print(f, Dict(), 4)
        end
        @info(order)
        sleep(2)
        if order !== nothing
            order["time"] = string(now())
            open_orders["SELL"][new_id] = order
            open(joinpath(bot_path, "open_orders.json"), "w") do f
                JSON.print(f, open_orders, 4)
            end
        else
            open_orders["SELL"][new_id] = Dict()
            open_orders["SELL"][new_id]["origQty"] = new_quant
            open_orders["SELL"][new_id]["price"]   = new_price
            open_orders["SELL"][new_id]["status"]  = "place order"
            open_orders["SELL"][new_id]["time"]    = string(now())
            open_orders["SELL"][new_id]["fail"]    = 0
            open_orders["SELL"][new_id]["if_fail"] = new_quant - 1/10^config["precision"]

            open(joinpath(bot_path, "open_orders.json"), "w") do f
                JSON.print(f, open_orders, 4)
            end
        end
    end
end
export check_why_there_is_no_sell


text:
using Pkg;Pkg.activate(".");Pkg.instantiate();

using TOML, DataFrames, CSV, ArgParse, Distributed, Dates

@everywhere import Base.@kwdef

# Used to search through command line for particular arguments
function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--paramFile","-p"
            help = "Specifies an inputfile to replace default parameters"
        "--verbosity", "-v"
            help = "Prints information to command line for debugging"
            action = :store_true
		"--cnFile", "-c"
			help = "Contains numeric array of copy numbers and header of which chromosomes"
			default = "birthLandscapeBrainCancer.txt"
		"--birthRateFile", "-b"
			help = "Contains numeric vector of birth rates"
			default = "myY.txt"
		"--outputfile","-o"
			help = "Specifies the filename for saving data from the simulation"
			default = "output_0.csv"
		"--initfile", "-i"
			help = "file that contains simulation information"
    end

    return parse_args(s)
end

# Struct containing input data needed for simulation with default values
@everywhere @kwdef mutable struct Input

	debugging::Int = 0				# prints info from ploidyMovement
	stepsize::Real = 1.0				# discretization of chromosome
	minChrom::Real = 1.0				# minimum chromosome allowed
	maxChrom::Real = 5.0				# maximum chromosome allowed
	deathRate::Float64 = 0.1			# universal death rate
	misRate::Float64 = 0.15				# universal missegregation rate
	finalDay::Real = 30.0				# end of simulation
	replating::Bool = false				# Whether we replate the cells
	startPop::Real = 1e3		# Starting population size
	maxPop::Real = 1e6				# Max population before replating
	compartmentMinimum::Bool = false	# Sets sizes < 1 to 0 if true

end

@everywhere function Input(paramFile::String)

	# Dictionary that contains info from the toml (INPUT) file
	params = TOML.tryparsefile(paramFile)

	# Get the parameters for the struct.
	debugging=get(params,"debugging",0)
	stepsize=get(params,"stepsize",1.0)
	minChrom=get(params,"minChrom",1.0)
	maxChrom=get(params,"maxChrom",5.0)
	deathRate=get(params,"deathRate",0.1)
	misRate=get(params,"misRate",0.15)
	finalDay=get(params,"finalDay",30.0)
	replating=get(params,"replating",false)
	startPop=get(params,"startPop",1e3)
	maxPop=get(params,"maxPop",1e6)
	compartmentMinimum=get(params,"compartmentMinimum",false)

	Input(
		debugging,
		stepsize,
		minChrom,
		maxChrom,
		deathRate,
		misRate,
		finalDay,
		replating,
		startPop,
		maxPop,
		compartmentMinimum
		)

end

function initialize()

	# Default values to be fed into ploidyMovement
	params = Input()

	# Grab command line args
	parsed_args = parse_commandline()

	# Replaces default parameters if input file is given
	if !isnothing(parsed_args["paramFile"]) && isfile(parsed_args["paramFile"])
		params = Input(parsed_args["paramFile"])
	end

	# Prints info if set to true
	verbosity = parsed_args["verbosity"]

	# Grab file containing copy number variation
	CNmatFilename = parsed_args["cnFile"]

	# Grab file containing birth rates
	birthFilename = parsed_args["birthRateFile"]

	# Grab output file name
	outputFile = parsed_args["outputfile"]

	# Grab initial condition file
	initfile = parsed_args["initfile"]

	# Error checking
	if !isfile(CNmatFilename)
		error("$CNmatFilename file cannot be found")
	end
	if !isfile(birthFilename)
		error("$birthFilename file cannot be found")
	end

	# Check to see if input file is given
    if verbosity
		println("parsed_args:")
		for (arg,val) in parsed_args
			println("  $arg => $val")
		end
	end

	return params, CNmatFilename, birthFilename, initfile, outputFile, verbosity
end

function main()

	# Grab input parameters and whether to print info to terminal
	params, CNmatFilename, birthFilename, initfile, outputFile, verbosity = initialize()

	date = today() # Get today's date to add to outputfile string

	# Printing stuff to terminal
	if verbosity
		println("Input:")
		for name in fieldnames(typeof(params))
			println("  $name => $(getfield(params,name))")
		end
	end

	# Read the birth rate file to be used for the interpolation
	X = readdlm(CNmatFilename, '\t')
	birthRates = readdlm(birthFilename, '\t')

	# FIXME!!!!!
	header, cn = X[1,:],Float64.(X[2:end,3:end])

	# Error if the elements in X or birthRates are not all subtypes of real
	if any((eltype(cn),eltype(birthRates)) .>: Real)
		error("cn and birthRates must contain only numeric values")
	end

	# readdlm gives a 2D array, we turn it into a vector (1d array) here.
	birthRates = dropdims(birthRates,dims=2)

	# Load initial condition from file if available:
	if !isfile(initfile)
		error("paralleldriver must have an initfile")
	else
		df = CSV.read(initfile,DataFrame)
	end

	# Grab the initial conditions and convert to integer for now
	u0mat = Int.(round.(convert(Matrix,df[:,[name for name in names(df)
	if occursin(r"Chr*",name)]])))

	# check if we update either of the two rates
	changeDeathRate = "deathRate" in names(df) ? true : false
	changeMisRate = "missegregationRate" in names(df) ? true : false

	if verbosity
		println("Running simulation in parallel with $(nworkers()) workers")
	end

	# Parallelization of the initial condition file (u0mat)
	pmap(enumerate(eachrow(u0mat))) do (i,initCN)

		# Define the outputfile string
		outputFile = "$(initfile)_$(date)_run-$i.csv"
		outputParams = "$(initfile)_$(date)_run-$(i)_params.csv"

		# Check to see if we replace the death or misRate from the initfile
		if changeDeathRate
			params.deathRate = df[i,"deathRate"]
			if verbosity
				println("new death rate: $(params.deathRate)")
			end
		end
		if changeMisRate
			params.misRate = df[i,"missegregationRate"]
			if verbosity
				println("new missegregation rate: $(params.misRate)")
			end
		end

		# run the simulation
		results, time = runPloidyMovement(params,cn,birthRates,initCN)

		# create header for the solution output
		outputHeader = permutedims(
			vcat(
				["Chr$chr" for chr in header if typeof(chr) == Int],
				time
				)
				)

		# concatnate the header with the results array
		output = vcat(outputHeader,results)

		# save to file
		writedlm( outputFile,  output)
		CSV.write( outputParams, [params])

	end

end

main()

max_stars_repo_name: MathOnco/ploidyEvolution

max_issues_repo_name: MathOnco/ploidyEvolution

max_forks_repo_name: MathOnco/ploidyEvolution


text:
function UpsampleBlock(inc, scale = 2)
    return Chain(
        Upsample(:nearest, scale = (scale, scale)),
        Conv((3, 3), inc => inc, Base.Fix2(leakyrelu, 0.2), pad = 1)
    )
end

struct DenseResidualBlock
    residual_beta
    blocks
end

function DenseResidualBlock(inc;
                            c = 32, residual_beta = 0.2)
    blocks = []
    for i in 0:4
        in_channels = inc + c * i
        out_channels = i<=3 ? c : inc
        act = i <= 3 ? Base.Fix2(leakyrelu, 0.2) : identity
        push!(blocks, Conv((3, 3), in_channels => out_channels, act, pad = 1))
    end

    return DenseResidualBlock(residual_beta, blocks)
end

@functor DenseResidualBlock

function (m::DenseResidualBlock)(x)
    new_inputs = x
    local out, new_inputs
    for block in m.blocks
        out = block(new_inputs)
        new_inputs = cat(new_inputs, out, dims = 3)
    end
    return m.residual_beta * out + x
end

mutable struct ResidualinResidualDenseBlock
    residual_beta
    rrdb
end

@functor ResidualinResidualDenseBlock

function ResidualinResidualDenseBlock(inc;
              residual_beta = 0.2)
    rrdb = Chain([DenseResidualBlock(inc) for _ in 1:3]...)
    ResidualinResidualDenseBlock(residual_beta, rrdb)
end

(m::ResidualinResidualDenseBlock)(x) = m.rrdb(x) * m.residual_beta + x

function esrgan_generator(inc = 3, nc = 64, nb = 23)
    initial = Conv((3, 3), inc => nc, pad = 1)
    residuals = Chain([ResidualinResidualDenseBlock(nc) for _ in 1:nb]...)
    conv = Conv((3, 3), nc => nc, pad = 1)
    upsamples = Chain(UpsampleBlock(nc), UpsampleBlock(nc))
    final = Chain(
        Conv((3, 3), nc => nc, Base.Fix2(leakyrelu, 0.2), pad = 1),
        Conv((3, 3), nc => inc, pad = 1)
    )
    return Chain(initial, SkipConnection(Chain(conv, residuals), +), upsamples, final)
end

function esrgan_discriminator(; in_c = 3, features = [64, 64, 128, 128, 256, 256, 512, 512])
    blocks = []
    for (idx, feature) in enumerate(features)
        s = 1 + (idx - 1) % 2
        push!(blocks, Conv((3, 3), in_c => feature, Base.Fix2(leakyrelu, 0.2), stride = s, pad = 1))
        in_c = feature
    end
    blocks = Chain(blocks...)
    classifier = Chain(
        AdaptiveMeanPool((6, 6)),
        MLUtils.flatten,
        Dense(512 * 6 * 6, 1024, Base.Fix2(leakyrelu, 0.2)),
        Dense(1024, 1)
    )
    return Chain(blocks, classifier)
end

function esrgan()
    return Chain(discriminator = esrgan_discriminator(), generator = esrgan_generator())
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Rotations from an Earth-Fixed, Earth-Centered (ECEF) reference frame to an
#   Earth-Fixed Inertial (ECI) reference frame.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export r_ecef_to_eci

"""
    r_ecef_to_eci([T,] ECEF, ECI, JD_UTC::Number [, eop_data])

Compute the rotation from an Earth-Centered, Earth-Fixed (`ECEF`) reference
frame to an Earth-Centered Inertial (`ECI`) reference frame at the Julian Day
[UTC] `JD_UTC`. The rotation description that will be used is given by `T`,
which can be `DCM` or `Quaternion`. The ECEF frame is selected by the input
`ECEF` and the `ECI` frame is selected by the input `ECI`. The possible values
are listed below. The model used to compute the rotation is specified by the
selection of the origin and destination frames. Currently, there are two models
supported: IAU-76/FK5 and IAU-2006 with 2010 conventions (CIO and equinox
approaches).

# Rotation description

The rotations that aligns the ECEF with ECI can be described by Direction Cosine
Matrices or Quaternions. This is selected by the parameter `T`. The possible
values are:

- `DCM`: The rotation will be described by a Direction Cosine Matrix.
- `Quaternion`: The rotation will be described by a Quaternion.

If no value is specified, then it falls back to `DCM`.

# Conversion model

The model that will be used to compute the rotation is automatically inferred
given the selection of the origin and destination frames. **Notice that mixing
IAU-76/FK5 and IAU-2006/2010 frames is not supported.**

# ECEF Frame

The ECEF frame is selected by the parameter `ECEF`. The possible values are:

- `ITRF()`: ECEF will be selected as the International Terrestrial Reference
    Frame (ITRF).
- `PEF()`: ECEF will be selected as the Pseudo-Earth Fixed (PEF) reference
    frame.
- `TIRS()`: ECEF will be selected as the Terrestrial Intermediate Reference
    System (TIRS).

# ECI Frame

The ECI frame is selected by the parameter `ECI`. The possible values are:

- `TEME()`: ECI will be selected as the True Equator Mean Equinox (TEME)
    reference frame.
- `TOD()`: ECI will be selected as the True of Date (TOD).
- `MOD()`: ECI will be selected as the Mean of Date (MOD).
- `J2000()`: ECI will be selected as the J2000 reference frame.
- `GCRF()`: ECI will be selected as the Geocentric Celestial Reference Frame
    (GCRF).
- `CIRS()`: ECI will be selected as the Celestial Intermediate Reference System
    (CIRS).
- `ERS()`: ECI will be selected as the Earth Reference System (ERS).
- `MOD06()`: ECI will be selected as the Mean of Date (MOD) according to the
    definition in IAU-2006/2010 theory.
- `MJ2000()`: ECI will be selected as the J2000 mean equatorial frame (MJ2000).

!!! note
    The frames `MOD()` and `MOD06()` are virtually the same. However, we
    selected different names to make clear which theory are being used since
    mixing transformation between frames from IAU-76/FK5 and IAU-2006/2010 must
    be performed with caution.

# EOP Data

The conversion between the frames depends on EOP Data (see
[`get_iers_eop`](@ref) and [`read_iers_eop`](@ref)). If IAU-76/FK5 model is
used, then the type of `eop_data` must be [`EOPData_IAU1980`](@ref). Otherwise,
if IAU-2006/2010 model is used, then the type of `eop_data` must be
[`EOPData_IAU2000A`](@ref). The following table shows the requirements for EOP
data given the selected frames.

|   Model                     |  ECEF  |   ECI    |    EOP Data     |
|:----------------------------|:-------|:---------|:----------------|
| IAU-76/FK5                  | `ITRF` | `GCRF`   | EOP IAU1980     |
| IAU-76/FK5                  | `ITRF` | `J2000`  | EOP IAU1980     |
| IAU-76/FK5                  | `ITRF` | `MOD`    | EOP IAU1980     |
| IAU-76/FK5                  | `ITRF` | `TOD`    | EOP IAU1980     |
| IAU-76/FK5                  | `ITRF` | `TEME`   | EOP IAU1980     |
| IAU-76/FK5                  | `PEF`  | `GCRF`   | EOP IAU1980     |
| IAU-76/FK5                  | `PEF`  | `J2000`  | Not required¹   |
| IAU-76/FK5                  | `PEF`  | `MOD`    | Not required¹   |
| IAU-76/FK5                  | `PEF`  | `TOD`    | Not required¹   |
| IAU-76/FK5                  | `PEF`  | `TEME`   | Not required¹   |
| IAU-2006/2010 CIO-based     | `ITRF` | `CIRS`   | EOP IAU2000A    |
| IAU-2006/2010 CIO-based     | `ITRF` | `GCRF`   | EOP IAU2000A    |
| IAU-2006/2010 CIO-based     | `TIRS` | `CIRS`   | Not required¹   |
| IAU-2006/2010 CIO-based     | `TIRS` | `GCRF`   | Not required¹ ² |
| IAU-2006/2010 Equinox-based | `ITRF` | `ERS`    | EOP IAU2000A    |
| IAU-2006/2010 Equinox-based | `ITRF` | `MOD06`  | EOP IAU2000A    |
| IAU-2006/2010 Equinox-based | `ITRF` | `MJ2000` | EOP IAU2000A    |
| IAU-2006/2010 Equinox-based | `TIRS` | `ERS`    | Not required¹ ³ |
| IAU-2006/2010 Equinox-based | `TIRS` | `MOD06`  | Not required¹ ³ |
| IAU-2006/2010 Equinox-based | `TIRS` | `MJ2000` | Not required¹ ³ |

`¹`: In this case, UTC will be assumed equal to UT1 to compute the Greenwich
Mean Sidereal Time. This is an approximation, but should be sufficiently
accurate for some applications. Notice that, if EOP Data is provided, UT1 will
be accurately computed.

`²`: In this case, the terms that account for the free core nutation and time
dependent effects of the Celestial Intermediate Pole (CIP) position with respect
to the GCRF will not be available, reducing the precision.

`³`: In this case, the terms that corrects the nutation in obliquity and in
longitude due to the free core nutation will not be available, reducing the
precision.

## MOD and TOD

In this function, if EOP corrections are not provided, then MOD and TOD frames
will be computed considering the original IAU-76/FK5 theory. Otherwise, the
corrected frame will be used.

# Returns

The rotation description represented by `T` that rotates the ECEF reference
frame into alignment with the ECI reference frame.

# Examples

```julia-repl
julia> eop_IAU1980 = get_iers_eop(Val(:IAU1980));

julia> r_ecef_to_eci(DCM, ITRF(), GCRF(), date_to_jd(1986, 06, 19, 21, 35, 0), eop_IAU1980)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.619267      0.78518     -0.00132979
 -0.78518      -0.619267     3.33509e-5
 -0.000797312   0.00106478   0.999999

julia> r_ecef_to_eci(ITRF(), GCRF(), date_to_jd(1986, 06, 19, 21, 35, 0), eop_IAU1980)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.619267      0.78518     -0.00132979
 -0.78518      -0.619267     3.33509e-5
 -0.000797312   0.00106478   0.999999

julia> r_ecef_to_eci(PEF(), J2000(), date_to_jd(1986, 06, 19, 21, 35, 0))
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.619271      0.785176    -0.00133066
 -0.785177     -0.619272     3.45854e-5
 -0.000796885   0.00106622   0.999999

julia> r_ecef_to_eci(PEF(), J2000(), date_to_jd(1986, 06, 19, 21, 35, 0), eop_IAU1980)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.619267      0.78518     -0.00133066
 -0.78518      -0.619267     3.45854e-5
 -0.000796879   0.00106623   0.999999

julia> r_ecef_to_eci(Quaternion, ITRF(), GCRF(), date_to_jd(1986, 06, 19, 21, 35, 0), eop_IAU1980)
Quaternion{Float64}:
  + 0.43631 - 0.000590997⋅i + 0.000305106⋅j + 0.000305106⋅k

julia> eop_IAU2000A = get_iers_eop(Val(:IAU2000A));

julia> r_ecef_to_eci(ITRF(), GCRF(), date_to_jd(1986, 06, 19, 21, 35, 0), eop_IAU2000A)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.619267      0.78518     -0.00132979
 -0.78518      -0.619267     3.33516e-5
 -0.000797311   0.00106478   0.999999

julia> r_ecef_to_eci(TIRS(), GCRF(), date_to_jd(1986, 06, 19, 21, 35, 0))
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.619271      0.785176    -0.00133066
 -0.785177     -0.619272     3.45884e-5
 -0.000796885   0.00106623   0.999999

julia> r_ecef_to_eci(Quaternion, ITRF(), GCRF(), date_to_jd(1986, 06, 19, 21, 35, 0), eop_IAU2000A)
Quaternion{Float64}:
  + 0.43631 - 0.000590997⋅i + 0.000305106⋅j + 0.000305106⋅k
```
"""
@inline function r_ecef_to_eci(
    T_ECEF::T_ECEFs,
    T_ECI::T_ECIs,
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    r_ecef_to_eci(DCM, T_ECEF, T_ECI, JD_UTC, eop_data)
end

@inline function r_ecef_to_eci(
    T_ECEF::T_ECEFs_IAU_2006,
    T_ECI::T_ECIs_IAU_2006,
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    r_ecef_to_eci(DCM, T_ECEF, T_ECI, JD_UTC, eop_data)
end

# Specializations for those cases that EOP Data is not needed.
@inline function r_ecef_to_eci(
    T_ECEF::Val{:PEF},
    T_ECI::Union{Val{:J2000}, Val{:MOD}, Val{:TOD}, Val{:TEME}},
    JD_UTC::Number
)
    r_ecef_to_eci(DCM, T_ECEF, T_ECI, JD_UTC)
end

@inline function r_ecef_to_eci(
        T_ECEF::Val{:TIRS},
        T_ECI::T_ECIs_IAU_2006,
        JD_UTC::Number
)
    r_ecef_to_eci(DCM, T_ECEF, T_ECI, JD_UTC)
end

################################################################################
#                                  IAU-76/FK5
################################################################################

#                                 ITRF => GCRF
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:GCRF},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    x_p      = eop_data.x(JD_UTC) * arcsec_to_rad
    y_p      = eop_data.y(JD_UTC) * arcsec_to_rad
    δΔϵ_1980 = eop_data.dEps(JD_UTC) * milliarcsec_to_rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC) * milliarcsec_to_rad

    # Compute the rotation.
    return r_itrf_to_gcrf_fk5(T, JD_UT1, JD_TT, x_p, y_p, δΔϵ_1980, δΔψ_1980)
end

#                                ITRF => J2000
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:J2000},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    arcsec_to_rad = π / 648000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    x_p = eop_data.x(JD_UTC) * arcsec_to_rad
    y_p = eop_data.y(JD_UTC) * arcsec_to_rad

    # Compute the rotation.
    return r_itrf_to_gcrf_fk5(T, JD_UT1, JD_TT, x_p, y_p, 0, 0)
end

#                                 ITRF => MOD
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:MOD},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    x_p      = eop_data.x(JD_UTC) * arcsec_to_rad
    y_p      = eop_data.y(JD_UTC) * arcsec_to_rad
    δΔϵ_1980 = eop_data.dEps(JD_UTC) * milliarcsec_to_rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC) * milliarcsec_to_rad

    # Compute the rotation.
    r_PEF_ITRF = r_itrf_to_pef_fk5(T, x_p, y_p)
    r_MOD_PEF  = r_pef_to_mod_fk5(T, JD_UT1, JD_TT, δΔϵ_1980, δΔψ_1980)

    return compose_rotation(r_PEF_ITRF, r_MOD_PEF)
end

#                                 ITRF => TOD
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:TOD},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    x_p      = eop_data.x(JD_UTC) * arcsec_to_rad
    y_p      = eop_data.y(JD_UTC) * arcsec_to_rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC) * milliarcsec_to_rad

    # Compute the rotation.
    r_PEF_ITRF = r_itrf_to_pef_fk5(T, x_p, y_p)
    r_TOD_PEF  = r_pef_to_tod_fk5(T, JD_UT1, JD_TT, δΔψ_1980)

    return compose_rotation(r_PEF_ITRF, r_TOD_PEF)
end

#                                 ITRF => TEME
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:TEME},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    arcsec_to_rad = π / 648000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    x_p = eop_data.x(JD_UTC) * arcsec_to_rad
    y_p = eop_data.y(JD_UTC) * arcsec_to_rad

    # Compute the rotation.
    r_PEF_ITRF = r_itrf_to_pef_fk5(T, x_p, y_p)
    r_TEME_PEF = r_pef_to_teme(T, JD_UT1)

    return compose_rotation(r_PEF_ITRF, r_TEME_PEF)
end

#                                 PEF => GCRF
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:PEF},
    ::Val{:GCRF},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    x_p      = eop_data.x(JD_UTC) * arcsec_to_rad
    y_p      = eop_data.y(JD_UTC) * arcsec_to_rad
    δΔϵ_1980 = eop_data.dEps(JD_UTC) * milliarcsec_to_rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC) * milliarcsec_to_rad

    # Compute the rotation.
    r_MOD_PEF  = r_pef_to_mod_fk5(T, JD_UT1, JD_TT, δΔϵ_1980, δΔψ_1980)
    r_GCRF_MOD = r_mod_to_gcrf_fk5(T, JD_TT)

    return compose_rotation(r_MOD_PEF, r_GCRF_MOD)
end

#                                PEF => J2000
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:PEF},
    ::Val{:J2000},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Compute the rotation.
    r_MOD_PEF  = r_pef_to_mod_fk5(T, JD_UT1, JD_TT, 0, 0)
    r_GCRF_MOD = r_mod_to_gcrf_fk5(T, JD_TT)

    return compose_rotation(r_MOD_PEF, r_GCRF_MOD)
end

function r_ecef_to_eci(T::T_ROT, ::Val{:PEF}, ::Val{:J2000}, JD_UTC::Number)
    # Since we do not have EOP Data, assume that JD_UTC is equal to JD_UT1.
    JD_UT1 = JD_UTC
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Compute the rotation.
    r_MOD_PEF  = r_pef_to_mod_fk5(T, JD_UT1, JD_TT, 0, 0)
    r_GCRF_MOD = r_mod_to_gcrf_fk5(T, JD_TT)

    return compose_rotation(r_MOD_PEF, r_GCRF_MOD)
end

#                                 PEF => MOD
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:PEF},
    ::Val{:MOD},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    x_p      = eop_data.x(JD_UTC) * arcsec_to_rad
    y_p      = eop_data.y(JD_UTC) * arcsec_to_rad
    δΔϵ_1980 = eop_data.dEps(JD_UTC) * milliarcsec_to_rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC) * milliarcsec_to_rad

    # Compute the rotation.
    return r_pef_to_mod_fk5(T, JD_UT1, JD_TT, δΔϵ_1980, δΔψ_1980)
end

function r_ecef_to_eci(T::T_ROT, ::Val{:PEF}, ::Val{:MOD}, JD_UTC::Number)
    # Since we do not have EOP Data, assume that JD_UTC is equal to JD_UT1.
    JD_UT1 = JD_UTC
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Compute the rotation.
    return r_pef_to_mod_fk5(T, JD_UT1, JD_TT, 0, 0)
end

#                                 PEF => TOD
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:PEF},
    ::Val{:TOD},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    x_p      = eop_data.x(JD_UTC) * arcsec_to_rad
    y_p      = eop_data.y(JD_UTC) * arcsec_to_rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC) * milliarcsec_to_rad

    # Compute the rotation.
    return r_pef_to_tod_fk5(T, JD_UT1, JD_TT, δΔψ_1980)
end

function r_ecef_to_eci(T::T_ROT, ::Val{:PEF}, ::Val{:TOD}, JD_UTC::Number)
    # Since we do not have EOP Data, assume that JD_UTC is equal to JD_UT1.
    JD_UT1 = JD_UTC
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Compute the rotation.
    return r_pef_to_tod_fk5(T, JD_UT1, JD_TT, 0)
end

#                                 PEF => TEME
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:PEF},
    ::Val{:TEME},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    # Get the time in UT1.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)

    # Compute the rotation.
    return r_pef_to_teme(T, JD_UT1)
end

function r_ecef_to_eci(T::T_ROT, ::Val{:PEF}, ::Val{:TEME}, JD_UTC::Number)
    # Since we do not have EOP Data, assume that JD_UTC is equal to JD_UT1.
    JD_UT1 = JD_UTC

    # Compute the rotation.
    return r_pef_to_teme(T, JD_UT1)
end

################################################################################
#                           IAU-2006/2010 CIO-based
################################################################################

#                                 ITRF => CIRS
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:CIRS},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    arcsec_to_rad = π / 648000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    x_p = eop_data.x(JD_UTC) * arcsec_to_rad
    y_p = eop_data.y(JD_UTC) * arcsec_to_rad

    # Compute the rotation.
    r_TIRS_ITRF = r_itrf_to_tirs_iau2006(T, JD_TT, x_p, y_p)
    r_CIRS_TIRS = r_tirs_to_cirs_iau2006(T, JD_UT1)

    return compose_rotation(r_TIRS_ITRF, r_CIRS_TIRS)
end

#                                 ITRF => GCRF
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:GCRF},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    x_p = eop_data.x(JD_UTC) * arcsec_to_rad
    y_p = eop_data.y(JD_UTC) * arcsec_to_rad
    dX  = eop_data.dX(JD_UTC) * milliarcsec_to_rad
    dY  = eop_data.dY(JD_UTC) * milliarcsec_to_rad

    # Compute the rotation.
    r_TIRS_ITRF = r_itrf_to_tirs_iau2006(T, JD_TT, x_p, y_p)
    r_CIRS_TIRS = r_tirs_to_cirs_iau2006(T, JD_UT1)
    r_GCRF_CIRS = r_cirs_to_gcrf_iau2006(T, JD_TT, dX, dY)

    return compose_rotation(r_TIRS_ITRF, r_CIRS_TIRS, r_GCRF_CIRS)
end

#                                 TIRS => CIRS
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:TIRS},
    ::Val{:CIRS},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)

    # Compute the rotation.
    return r_tirs_to_cirs_iau2006(T, JD_UT1)
end

function r_ecef_to_eci(T::T_ROT, ::Val{:TIRS}, ::Val{:CIRS}, JD_UTC::Number)
    # Since we do not have EOP Data, assume that JD_UTC is equal to JD_UT1.
    JD_UT1 = JD_UTC

    # Compute the rotation.
    return r_tirs_to_cirs_iau2006(T, JD_UT1)
end

#                                 TIRS => GCRF
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:TIRS},
    ::Val{:GCRF},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    milliarcsec_to_rad = π / 648000000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    dX = eop_data.dX(JD_UTC) * milliarcsec_to_rad
    dY = eop_data.dY(JD_UTC) * milliarcsec_to_rad

    # Compute the rotation.
    r_CIRS_TIRS = r_tirs_to_cirs_iau2006(T, JD_UT1)
    r_GCRF_CIRS = r_cirs_to_gcrf_iau2006(T, JD_TT, dX, dY)

    return compose_rotation(r_CIRS_TIRS, r_GCRF_CIRS)
end

function r_ecef_to_eci(T::T_ROT, ::Val{:TIRS}, ::Val{:GCRF}, JD_UTC::Number)
    # Since we do not have EOP Data, assume that JD_UTC is equal to JD_UT1.
    JD_UT1 = JD_UTC

    # Get the time in TT.
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Compute the rotation.
    r_CIRS_TIRS = r_tirs_to_cirs_iau2006(T, JD_UT1)
    r_GCRF_CIRS = r_cirs_to_gcrf_iau2006(T, JD_TT)

    return compose_rotation(r_CIRS_TIRS, r_GCRF_CIRS)
end

################################################################################
#                         IAU-2006/2010 equinox-based
################################################################################

#                                 ITRF => ERS
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:ERS},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    x_p = eop_data.x(JD_UTC) * arcsec_to_rad
    y_p = eop_data.y(JD_UTC) * arcsec_to_rad

    # Obtain the correction of the nutation in longitude.
    ~, δΔΨ_2000 = deps_dpsi(eop_data, JD_UTC)
    δΔΨ_2000 *= milliarcsec_to_rad

    # Compute the rotation.
    r_TIRS_ITRF = r_itrf_to_tirs_iau2006(T, JD_TT, x_p, y_p)
    r_ERS_TIRS = r_tirs_to_ers_iau2006(T, JD_UT1, JD_TT, δΔΨ_2000)

    return compose_rotation(r_TIRS_ITRF, r_ERS_TIRS)
end

#                                 ITRF => MOD
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:MOD06},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    x_p = eop_data.x(JD_UTC) * arcsec_to_rad
    y_p = eop_data.y(JD_UTC) * arcsec_to_rad

    # Obtain the correction of the nutation in obliquity and longitude.
    δΔϵ_2000, δΔΨ_2000 = deps_dpsi(eop_data, JD_UTC)
    δΔϵ_2000 *= milliarcsec_to_rad
    δΔΨ_2000 *= milliarcsec_to_rad

    # Compute the rotation.
    r_TIRS_ITRF = r_itrf_to_tirs_iau2006(T, JD_TT, x_p, y_p)
    r_MOD_TIRS = r_tirs_to_mod_iau2006(T, JD_UT1, JD_TT, δΔϵ_2000, δΔΨ_2000)

    return compose_rotation(r_TIRS_ITRF, r_MOD_TIRS)
end

#                                ITRF => MJ2000
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:ITRF},
    ::Val{:MJ2000},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    arcsec_to_rad = π / 648000
    milliarcsec_to_rad = arcsec_to_rad / 1000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    x_p = eop_data.x(JD_UTC) * arcsec_to_rad
    y_p = eop_data.y(JD_UTC) * arcsec_to_rad

    # Obtain the correction of the nutation in obliquity and longitude.
    δΔϵ_2000, δΔΨ_2000 = deps_dpsi(eop_data, JD_UTC)
    δΔϵ_2000 *= milliarcsec_to_rad
    δΔΨ_2000 *= milliarcsec_to_rad

    # Compute the rotation.
    r_TIRS_ITRF  = r_itrf_to_tirs_iau2006(T, JD_TT, x_p, y_p)
    r_MOD_TIRS   = r_tirs_to_mod_iau2006(T, JD_UT1, JD_TT, δΔϵ_2000, δΔΨ_2000)
    r_MJ2000_MOD = r_mod_to_mj2000_iau2006(T, JD_TT)

    return compose_rotation(r_TIRS_ITRF, r_MOD_TIRS, r_MJ2000_MOD)
end

# NOTE: We do not implement the conversion ITRF => GCRF using equinox-based
# IAU-2006/2010 theory because the CIO-based approach is faster and more
# precise.

#                                  TIRS => ERS
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:TIRS},
    ::Val{:ERS},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    milliarcsec_to_rad = π / 648000000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Obtain the correction of the nutation in longitude.
    ~, δΔΨ_2000 = deps_dpsi(eop_data, JD_UTC)
    δΔΨ_2000 *= milliarcsec_to_rad

    # Compute the rotation.
    return r_tirs_to_ers_iau2006(T, JD_UT1, JD_TT, δΔΨ_2000)
end

function r_ecef_to_eci(T::T_ROT, ::Val{:TIRS}, ::Val{:ERS}, JD_UTC::Number)
    # Since we do not have EOP Data, assume that JD_UTC is equal to JD_UT1.
    JD_UT1 = JD_UTC
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Compute the rotation.
    return r_tirs_to_ers_iau2006(T, JD_UT1, JD_TT)
end

#                                 TIRS => MOD
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:TIRS},
    ::Val{:MOD06},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    milliarcsec_to_rad = π / 648000000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Obtain the correction of the nutation in obliquity and longitude.
    δΔϵ_2000, δΔΨ_2000 = deps_dpsi(eop_data, JD_UTC)
    δΔϵ_2000 *= milliarcsec_to_rad
    δΔΨ_2000 *= milliarcsec_to_rad

    # Compute the rotation.
    return r_tirs_to_mod_iau2006(T, JD_UT1, JD_TT, δΔϵ_2000, δΔΨ_2000)
end

function r_ecef_to_eci(T::T_ROT, ::Val{:TIRS}, ::Val{:MOD06}, JD_UTC::Number)
    # Since we do not have EOP Data, assume that JD_UTC is equal to JD_UT1.
    JD_UT1 = JD_UTC
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Compute the rotation.
    return r_tirs_to_mod_iau2006(T, JD_UT1, JD_TT)
end

#                                TIRS => MJ2000
# ==============================================================================

function r_ecef_to_eci(
    T::T_ROT,
    ::Val{:TIRS},
    ::Val{:MJ2000},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    milliarcsec_to_rad = π / 648000000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Obtain the correction of the nutation in obliquity and longitude.
    δΔϵ_2000, δΔΨ_2000 = deps_dpsi(eop_data, JD_UTC)
    δΔϵ_2000 *= milliarcsec_to_rad
    δΔΨ_2000 *= milliarcsec_to_rad

    # Compute the rotation.
    r_MOD_TIRS   = r_tirs_to_mod_iau2006(T, JD_UT1, JD_TT, δΔϵ_2000, δΔΨ_2000)
    r_MJ2000_MOD = r_mod_to_mj2000_iau2006(T, JD_TT)

    return compose_rotation(r_MOD_TIRS, r_MJ2000_MOD)
end

function r_ecef_to_eci(T::T_ROT, ::Val{:TIRS}, ::Val{:MJ2000}, JD_UTC::Number)
    # Since we do not have EOP Data, assume that JD_UTC is equal to JD_UT1.
    JD_UT1 = JD_UTC
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Compute the rotation.
    r_MOD_TIRS   = r_tirs_to_mod_iau2006(T, JD_UT1, JD_TT)
    r_MJ2000_MOD = r_mod_to_mj2000_iau2006(T, JD_TT)

    return compose_rotation(r_MOD_TIRS, r_MJ2000_MOD)
end

# NOTE: We do not implement the conversion TIRS => GCRF using equinox-based
# IAU-2006/2010 theory because the CIO-based approach is faster and more
# precise.


text:
include("turbines.jl")
import YAML
import Base.Filesystem

abstract type AbstractFileType end

struct InputReader
    data_path:: String
    read_path:: String
    write_path:: String
end

struct YAMLFileType <: AbstractFileType
end

function read_file()
    data = YAML.load_file(open(InputReader.read_path))
    println(data)
end

function write_to_file()
    YAML.write_file(InputReader.write_path, InputReader.data_path)
end

file_type = YAMLFileType
reader = InputReader("D:\\source\\WindFarmCode\\Julia\\FlowFarm.jl\\test\\format_file.data", "D:\\source\\WindFarmCode\\Julia\\FlowFarm.jl\\test\\wind_turbine.yml","D:\\source\\WindFarmCode\\Julia\\FlowFarm.jl\\test\\wind_turbine.yml")
w = write_to_file()
r = read_file()


text:
filename() = "day19.txt"

# wraps patterns in parentheses
datum(x::String) = "($x)"
datum(x::Vector{String}) = datum(join(datum.(x)))
datum(x::Vector{Vector{String}}) = datum(join(datum.(x), "|"))

"""
The rules are parsed into a Dict that looks like this:

julia> read_data("day19_sample.txt")[1]
Dict{String,String} with 6 entries:
  "(2)" => "(((4)(4))|((5)(5)))"
  "(5)" => "(b)"
  "(1)" => "(((2)(3))|((3)(2)))"
  "(3)" => "(((4)(5))|((5)(4)))"
  "(4)" => "(a)"
  "(0)" => "((4)(1)(5))"
"""
function read_data(file)
    blocks = split(String(read(file)), "\n\n")

    rules = Dict(let v = string.(split(r, r":* "))
                    k = v[1]
                    if startswith(v[2], "\"")  # final term
                        datum(k) => datum(v[2][2:end-1])
                    elseif "|" ∈ v[2:end]      # multiple choices
                        d = findfirst(==("|"), v)
                        datum(k) => datum([v[2:d-1], v[d+1:end]])
                    else                       # single choice
                        datum(k) => datum(v[2:end])
                    end
                 end for r in split(blocks[1], "\n"))

    messages = split(blocks[2], "\n")

    return rules, messages
end

"""
Given the rules, resolve the variables into terminal strings "a" and "b"
and then return a `Regex` object.
"""
function make_regex(rules)
    str = rules["(0)"]
    while match(r"\([0-9]+\)", str) !== nothing
        for (k,v) in rules
            str = replace(str, k => v)
        end
    end
    for i in 1:10  # trim parentheses to reduce size of regex
        str = replace(str, r"\(([ab]*)\)" => s"\1")
    end
    str = "^" * str * "\$"
    # @show str
    return Regex(str)
end

function part1()
    rules, messages = read_data(filename())
    re = make_regex(rules)
    return count(m -> !isnothing(match(re, m)), messages)
end

function part2()
    rules, messages = read_data(filename())
    rules["(8)"] = "((42)+)"
    return sum(1:8) do i
        rules["(11)"] = "((42){$i}(31){$i})"
        re = make_regex(rules)
        count(m -> !isnothing(match(re, m)), messages)
    end
end

# For some reasons, regex recursion did not work: "((42)(?R)?(31))"
function part2b()
    rules, messages = read_data(filename())
    rules["(8)"] = "((42)+)"
    rules["(11)"] = "((42)(?1)?(31))"
    re = make_regex(rules)
    return count(m -> !isnothing(match(re, m)), messages)
end


text:
include("defs.jl")  # testing setup

my_tests = [
			"unit_tests.jl",
            "parsing.jl",
            "syntax.jl",
            "firstorder.jl",
            "indexing.jl",
            "types.jl",
            "loops.jl",
            "functions.jl",
            "more.jl"
           ]

println("Running tests:")

for my_test in my_tests
    println("  * $(my_test) *")
    include(my_test)
end

println("Finished")


text:
# T: coeff type
# Degree: degree
# Dim1: dimension of the support + 1
mutable struct LagrangeRefSpace{T,Degree,Dim1,NF} <: RefSpace{T,NF} end

numfunctions(s::LagrangeRefSpace{T,D,2}) where {T,D} = D+1
numfunctions(s::LagrangeRefSpace{T,0,3}) where {T} = 1
numfunctions(s::LagrangeRefSpace{T,1,3}) where {T} = 3

valuetype(ref::LagrangeRefSpace{T}, charttype) where {T} =
        SVector{numfunctions(ref), Tuple{T,T}}

# Evaluate constant lagrange elements on anything
(ϕ::LagrangeRefSpace{T,0})(tp) where {T} = SVector(((value=one(T), derivative=zero(T)),))

# Evaluate linear Lagrange elements on a segment
function (f::LagrangeRefSpace{T,1,2})(mp) where T
    u = mp.bary[1]
    j = jacobian(mp)
    SVector(
        (value=  u, derivative=-1/j),
        (value=1-u, derivative= 1/j))
end

# Evaluete linear lagrange elements on a triangle
function (f::LagrangeRefSpace{T,1,3})(t) where T
    u,v,w, = barycentric(t)
    SVector(
        (value=u,),
        (value=v,),
        (value=w,))
end

"""
    f(tangent_space, Val{:withcurl})

Compute the values of the shape functions together with their curl.
"""
function (f::LagrangeRefSpace{T,1,3})(t, ::Type{Val{:withcurl}}) where T
    # Evaluete linear Lagrange elements on a triange, together with their curl
    j = jacobian(t)
    u,v,w, = barycentric(t)
    p = t.patch
    SVector(
        (value=u, curl=(p[3]-p[2])/j),
        (value=v, curl=(p[1]-p[3])/j),
        (value=w, curl=(p[2]-p[1])/j)
    )
end

# Evaluate constant Lagrange elements on a triangle, with their curls
function (f::LagrangeRefSpace{T,0,3})(t, ::Type{Val{:withcurl}}) where T
    i = one(T)
    z = zero(cartesian(t))
    (
        (value=i, curl=z,),
    )
end

function curl(ref::LagrangeRefSpace, sh, el)
    sh1 = Shape(sh.cellid, mod1(sh.refid+1,3), -sh.coeff)
    sh2 = Shape(sh.cellid, mod1(sh.refid+2,3), +sh.coeff)
    return [sh1, sh2]
end

function gradient(ref::LagrangeRefSpace{T,1,4}, sh, tet) where {T}
    this_vert = tet.vertices[sh.refid]
    # other_verts = deleteat(tet.vertices, sh.refid)
    # opp_face = simplex(other_verts...)
    opp_face = faces(tet)[sh.refid]
    ctr_opp_face = center(opp_face)
    n = normal(ctr_opp_face)
    h = -dot(this_vert - cartesian(ctr_opp_face), n)
    @assert h > 0
    gradval = -h*n
    output = Vector{Shape{T}}()
    for (i,edge) in enumerate(CompScienceMeshes.edges(tet))
        ctr_edge = center(edge)
        tgt = tangents(ctr_edge,1)
        tgt = normalize(tgt)
        lgt = volume(edge)
        cff = -lgt * dot(tgt, gradval)
        isapprox(cff, 0, atol=sqrt(eps(T))) && continue
        push!(output, Shape(sh.cellid, i, cff))
    end
    return output
end

function strace(x::LagrangeRefSpace, cell, localid, face)

    Q = zeros(scalartype(x),2,3)

    p1 = neighborhood(face, 1)
    p2 = neighborhood(face, 0)

    u1 = carttobary(cell, cartesian(p1))
    u2 = carttobary(cell, cartesian(p2))

    P1 = neighborhood(cell, u1)
    P2 = neighborhood(cell, u2)

    vals1 = x(P1)
    vals2 = x(P2)

    for j in 1:numfunctions(x)
        Q[1,j] = vals1[j].value
        Q[2,j] = vals2[j].value
    end

    Q
end

function restrict(refs::LagrangeRefSpace{T,0}, dom1, dom2) where T
    #Q = eye(T, numfunctions(refs))
    Q = Matrix{T}(I, numfunctions(refs), numfunctions(refs))
end

function restrict(f::LagrangeRefSpace{T,1}, dom1, dom2) where T

    D = numfunctions(f)
    Q = zeros(T, D, D)

    # for each point of the new domain
    for i in 1:D
        v = dom2.vertices[i]

        # find the barycentric coordinates in dom1
        uvn = carttobary(dom1, v)

        # evaluate the shape functions in this point
        x = neighborhood(dom1, uvn)
        fx = f(x)

        for j in 1:D
            Q[j,i] = fx[j][1]
        end
    end

    return Q
end