text:
using FluxRC, KitBase, Plots

cd(@__DIR__)
ps = UnstructPSpace("square.msh")

N = 2
Np = (N + 1) * (N + 2) ÷ 2
pl, wl = tri_quadrature(N)
pf, wf = triface_quadrature(N)

spg = zeros(size(ps.cellid, 1), Np, 2)
for i in axes(spg, 1), j in axes(spg, 2)
    id1, id2, id3 = ps.cellid[i, :]
    spg[i, j, :] .=
        rs_xy(pl[j, :], ps.points[id1, 1:2], ps.points[id2, 1:2], ps.points[id3, 1:2])
end

nface = size(ps.faceType, 1)
fpg = zeros(size(ps.faceType, 1), N + 1, 2)
for i in axes(fpg, 1), j in axes(fpg, 2)
    idc = ifelse(ps.faceCells[i, 1] != -1, ps.faceCells[i, 1], ps.faceCells[i, 2])
    id1, id2, id3 = ps.cellid[idc, :]

    if !(id3 in ps.facePoints[i, :])
        idf = 1
    elseif !(id1 in ps.facePoints[i, :])
        idf = 2
    elseif !(id2 in ps.facePoints[i, :])
        idf = 3
    end

    fpg[i, j, :] .=
        rs_xy(pf[idf, j, :], ps.points[id1, 1:2], ps.points[id2, 1:2], ps.points[id3, 1:2])
end

a = 1.0
u = zeros(size(ps.cellid, 1), Np)
for i in axes(u, 1), j in axes(u, 2)
    u[i, j] = exp(-300 * ((spg[i, j, 1] - 0.5)^2 + (spg[i, j, 2] - 0.5)^2))
end

f = zeros(size(ps.cellid, 1), Np, 2)
for i in axes(f, 1)
    xr, yr = ps.points[ps.cellid[i, 2], 1:2] - ps.points[ps.cellid[i, 1], 1:2]
    xs, ys = ps.points[ps.cellid[i, 3], 1:2] - ps.points[ps.cellid[i, 1], 1:2]
    J = xr * ys - xs * yr
    for j in axes(f, 2)
        fg = a * u[i, j]
        gg = a * u[i, j]
        f[i, j, :] .= [ys * fg - xs * gg, -yr * fg + xr * gg] ./ J
    end
end

V = vandermonde_matrix(N, pl[:, 1], pl[:, 2])
Vr, Vs = ∂vandermonde_matrix(N, pl[:, 1], pl[:, 2])

∂l = zeros(Np, Np, 2)
for i = 1:Np
    ∂l[i, :, 1] .= V' \ Vr[i, :]
    ∂l[i, :, 2] .= V' \ Vs[i, :]
end

du = zeros(size(ps.cellid, 1), Np)
for i in axes(du, 1), j in axes(du, 2)
    du[i, j] = -sum(f[i, :, 1] .* ∂l[j, :, 1]) - sum(f[i, :, 2] .* ∂l[j, :, 2])
end

ψf = zeros(3, N + 1, Np)
for i = 1:3
    ψf[i, :, :] .= vandermonde_matrix(N, pf[i, :, 1], pf[i, :, 2])
end

lf = zeros(3, N + 1, Np)
for i = 1:3, j = 1:N+1
    lf[i, j, :] .= V' \ ψf[i, j, :]
end

# correction field
ϕ = zeros(3, N + 1, Np)
σ = zeros(3, N + 1, Np)

for k = 1:Np
    for j = 1:N+1
        for i = 1:3
            σ[i, j, k] = wf[i, j] * ψf[i, j, k]
        end
    end
end

for f = 1:3, j = 1:N+1, i = 1:Np
    ϕ[f, j, i] = sum(σ[f, j, :] .* V[i, :])
end

function dudt!(du, u, p, t)
    du .= 0.0

    a, deg, nface = p

    ncell = size(u, 1)
    nsp = size(u, 2)

    f = zeros(ncell, nsp, 2)
    for i in axes(f, 1)
        xr, yr = ps.points[ps.cellid[i, 2], 1:2] - ps.points[ps.cellid[i, 1], 1:2]
        xs, ys = ps.points[ps.cellid[i, 3], 1:2] - ps.points[ps.cellid[i, 1], 1:2]
        J = xr * ys - xs * yr
        for j in axes(f, 2)
            fg = a * u[i, j]
            gg = a * u[i, j]
            f[i, j, :] .= [ys * fg - xs * gg, -yr * fg + xr * gg] ./ J
        end
    end

    u_face = zeros(ncell, 3, deg + 1)
    f_face = zeros(ncell, 3, deg + 1, 2)
    for i = 1:ncell, j = 1:3, k = 1:deg+1
        u_face[i, j, k] = sum(u[i, :] .* lf[j, k, :])
        f_face[i, j, k, 1] = sum(f[i, :, 1] .* lf[j, k, :])
        f_face[i, j, k, 2] = sum(f[i, :, 2] .* lf[j, k, :])
    end

    n = [[0.0, -1.0], [1 / √2, 1 / √2], [-1.0, 0.0]]

    fn_face = zeros(ncell, 3, deg + 1)
    for i = 1:ncell, j = 1:3, k = 1:deg+1
        fn_face[i, j, k] = sum(f_face[i, j, k, :] .* n[j])
    end

    au = zeros(nface, deg + 1, 2)
    f_interaction = zeros(nface, deg + 1, 2)
    for i = 1:nface
        c1, c2 = ps.faceCells[i, :]
        if -1 in (c1, c2)
            continue
        end

        pids = ps.facePoints[i, :]
        if ps.cellid[c1, 1] ∉ pids
            fid1 = 2
        elseif ps.cellid[c1, 2] ∉ pids
            fid1 = 3
        elseif ps.cellid[c1, 3] ∉ pids
            fid1 = 1
        end

        if ps.cellid[c2, 1] ∉ pids
            fid2 = 2
        elseif ps.cellid[c2, 2] ∉ pids
            fid2 = 3
        elseif ps.cellid[c2, 3] ∉ pids
            fid2 = 1
        end

        for j = 1:deg+1, k = 1:2
            au[i, j, k] =
                (f_face[c1, fid1, j, k] - f_face[c2, fid2, j, k]) /
                (u_face[c1, fid1, j] - u_face[c2, fid2, j] + 1e-6)

            f_interaction[i, j, k] =
                0.5 * (f_face[c1, fid1, j, k] + f_face[c2, fid2, j, k]) -
                0.5 * abs(au[i, j, k]) * (u_face[c1, fid1, j] - u_face[c2, fid2, j])
        end
    end

    fn_interaction = zeros(ncell, 3, deg + 1)
    for i = 1:ncell
        for j = 1:3, k = 1:deg+1
            fn_interaction[i, j, k] = sum(f_interaction[ps.cellFaces[i, j], k, :] .* n[j])
        end
    end

    rhs1 = zeros(ncell, nsp)
    for i in axes(rhs1, 1), j in axes(rhs1, 2)
        rhs1[i, j] = -sum(f[i, :, 1] .* ∂l[j, :, 1]) - sum(f[i, :, 2] .* ∂l[j, :, 2])
    end

    for i = 1:ncell
        if ps.cellType[i] != 1
            for j = 1:nsp
                du[i, j] =
                    rhs1[i, j] -
                    sum(fn_interaction[i, :, :] .- fn_face[i, :, :] .* ϕ[:, :, j])
            end
        end
    end

end

du = zero(u)
dudt!(du, u, (a, N, nface), 0.0)

setdiff([1, 2], [2])

ps.cellFaces[1, :]

ps.cellid[1, :]

ps.facePoints[3, :]

u[idx, :] .* lf[1, 2, :] |> sum
u[idx, :] .* lf[2, 2, :] |> sum
u[idx, :] .* lf[3, 2, :] |> sum

u[idx, :] .* ∂l[1, :, 2] |> sum
u[idx, :] .* ∂l[2, :, 2] |> sum
u[idx, :] .* ∂l[3, :, 2] |> sum
u[idx, :] .* ∂l[4, :, 2] |> sum
u[idx, :] .* ∂l[5, :, 2] |> sum
u[idx, :] .* ∂l[6, :, 2] |> sum

spg[idx, :, :]

u[idx, 4]
u[idx, 6]

idx = 468

id1, id2, id3 = ps.cellFaces[idx, :]
scatter(spg[idx, :, 1], spg[idx, :, 2], legend = :none)
scatter!(fpg[id1, :, 1], fpg[id1, :, 2])
scatter!(fpg[id2, :, 1], fpg[id2, :, 2])
scatter!(fpg[id3, :, 1], fpg[id3, :, 2])

write_vtk(ps.points, ps.cellid, u[:, 1])

scatter(spg[1:10, :, 1], spg[1:10, :, 2], legend = :none, ratio = 1 / 3)


text:
cd(@__DIR__); include("setups/simpletree.jl");
pyplot(dpi = 200)

## (a)
println("l = $(ieb1[end] - 1) and λₗ = $(𝛌[ieb1[end] - 1])")
plot(size = (200, 500), framestyle = :none)
gplot!(W, X, width = 1)
scatter_gplot!(X; marker = 𝚽[:, ieb1[end]])
plt = plot!(cbar = false, clims = (-0.3, 0.3))
savefig(plt, "../figs/SimpleTree_eigenvector$(ieb1[end]-1)_pink.png")

## (b)
println("l = $(ieb2[end] - 1) and λₗ = $(𝛌[ieb2[end] - 1])")
plot(size = (200, 500), framestyle = :none)
gplot!(W, X, width = 1)
scatter_gplot!(X; marker = 𝚽[:, ieb2[end]])
plt = plot!(cbar = false, clims = (-0.3, 0.3))
savefig(plt, "../figs/SimpleTree_eigenvector$(ieb2[end]-1)_orange.png")

## (c)
println("l = $(ieb3[end] - 1) and λₗ = $(𝛌[ieb3[end] - 1])")
plot(size = (200, 500), framestyle = :none)
gplot!(W, X, width = 1)
scatter_gplot!(X; marker = 𝚽[:, ieb3[end]])
plt = plot!(cbar = false, clims = (-0.3, 0.3))
savefig(plt, "../figs/SimpleTree_eigenvector$(ieb3[end]-1)_green.png")

## (d)
println("l = $(ieb4[end] - 1) and λₗ = $(𝛌[ieb4[end] - 1])")
plot(size = (200, 500), framestyle = :none)
gplot!(W, X, width = 1)
scatter_gplot!(X; marker = 𝚽[:, ieb4[end]])
plt = plot!(cbar = false, clims = (-0.3, 0.3))
savefig(plt, "../figs/SimpleTree_eigenvector$(ieb4[end]-1)_yellow.png")

ext: jl


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

JLLWrappers.@generate_main_file_header("boostpython")
JLLWrappers.@generate_main_file("boostpython", UUID("398de629-0a17-50a6-9837-8b3a70a53854"))
end  # module boostpython_jll


text:
module Oceananigans

if VERSION < v"1.1"
    @error "Oceananigans requires Julia v1.1 or newer."
end

export
    # Architectures
    CPU, GPU,

    # Logging
    ModelLogger, Diagnostic, Setup,

    # Grids
    Periodic, Bounded, Flat,
    RegularCartesianGrid, VerticallyStretchedCartesianGrid,

    # Boundary conditions
    BoundaryCondition,
    Flux, Value, Gradient,
    FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition,
    CoordinateBoundaryConditions, FieldBoundaryConditions,
    UVelocityBoundaryConditions, VVelocityBoundaryConditions, WVelocityBoundaryConditions,
    TracerBoundaryConditions, PressureBoundaryConditions,
    BoundaryFunction,

    # Fields and field manipulation
    Field, CellField, XFaceField, YFaceField, ZFaceField,
    interior, set!,

    # Forcing functions
    ModelForcing, SimpleForcing, ParameterizedForcing,

    # Coriolis forces
    FPlane, BetaPlane, NonTraditionalFPlane,

    # Buoyancy and equations of state
    BuoyancyTracer, SeawaterBuoyancy,
    LinearEquationOfState, RoquetIdealizedNonlinearEquationOfState, TEOS10,

    # Surface waves via Craik-Leibovich equations
    SurfaceWaves,

    # Time stepping
    time_step!,
    TimeStepWizard, update_Δt!,

    # Models
    IncompressibleModel, NonDimensionalModel, Clock,

    # Simulations
    Simulation, run!,
    iteration_limit_exceeded, stop_time_exceeded, wall_time_limit_exceeded,

    # Utilities
    prettytime, pretty_filesize,

    # Turbulence closures
    ConstantIsotropicDiffusivity, ConstantAnisotropicDiffusivity,
    AnisotropicBiharmonicDiffusivity,
    ConstantSmagorinsky, AnisotropicMinimumDissipation

# Standard library modules
using Printf
using Logging
using Statistics
using LinearAlgebra

# Third-party modules
using CUDA
using Adapt
using OffsetArrays
using FFTW
using JLD2
using NCDatasets

import GPUifyLoops

using Base: @propagate_inbounds
using Statistics: mean
using GPUifyLoops: @launch, @loop, @unroll

import Base:
    +, -, *, /,
    size, length, eltype,
    iterate, similar, show,
    getindex, lastindex, setindex!,
    push!

#####
##### Abstract types
#####

"""
    AbstractDiagnostic

Abstract supertype for diagnostics that compute information from the current
model state.
"""
abstract type AbstractDiagnostic end

"""
    AbstractOutputWriter

Abstract supertype for output writers that write data to disk.
"""
abstract type AbstractOutputWriter end

#####
##### Place-holder functions
#####

function TimeStepper end
function run_diagnostic end
function write_output end

#####
##### Include all the submodules
#####

include("Architectures.jl")

using Oceananigans.Architectures: @hascuda
@hascuda begin
    println("CUDA-enabled GPU(s) detected:")
    for (gpu, dev) in enumerate(CUDA.devices())
        println(dev)
    end
end

include("Utils/Utils.jl")
include("Logger.jl")
include("Grids/Grids.jl")
include("Operators/Operators.jl")
include("BoundaryConditions/BoundaryConditions.jl")
include("Fields/Fields.jl")
include("Coriolis/Coriolis.jl")
include("Buoyancy/Buoyancy.jl")
include("SurfaceWaves.jl")
include("TurbulenceClosures/TurbulenceClosures.jl")
include("Solvers/Solvers.jl")
include("Forcing/Forcing.jl")
include("Models/Models.jl")
include("TimeSteppers/TimeSteppers.jl")
include("Diagnostics/Diagnostics.jl")
include("OutputWriters/OutputWriters.jl")
include("Simulations/Simulations.jl")
include("AbstractOperations/AbstractOperations.jl")

#####
##### Re-export stuff from submodules
#####

using .Architectures
using .Utils
using .Grids
using .BoundaryConditions
using .Fields
using .Coriolis
using .Buoyancy
using .SurfaceWaves
using .TurbulenceClosures
using .Solvers
using .Forcing
using .Models
using .TimeSteppers
using .Simulations

end # module


text:
"""
    $SIGNATURES

Calculate the ellipse from covariance matrix and return as lambda function.

Notes
- https://cookierobotics.com/007/

Related

plotCovEllipseLayer
"""
function covEllipseParameterized(distr::MvNormal; meanOffset::Bool=true)
  # println(round.(cov(distr), digits=3) )
  a = cov(distr)[1,1]
  b = cov(distr)[1,2] + cov(distr)[2,1]
  b *= 0.5
  c = cov(distr)[2,2]

  λ1 = 0.5*(a+c) + sqrt(0.25*(a-c)^2 + b^2)
  λ2 = 0.5*(a+c) - sqrt(0.25*(a-c)^2 + b^2)
  θ = if b ≈ 0 && a >= c
    0
  elseif b ≈ 0 && a < c
    pi/2
  else
    atan(λ1-a, b)
  end
  sλ1 = sqrt(λ1)
  sλ2 = sqrt(λ2)
  sθ =  sin(θ)
  cθ =  cos(θ)
  ox = meanOffset ? distr.μ[1] : 0
  oy = meanOffset ? distr.μ[2] : 0
  (t) -> [sλ1*cθ*cos(t)-sλ2*sθ*sin(t)+ox; sλ1*sθ*cos(t)+sλ2*cθ*sin(t)+oy]
end

covEllipseParameterized(pts::Array{Float64,2}; meanOffset::Bool=true) = covEllipseParameterized( fit(MvNormal, pts), meanOffset=meanOffset )
covEllipseParameterized(X::Union{<:BallTreeDensity,<:ManifoldKernelDensity}; meanOffset::Bool=true) = covEllipseParameterized( getPoints(X), meanOffset=meanOffset )
covEllipseParameterized(dfg::AbstractDFG, sym::Symbol; meanOffset::Bool=true, solveKey::Symbol=:default) = covEllipseParameterized( getBelief(dfg, sym, solveKey), meanOffset=meanOffset, solveKey=solveKey )

"""
    $SIGNATURES

Plotting tool to draw Gadfly layers of ellipses of 2D covariance fitted to the belief of factor graph variable nonparametric points.

Related

covEllipseParameterized, plotSLAM2DPosesLandms
"""
function plotCovEllipseLayer( dfg::AbstractDFG,
                              vsym::Symbol;
                              solveKey::Symbol=:default,
                              drawPoints::Bool=true,
                              ellipseColor::AbstractString="gray30",
                              pointsColor::AbstractString="gray30",
                              drawEllipse::Bool=true  )
  #
  PL = []

  if !(drawEllipse || drawPoints)
    return PL
  end

  # points to work from
  bel = getBelief(dfg, vsym, solveKey)
  pp__ = getPoints(bel)
  pp_ = if bel.manifold isa SpecialEuclidean
    # assume Pose2
    (x->x.parts[1]).(pp__)
  else
    # assume Point2
    pp__
  end

  @cast pp[i,j] := pp_[j][i]

  if drawEllipse
    # get ellipse function
    eX = covEllipseParameterized(pp[1:2,:], meanOffset=false)
    vEl = eX.(0:0.02:2pi)
    el = [(x->x[1]).(vEl) (x->x[2]).(vEl)]
    # add suggested PPE mean offset
    el[:,1] .+= getVariablePPE(dfg, vsym, solveKey).suggested[1]
    el[:,2] .+= getVariablePPE(dfg, vsym, solveKey).suggested[2]

    # add the ellipse layers
    plelX2 = Gadfly.layer(x=el[:,1], y=el[:,2], Geom.path, Theme(default_color=parse(Colorant, ellipseColor)))
    push!(PL, plelX2[1])
  end

  # add the points layer if needed
  if drawPoints
    plelX1 = Gadfly.layer(x=pp[1,:],
                          y=pp[2,:],
                          Geom.point,
                          Theme(default_color=parse(Colorant, pointsColor),
                          point_size=1pt))
    #
    push!(PL, plelX1[1])
  end

  return PL
end


text:
### A Pluto.jl notebook ###
# v0.17.5

using Markdown
using InteractiveUtils

# ╔═╡ 9188cab0-7304-11ec-24ed-8d1d5505988b
# hideall
let
    docs_dir = dirname(dirname(@__DIR__))
    pkg_dir = dirname(docs_dir)

    using Pkg: Pkg
    Pkg.activate(docs_dir)
    Pkg.develop(; path=pkg_dir)
    Pkg.instantiate()

    # Putting the include here to avoid Pluto getting confused about cell order.
    include(joinpath(docs_dir, "src", "tutorials_utils.jl"))
end;

# ╔═╡ 38fef837-ff32-4f38-83b4-cf1ab07d2b1e
using CSV

# ╔═╡ 79253973-b0cd-4b57-92dc-b708577f3a26
using DataFrames

# ╔═╡ e06ca210-b517-41f2-986c-c20c95ab7fb4
using TuringGLM

# ╔═╡ 635eae65-3d2c-4f4b-8987-6d93534bbcfd
md"""
Currently, TuringGLM only supports hierarchical models with a single random-intercept.
This is done by using the `(1 | group)` inside the `@formula` macro.

For our **Hierarchical Model** example, let's use a famous dataset called `cheese` (Boatwright, McCulloch & Rossi, 1999), which is data from cheese ratings.
A group of 10 rural and 10 urban raters rated 4 types of different cheeses (A, B, C and D) in two samples.
So we have $4 \cdot 20 \cdot 2 = 160$ observations and 4 variables:

* `cheese`: type of cheese from `A` to `D`
* `rater`: id of the rater from `1` to `10`
* `background`: type of rater, either `rural` or `urban`
* `y`: rating of the cheese
"""

# ╔═╡ 658a690a-bbec-4b80-9fb1-accf65e15832
url = "https://github.com/TuringLang/TuringGLM.jl/raw/main/data/cheese.csv";

# ╔═╡ 64459b91-4e2b-4e52-92ac-bda2f2377a0f
cheese = CSV.read(download(url), DataFrame)

# ╔═╡ 0f38a9fb-20e8-4ff0-8617-add48850892e
md"""
Using `y` as dependent variable and `background` is independent variable with a varying-intercept per `cheese` type:
"""

# ╔═╡ e2cc4deb-cb7a-45df-86b5-6bb47a7ed35e
fm = @formula(y ~ (1 | cheese) + background)

# ╔═╡ 934b7ec9-d4b3-435c-876f-49215dca8809
md"""
We instantiate our model with `turing_model` without specifying any model, thus the default model will be used: `Gaussian()`
"""

# ╔═╡ 5fc736a2-068e-4deb-9a77-bd21d93f6f32
model = turing_model(fm, cheese);

# ╔═╡ f34064fe-fe17-45e6-af79-0a727d408ac0
chn = sample(model, NUTS(), 2_000);

# ╔═╡ ae013cfb-1634-4341-88d0-ca68bc29dab7
# hide
plot_chains(chn)

# ╔═╡ 30622375-7018-45f6-ba84-0171c40b9c3e
md"""
## References

Boatwright, P., McCulloch, R., & Rossi, P. (1999).
Account-level modeling for trade promotion: An application of a constrained parameter hierarchical model.
Journal of the American Statistical Association, 94(448), 1063–1073.
"""

# ╔═╡ Cell order:
# ╠═9188cab0-7304-11ec-24ed-8d1d5505988b
# ╠═635eae65-3d2c-4f4b-8987-6d93534bbcfd
# ╠═38fef837-ff32-4f38-83b4-cf1ab07d2b1e
# ╠═79253973-b0cd-4b57-92dc-b708577f3a26
# ╠═658a690a-bbec-4b80-9fb1-accf65e15832
# ╠═64459b91-4e2b-4e52-92ac-bda2f2377a0f
# ╠═e06ca210-b517-41f2-986c-c20c95ab7fb4
# ╠═0f38a9fb-20e8-4ff0-8617-add48850892e
# ╠═e2cc4deb-cb7a-45df-86b5-6bb47a7ed35e
# ╠═934b7ec9-d4b3-435c-876f-49215dca8809
# ╠═5fc736a2-068e-4deb-9a77-bd21d93f6f32
# ╠═f34064fe-fe17-45e6-af79-0a727d408ac0
# ╠═ae013cfb-1634-4341-88d0-ca68bc29dab7
# ╠═30622375-7018-45f6-ba84-0171c40b9c3e


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct VariableReserve{T <: ReserveDirection} <: Reserve{T}
        name::String
        available::Bool
        time_frame::Float64
        requirement::Float64
        ext::Dict{String, Any}
        time_series_container::InfrastructureSystems.TimeSeriesContainer
        internal::InfrastructureSystemsInternal
    end

Data Structure for the procurement products for system simulations.

# Arguments
- `name::String`
- `available::Bool`
- `time_frame::Float64`: the relative saturation time_frame, validation range: `(0, nothing)`, action if invalid: `error`
- `requirement::Float64`: the required quantity of the product should be scaled by a TimeSeriesData
- `ext::Dict{String, Any}`
- `time_series_container::InfrastructureSystems.TimeSeriesContainer`: internal time_series storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct VariableReserve{T <: ReserveDirection} <: Reserve{T}
    name::String
    available::Bool
    "the relative saturation time_frame"
    time_frame::Float64
    "the required quantity of the product should be scaled by a TimeSeriesData"
    requirement::Float64
    ext::Dict{String, Any}
    "internal time_series storage"
    time_series_container::InfrastructureSystems.TimeSeriesContainer
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function VariableReserve{T}(name, available, time_frame, requirement, ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), ) where T <: ReserveDirection
    VariableReserve{T}(name, available, time_frame, requirement, ext, time_series_container, InfrastructureSystemsInternal(), )
end

function VariableReserve{T}(; name, available, time_frame, requirement, ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), internal=InfrastructureSystemsInternal(), ) where T <: ReserveDirection
    VariableReserve{T}(name, available, time_frame, requirement, ext, time_series_container, internal, )
end

# Constructor for demo purposes; non-functional.
function VariableReserve{T}(::Nothing) where T <: ReserveDirection
    VariableReserve{T}(;
        name="init",
        available=false,
        time_frame=0.0,
        requirement=0.0,
        ext=Dict{String, Any}(),
        time_series_container=InfrastructureSystems.TimeSeriesContainer(),
    )
end

"""Get [`VariableReserve`](@ref) `name`."""
get_name(value::VariableReserve) = value.name
"""Get [`VariableReserve`](@ref) `available`."""
get_available(value::VariableReserve) = value.available
"""Get [`VariableReserve`](@ref) `time_frame`."""
get_time_frame(value::VariableReserve) = value.time_frame
"""Get [`VariableReserve`](@ref) `requirement`."""
get_requirement(value::VariableReserve) = value.requirement
"""Get [`VariableReserve`](@ref) `ext`."""
get_ext(value::VariableReserve) = value.ext
"""Get [`VariableReserve`](@ref) `time_series_container`."""
get_time_series_container(value::VariableReserve) = value.time_series_container
"""Get [`VariableReserve`](@ref) `internal`."""
get_internal(value::VariableReserve) = value.internal

"""Set [`VariableReserve`](@ref) `available`."""
set_available!(value::VariableReserve, val) = value.available = val
"""Set [`VariableReserve`](@ref) `time_frame`."""
set_time_frame!(value::VariableReserve, val) = value.time_frame = val
"""Set [`VariableReserve`](@ref) `requirement`."""
set_requirement!(value::VariableReserve, val) = value.requirement = val
"""Set [`VariableReserve`](@ref) `ext`."""
set_ext!(value::VariableReserve, val) = value.ext = val
"""Set [`VariableReserve`](@ref) `time_series_container`."""
set_time_series_container!(value::VariableReserve, val) = value.time_series_container = val


text:
function demo_steady_solution()
  cur_T = 18.067

  cur_reactor = demo_steady_deck(cur_T)

  cur_reactor.n_bar = 0.8746
  cur_reactor.I_P = 21.627

  cur_reactor.R_0 = 7.500
  cur_reactor.B_0 = 5.627

  cur_reactor.P_F = 3255

  cur_reactor.f_BS = 0.611
  cur_reactor.f_CD = 1.00 - cur_reactor.f_BS
  cur_reactor.f_IN = 0.00

  cur_reactor.beta_N = cur_reactor.max_beta_N
  cur_reactor.q_95 = 4.405
  cur_reactor.P_W = 1.911

  cur_reactor.volume = 2217

  cur_reactor
end


text:
module GLMdiagnostics

import
    GLM,
    Gadfly,
    Distributions,
    Random,
    StatsModels,
    StatsBase,
    Statistics

# Separate .jl files for each function (have to manually add file, otherwise error)

# Function #1: Plot quantile-quantile plot (qqplot)
include("plot_qq.jl")
export plot_qq

# Function #2: Plot fitted versus residuals plot
include("plot_fit_res.jl")
export plot_fit_res

# Function #3: Plot scale-location plot
include("plot_scale_loc.jl")
export plot_scale_loc

# Function #4: Plot all diagnostics plots together
include("plot_diagnostics.jl")
export plot_diagnostics

end

max_stars_repo_name: guysutton/GLMdiagnostics

max_issues_repo_name: guysutton/GLMdiagnostics

max_forks_repo_name: guysutton/GLMdiagnostics


text:
"""
Load all summary tables containing average time series into a dictonary.
"""
function get_summary_tables(path, opt)

    resolutions = opt["resolutions"]

    watersheds = opt["stat_sel"]

    dict_wsh = Dict()

    for watershed in watersheds

        dict_res = Dict()

        for resolution in resolutions

            file = joinpath(path, watershed, "res_$(resolution).csv")

            df = CSV.File(file) |> DataFrame

            df[:time] = DateTime.(df[:time], DateFormat("yyyy-mm-ddTHH:MM:SS"))

            dict_res[resolution] = df

        end

        dict_wsh[watershed] = dict_res

    end

    return dict_wsh

end

"""
Extract one variable from the processed vic results.
"""
function get_variable(res_vic::VicRes, var_name; average = true)

    ivar = find(res_vic.var_names .== var_name)

    if average == true
        return res_vic.time[366:end], res_vic.data_mean[366:end,ivar[1]]
    else
        return res_vic.time[366:end], res_vic.data_all[366:end,ivar[1],:]
    end

end


text:
# test/runtests.jl

using MIRTio
using Test: @test, @testset, detect_ambiguities

@testset "pfile/header" begin
    include("pfile/header.jl")
end

@testset "read_rdb_hdr" begin
    include("pfile/read_rdb_hdr.jl")
end

@testset "hdf5utils" begin
    include("fastmri/hdf5utils.jl")
end

@testset "pfile/pfile" begin
    include("pfile/pfile.jl")
end

@test length(detect_ambiguities(MIRTio)) == 0


text:
# This file was generated, do not modify it. # hide
vv = inverse_transform(stand, w)
sum(abs.(vv .- v))


text:
struct TestQs <: QueryStrategy end

qs_score(qs::TestQs, data::Array{T, 2}, pools::Dict{Symbol, Vector{Int}}) where T <: Real = collect(1:size(data, 2))


text:
#=
Functionalities supporting the schema of `X` when `X` is a `Tables.jl`
compatible table.
=#
struct Schema{names, scitypes, types}
    storednames::Union{Nothing, Vector{Symbol}}
    storedscitypes::Union{Nothing, Vector{Type}}
    storedtypes::Union{Nothing, Vector{Type}}
end

"""
    Schema(names, scitypes, types)

Constructor for a `Schema` object.
"""
function Schema(
    names::Tuple{Vararg{Symbol}},
    scitypes::Type{S},
    types::Type{T},
) where {T<:Tuple, S<:Tuple}
    return Schema{names, S, T}()
end

function Schema{names, scitypes, types}() where {names, scitypes, types}
    return Schema{names, scitypes, types}(nothing, nothing, nothing)
end

function Schema(
    names::Tuple{Vararg{Symbol}},
    scitypes::Type{S},
    ::Nothing
) where {S<:Tuple}
    return Schema{names, S, nothing}(nothing, nothing, nothing)
end

# whether names/types are stored or not
function stored(::Schema{names, scitypes, types}) where {names, scitypes, types}
    return (names === nothing && scitypes === nothing && types === nothing)
end

stored(::Nothing) = false
tuple_of_symbols(x::NTuple{N, Symbol}) where {N} = x
tuple_of_symbols(x) = Tuple(map(Symbol, x))

@inline function Schema(names, scitypes, types; stored::Bool=false)
    if stored || length(names) > SCHEMA_SPECIALIZATION_THRESHOLD
        return Schema{nothing, nothing, nothing}(
            [Symbol(x) for x in names],
            Type[T for T in scitypes],
            Type[T for T in types]
        )
    else
        return Schema{tuple_of_symbols(names), Tuple{scitypes...}, Tuple{types...}}()
    end
end

if VERSION < v"1.1"
    fieldtypes(t) = Tuple(fieldtype(t, i) for i = 1:fieldcount(t))
end

# Note that `getproperty(::Schema, :name)`, `getproperty(::Schema, :scitypes)`
# and `getproperty(::Schema, :types)` cannot return `nothing` even though the
# definition below allows for this case. This is because the nature of the `schema`
# function defined below.
function Base.getproperty(sch::Schema{names, scitypes, types},
                          field::Symbol) where {names, scitypes, types}
    if field === :names
        return names === nothing ? getfield(sch, :storednames) : names
    elseif field === :scitypes
        return scitypes === nothing ?
            (S = getfield(sch, :storedscitypes); S !== nothing ? S : nothing) : fieldtypes(scitypes)
    elseif field === :types
        return types === nothing ?
            (T = getfield(sch, :storedtypes); T !== nothing ? T : nothing) : fieldtypes(types)
    else
        throw(ArgumentError("unsupported property for ScientificTypes.Schema"))
    end
end

Base.propertynames(sch::Schema) = (:names, :scitypes, :types)

# make `Schema` a Tables.jl compatible object.
Tables.istable(::Type{<:Schema}) = true
Tables.columnaccess(::Schema) = true
Tables.columns(s::Schema) = s
Tables.columnnames(sch::Schema) = propertynames(sch)
Tables.schema(sch::Schema) = Tables.Schema(Tables.columnnames(sch), Tuple{Symbol, Type, Type})
Tables.getcolumn(sch::Schema, i::Int) = Tables.getcolumn(sch, Tables.columnnames(sch)[i])

function Tables.getcolumn(sch::Schema, nm::Symbol)
    col = getproperty(sch, nm)
    if col === nothing
        N = length(getproperty(sch, :names))
        return fill(nothing, N)
    end
    return col
end

"""
    schema(X)

Inspect the column types and scitypes of a tabular object.
returns `nothing` if the column types and/or scitypes can't be inspected.

## Example

```
X = (ncalls=[1, 2, 4], mean_delay=[2.0, 5.7, 6.0])
schema(X)
```
"""
schema(X) = schema(X, vtrait(X))

# Fallback
function schema(X, ::Val{:other})
    throw(
        ArgumentError("Cannot inspect the internal scitypes of a non-tabular object. ")
    )
end

function schema(X, ::Val{:table})
    if Tables.columnaccess(X)
        return cols_schema(X)
    else
        return rows_schema(X)
    end
end

function cols_schema(X)
    cols = Tables.columns(X)
    sch = Tables.schema(cols)
    sch === nothing && return nothing
    return _cols_schema(cols, sch)
end

function _cols_schema(cols, sch::Tables.Schema{names, types}) where {names, types}
    N = length(names)
    if N <= COLS_SPECIALIZATION_THRESHOLD
        return __cols_schema(cols, sch)
    else
        stypes = if types === nothing
            Type[
                elscitype(Tables.getcolumn(cols, names[i]))
                for i in Base.OneTo(N)
            ]
        else
            Type[
                elscitype(Tables.getcolumn(cols, fieldtype(types, i), i, names[i]))
                for i in Base.OneTo(N)
            ]
        end
        return Schema(names, Tuple{stypes...}, types)
    end
end

@inline function __cols_schema(cols, sch::Tables.Schema{names, types}) where {names, types}
    N = length(names)
    if @generated
        stypes = if types === nothing
            (
                :(elscitype(Tables.getcolumn(cols, $(Meta.QuoteNode(names[i])))))
                for i in Base.OneTo(N)
            )
        else
            (
                quote
                    elscitype(
                        Tables.getcolumn(
                            cols, $(fieldtype(types, i)), $i,  $(Meta.QuoteNode(names[i]))
                        )
                    )
                end
                for i in Base.OneTo(N)
            )
        end

        return :(Schema(names, Tuple{$(stypes...)}, types))

    else

        stypes = if types === nothing
            (
                elscitype(Tables.getcolumn(cols, names[i]))
                for i in Base.OneTo(N)
            )
        else
            (
                elscitype(Tables.getcolumn(cols, fieldtype(types, i), i, names[i]))
                for i in Base.OneTo(N)
            )
        end

        return Schema(names, Tuple{stypes...}, types)

    end

end

# For extremely wide column oriented tables.
@inline function _cols_schema(cols, sch::Tables.Schema{nothing, nothing})
    names = sch.names
    types = sch.types
    len = length(names)
    stypes = Vector{Type}(undef, len)
    @simd for ind in eachindex(names)
        name = @inbounds(names[ind])
        @inbounds(
            stypes[ind] = elscitype(Tables.getcolumn(cols, name))
        )
    end
    return Schema(names, stypes, types)
end

struct RowsWrapper{T}
    rows::T
end

Tables.istable(::Type{<:RowsWrapper}) = true
Tables.rowaccess(::RowsWrapper) = true
Tables.rows(x::RowsWrapper) = x.rows
Tables.columns(x::RowsWrapper) = Tables.columns(x.rows)

function rows_schema(X)
    rows = Tables.rows(X)
    sch = Tables.schema(rows)
    sch === nothing && return nothing
    return _rows_schema(rows, sch)
end

function _rows_schema(rows, sch::Tables.Schema{names, types}) where {names, types}
    rows_ = RowsWrapper(rows)
    cols = if length(names) <= ROWS_SPECIALIZATION_THRESHOLD
            Tables.columns(rows_)
    else
        Tables.dictcolumntable(rows_)
    end
    return _cols_schema(cols, sch)
end

# for extremely wide row oriented tables
function _rows_schema(rows, sch::Tables.Schema{nothing, nothing})
    cols = Tables.dictcolumntable(RowsWrapper(rows))
    return _cols_schema(cols, sch)
end

function Base.show(io::IO, ::MIME"text/plain", s::Schema)
    data = Tables.matrix(s)
    header = (["names", "scitypes", "types"],)
    pretty_table(io, data, header=header;
                 header_crayon=Crayon(bold=false),
                 alignment=:l)
end


text:
# module MySQL_Octo.SQL

function execute(adapter::QueryBuilders.AbstractAdapter, commands::String)
    print_with_color(:bold, commands, "\n")
    res = mysql_execute(adapter.handle, commands; opformat=MYSQL_TUPLES)
    print_with_color(:normal, join(res, ' '), "\n")
end


text:
function FieldEcksteinHechlerPropagator(arg0::FieldOrbit, arg1::AttitudeProvider, arg10::PropagationType, arg2::RealFieldElement, arg3::jdouble, arg4::RealFieldElement, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return FieldEcksteinHechlerPropagator((FieldOrbit, AttitudeProvider, PropagationType, RealFieldElement, jdouble, RealFieldElement, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function FieldEcksteinHechlerPropagator(arg0::FieldOrbit, arg1::AttitudeProvider, arg2::RealFieldElement, arg3::UnnormalizedSphericalHarmonicsProvider)
    return FieldEcksteinHechlerPropagator((FieldOrbit, AttitudeProvider, RealFieldElement, UnnormalizedSphericalHarmonicsProvider), arg0, arg1, arg2, arg3)
end

function FieldEcksteinHechlerPropagator(arg0::FieldOrbit, arg1::AttitudeProvider, arg2::RealFieldElement, arg3::UnnormalizedSphericalHarmonicsProvider, arg4::PropagationType)
    return FieldEcksteinHechlerPropagator((FieldOrbit, AttitudeProvider, RealFieldElement, UnnormalizedSphericalHarmonicsProvider, PropagationType), arg0, arg1, arg2, arg3, arg4)
end

function FieldEcksteinHechlerPropagator(arg0::FieldOrbit, arg1::AttitudeProvider, arg2::RealFieldElement, arg3::UnnormalizedSphericalHarmonicsProvider, arg4::UnnormalizedSphericalHarmonicsProvider_UnnormalizedSphericalHarmonics)
    return FieldEcksteinHechlerPropagator((FieldOrbit, AttitudeProvider, RealFieldElement, UnnormalizedSphericalHarmonicsProvider, UnnormalizedSphericalHarmonicsProvider_UnnormalizedSphericalHarmonics), arg0, arg1, arg2, arg3, arg4)
end

function FieldEcksteinHechlerPropagator(arg0::FieldOrbit, arg1::AttitudeProvider, arg2::RealFieldElement, arg3::UnnormalizedSphericalHarmonicsProvider, arg4::UnnormalizedSphericalHarmonicsProvider_UnnormalizedSphericalHarmonics, arg5::PropagationType)
    return FieldEcksteinHechlerPropagator((FieldOrbit, AttitudeProvider, RealFieldElement, UnnormalizedSphericalHarmonicsProvider, UnnormalizedSphericalHarmonicsProvider_UnnormalizedSphericalHarmonics, PropagationType), arg0, arg1, arg2, arg3, arg4, arg5)
end

function FieldEcksteinHechlerPropagator(arg0::FieldOrbit, arg1::AttitudeProvider, arg2::RealFieldElement, arg3::jdouble, arg4::RealFieldElement, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return FieldEcksteinHechlerPropagator((FieldOrbit, AttitudeProvider, RealFieldElement, jdouble, RealFieldElement, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function FieldEcksteinHechlerPropagator(arg0::FieldOrbit, arg1::AttitudeProvider, arg2::UnnormalizedSphericalHarmonicsProvider)
    return FieldEcksteinHechlerPropagator((FieldOrbit, AttitudeProvider, UnnormalizedSphericalHarmonicsProvider), arg0, arg1, arg2)
end

function FieldEcksteinHechlerPropagator(arg0::FieldOrbit, arg1::AttitudeProvider, arg2::jdouble, arg3::RealFieldElement, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble)
    return FieldEcksteinHechlerPropagator((FieldOrbit, AttitudeProvider, jdouble, RealFieldElement, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8)
end

function FieldEcksteinHechlerPropagator(arg0::FieldOrbit, arg1::RealFieldElement, arg2::UnnormalizedSphericalHarmonicsProvider)
    return FieldEcksteinHechlerPropagator((FieldOrbit, RealFieldElement, UnnormalizedSphericalHarmonicsProvider), arg0, arg1, arg2)
end

function FieldEcksteinHechlerPropagator(arg0::FieldOrbit, arg1::RealFieldElement, arg2::jdouble, arg3::RealFieldElement, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble)
    return FieldEcksteinHechlerPropagator((FieldOrbit, RealFieldElement, jdouble, RealFieldElement, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8)
end

function FieldEcksteinHechlerPropagator(arg0::FieldOrbit, arg1::UnnormalizedSphericalHarmonicsProvider)
    return FieldEcksteinHechlerPropagator((FieldOrbit, UnnormalizedSphericalHarmonicsProvider), arg0, arg1)
end

function FieldEcksteinHechlerPropagator(arg0::FieldOrbit, arg1::UnnormalizedSphericalHarmonicsProvider, arg2::PropagationType)
    return FieldEcksteinHechlerPropagator((FieldOrbit, UnnormalizedSphericalHarmonicsProvider, PropagationType), arg0, arg1, arg2)
end

function FieldEcksteinHechlerPropagator(arg0::FieldOrbit, arg1::jdouble, arg2::RealFieldElement, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble)
    return FieldEcksteinHechlerPropagator((FieldOrbit, jdouble, RealFieldElement, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function propagate_orbit(obj::FieldEcksteinHechlerPropagator, arg0::FieldAbsoluteDate)
    return jcall(obj, "propagateOrbit", FieldCartesianOrbit, (FieldAbsoluteDate,), arg0)
end

function reset_initial_state(obj::FieldEcksteinHechlerPropagator, arg0::FieldSpacecraftState)
    return jcall(obj, "resetInitialState", void, (FieldSpacecraftState,), arg0)
end

function reset_initial_state(obj::FieldEcksteinHechlerPropagator, arg0::FieldSpacecraftState, arg1::PropagationType)
    return jcall(obj, "resetInitialState", void, (FieldSpacecraftState, PropagationType), arg0, arg1)
end


text:
@testset "ExponentialSmoothing" begin
    @test has_fit_methods(ExponentialSmoothing)

    nile = CSV.File(StateSpaceModels.NILE) |> DataFrame
    model = ExponentialSmoothing(nile.flow)
    fit!(model)
    # statsmodels gives -638.031
    @test loglike(model) ≈ -638.025 atol = 1e-5 rtol = 1e-5
    filt = kalman_filter(model)
    smoother = kalman_smoother(model)
    forec = forecast(model, 10)
    @test monotone_forecast_variance(forec)

    finland_fatalities = CSV.File(StateSpaceModels.VEHICLE_FATALITIES) |> DataFrame
    log_finland_fatalities = log.(finland_fatalities.ff)
    model = ExponentialSmoothing(log_finland_fatalities; trend = true)
    fit!(model)
    @test loglike(model) ≈ 32.258 atol = 1e-5 rtol = 1e-5
    forec = forecast(model, 10)
    @test monotone_forecast_variance(forec)

    finland_fatalities = CSV.File(StateSpaceModels.VEHICLE_FATALITIES) |> DataFrame
    log_finland_fatalities = log.(finland_fatalities.ff)
    model = ExponentialSmoothing(log_finland_fatalities; trend = true, damped_trend = true)
    fit!(model)
    @test loglike(model) ≈ 21.000 atol = 1e-3 rtol = 1e-3
    forec = forecast(model, 10)
    @test monotone_forecast_variance(forec)

    air_passengers = CSV.File(StateSpaceModels.AIR_PASSENGERS) |> DataFrame
    log_air_passengers = log.(air_passengers.passengers)
    model = ExponentialSmoothing(log_air_passengers; trend = true, seasonal = 12)
    # statsmodels gives -273.215
    fit!(model)
    @test loglike(model) ≈ 275.313 atol = 1e-3 rtol = 1e-3
    forec = forecast(model, 10)
    @test monotone_forecast_variance(forec)

    # Automatic ETS
    model = auto_ets(nile.flow)
    @test model.trend == false

    model = auto_ets(log_finland_fatalities)
    @test model.trend == true
    @test model.damped_trend == false

    model = auto_ets(log_air_passengers; seasonal = 12)
    @test model.trend == true
    @test model.damped_trend == false
end


text:
# MIT License
# Copyright (c) 2017: Xavier Gandibleux, Anthony Przybylski, Gauthier Soleilhac, and contributors.
module vOptSpecific

const libLAPpath = joinpath(@__DIR__, "..", "deps", "libLAP.so")
const libUMFLPpath = joinpath(@__DIR__, "..", "deps", "libUMFLP.so")
const libcomboPath = joinpath(@__DIR__,"..","deps","libcombo.so")

include("LAP.jl")
include("UMFLP.jl")
include("jMOScheduling.jl")
include("Bi01KP/Bi01KP.jl")
include("KP.jl")

@static if VERSION > v"0.7-"
	using DelimitedFiles
	using SparseArrays
	_unsafe_wrap(Atype, p, dims ; own=false) = unsafe_wrap(Atype, p, dims, own=own)
end

@static if VERSION < v"0.7-"
	_unsafe_wrap(Atype, p, dims ; own=false) = unsafe_wrap(Atype, p, dims, own)
end

export vSolve,
    load2LAP, set2LAP, LAP_Przybylski2008,
    load2OSP, set2OSP, OSP_VanWassenhove1980, generateHardInstance,
    load2UMFLP, set2UMFLP, UMFLP_Delmee2017,
    load2UKP, set2UKP, UKP_Jorge2010

__init__() = (!isfile(libLAPpath) || !isfile(libUMFLPpath)) || !isfile(libcomboPath) && warn("""Some libraries are missing, run Pkg.build("vOptSpecific")""")

end # module

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific


text:
function init_manager(manager::ExperimentManager)
    Mongoc.delete_many(manager.instances, Mongoc.BSON())
    Mongoc.delete_many(manager.decompositions, Mongoc.BSON())
    load_instance_by_paths(manager, "./data/case9.m", "./data/case9_sdp_ctr.txt", "./data/case9_sdp_mat.txt", "./data/case9.csv")
    load_instance_by_paths(manager, "./data/case9mod.m", "./data/case9mod_sdp_ctr.txt", "./data/case9mod_sdp_mat.txt", "./data/case9mod.csv")
    #load_instance_by_paths(manager, "./data/case1354pegase.m", "./data/case1354pegase_sdp_ctr.txt", "./data/case1354pegase_sdp_mat.txt", "./data/case1354pegase.csv")
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
frule(Δargs, ::Type{MersenneTwister}, args...) = MersenneTwister(args...), Zero()

function rrule(::Type{MersenneTwister}, args...)
    function MersenneTwister_pullback(ΔΩ)
        return (NO_FIELDS, map(_ -> Zero(), args)...)
    end
    return MersenneTwister(args...), MersenneTwister_pullback
end


text:
using EngEconomics
using Test
using Dates

@testset "Simple Interest" begin
	@test simpleInterest(4500, 0.08, 1) == 4860
	@test simpleInterest(4500, 0.08, 4) == 5940
	@test isapprox(simpleInterest(1000, 0.08, 80/360), (1000 + 17.78); atol=0.01)
end

@testset "Fractional Payment Period" begin
	apr1 = Date(Dates.Month(4), Dates.Day(1))
	jun20 = Date(Dates.Month(6), Dates.Day(20))
	@test findFractionalPeriod2(apr1, jun20) == 80/360
end


text:
@testset "c\"extern ...\"" begin
	@eval c``

	asprintf = @eval c"extern int asprintf(char const **str, char *fmt, ...);"
	f(x) = asprintf(x, "%s i%c %d great test of CBinding.jl v%3.1f%c\n", "This", 's', 0x01, 0.1, '!')
	str = Ref(c"char const *"(C_NULL))
	expect = "This is 1 great test of CBinding.jl v0.1!\n"
	bytes = 0
	for i in 1:100
		b = @allocated len = f(str)
		bytes += i > 1 ? b : 0
		@test len == length(expect)
		@test unsafe_string(str[]) == expect
		Libc.free(str[])
		str[] = C_NULL
	end
	@test bytes == 0
	@test_throws MethodError asprintf(1234)
	@test_throws MethodError asprintf(1234, "still wrong")

	# TODO: add test for accessing libjulia global functions and variables
end


text:
# visualization related functions

# can be generalized with a few constants
function _saveSolutionToMesh(mesh::PumiMesh, u::AbstractVector,
                             mnew_ptr=mesh.mnew_ptr, fnew_ptr=mesh.fnew_ptr)
# saves the solution in the vector u to the mesh.f_ptr, and then
# callcs transferFieldToSubmesh to transfer the field fo
# fnew_ptr on mesh mnew_ptr, if needed
# it uses mesh.elementNodeOffsets to access the pumi field values in the
# right order of the given element
# note that this performs a loop over elements, so some values get
# written several times.  This is why is is necessary to write to the
# field in the right order
# because this function used apf.getNumberJ to get dof numbers, and set
# values, it doesn't really need to use the nodemap because values
# are set/get *consistently*, even if not in the same order for every
# element depending on the orientation
 # dofnums = zeros(Int, mesh.numDofPerNode)
#  u_vals = zeros(mesh.numDofPerNode)

  num_entities = [3, 3, 1] # number of vertices, edges, faces
  q_vals = zeros(mesh.numDofPerNode)
  node_entities = apf.ElementNodeEntities(mesh.m_ptr, mesh.mshape_ptr, mesh.dim)

  for el=1:mesh.numEl
    el_i = mesh.elements[el]
    apf.getNodeEntities(node_entities, el_i)
    #node_entities = apf.getNodeEntities(mesh.m_ptr, mesh.mshape_ptr, el_i)
    col = 1 # current node of the element
    for i=1:node_entities.nentities
      entity = node_entities.entities[i]
      for k=1:node_entities.nodecounts[i]
#    for i=1:3  # loop over verts, edges, faces
#      for j=1:num_entities[i]  # loop over all entities of this type
#	for k=1:mesh.numNodesPerType[i]  # loop over nodes on this entity
#          entity = node_entities[col]  # current entity
	  offset_k = mesh.elementNodeOffsets[col, el] # offset for current node
	  new_node = abs(offset_k - k) - 1

          # get solution values
	  for p=1:mesh.numDofPerNode  # loop over all dofs
	    dofnum_p = apf.getNumberJ(mesh.dofnums_Nptr, entity, new_node, p-1)
	    q_vals[p] = u[dofnum_p]
	  end
          # save to mesh
          apf.setComponents(mesh.f_ptr, entity, abs(offset_k - k) - 1, q_vals)

	  col += 1
	end  # end loop over nodes on current entity
      end  # end loop over entities
#    end
  end  # end loop over elements

  transferFieldToSubmesh(mesh, u, mnew_ptr, fnew_ptr)

  return nothing
end  # end function saveSolutionToMesh

function saveNodalSolution(mesh::PumiMesh, u::AbstractArray{Float64, 3})
# u must be a nodal solution

  down_verts = Array{Ptr{Void}}(12)
  for el = 1:mesh.numEl
    el_i = mesh.elements[el]
    nverts = apf.getDownward(mesh.m_ptr, el_i, 0, down_verts)

    for j=1:nverts
      vals_j = sview(u, :, j, el)
      vert_j = down_verts[j]
      apf.setComponents(mesh.fnew_ptr, vert_j, 0, vals_j)
    end
  end

  return nothing
end

function saveSolutionToMesh(mesh::PumiMesh, u::AbstractVector)

  checkMeshPtr(mesh)

  if mesh.isDG
    # all DG meshes interpolate directly
    interpolateToMesh(mesh, u)
  else
    _saveSolutionToMesh(mesh, u)
  end

  if mesh.isDG && mesh.mexact_ptr != C_NULL
    _saveSolutionToMesh(mesh, u, mesh.mexact_ptr, mesh.fexact_ptr)
  end
end

"""
  Function for doing average reduction with [`interpolateToMesh`](@ref).
"""
function avgReduction(u_new::Number, jac_new::Number, u_old::Number)

  return u_new*jac_new + u_old
end

function maxReduction(u_new::Number, jac_new::Number, u_old::Number)

  return max(u_new, u_old)
end

function minReduction(u_new::Number, jac_new::Number, u_old::Number)

  return min(u_new, u_old)
end

"""
  This function interpolates a solution vector onto the coordinate field of
  the mesh.  This field will be printed to Paraview files if they are written.

  By default, this function does a mapping jacobian determinant weighted
  average of all contributions to the node.

  Currently values are averaged only at the vertices, and not mid edge nodes
  for quadratic meshes.

  **Inputs**

   * mesh: mesh object
   * u: solution vector, length mesh.numDof.  Can have any element type that
        supports the real() function.  Only the part returned by real() will be
        written to the field (must be convertable to a C double)
   * reduce_op: this function is applied to combine values at each node.
                Because the field is C0 continuous but the solution stroed
                in `u` may not be, it is necessary to combine the potentially
                different values for nodes shared by more than one element.
                `reduce_op` should have the signature

                `function reduce_op_name(u_new::Number, jac_new::Number, u_old::Number)`

                where `u_old` is the existing value at the node, `u_new` is the
                value to be combined with `u_old`, and `jac_new` is the determinant
                of the mapping jacobian at the node.  The function should compute
                the new value at the node and return it.  `u_old` will be zero
                the first time the function is called at a given node.
                No guarantees are made about the order in which values are
                incorporated to each node, so `reduce_op` should be associated.
                In the special case when `reduce_op` is [`avgReduction`](@ref),
                a second loop will be performed to divide the value at each node
                by the sum of the mapping jacobian determinants of all values
                that contribute to it (providing an average weighted by the
                determinant).

                Default value: `avgReduction`
"""
function interpolateToMesh(mesh::PumiMesh{T}, u::AbstractVector, reduce_op::Function=avgReduction) where T
# interpolate solution stored in u to the field mesh.fnew_ptr which resides on
# mesh.mnew_ptr
# u is the vector containing the solution, even though for DG the vector and
# the 3D array forms of the solution have identical memory layout
# change this in the future with ReshapedArrays?

# the interpolation operator *must* interpolate the solution in a particular
# order.  The order must be the order of the nodes returned by
# apf.getNodeEntities(mnew_ptr)

  @assert mesh.m_ptr == mesh.mnew_ptr  # needed for averaging step
  @assert size(mesh.interp_op, 1) == mesh.coord_numNodesPerElement

 # println(mesh.f, "entered interpolateToMesh")
 # flush(mesh.f)

  myrank = mesh.myrank
  interp = mesh.interp_op
  u_el = zeros(Float64, mesh.numNodesPerElement, mesh.numDofPerNode)
  u_verts = zeros(Float64, size(interp, 1), mesh.numDofPerNode)
  jac_verts = zeros(T, size(interp, 1))
  u_node = zeros(Float64, mesh.numDofPerNode)  # hold new node values
  u_node2 = zeros(u_node)  # hold existing node values
  node_entities = Array{Ptr{Void}}(mesh.coord_numNodesPerElement)
  dofs = mesh.dofs
  numTypePerElement = mesh.numTypePerElement
  numEntitiesPerType = mesh.numEntitiesPerType
  fshape_ptr = mesh.fnewshape_ptr
  vshare = mesh.vert_sharing
  # this is ugly and inefficient
  match_data = zeros(mesh.numDofPerNode + 1, mesh.numVert)

  matches_partnums = zeros(Cint, 400 + 8)  # upper bound
  matches_entities = Array{Ptr{Void}}(400 + 8)

  # count nodes on solution field
  numNodesPerType = Array{Int}(mesh.dim + 1)
  numNodesPerType[1] = apf.countNodesOn(fshape_ptr, 0)
  numNodesPerType[2] = apf.countNodesOn(fshape_ptr, 1)
  numNodesPerType[3] = apf.countNodesOn(fshape_ptr, 2)
  if mesh.dim == 3
    numNodesPerType[4] = apf.countNodesOn(fshape_ptr, 4) # tetrahedron
  end
  node_entities_mnew = apf.ElementNodeEntities(mesh.mnew_ptr, fshape_ptr, mesh.dim)

  node_entities_m = apf.ElementNodeEntities(mesh.m_ptr, fshape_ptr, mesh.dim)

  apf.zeroField(mesh.fnew_ptr)

  # data for MPI
  # values solution values in first n indices of inner array, followed
  # by the weighting factor
  peer_vals_send = Array{Array{Float64, 2}}(vshare.npeers)
  peer_vals_recv = Array{Array{Float64, 2}}(vshare.npeers)
  for i=1:vshare.npeers
    peer_vals_send[i] = zeros(Float64, mesh.numDofPerNode + 1, vshare.counts[i])
    peer_vals_recv[i] = zeros(peer_vals_send[i])

  end

  coords = zeros(3)  # debugging
  for el=1:mesh.numEl
#    println(mesh.f, "element ", el)
    el_i = mesh.elements[el]

    # get the solution values out of u
    for j=1:mesh.numNodesPerElement
      for k=1:mesh.numDofPerNode
        dof_k = dofs[k, j, el]
        u_el[j, k] = real(u[dof_k])
      end
    end

    # interpolate solution
    smallmatmat!(interp, u_el, u_verts)
#    println(mesh.f, "u_el = \n", u_el)
#    println(mesh.f, "u_verts = \n", u_verts)

    # interpolate jacobian
    jac_el = sview(mesh.jac, :, el)
    smallmatvec!(interp, jac_el, jac_verts)

    # save values to mesh
    # assumes the solution fieldshape is the same as the coordinate fieldshape
    apf.getNodeEntities(node_entities_mnew, el_i)
#    apf.getNodeEntities(mesh.mnew_ptr, fshape_ptr, el_i, node_entities)
    col = 1  # current node of the element

    for i=1:node_entities_mnew.nentities
      entity = node_entities_mnew.entities[i]
      edim = node_entities_mnew.dimensions[i]
      for k=1:node_entities_mnew.nodecounts[i]
#    for i=1:(mesh.dim+1)
#      for j=1:numTypePerElement[i]
#        for k=1:numNodesPerType[i]
#          entity = node_entities[col]
#          @assert i == 1
          vertnum = apf.getNumberJ(mesh.entity_Nptrs[edim+1], entity, 0, 0) + 1
          getCoords(mesh, entity, 0, coords)

          # pack array to send to other processes
          if edim == 0 && haskey(vshare.rev_mapping, vertnum)
#            println(mesh.f, "\nvert number ", vertnum)
            pair = vshare.rev_mapping[vertnum]
#            println(mesh.f, "peers = \n", pair.first)
#            println(mesh.f, "shared vert idxs = \n", pair.second)
            for first_idx = 1:length(pair.first)
              peer_idx = findfirst( vshare.peer_nums, pair.first[first_idx])
#              println(mesh.f, "peer_idx = ", peer_idx)
              shared_vert_idx = pair.second[first_idx]
              arr_peer = peer_vals_send[peer_idx]
              # accumulate weighted solution values, weighting factor
#              println(mesh.f, "u_verts = \n", u_verts)

              for p=1:mesh.numDofPerNode
                arr_peer[p, shared_vert_idx] = reduce_op(u_verts[col, p], real(jac_verts[col]), arr_peer[p, shared_vert_idx])
              end
              arr_peer[mesh.numDofPerNode + 1, shared_vert_idx] += real(jac_verts[col])
            end  # end loop over peers
          end  # end if haskey

          # skip elementNodeOffsets - maximum of 1 node per entity
          apf.getComponents(mesh.fnew_ptr, entity, 0, u_node2)
          for p=1:mesh.numDofPerNode
            # multiply by element volume as a weighting factor
#            println(mesh.f, "adding contribution ", jac_verts[col]*u_verts[col, p], " to vert at ", coords[1], ", ", coords[2], ", ", coords[3])
#            println(mesh.f, "jac = ", jac_verts[col])
#            println(mesh.f, "u_verts = ", u_verts[col, p])

            u_node[p] = reduce_op(u_verts[col, p], real(jac_verts[col]), u_node2[p])
#            u_node[p] = real(jac_verts[col])*u_verts[col, p] + u_node2[p]
          end

          # check for local matches
          nmatches = apf.countMatches(mesh.m_ptr, entity)
          if edim == 0 && nmatches > 0
            @assert nmatches <= length(matches_partnums)
            apf.getMatches(matches_partnums, matches_entities)
            for match_idx = 1:nmatches
              if matches_partnums[match_idx] == myrank
                # save the data to the other match
                matched_vertnum = apf.getNumberJ(mesh.vert_Nptr, matches_entities[match_idx], 0, 0) + 1
                for p=1:mesh.numDofPerNode
                  match_data[p, matched_vertnum] = reduce_op(u_verts[col, p], real(jac_verts[col]), match_data[p, matched_vertnum])
                end
                match_data[mesh.numDofPerNode + 1, matched_vertnum] += real(jac_verts[col])
              end  # end if match is on this part
            end  # end loop over matches
          end  # end if nmatches > 0

#          println(mesh.f, "u_node is now ", u_node[1])

          apf.setComponents(mesh.fnew_ptr, entity, 0, u_node)
          col += 1
        end
     # end
    end  # end loop over entity dimensions

  end  # end loop over elements

  # send the data
  send_reqs = Array{MPI.Request}(vshare.npeers)
  recv_reqs = Array{MPI.Request}(vshare.npeers)
  for i=1:vshare.npeers
    # use the tag 20 in case any other parts of the program have active
    # communications
    send_reqs[i] = MPI.Isend(peer_vals_send[i], vshare.peer_nums[i], 20, mesh.comm)
    recv_reqs[i] = MPI.Irecv!(peer_vals_recv[i], vshare.peer_nums[i], 20, mesh.comm)
  end

  # wait for all recieves to finish
  for i=1:vshare.npeers
    MPI.Wait!(recv_reqs[i])
  end

  # add in the contributions from each peer to the weighted sum
  #TODO: use MPI.WaitAny! to combine this with the above loop
  for peer=1:vshare.npeers
    peer_vals_p = peer_vals_recv[peer]
    vertnums_p = vshare.vert_nums[peer]

    for i=1:length(vertnums_p)
      vert_i = mesh.verts[vertnums_p[i]]
      getCoords(mesh, vert_i, 0, coords)
      apf.getComponents(mesh.fnew_ptr, vert_i, 0, u_node)
      weight_i = peer_vals_p[mesh.numDofPerNode + 1, i]
      for p=1:mesh.numDofPerNode
#      println(mesh.f, "adding parallel contribution ", peer_vals_p[p, i], " to vert at ", coords[1], ", ", coords[2], ", ", coords[3])
        # set jac_det to 1 here because both the old and new values have
        # already been scaled
        u_node[p] = reduce_op(peer_vals_p[p, i], 1.0, u_node[p])
      end

      # check for local matches
      apf.setComponents(mesh.fnew_ptr, vert_i, 0, u_node)
    end
  end

  # add in the local match data
  for i=1:mesh.numVert
    vert_i = mesh.verts[i]
    getCoords(mesh, vert_i, 0, coords)
    apf.getComponents(mesh.fnew_ptr, vert_i, 0, u_node)
    for p=1:mesh.numDofPerNode
#      println(mesh.f, "adding local match contribution ", match_data[p, i], " to vert at ", coords[1], ", ", coords[2], ", ", coords[3])
      u_node[p] = reduce_op(match_data[p, i], 1.0, u_node[p])
    end
    apf.setComponents(mesh.fnew_ptr, vert_i, 0, u_node)
  end

  if reduce_op == avgReduction
    # divide by the total volume of elements that contributed to each node
    # so the result is the average value
    up_els = Array{Ptr{Void}}(400)  # equivalent of apf::Up
        for (entity_i, dim) in apf.FieldEntityIt(mesh.m_ptr, fshape_ptr)
          entity_num = apf.getNumberJ(mesh.entity_Nptrs[dim+1], entity_i, 0, 0) + 1
          nel = apf.countAdjacent(mesh.mnew_ptr, entity_i, mesh.dim)
          apf.getAdjacent(up_els)
          # compute the sum of the volumes of the elements
          jac_sum = 0.0
          for j=1:nel
            el_j = up_els[j]
            elnum_j = apf.getNumberJ(mesh.el_Nptr, el_j, 0, 0) + 1

            jac_el = sview(mesh.jac, :, elnum_j)
            # searching apf.getNodeEntities for the index guarantees we can
            # identify the right interpolated point for all elements,
            # independent of topology
            apf.getNodeEntities(node_entities_m, el_j)
#            apf.getNodeEntities(mesh.m_ptr, fshape_ptr, el_j, node_entities)
            entity_idx = findfirst(node_entities_m.entities, entity_i)
            node_idx = 1
            for p=1:(entity_idx-1)
              node_idx += node_entities_m.nodecounts[p]
            end
#            node_idx = findfirst(node_entities, entity_i)

            # interpolate the jacobian to this point
            jac_entity = 0.0
            for p=1:size(interp, 2)
              jac_entity += interp[node_idx, p]*real(jac_el[p])
            end
            jac_sum += jac_entity
          end  # end loop j

          # add in the parallel contributions
          if dim == 0 && haskey(vshare.rev_mapping, entity_num)
            pair = vshare.rev_mapping[entity_num]
            for i=1:length(pair.first)  # loop over peers that share this vert
              peer_idx = findfirst(vshare.peer_nums, pair.first[i])
              vert_idx = pair.second[i]
              jac_entity = peer_vals_recv[peer_idx][mesh.numDofPerNode + 1, vert_idx]
              jac_sum += jac_entity
            end
          end  # end if haskey

          if dim == 0
            jac_entity = match_data[mesh.numDofPerNode + 1, entity_num]
            jac_sum += jac_entity
          end

          apf.getComponents(mesh.fnew_ptr, entity_i, 0, u_node)
          fac = 1/jac_sum
          for p=1:mesh.numDofPerNode
            u_node[p] = fac*u_node[p]
          end
          apf.setComponents(mesh.fnew_ptr, entity_i, 0, u_node)
        end  # end iterator
  end  # end if reduce_op == avgReduce

  # wait for sends to finish before exiting
  for i=1:vshare.npeers
    MPI.Wait!(send_reqs[i])
  end

#  println(mesh.f, "finished interpolateToMesh")
#  flush(mesh.f)

  return nothing
end  # end function

"""
  This function is the counterpart to [`interpolateToMesh`](@ref).  It retrieves
  a solution saved to the mesh and puts in into a vector

  **Inputs**

   * mesh
   * fnew_ptr: the apf::Field* to retrieve from. Must have fieldshape
               mesh.fnewshape_ptr (typically the coordinate fieldshape)

  **Inputs/Outputs**

   * u_vec: to be overwritten
"""
function retrieveSolutionFromMesh_interp(mesh::PumiMeshDG, u_vec::AbstractVector, fnew_ptr=mesh.fnew_ptr)

  @assert mesh.m_ptr == mesh.mnew_ptr  # needed for averaging step
  @assert size(mesh.interp_op, 1) == mesh.coord_numNodesPerElement

  myrank = mesh.myrank
  interp = mesh.interp_op2
  u_el = zeros(Float64, mesh.numNodesPerElement, mesh.numDofPerNode)
  u_verts = zeros(Float64, size(interp, 2), mesh.numDofPerNode)
  u_node = zeros(Float64, mesh.numDofPerNode)  # hold new node values
  node_entities = Array{Ptr{Void}}(mesh.coord_numNodesPerElement)
  dofs = mesh.dofs
  numTypePerElement = mesh.numTypePerElement
#  numEntitiesPerType = mesh.numEntitiesPerType
  fshape_ptr = mesh.fnewshape_ptr

  # count nodes on solution field
  numNodesPerType = Array{Int}(mesh.dim + 1)
  numNodesPerType[1] = apf.countNodesOn(fshape_ptr, 0)
  numNodesPerType[2] = apf.countNodesOn(fshape_ptr, 1)
  numNodesPerType[3] = apf.countNodesOn(fshape_ptr, 2)
  if mesh.dim == 3
    numNodesPerType[4] = apf.countNodesOn(fshape_ptr, 4) # tetrahedron
  end

  for i=1:mesh.numEl
    el_i = mesh.elements[i]

    # save values to mesh
    # assumes the solution fieldshape is the same as the coordinate fieldshape
    apf.getNodeEntities(mesh.mnew_ptr, fshape_ptr, el_i, node_entities)
    col = 1  # current node of the element
    for d=1:(mesh.dim+1)
      for j=1:numTypePerElement[d]
        for k=1:numNodesPerType[d]
          entity = node_entities[col]

          # skip elementNodeOffsets - maximum of 1 node per entity
          apf.getComponents(fnew_ptr, entity, 0, u_node)
          for p=1:mesh.numDofPerNode
            u_verts[col, p] = u_node[p]
          end

          col += 1

        end  # end k
      end # end j
    end  # end d

    # interpolate solution back to the solution nodes
    smallmatmat!(interp, u_verts, u_el)

    for j=1:mesh.numNodesPerElement
      for k=1:mesh.numDofPerNode
        u_vec[mesh.dofs[k, j, i]] = u_el[j, k]
      end
    end

  end  # end loop i

  return nothing
end

function writeVisFiles(mesh::PumiMeshDG, fname::AbstractString; writeall::Bool=false)
  # writes vtk files

  checkMeshPtr(mesh)

  if mesh.myrank == 0
    println("writing visualization files ", fname)
  end
  f_arr = collect(mesh.fields)
  n_arr = collect(mesh.numberings)
  gn_arr = Vector{Ptr{Void}}(0)
  apf.writeVtkFiles(fname, mesh.mnew_ptr, f_arr, n_arr, gn_arr, writeall=writeall)

  if mesh.mexact_ptr != C_NULL
    fname_exact = fname*"_exact"
    apf.writeVtkFiles(fname_exact, mesh.mexact_ptr, writeall=writeall)
  end

  return nothing
end

function writeVisFiles(mesh::PumiMesh2CG, fname::AbstractString; writeall::Bool=false)
  # writes vtk files

  checkMeshPtr(mesh)

  if mesh.order <= 2
    f_arr = collect(mesh.fields)
    n_arr = collect(mesh.numberings)
    gn_arr = Vector{Ptr{Void}}(0)

    apf.writeVtkFiles(fname, mesh.m_ptr, f_arr, n_arr, gn_arr,
                      writeall=writeall)
  else
    apf.writeVtkFiles(fname, mesh.mnew_ptr, writeall=writeall)
  end

  return nothing
end


text:
import Pkg

Pkg.add("Mustache")
Pkg.add("DataFrames")


text:
# mri-sense.jl test

const DEBUG = false

using ImagePhantoms: ellipse_parameters, SheppLoganBrainWeb, Ellipse
using ImagePhantoms: phantom, spectrum
using ImagePhantoms: mri_smap_fit, mri_spectra
using FFTW: fft, fftshift
using Unitful: mm
using Random: seed!
using ImageGeoms: embed
using LazyGrids: ndgrid
using LinearAlgebra: norm
using Test: @test, @testset, @test_throws, @inferred

if DEBUG
    include("helper.jl")
    using MIRTjim: jim, prompt; jim(:prompt, true)
else
    jim(args...; kwargs...) = nothing
end
jimf = (args...; kwargs...) -> jim(args...; prompt=false, kwargs...)

# image geometry
fovs = (256mm, 250mm)
nx, ny = (64, 50) .* 1
dx, dy = fovs ./ (nx,ny)
x = (-(nx÷2):(nx÷2-1)) * dx
y = (-(ny÷2):(ny÷2-1)) * dy

# define object
oa = ellipse_parameters(SheppLoganBrainWeb() ; disjoint=true, fovs)
seed!(0)
oa[:,end] = [1; rand(ComplexF32,9)] # random phases
oa = Ellipse(oa)
oversample = 3
image0 = phantom(x, y, oa, oversample)
cfun = z -> cat(dims = ndims(z)+1, real(z), imag(z))
jim(x, y, cfun(image0), "image0")

@test spectrum(oa[1]) isa Function

# define mask
mask = trues(nx,ny)
mask[:,1:2] .= false
mask[[1:8;end-8:end],:] .= false
@assert mask .* image0 == image0
jim(x, y, mask, "mask")

# sensitivity maps (smaps) (not normalized)
ncoil = 2
smap = cat(dims=3,
    7 * exp.(im * 2pi * ((0:nx-1) .- 27) / nx * 3) * ones(ny)',
    5 * ones(nx) * exp.(-im * 2pi * ((0:ny-1) .+ 39) / ny * 1)',
)
smap .*= mask
stacker = x -> [(@view x[:,:,i]) for i=1:size(x,3)]
smaps = stacker(smap)
jim(x, y, cfun(smaps), "Sensitivity maps")

# fit each smap
deltas = (dx, dy)
kmax = 7
fit = @NOTinferred mri_smap_fit(smaps, embed; mask, kmax, deltas)
@test fit isa NamedTuple
@test fit.nrmse ≤ 1e-6

jim(
 jim(x, y, cfun(smaps), "orig"; prompt=false),
 jim(x, y, cfun(fit.smaps), "fit"; prompt=false),
 jim(x, y, cfun(1e6*(fit.smaps - smaps)), "err * 1e6"; prompt=false),
)

coefs = map(x -> reshape(x,15,15), fit.coefs)
jim(-kmax:kmax, -kmax:kmax, cfun(coefs), "coefs")

fx = (-(nx÷2):(nx÷2-1)) / (nx*dx) # crucial to match smap_basis internals!
fy = (-(ny÷2):(ny÷2-1)) / (ny*dy)
gx, gy = ndgrid(fx, fy)

# test one object with no smap
ob = oa[3]
coil = 2
image1 = phantom(x, y, [ob], oversample)

#=

fun0 = @inferred spectrum(ob)
kspace0 = fun0.(fx, fy')
kspace1 = spectrum(fx, fy, [ob])
#@test kspace0 ≈ kspace1 # https://github.com/PainterQubits/Unitful.jl/pull/468
@test kspace0 / oneunit(eltype(kspace0)) ≈ kspace1 / oneunit(eltype(kspace1))

kspace1 = myfft(image1) * dx * dy

jim(
 jimf(fx, fy, cfun(kspace0), "analytical"),
 jimf(fx, fy, cfun(kspace1), "fft"),
 jimf(fx, fy, cfun(kspace1 - kspace0), "error"),
)

#@test norm
=#

# test one object with one smap

image2 = image1 .* smaps[coil] # digital
jim(
 jimf(x, y, cfun(image1), "complex image"),
 jimf(x, y, cfun(image2), "with smap"),
 jimf(x, y, cfun(smaps[coil]), "smaps"),
)

kspace2 = myfft(image2) * dx * dy
fun = spectrum(ob, fit, coil)
kspace3 = fun.(fx, fy')
jim(
 jimf(fx, fy, cfun(kspace2), "fft"),
 jimf(fx, fy, cfun(kspace3), "analytical"),
 jimf(fx, fy, cfun(kspace3 - kspace2), "err"),
)

@test norm(kspace3 - kspace2) / norm(kspace2) ≤ 0.09

# test multiple objects with all smaps

#kspace0 = spectrum(oa).(vec(gx), vec(gy))
#kspace0 = reshape(spectrum(oa).(vec(gx), vec(gy)), nx, ny)
#kspace0 = spectrum(oa).(fx, fy') # without smaps
#kspace1 = mri_spectra(fx, fy, oa, fit) # no!
kspace1 = mri_spectra(vec(gx), vec(gy), oa, fit)
kspace1 = [reshape(k, nx, ny) for k in kspace1]
#@test kspace0 / oneunit(eltype(kspace0)) ≈ kspace1 / oneunit(eltype(kspace1[1]))

image2 = [image0 .* s for s in smaps] # digital
kspace2 = myfft.(image2) * (dx * dy)

#p0 = jim(fx, fy, cfun(kspace0), "kspace0: before smap"; prompt=false)
jim(
 jimf(x, y, image0, "image0"),
 jimf(fx, fy, cfun(kspace1), "anal with smaps"),
 jimf(fx, fy, cfun(kspace2), "fft with smaps"),
 jimf(fx, fy, cfun(kspace2 - kspace1), "error"),
)

norma = za -> sqrt(sum(abs2, norm.(za))) # for array of arrays
@test norma(kspace2 - kspace1) / norma(kspace2) ≤ 0.03


text:
"""
Orthogonalize!(V, hv, vv, orth)

Orthogonalize the Krylov vectors using your (my) choice of
methods. Anything other than classical Gram-Schmidt twice (cgs2) is
likely to become an undocumented option. Methods other than cgs2 are
for CI for the linear solver.
"""
function Orthogonalize!(V, hv, vv, orth = "mgs1")
    orthopts = ["mgs1", "mgs2", "cgs1", "cgs2"]
    orth in orthopts || error("Impossible orth spec in Orthogonalize!")
    if orth == "mgs1"
        mgs!(V, hv, vv)
    elseif orth == "mgs2"
        mgs!(V, hv, vv, "twice")
    elseif orth == "cgs1"
        cgs!(V, hv, vv, "once")
    else
        cgs!(V, hv, vv, "twice")
    end
end

"""
cgs!(V, hv, vv, orth="twice")

Classical Gram-Schmidt.
"""
function cgs!(V, hv, vv, orth)
    T = eltype(V)
    k = length(hv)
    @views rk = hv[1:k-1]
    onep = T(1.0)
    zerop = T(0.0)
    pk = zeros(T, k - 1)
    qk = vv
    Qkm = V
    #
    # As the inventors of the BLAS said, "RTFM"!
    #
    # Orthogonalize!
    BLAS.gemv!('T', onep, Qkm, qk, onep, rk)
    BLAS.gemv!('N', -onep, Qkm, rk, onep, qk)
    if orth == "twice"
        #
        # Orthogonalize again!
        BLAS.gemv!('T', onep, Qkm, qk, zerop, pk)
        BLAS.gemv!('N', -onep, Qkm, pk, onep, qk)
        rk .+= pk
    end
    #
    # Keep track of what you did.
    nqk = norm(qk)
    qk ./= nqk
    hv[k] = nqk
end

#
function cgs!(V::SubArray{Float16,2}, hv, vv, orth)
    #
    #   no BLAS
    #
    k = length(hv)
    T = eltype(V)
    onep = T(1.0)
    zerop = T(0.0)
    @views rk = hv[1:k-1]
    pk = zeros(T, size(rk))
    qk = vv
    Qkm = V
    # Orthogonalize
    rk .+= Qkm' * qk
    qk .-= Qkm * rk
    if orth == "twice"
        # Orthogonalize again
        pk .= Qkm' * qk
        qk .-= Qkm * pk
        rk .+= pk
    end
    # Keep track of what you did.
    nqk = norm(qk)
    nqk != 0.0 || println("breakdown")
    nqk == 0.0 || qk ./= nqk
    hv[k] = nqk
end

"""
mgs!(V, hv, vv, orth)
"""
function mgs!(V, hv, vv, orth = "once")
    k = length(hv) - 1
    normin = norm(vv)
    #p=copy(vv)
    @views for j = 1:k
        p = vec(V[:, j])
        hv[j] = p' * vv
        vv .-= hv[j] * p
    end
    hv[k+1] = norm(vv)
    if (normin + 0.001 * hv[k+1] == normin) && (orth == "twice")
        @views for j = 1:k
            p = vec(V[:, j])
            hr = p' * vv
            hv[j] += hr
            vv .-= hr * p
        end
        hv[k+1] = norm(vv)
    end
    nv = hv[k+1]
    #
    # Watch out for happy breakdown
    #
    #if hv[k+1] != 0
    #@views vv .= vv/hv[k+1]
    nv != 0 || println("breakdown")
    if nv != 0
        vv ./= nv
    end
end


text:
# Default implementation forcing specific definitions of interface methods

struct MissingImplementationError <: Exception
    method::Symbol
    whatfor::String
end

Base.showerror(io::IO, e::MissingImplementationError) =
    print(io, "$(e.method) need to be implemented for $(e.whatfor).")

"""
    bisect(::AbstractSearch, elem)

Return two new elements built by bisecting `elem`.

Returned elements must have the same type as the bisected element.

Must be implemented by the user for custom search types, default implementation
throw a `MissingImplementationError`.
"""
bisect(::AbstractSearch, elem) = throw(MissingImplementationError(:bisect, "custom search types"))

"""
    get_leaf_id!(::AbstractSearch, wt::BPTree)

Return the `id` of the next leaf that must be processed and remove this leaf
from the list of working leaves of the working tree `wt`.
"""
get_leaf_id!(::AbstractSearch, wt) = throw(MissingImplementationError(:get_leaf_id, "custom search types"))

"""
    insert_leaf!(::AbstractSearch, wt::BPTree, leaf::BPLeaf)

Insert a leaf in the list of working leaves in the tree `wt`.

Must be implemented by the user for custom search types, default implementation
throw a `MissingImplementationError`.
"""
insert_leaf!(::AbstractSearch, wt) = throw(MissingImplementationError(:insert_leaf!, "custom search types"))

"""
    keyfunc(::AbstractKeySearch, elem)

Return the key associated with element `elem` for a `AbstractKeySearch`.

The `AbstractKeySearch` processes elements with the largest key first.

Must be implemented by the user, default implementation throw a
`MissingImplementationError`.
"""
keyfunc(::AbstractKeySearch, elem) = throw(MissingImplementationError(:keyfunc, "AbstractKeySearch"))

"""
    process(::AbstractSearch, elem)`

Return a symbol representing the action to perform with the element `elem` and
an object (of the same type) representing the state of the element after
processing (may return `elem` unchanged).

# Valid symbols returned by the process function
  - `:store`: the element is considered as final and is stored, it will not be
        further processed
  - `:bisect`: the element is bisected and each of the two resulting part will
        be processed
  - `:discard`: the element is discarded from the tree, allowing to free memory

Must be implemented by the user for custom search types, default implementation
throw a `MissingImplementationError`.
"""
process(::AbstractSearch, wt) = throw(MissingImplementationError(:process, "custom search types"))


text:
# Test `SocketCounterRecord`
@testset "SocketCounterRecord" begin
    a = Counters.SocketCounterRecord{6}()

    # Populate some fields
    a.aggregate_core[:dtlb_store_misses] = [10, 11]
    a.imc_channels[1][:tag_hits] = [20, 21]

    # Create another one and make sure `merge!` works.
    b = Counters.SocketCounterRecord{6}()
    b.aggregate_core[:dtlb_load_misses] = [20]
    b.imc_channels[1][:tag_hits] = [10]
    b.imc_channels[2][:tag_misses] = [100]

    merge!(a, b)
    @test a.aggregate_core[:dtlb_store_misses] == [10, 11]
    @test a.aggregate_core[:dtlb_load_misses] == [20]
    # Should be replaced by `b`
    @test a.imc_channels[1][:tag_hits] == [10]
    @test a.imc_channels[2][:tag_misses] == [100]

    # Create another dummy channel to make sure the non-mutating `merge` works.
    c = Counters.SocketCounterRecord{6}()
    c.aggregate_core[:entry] = [0]
    c.imc_channels[3][:hello] = [-1]

    d = merge(a, c)
    @test d.aggregate_core[:dtlb_store_misses] == [10, 11]
    @test d.aggregate_core[:dtlb_load_misses] == [20]
    @test d.aggregate_core[:entry] == [0]
    @test d.imc_channels[1][:tag_hits] == [10]
    @test d.imc_channels[2][:tag_misses] == [100]
    @test d.imc_channels[3][:hello] == [-1]
end

@testset "Testing PCM Hooks" begin
    measurements = (
        counters = Counters.CoreMonitorWrapper(),
    )

    # Make sure we get the resulting data back in the way we expect.
    trace, processed = SystemSnoop.snoop(measurements) do snooper
        for i in 1:10
            sleep(0.1)
            SystemSnoop.measure!(snooper)
        end
        return SystemSnoop.postprocess(snooper)
    end

    # Get out the socket data from the returned NamedTuple
    data = processed.socket_1
    @test isa(data, Counters.SocketCounterRecord)

    events = Counters.default_events()
    for event in events
        # Make sure the events were recorded and that they are the correct length.
        @test haskey(data.aggregate_core, event.name)
        @test length(data.aggregate_core[event.name]) == 10
    end

    # Now try Uncore Counters
    counters = Counters.DEFAULT_NT
    measurements = (
        counters = Counters.Uncore{2,2,6}(counters),
    )

    # Make sure we get the resulting data back in the way we expect.
    trace, processed = SystemSnoop.snoop(measurements) do snooper
        for i in 1:10
            sleep(0.1)
            SystemSnoop.measure!(snooper)
        end
        return SystemSnoop.postprocess(snooper)
    end
    nsockets = 2

    for i in 1:nsockets
        # Get out the socket data from the returned NamedTuple
        data = getproperty(processed, Symbol("socket_$(i-1)"))
        @test isa(data, Counters.SocketCounterRecord)

        L = length(data.imc_channels)
        for i in 1:L
            for name in keys(counters)
                @test haskey(data.imc_channels[i], name)
                @test length(data.imc_channels[i][name]) == 10
            end
        end
    end
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
@testset "invalid" begin
    c = lxmock(raw"""\begin{:sectiom}""")
    @test_throws ErrorException lx_begin(c,0)
end

@testset "section" begin
    F.LOCAL_VARS["sections"] = Pair{String,String}[] => (Pair{String,String},)
    F.GLOBAL_VARS["section_counter"] = 1 => (Int,)
    counter = F.globvar("section_counter")
    c = lxmock(raw"""\begin{:section, title="aaa"}""")
    h = lx_begin(c,0)
    @test F.globvar("section_counter") == counter + 1
    @test isapproxstr(strip(h), """
        ~~~

                  aaa

        ~~~""")
    c = lxmock(raw"""\end{:section}""")
    h = lx_end(c,0)
    @test isapproxstr(h, """~~~

          ~~~""")
end

@testset "columns" begin
    c = lxmock(raw"""\begin{:columns}""")
    h = lx_begin(c,0)
    @test isapproxstr(h, """~~~

        ~~~""")
    c = lxmock(raw"""\end{:columns}""")
    h = lx_end(c,0)
    @test isapproxstr(h, """~~~

        ~~~""")
end


text:
function mcs_m!(g::AbstractGraph)
    origin_edges = edges(g)
    weights = zeros(Int, nv(g))
    order = zeros(Int, nv(g))
    unumbered = collect(vertices(g))
    for i in nv(g):-1:1
        old_weights = copy(weights)
        z = argmax(weights)
        order[i] = z
        filter!(x->x != z, unumbered)
        weights[z] = -1
        for y in unumbered
            lower_weight = filter(x -> old_weights[x] < old_weights[y], unumbered)
            sg, vmap = induced_subgraph(g, vcat(lower_weight, [z, y]))
            if has_path(sg, findall(==(y), vmap)[1], findall(==(z), vmap)[1])
                weights[y] += 1
                add_edge!(g, y, z)
            end
        end
    end
    data = Dict("order" => order,
                "added_edges" => setdiff(edges(g), origin_edges))
    return data
end

function mcs_m(g::AbstractGraph; kwargs...)
    cg = SimpleGraph(g)
    data = MCS_M!(cg)
    return cg, data
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
"""
    struct InfoLine([  block_name="", data_type=Float32, n_dim=Int32(0),
                                is_present=zeros(Int32, 6) ])

Contains the data of a single entry in the `INFO` block of a Gadget snapshot.

# Fields
| Name                                 | Meaning                                                                                |
| :----------------------------------  | :------------------------------------------------------------------------------------- |
| `block_name::String`                 | name of the data block, e.g. "POS"                                                     |
| `data_type::DataType`                | datatype of the block, e.g. Float32 for single precision, Float64 for double           |
| `n_dim::Int32`                       | number of dimensions of the block, usually 1 or 3                                      |
| `is_present::Vector{Int32}`          | array of flags for which particle type this block is present,                          |
|                                      |  e.g. gas only:  [ 1, 0, 0, 0, 0, 0 ],                                                 |
|                                      |  or gas + BHs: [ 1, 0, 0, 0, 0, 1 ]                                                    |

"""
struct InfoLine
    block_name::String              # name of the data block, e.g. "POS"
    data_type::DataType             # datatype of the block, e.g. Float32 for single precision, Float64 for double
    n_dim::Int32                    # number of dimensions of the block, usually 1 or 3
    is_present::Vector{Int32}       # array of flags for which particle type this block is present,
                                    # e.g. gas only:  [ 1, 0, 0, 0, 0, 0 ]
                                    # e.g. gas + BHs: [ 1, 0, 0, 0, 0, 1 ]

    function InfoLine(block_name="", data_type=Float32, n_dim=Int32(0),
                        is_present=zeros(Int32, 6))

        new(block_name, data_type, n_dim, is_present)
    end
end

"""
    read_info(filename; verbose::Bool=false)

Reads the info block of a snapshot and returns the information in an array of `InfoLine` types.
If `verbose=true` the blocknames are also printed to console.
"""
function read_info(filename::String; verbose::Bool=false)

    if !isfile(filename)
        filename = select_file(filename, 0)
    end

    f = open(filename)
    blocksize = read(f, Int32)

    while eof(f) != true

        # read block name
        blockname = read_bockname(f)

        p = position(f)

        seek(f,p+8)

        skipsize = read(f, UInt32)

        skipsize = check_blocksize(f, p, skipsize)

        if blockname == "INFO"
            seek(f,p+12)
            n_blocks = Int(skipsize/40) # one info line is 40 bytes
            arr_info = Array{InfoLine,1}(undef,n_blocks)

            for i = 1:n_blocks
                arr_info[i] = read_info_line(f)
            end # for

            close(f)

            if verbose == true
                println("Found Info block.\nEntries are:\n")
                for i ∈ 1:length(arr_info)
                    println(i, " - ", arr_info[i].block_name)
                end # for
            end # verbose

            return arr_info

        else
            seek(f,p+skipsize+20)
        end # if blockname == "INFO"

    end # while eof(f) != true

    println("No info block present!")

    return 1
end

"""
    read_info_line(f::IOStream)

Helper function to read the binary data into a `InfoLine` struct.
"""
function read_info_line(f::IOStream)

    block_name = read_bockname(f)

    # the datatype is stored as a 8 char word.
    letters = read!(f, Array{Int8,1}(undef,8))
    letters = Char.(letters)
    letters = string.(letters)

    data_type = ""
    for i = 1:length(letters)
        data_type *= letters[i]
    end

    # erase whitespace and make lower case for comparison
    data_type = lowercase(strip(data_type))

    # assign data types
    if data_type == "float" ||
       data_type == "floatn"
            dt = Float32
    elseif data_type == "long"
            dt = UInt32
    elseif data_type == "llong"
            dt = UInt64
    elseif data_type == "double" ||
           data_type == "doublen"
            dt = Float64
    end

    # array dimensions are stored as Int32
    n_dim = read(f,Int32)

    # read the information for which particles the block is relevant.
    # e.g. for gas particles: [ 1, 0, 0, 0, 0, 0 ]
    is_present = read!(f, Array{Int32,1}(undef,6))

    # construct the info line struct and return it.
    return InfoLine(block_name, dt, n_dim, is_present)
end


text:
module GLSBenchmarks

using GeneralizedLeastSquares: coef, LinRegQR, LinRegCholesky, LinRegBunchKaufman, LinRegSweep
using LinearAlgebra
using BenchmarkTools
using Statistics
using UnicodePlots
using OrderedCollections: OrderedDict

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 5

function sim(n, p, w=nothing)
    title="n=$n, p=$p, w=$(typeof(w))"
    x = randn(n, p)
    errs = isnothing(w) ? randn(n) : cholesky(w).U * randn(n)
    y = x * (1:p) + errs

    args = isnothing(w) ? (x, y) : (x, y, w)

    out = OrderedDict{String, BenchmarkTools.Trial}()

    out["LinRegSweep"] =        @benchmark LinRegSweep($args...)
    out["LinRegQR"] =           @benchmark LinRegQR($args...)
    out["LinRegCholesky"] =     @benchmark LinRegCholesky($args...)
    out["LinRegBunchKaufman"] = @benchmark LinRegBunchKaufman($args...)

    data = round.(map(x -> median(x).time, values(out)) ./ 10^9, digits=4)
    b = barplot(collect(keys(out)), data; width=70, title)
    show(b)
    println()
end

sim(10^6, 50)
sim(10^6, 50, Diagonal(rand(10^6)))

sim(1000, 500)
sim(1000, 500, Diagonal(rand(1000)))

end


text:
# fill and reverse Cint array with padding etc. for cudnn calls
function cdsize(w, nd)
    if isa(w,Number)
        fill(Cint(w),nd)
    elseif length(w)==nd
        [ Cint(w[nd-i+1]) for i=1:nd ]
    else
        throw(DimensionMismatch("$w $nd"))
    end
end

# convert padding etc. size to an Int array of the right dimension
function psize(p, x)
    nd = ndims(x)-2
    if isa(p,Number)
        fill(Int(p),nd)
    elseif length(p)==nd
        collect(Int,p)
    else
        throw(DimensionMismatch("psize: $p $nd"))
    end
end

DT(::CuArray{Float32})=UInt32(0)
DT(::CuArray{Float64})=UInt32(1)
DT(::CuArray{Float16})=UInt32(2)

function cdims(w,x; padding=0, stride=1, o...)
    N = ndims(x)
    ntuple(N) do i
        if i < N-1
            pi = (if isa(padding,Number); padding; else padding[i]; end)
            si = (if isa(stride,Number); stride; else stride[i]; end)
            1 + div(size(x,i) - size(w,i) + 2*pi, si)
        elseif i == N-1
            size(w,N)
        else # i == N
            size(x,N)
        end
    end
end

function pdims(x; window=2, padding=0, stride=window, o...)
    N = ndims(x)
    ntuple(N) do i
        if i < N-1
            wi = (if isa(window,Number); window; else window[i]; end)
            pi = (if isa(padding,Number); padding; else padding[i]; end)
            si = (if isa(stride,Number); stride; else stride[i]; end)
            1 + div(size(x,i) + 2*pi - wi, si)
        else
            size(x,i)
        end
    end
end

function dcdims(w,x; padding=0, stride=1, o...)
    N = ndims(x)
    ntuple(N) do i
        if i < N-1
            pi = (if isa(padding,Number); padding; else padding[i]; end)
            si = (if isa(stride,Number); stride; else stride[i]; end)
            si*(size(x,i)-1) + size(w,i) - 2*pi
        elseif i == N-1
            size(w,N)
        else
            size(x,N)
        end
    end
end

function updims(x; window=2, padding=0, stride=window, o...)
    window = psize(window,x)
    stride = psize(stride,x)
    padding = psize(padding,x)
    N = ndims(x)
    ntuple(N) do i
        if i < N-1
            (size(x,i)-1)*stride[i]+window[i]-2*padding[i]
        else
            size(x,i)
        end
    end
end

# convolution padding size that preserves the input size when filter size is odd and stride=1
padsize(w)=ntuple(i->div(size(w,i)-1,2), ndims(w)-2)


text:
"""
    MinTimeFinal <: AbstractSequencingAlgorithm

The sequence construction proceeds in two stages:

1. Compute an ordering of the zones to be explored based on the spatial distances between
   the polygons obtained by taking the convex hull of the stops of each zone.

2. Once an ordering of the zones has been obtained, we compute the ordering of the stops within each zone.
   The optimization is such that the first stop of the current zone is optimized (to minimize travel time) together
   with all the stops of the next zone.

The function used in step 1 is `minzone`.
The function used in step 2 is `minfunc`.
"""
struct MinTimeFinal{F, Z} <: AbstractSequencingAlgorithm
    minfunc::F
    minzone::Z
end

MinTimeFinal() = MinTimeFinal(_min_length_path_mip, _min_zone)

_travel_times(T::AbstractMatrix) = T
_travel_times(T::TravelTimes) = T.travel_times

function solve(alg::MinTimeFinal, route::AbstractRoute, times;
               optimizer, initial_stop=station(route))

    T = _travel_times(times)
    S = stops(route)
    num_stops = length(S)
    minfunc = alg.minfunc
    minzone = alg.minzone

    # salida
    route_id = name(route)
    zones_list = minzone(route)

    seq = Vector{BareStop}()
    out = Sequence(route_id, seq)

    x = initial_stop |> name
    push!(out, BareStop(x))

    for Z in zones_list
        # FIXME use split_by_zone
        stops_per_zone = stops(route)[findall(x->zone(x)==Z, stops(route))]

        if length(stops_per_zone) == 1
            push!(out, BareStop(name(stops_per_zone[1])))

        else
            insert!(stops_per_zone, 1, stops(route)[findfirst(x->x.id==initial_stop.id, stops(route))])
            stops_str = [name(stops) for stops in stops_per_zone]

            _, result = minfunc(T, stops_per_zone, optimizer=optimizer)

            @inbounds for i in 2:length(result)
                push!(out, BareStop(name(result[i])))
            end

            # la parada inicial de la zona siguiente es aquella donde termine en la zona actual
            initial_stop = result[end] # stops(route)[findfirst(x->x.id==result[length(result)].id, stops(route))]
        end
    end
    return out
end


text:
using ModelConstructors, FileIO, Random, SMC, JLD2

### Construct a generic model and populate it with parameters
reg = GenericModel()
fn = dirname(@__FILE__)
reg <= Setting(:dataroot, "$(fn)/save/input_data/")
reg <= Setting(:saveroot, "$(fn)/save/")

reg <= parameter(:α1, 0., (-1e5, 1e5), (-1e5, 1e5), Untransformed(), Normal(0, 10), fixed = false)
reg <= parameter(:β1, 0., (-1e5, 1e5), (-1e5, 1e5), Untransformed(), Normal(0, 10), fixed = false)

### Estimate with SMC

## Make data
X = rand(100)
β = 1.
α = 1.
y = β*X + α
N = 1
data = hcat(y, X)
jldopen(dataroot(reg)*"reg_data.jld2", "w") do file
    file["data"] = data
end

## Construct likelihood function:
function likelihood_fnct(p, d)
    α = p[1]
    β = p[2]
    Σ = 1 #ones(N,N)
    det_Σ = det(Σ)
    inv_Σ = inv(Σ)
    term1 = -N / 2 * log(2 * π) - 1 /2 * log(det_Σ)
    logprob = 0.
    errors = d[:, 1] .- α .- β .* d[:, 2]
    for t in 1:size(d,1)
        logprob += term1 - 1/2 * dot(errors, inv_Σ * errors)
    end
    return logprob
end

Random.seed!(1793)
println("Starting to estimate Regression with SMC . . .")
@everywhere using SMC
smc(likelihood_fnct, reg.parameters, data, n_parts = 10, use_fixed_schedule = false, tempering_target = 0.97)


text:
# Support for the DAVIS240 camera.
# DAVIS240 has three different event types.

const DAVIS240_SIZE = (180,240)

""" DVS event (from event-based sensor) """
struct DAVIS240_DVS
    y::Int16
    x::Int16
    ext::Bool
    pol::Bool
end

""" APS packet (from frame-based sensor) """
struct DAVIS240_APS
    y::Int16
    x::Int16
    kind::Int8
    sample::Int16
end

""" IMU packet """
struct DAVIS240_IMU
    channel::Int8
    sample::Int8
end

DAVIS240_Any = Union{DAVIS240_DVS, DAVIS240_APS, DAVIS240_IMU}

""" Extracts the event subtype: DVS, APS or IMU """
event_subtype(data::UInt32) = if ((data >> 31) & 0x1 == 0) DAVIS240_DVS
    elseif ((data >> 10) & 0x03) == 0x03 DAVIS240_IMU
    else DAVIS240_APS
end

# Events can be interpreted as a specific class
# if they have the corresponding subtype:
isevent(T::Type{<:DAVIS240_Any}, e::Event{UInt32}) =
    event_subtype(e.address) === T

event_pol(e::DAVIS240_DVS) = e.pol
event_coord(e::Union{DAVIS240_DVS,DAVIS240_APS}) = (e.x, e.y)
event_location(e::Union{DAVIS240_DVS,DAVIS240_APS}, imsize) = (DAVIS240_SIZE[1]-e.y, e.x+1)

image_size(e::Type{<:Union{DAVIS240_DVS,DAVIS240_APS}}) = DAVIS240_SIZE

function Base.convert(::Type{DAVIS240_DVS}, a::UInt32)
    y = (a >> 22) & 0x01FF
    x = (a >> 12) & 0x01FF
    ext = (a >> 10) & 0x1
    pol = (a >> 11) & 0x1
    DAVIS240_DVS(y, x, ext, pol)
end

function Base.convert(::Type{DAVIS240_APS}, a::UInt32)
    y = (a >> 22) & 0x01FF
    x = (a >> 12) & 0x01FF
    kind = (a >> 10) & 0x03
    sample = a & 0x03FF
    DAVIS240_APS(y, x, kind, sample)
end

function Base.convert(::Type{DAVIS240_IMU}, a::UInt32)
    channel = (a >> 28) & 0x07
    sample = (a >> 12) & 0x7fff
    DAVIS240_IMU(channel, sample)
end


text:
loadfish(filename)::Vector{Int} = parse.(Int, split(readchomp(filename), ','))
examplefish = loadfish("example.txt")
inputfish = loadfish("input.txt")

descendants = let
    cache = Dict{Int,Int}()
    function(age)
        get!(cache, age) do
            kids = age - 9:-7:0
            if length(kids) == 0
                0
            else
                length(kids) + sum(descendants.(kids))
            end
        end
    end
end

part2(fish, days) = length(fish) + sum(descendants.(days + 8 .- fish))

using Test
@testset begin
    @test part2(examplefish, 18) == 26
    @test part2(examplefish, 80) == 5934
    @test part2(examplefish, 256) == 26984457539
    println("The answer to part 1 is: ", part2(inputfish, 80))
    println("The answer to part 2 is: ", part2(inputfish, 256))
end

# Part 2 can be done very efficiently by using matrix exponentiation.
function matrixpart2(fish, age)
    F = collect(count(==(n), fish) for n in 0:8)
    P = [
        0 1 0 0 0 0 0 0 0
        0 0 1 0 0 0 0 0 0
        0 0 0 1 0 0 0 0 0
        0 0 0 0 1 0 0 0 0
        0 0 0 0 0 1 0 0 0
        0 0 0 0 0 0 1 0 0
        1 0 0 0 0 0 0 1 0
        0 0 0 0 0 0 0 0 1
        1 0 0 0 0 0 0 0 0
    ]
    sum(P ^ age * F)
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code