text:
module PrinzSingleCompartment

using Test, OrdinaryDiffEq

include(normpath(@__DIR__, "..", "demo", "prinz_kinetics.jl"))

@testset "Prinz Single-Compartment Neuron" begin

channels = [NaV(100mS/cm^2),
            CaT(0mS/cm^2),
            CaS(4mS/cm^2),
            KA(0mS/cm^2),
            KCa(15mS/cm^2),
            Kdr(50mS/cm^2),
            H(.02mS/cm^2),
            leak(.03mS/cm^2)]

@named neuron = Soma(channels, gradients, area = area, V0 = -50mV, aux = [calcium_conversion]);

@test length.([equations(neuron.sys),
               states(neuron.sys),
               parameters(neuron.sys)]) == [39,40,19]

t = 2000
simul_sys = Simulation(neuron, time = t*ms, system = true)

@test length.([equations(simul_sys),
               states(simul_sys),
               parameters(simul_sys)]) == [15,15,17]

# Prinz STG neuron hand-written reference implementation

# Steady-state functions
# Input vectors
V0 = -50.0
u0 = zeros(Float64, 13)
p  = zeros(Float64,17)

# Initial state
u0[1]  = 0.05 # Caᵢ
u0[2]  = V0 # Vₘ
u0[3]  = (1.0 + exp(-4.82042 - (0.189036V0)))^-1 # NaV₊m
u0[4]  = (1.0 + exp(9.44015 + 0.19305V0))^-1 # NaV₊h
u0[5]  = (1.0 + exp(-3.76389 - (0.138889V0)))^-1 # CaT₊m
u0[6]  = (1.0 + exp(5.83636 + 0.181818V0))^-1 # CaT₊h
u0[7]  = (1.0 + exp(-4.07407 - (0.123457V0)))^-1 # CaS₊m
u0[8]  = (1.0 + exp(9.67742 + 0.16129V0))^-1 # CaS₊h
u0[9]  = (1.0 + exp(-3.12644 - (0.114943V0)))^-1 # KA₊m
u0[10] = (1.0 + exp(11.6122 + 0.204082V0))^-1 # KA₊h
u0[11] = 0.05*((3.0 + 0.05)^-1)*((1.0 + exp(-2.24603 - (0.0793651V0)))^-1) # KCa₊m
u0[12] = (1.0 + exp(-1.04237 - (0.0847458V0)))^-1 # Kdr₊m
u0[13] = (1.0 + exp(13.6364 + 0.181818V0))^-1 # H₊m

# Parameters
p[1]  = 0.001 # cₘ
p[2]  = 0.000628 # aₘ
p[3]  = 200 # τCa
p[4]  = 0.05 # Ca∞
p[5]  = 14.96 # f
p[6]  = 50.0 # ENa
p[7]  = -80.0 # EK
p[8]  = -20.0 # EH
p[9]  = -50.0 # El
p[10] = 100.0 # NaV₊gbar
p[11] = 0.0 # CaT₊gbar
p[12] = 4.0 # CaS₊gbar
p[13] = 0.0 # KA₊gbar
p[14] = 15.0 # KCa₊gbar
p[15] = 50.0 # Kdr₊gbar
p[16] = 0.02 # H₊gbar
p[17] = 0.03 # leak₊g

function prinz_neuron!(du, u, p, t)
    @inbounds let Caᵢ         = u[1],
                  Vₘ          = u[2],
                  NaV₊m       = u[3],
                  NaV₊h       = u[4],
                  CaT₊m       = u[5],
                  CaT₊h       = u[6],
                  CaS₊m       = u[7],
                  CaS₊h       = u[8],
                  KA₊m        = u[9],
                  KA₊h        = u[10],
                  KCa₊m       = u[11],
                  Kdr₊m       = u[12],
                  H₊m         = u[13],

                  cₘ          = p[1],
                  aₘ          = p[2],
                  τCa         = p[3],
                  Ca∞         = p[4],
                  f           = p[5],
                  ENa         = p[6],
                  EK          = p[7],
                  EH          = p[8],
                  El          = p[9],
                  NaV₊gbar    = p[10],
                  CaT₊gbar    = p[11],
                  CaS₊gbar    = p[12],
                  KA₊gbar     = p[13],
                  KCa₊gbar    = p[14],
                  Kdr₊gbar    = p[15],
                  H₊gbar      = p[16],
                  leak₊g      = p[17],
                  ECa = (Vₘ + 12.199306599503455 * log(0.0003333333333333333 * max(Caᵢ, 0.001)))
        let
            du[1] =
                (
                    Ca∞ - Caᵢ - 1000 * f *
                    (CaS₊gbar * aₘ * CaS₊h * CaS₊m^3 *
                    ECa +
                    CaT₊gbar * aₘ * CaT₊h * CaT₊m^3 *
                    ECa)
                   ) * inv(τCa)

            du[2] =
                (
                    -aₘ * leak₊g * (Vₘ - El) +
                    -H₊gbar * aₘ * H₊m * (Vₘ - EH) +
                    -KCa₊gbar * aₘ * KCa₊m^4 * (Vₘ - EK) +
                    -Kdr₊gbar * aₘ * Kdr₊m^4 * (Vₘ - EK) +
                    -CaS₊gbar * aₘ * CaS₊h * CaS₊m^3 * ECa +
                    -CaT₊gbar * aₘ * CaT₊h * CaT₊m^3 * ECa +
                    -KA₊gbar * aₘ * KA₊h * KA₊m^3 * (Vₘ - EK) +
                    -NaV₊gbar * aₘ * NaV₊h * NaV₊m^3 * (Vₘ - ENa)
                ) * inv(aₘ) * inv(cₘ)

            du[3] =
                (-NaV₊m + inv(1.0 + exp(-4.8204158790170135 - 0.1890359168241966Vₘ))) *
                inv(2.64 - 2.52 * inv(1 + exp(-4.8 - 0.04Vₘ)))

            du[4] =
                0.7462686567164178 * (1.0 + exp(-6.29 - 0.1Vₘ)) *
                (-NaV₊h + inv(1.0 + exp(9.44015444015444 + 0.19305019305019305Vₘ))) *
                inv(1.5 + inv(1.0 + exp(9.694444444444445 + 0.2777777777777778Vₘ)))

            du[5] =
                (-CaT₊m + inv(1.0 + exp(-3.7638888888888893 + -0.1388888888888889Vₘ))) *
                inv(43.4 - 42.6 * inv(1.0 + exp(-3.321951219512195 - 0.04878048780487805Vₘ)))

            du[6] =
                (-CaT₊h + inv(1.0 + exp(5.836363636363637 + 0.18181818181818182Vₘ))) *
                inv(210.0 - 179.6 * inv(1.0 + exp(-3.2544378698224854 - 0.0591715976331361Vₘ)))

            du[7] =
                (-CaS₊m + inv(1.0 + exp(-4.074074074074074 - 0.1234567901234568Vₘ))) *
                inv(2.8 + 14.0 * inv(exp(-5.384615384615385 - 0.07692307692307693Vₘ) + exp(2.7 + 0.1Vₘ)))

            du[8] =
                (-CaS₊h + inv(1.0 + exp(9.67741935483871 + 0.16129032258064516Vₘ))) *
                inv(120.0 + 300.0 * inv(exp(-4.0625 - 0.0625Vₘ) +
                                        exp(6.111111111111111 + 0.1111111111111111Vₘ)))

            du[9] =
                (-KA₊m + inv(1.0 + exp(-3.1264367816091956 - 0.1149425287356322Vₘ))) *
                inv(23.2 - 20.8 * inv(1.0 + exp(-2.164473684210526 - 0.06578947368421052Vₘ)))

            du[10] =
                (-KA₊h + inv(1.0 + exp(11.612244897959181 + 0.2040816326530612Vₘ))) *
                inv(77.2 - 58.4 * inv(1.0 + exp(-1.4679245283018867 - 0.03773584905660377Vₘ)))
            du[11] =
                (-KCa₊m + Caᵢ * inv(3.0 + Caᵢ) *
                inv(1.0 + exp(-2.246031746031746 - 0.07936507936507936Vₘ))) *
                inv(180.6 + -150.2 * inv(1.0 + exp(-2.026431718061674 - 0.04405286343612335Vₘ)))

            du[12] =
                (-Kdr₊m + inv(1.0 + exp(-1.0423728813559323 - 0.0847457627118644Vₘ))) *
                inv(14.4 - 12.8 * inv(1.0 + exp(-1.4739583333333335 - 0.052083333333333336Vₘ)))

            du[13] =
                1 // 2 *
                (exp(-1.8671328671328669 + 0.06993006993006992Vₘ) +
                 exp(-14.629310344827585 - 0.08620689655172414Vₘ)) *
                (-H₊m + inv(1.0 + exp(13.636363636363637 + 0.18181818181818182Vₘ)))
            nothing
        end
    end
end

byhand_prob = ODEProblem{true}(prinz_neuron!, u0, (0.,2000.), p)
mtk_prob = ODAEProblem(simul_sys, [], (0., t), [])
byhand_sol = solve(byhand_prob, Rosenbrock23(), reltol=1e-8, abstol=1e-8, saveat=0.025);
current_mtk_sol = solve(mtk_prob, Rosenbrock23(), reltol=1e-8, abstol=1e-8, saveat=0.025);

tsteps = 0.0:0.025:2000.0
byhand_out = Array(byhand_sol(tsteps, idxs=2))
current_mtk_out = current_mtk_sol(tsteps)[Vₘ]

@test isapprox(byhand_out, current_mtk_out, rtol=0.001)

end # testset
end # module


text:
function ueberpruefeSatzvervollstaendigung2(c,d)

if cmp(c,"senkrecht") == 0 #Null bedeutet richtig, stimmt überein
    if cmp(d,"Normalenvektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"normalenvektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"normalen vektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"Normalen Vektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"normalen Vektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"Normalen vektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    else
        println("Eure Antwort für Wort1 ist korrekt, nicht aber die für Wort2. Lest euch noch einmal den Satz durch und überprüft eure Rechtschreibung.");
    end
elseif cmp(c,"Senkrecht") == 0 #Null bedeutet richtig, stimmt überein
    if cmp(d,"Normalenvektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"normalenvektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"normalen vektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"Normalen Vektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"normalen Vektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"Normalen vektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    else
        println("Eure Antwort für Wort1 ist korrekt, nicht aber die für Wort2. Lest euch noch einmal den Satz durch und überprüft eure Rechtschreibung.");
    end
elseif cmp(c,"orthogonal") == 0 #Null bedeutet richtig, stimmt überein
    if cmp(d,"Normalenvektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"normalenvektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"normalen vektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"Normalen Vektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"normalen Vektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"Normalen vektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    else
        println("Eure Antwort für Wort1 ist korrekt, nicht aber die für Wort2. Lest euch noch einmal den Satz durch und überprüft eure Rechtschreibung.");
    end
elseif cmp(c,"Orthogonal") == 0 #Null bedeutet richtig, stimmt überein
    if cmp(d,"Normalenvektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"normalenvektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"normalen vektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"Normalen Vektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"normalen Vektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    elseif cmp(d,"Normalen vektor") == 0
        println("Eure Antwort für Wort1 und Wort2 sind beide korrekt!");
    else
        println("Eure Antwort für Wort1 ist korrekt, nicht aber die für Wort2. Lest euch noch einmal den Satz durch und überprüft eure Rechtschreibung.");
    end
else
    if cmp(d,"Normalenvektor") == 0
        println("Eure Antwort für Wort2 ist korrekt, nicht aber für Wort1.");
    elseif cmp(d,"normalenvektor") == 0
        println("Eure Antwort für Wort2 ist korrekt, nicht aber für Wort1.");
    elseif cmp(d,"normalen vektor") == 0
        println("Eure Antwort für Wort2 ist korrekt, nicht aber für Wort1.");
    elseif cmp(d,"Normalen Vektor") == 0
        println("Eure Antwort für Wort2 ist korrekt, nicht aber für Wort1.");
    elseif cmp(d,"normalen Vektor") == 0
        println("Eure Antwort für Wort2 ist korrekt, nicht aber für Wort1.");
    elseif cmp(d,"Normalen vektor") == 0
        println("Eure Antwort für Wort2 ist korrekt, nicht aber für Wort1.");
    else
        println("Eure Antworten für Wort1 und Wort2 sind noch nicht korrekt. Lest euch noch einmal den Satz durch und überprüft eure Rechtschreibung.");
    end
end

end

ext: jl


text:
"""
Constant to be used in order to indicate a hex-like extrusion axis.
"""
const HEX_AXIS = 1

"""
Constant to be used in order to indicate a tet-like extrusion axis.
"""
const TET_AXIS = 2

struct NFace{D}
  anchor::Point{D,Int}
  extrusion::Point{D,Int}
end

struct DFace{D} <: GridapType
  extrusion::Point{D,Int}
  nfaces::Vector{NFace{D}}
  dimranges::Vector{UnitRange{Int}}
  dims::Vector{Int}
  nf_nfs::Vector{Vector{Int}}
  nf_dimranges::Vector{Vector{UnitRange{Int}}}
  nf_dims::Vector{Vector{Int}}
end

"""
    struct ExtrusionPolytope{D} <: Polytope{D}
      extrusion::Point{D,Int}
      # + private fields
    end

Concrete type for polytopes that can be
represented with an "extrusion" tuple. The underlying extrusion is available
in the field `extrusion`. Instances of this type can be obtained with the constructors

- [`Polytope(extrusion::Int...)`](@ref)
- [`ExtrusionPolytope(extrusion::Int...)`](@ref)

"""
struct ExtrusionPolytope{D} <: Polytope{D}
  extrusion::Point{D,Int}
  dface::DFace{D}
  vertex_coords::Vector{Point{D,Float64}}
  face_normals::Vector{Point{D,Float64}}
  face_orientations::Vector{Int}
  vertex_perms::Vector{Vector{Int}}
  face_vertex_perms::Vector{Vector{Vector{Int}}}
  m_n_to_mface_to_nface::Matrix{Vector{Vector{Int}}}
  function ExtrusionPolytope(dface::DFace{D}) where D
    vertex_coords = _vertices_coordinates(Float64,dface)
    face_normals, face_orientations = _face_normals(Float64,dface)
    vertex_perms = _precompute_vertex_perms_if_possible(dface)
    face_vertex_perms = _admissible_face_vertex_permutations(dface)
    m_n_to_mface_to_nface = _precompute_m_n_to_mface_to_nface(dface)
    new{D}(
      dface.extrusion,dface,vertex_coords,face_normals,
      face_orientations,vertex_perms,face_vertex_perms,
      m_n_to_mface_to_nface)
  end
end

# Constructors

"""
    Polytope(extrusion::Int...)

Equivalent to `ExtrusionPolytope(extrusion...)`
"""
Polytope(extrusion::Int...) = Polytope(extrusion)

function Polytope(extrusion::NTuple{D,Int}) where D
  ExtrusionPolytope(extrusion)
end

"""
    ExtrusionPolytope(extrusion::Int...)

Generates an `ExtrusionPolytope` from the tuple `extrusion`.
The values in `extrusion` are either equal to the constant
[`HEX_AXIS`](@ref) or the constant [`TET_AXIS`](@ref).

# Examples

Creating a quadrilateral, a triangle, and a wedge

```jldoctest
using Gridap.ReferenceFEs

quad = ExtrusionPolytope(HEX_AXIS,HEX_AXIS)

tri = ExtrusionPolytope(TET_AXIS,TET_AXIS)

wedge = ExtrusionPolytope(TET_AXIS,TET_AXIS,HEX_AXIS)

println(quad == QUAD)
println(tri == TRI)
println(wedge == WEDGE)

# output
true
true
true

```

"""
function ExtrusionPolytope(extrusion::Int...)
  ExtrusionPolytope(extrusion)
end

function ExtrusionPolytope(extrusion::NTuple{N,Int}) where N
  ExtrusionPolytope(Point{N,Int}(extrusion))
end

function ExtrusionPolytope(extrusion::Point{N,Int}) where N
  dface = DFace(extrusion)
  ExtrusionPolytope(dface)
end

# Implementation of the interface

get_faces(p::ExtrusionPolytope) = p.dface.nf_nfs

function get_faces(p::ExtrusionPolytope,dimfrom::Integer,dimto::Integer)
  p.m_n_to_mface_to_nface[dimfrom+1,dimto+1]
end

get_dimranges(p::ExtrusionPolytope) = p.dface.dimranges

get_facedims(p::ExtrusionPolytope) = p.dface.dims

function get_face_dimranges(p::ExtrusionPolytope)
  p.dface.nf_dimranges
end

function get_face_dimranges(p::ExtrusionPolytope,d::Integer)
  r = get_dimrange(p,d)
  p.dface.nf_dimranges[r]
end

function Polytope{D}(p::ExtrusionPolytope,Dfaceid::Integer) where D
  offset = get_offset(p,D)
  id = Dfaceid + offset
  dface = DFace{D}(p.dface,id)
  ExtrusionPolytope(dface)
end

function Polytope{0}(p::ExtrusionPolytope,Dfaceid::Integer)
  @assert Dfaceid in 1:num_vertices(p)
  VERTEX
end

function Polytope{D}(p::ExtrusionPolytope{D},Dfaceid::Integer) where D
  @assert Dfaceid == 1
  p
end

function (==)(a::ExtrusionPolytope{D},b::ExtrusionPolytope{D}) where D
  #The first axis is irrelevant here
  ea = Point(Tuple(a.extrusion)[2:end])
  eb = Point(Tuple(b.extrusion)[2:end])
  ea == eb
end

function (==)(a::ExtrusionPolytope{1},b::ExtrusionPolytope{1})
  true
end

function (==)(a::ExtrusionPolytope{0},b::ExtrusionPolytope{0})
  true
end

function get_vertex_coordinates(p::ExtrusionPolytope)
  p.vertex_coords
end

function get_edge_tangent(p::ExtrusionPolytope)
  _edge_tangents(Float64,p.dface)
end

function get_facet_normal(p::ExtrusionPolytope)
  p.face_normals
end

function get_facet_orientations(p::ExtrusionPolytope)
  p.face_orientations
end

function get_vertex_permutations(p::ExtrusionPolytope)
  p.vertex_perms
end

function get_face_vertex_permutations(p::ExtrusionPolytope,d::Integer)
  r = get_dimrange(p,d)
  perms = get_face_vertex_permutations(p)
  perms[r]
end

function get_face_vertex_permutations(p::ExtrusionPolytope)
  p.face_vertex_perms
end

function is_simplex(p::ExtrusionPolytope)
  all(Tuple(p.extrusion) .== TET_AXIS)
end

function is_n_cube(p::ExtrusionPolytope)
  all(Tuple(p.extrusion) .== HEX_AXIS)
end

function is_simplex(p::ExtrusionPolytope{0})
  true
end

function is_n_cube(p::ExtrusionPolytope{0})
  true
end

function is_simplex(p::ExtrusionPolytope{1})
  true
end

function is_n_cube(p::ExtrusionPolytope{1})
  true
end

function simplexify(p::ExtrusionPolytope)
  @notimplemented
end

function simplexify(p::ExtrusionPolytope{0})
  [[1,],], VERTEX
end

function simplexify(p::ExtrusionPolytope{1})
  [[1,2],], SEGMENT
end

function simplexify(p::ExtrusionPolytope{2})
  if p == QUAD
    ([[1,2,3],[2,3,4]], TRI)
  elseif p == TRI
    ([[1,2,3],], TRI)
  else
     @notimplemented
  end
end

function simplexify(p::ExtrusionPolytope{3})
  if p == HEX
    simplices = [
      [1,2,3,7], [1,2,5,7], [2,3,4,7],
      [2,4,7,8], [2,5,6,7], [2,6,7,8]]
    (simplices, TET)
  elseif p == TET
    simplices = [[1,2,3,4],]
    (simplices, TET)
  else
     @notimplemented
  end
end

function Base.show(io::IO,p::ExtrusionPolytope)
  if p == VERTEX
    s = "VERTEX"
  elseif p==SEGMENT
    s = "SEGMENT"
  elseif p==QUAD
    s = "QUAD"
  elseif p==TRI
    s = "TRI"
  elseif p==TET
    s = "TET"
  elseif p== HEX
    s = "HEX"
  elseif p==WEDGE
    s = "WEDGE"
  elseif p==PYRAMID
    s = "PYRAMID"
  else
    s = nameof(typeof(p))
  end
  print(io,s)
end

"""
    get_extrusion(p::ExtrusionPolytope)

Equivalent to `p.extrusion`.
"""
get_extrusion(p::ExtrusionPolytope) = p.extrusion

# Helpers for the construction of a DFace

function DFace(extrusion::Point{D,Int}) where D
  zerop = zero(Point{D,Int})
  pol_nfs, pol_dim = _polytopenfaces(zerop, extrusion)
  nfs_id = Dict{NFace{D},Int}(nf => i for (i, nf) in enumerate(pol_nfs))
  nf_nfs, nf_dimranges = _polytopemesh(pol_nfs, nfs_id)
  dims = _nface_to_nfacedim(length(pol_nfs), pol_dim)
  nf_dims = _nf_dims(nf_nfs,nf_dimranges)
  DFace{D}(extrusion, pol_nfs, pol_dim, dims, nf_nfs, nf_dimranges, nf_dims)
end

# Provides for all n-faces of a polytope the d-faces for 0 <= d <n on its
# boundary (e.g., given a face, it provides the set of edges and corners on its
# boundary) using the global n-face numbering of the base polytope
function _polytopemesh(nfaces::Vector{NFace{D}}, nfaceid::Dict{NFace{D},Int}) where D
  num_nfs = length(nfaces)
  nfnfs = Vector{Vector{Int}}(undef, num_nfs)
  nfnfs_dim = Vector{Vector{UnitRange{Int}}}(undef, num_nfs)
  for (inf, nf) in enumerate(nfaces)
    nf_nfs, dimnfs = _polytopenfaces(nf.anchor, nf.extrusion)
    nfnfs[inf] = [ nfaceid[nf] for nf in nf_nfs]
    nfnfs_dim[inf] = dimnfs
  end

  return (nfnfs, nfnfs_dim)
end

function _nface_to_nfacedim(nnfaces::Int,dimranges::Vector{UnitRange{Int}})
  v = zeros(Int,nnfaces)
  for (d,r) in enumerate(dimranges)
    v[r] .= d-1
  end
  v
end

function _nf_dims(nf_nfs,nf_dimranges)
  nf_dims = Vector{Int}[]
  for i in 1:length(nf_nfs)
    nfs = nf_nfs[i]
    dimranges = nf_dimranges[i]
    dims = _nface_to_nfacedim(length(nfs),dimranges)
    push!(nf_dims,dims)
  end
  nf_dims
end

# Generates the array of n-faces of a polytope
function _polytopenfaces(anchor, extrusion)

  D = num_components(extrusion)
  zerop = zero(Point{D,Int})
  nf_nfs = Vector{NFace{D}}(undef,0)
  _nfaceboundary!(anchor, zerop, extrusion, true, nf_nfs)

  for i in 1:D
    sort!(nf_nfs, by = x -> x.anchor[i])
  end
  for i in 1:D
    sort!(nf_nfs, by = x -> x.extrusion[i])
  end
  sort!(nf_nfs, by = x -> count(p->p!=0,x.extrusion))

  numnfs = length(nf_nfs)
  nfsdim = [_nfdim(nf_nfs[i].extrusion) for i = 1:numnfs]
  dnf = _nfdim(extrusion)
  dimnfs = Array{UnitRange{Int},1}(undef, dnf + 1)
  dim = 0
  i = 1
  for iface = 1:numnfs
    if (nfsdim[iface] > dim)
      # global dim; # global i
      dim += 1
      dimnfs[dim] = i:iface-1
      i = iface
    end
  end
  dimnfs[dnf+1] = numnfs:numnfs

  (nf_nfs, dimnfs)
end

function _nfdim(extrusion)
  z = zero(eltype(extrusion))
  for e in extrusion
    if e > 0
      z +=1
    end
  end
  z
end

# Generates the list of n-face of a polytope the d-faces for 0 <= d <n on its
# boundary
function _nfaceboundary!(anchor, extrusion, extend, isanchor, list)

  D = num_components(extrusion)
  newext = extend
  push!(list,NFace(anchor,extrusion))

  for i = 1:D
    curex = newext[i]
    if (curex > 0) # Perform extension

      newext = _newext(newext,i)
      edim = _edim(newext,i)
      tetp = _tetp(anchor,i)

      if (curex == HEX_AXIS) # Quad extension
        list = _nfaceboundary!(anchor + edim, extrusion, newext, false, list)
      elseif isanchor
        list = _nfaceboundary!(tetp + edim, extrusion, newext, false, list)
      end
      list = _nfaceboundary!( anchor, extrusion + edim * curex, newext, false, list)

    end
  end

  list
end

function _newext(newext,i)
  m = zero(Mutable(newext))
  D = num_components(newext)
  for j in 1:D
    m[j] = j == i ? 0 : newext[j]
  end
  Point(m)
end

function _edim(newext,i)
  m = zero(Mutable(newext))
  D = num_components(newext)
  for j in 1:D
    m[j] = j == i ? 1 : 0
  end
  Point(m)
end

function _tetp(anchor,i)
  m = zero(Mutable(anchor))
  D = num_components(anchor)
  for j in 1:D
    m[j] = j >= i ? anchor[j] : 0
  end
  Point(m)
end

# Helper API on a DFace

function _dimfrom_fs_dimto_fs(p::DFace{D}, dim_from::Int, dim_to::Int) where D
  @assert dim_from <= D
  if dim_to == dim_from
    return [ [i,] for i in 1:_num_nfaces(p,dim_from)]
  end
  @assert dim_to < dim_from
  dim_from += 1
  dim_to += 1
  dffs_r = p.nf_dimranges[end][dim_from]
  dffs_dtfs = Vector{Vector{Int}}(undef, dffs_r[end] - dffs_r[1] + 1)
  offs = p.nf_dimranges[end][dim_to][1] - 1
  for (i_dff, dff) in enumerate(dffs_r)
    dff_nfs = p.nf_nfs[dff]
    dff_dtfs_r = p.nf_dimranges[dff][dim_to]
    dff_dtfs = dff_nfs[dff_dtfs_r]
    dffs_dtfs[i_dff] = dff_dtfs .- offs
  end
  dffs_dtfs
end

function _precompute_m_n_to_mface_to_nface(p::DFace{D}) where D
  m_n_to_mface_to_nface = Matrix{Vector{Vector{Int}}}(undef,D+1,D+1)

  for m in 0:D
    for n in 0:D
      if n <= m
       mface_to_nface = _dimfrom_fs_dimto_fs(p,m,n)
      else
       mface_to_nface = _dimfrom_fs_dimto_fs_dual(p,m,n)
      end
      m_n_to_mface_to_nface[m+1,n+1] = mface_to_nface
    end
  end

  m_n_to_mface_to_nface

end

function _dimfrom_fs_dimto_fs_dual(p::DFace,dimfrom,dimto)
  @assert dimfrom < dimto
  tface_to_ffaces = _dimfrom_fs_dimto_fs(p,dimto,dimfrom)
  nffaces = _num_nfaces(p,dimfrom)
  fface_to_tfaces = [Int[] for in in 1:nffaces]
  for (tface,ffaces) in enumerate(tface_to_ffaces)
    for fface in ffaces
      push!(fface_to_tfaces[fface],tface)
    end
  end
  fface_to_tfaces
end

# iface in global numeration
function DFace{d}(p::DFace{D},iface::Int) where {d,D}
  @assert d == p.dims[iface]
  nf = p.nfaces[iface]
  r_ext = _eliminate_zeros(Val{d}(),nf.extrusion)
  DFace(r_ext)
end

function DFace{D}(p::DFace{D},iface::Int) where D
  @assert D == p.dims[iface]
  p
end

function _eliminate_zeros(::Val{d},a) where d
  b = zero(Mutable(Point{d,Int}))
  D = num_components(a)
  k = 1
  for i in 1:D
    m = a[i]
    if (m != 0)
      b[k] = m
      k += 1
    end
  end
  Point(b)
end

# It generates the list of coordinates of all vertices in the polytope. It is
# assumed that the polytope has the bounding box [0,1]**dim
function _vertices_coordinates(::Type{T},p::DFace{D}) where {D,T}
  vs = first(_dimfrom_fs_dimto_fs(p, D, 0))
  vcs = zeros(Point{D,T},length(vs))
  for i = 1:length(vs)
    vx = p.nfaces[vs[i]]
    vc = Tuple(vx.anchor)
    vcs[i] = vc
  end
  vcs
end

# Returns number of nfaces of dimension dim
function _num_nfaces(polytope::DFace, dim::Integer)
  n = length(polytope.nf_dimranges)
  k = 0
  for nface = 1:n
    d = length(polytope.nf_dimranges[nface]) - 1
    if d == dim
      k += 1
    end
  end
  k
end

function _nfaces_vertices(::Type{T},p::DFace,d::Integer) where T
  nc = _num_nfaces(p,d)
  verts = _vertices_coordinates(T,p)
  faces_vs = _dimfrom_fs_dimto_fs(p,d,0)
  cfvs = collect(lazy_map(Broadcasting(Reindex(verts)),faces_vs))
end

# Return the n-faces vertices coordinates array for a given n-face dimension
function _nfaces_vertices(::Type{T},p::ExtrusionPolytope,d::Integer) where T
  _nfaces_vertices(T,p.dface,d)
end

# It generates the outwards normals of the facets of a polytope. It returns two
# arrays, the first one being the outward normal and the second one the orientation.
function _face_normals(::Type{T},p::DFace{D}) where {D,T}
  nf_vs = _dimfrom_fs_dimto_fs(p, D - 1, 0)
  vs = _vertices_coordinates(T,p)
  n = length(p.nf_dimranges[end][end-1])
  f_ns = zeros(Point{D,T},n)
  f_os = zeros(Int,n)
  for i_f = 1:n
    f_n, f_o = _facet_normal(T,p, nf_vs, vs, i_f)
    f_ns[i_f] = f_n
    f_os[i_f] = f_o
  end
  f_ns, f_os
end

function _face_normals(::Type{T},p::DFace{0}) where T
  (VectorValue{0,T}[], Int[])
end

function _facet_normal(::Type{T},p::DFace{D}, nf_vs, vs, i_f) where {D,T}
  if length(p.extrusion) > 1
    n1 = length(nf_vs[i_f])-1
    n2 = D
    v = zeros(T,n1,n2)
    for i in 2:length(nf_vs[i_f])
      vi = vs[nf_vs[i_f][i]] - vs[nf_vs[i_f][1]]
      for d in 1:num_components(vi)
        v[i-1,d] = vi[d]
      end
    end
    _n = _nullspace(v)
    n = Point{D,T}(_n)
    n = n * 1 / sqrt(dot(n, n))
    ext_v = _vertex_not_in_facet(p, i_f, nf_vs)
    v3 = vs[nf_vs[i_f][1]] - vs[ext_v]
    f_or = 1
    if dot(v3, n) < 0.0
      n *= -1
      f_or = -1
    end
  elseif (length(p.extrusion) == 1)
    ext_v = _vertex_not_in_facet(p, i_f, nf_vs)
    n = vs[nf_vs[i_f][1]] - vs[ext_v]
    n = n .* 1 / dot(n, n)
    f_or = 1
  else
    @unreachable "O-dim polytopes do not have properly define outward facet normals"
  end
  n, f_or
end

function _nullspace(v)
  if size(v,1) == 1
    return [-v[2],v[1]]
  else
    return nullspace(v)
  end
end

function _vertex_not_in_facet(p::DFace, i_f, nf_vs)
  for i in p.nf_dimranges[end][1]
    is_in_f = false
    for j in nf_vs[i_f]
      if i == j
        is_in_f = true
        break
      end
    end
    if !is_in_f
      return i
      break
    end
  end
  -1
end

# It generates the tangent vectors for polytope edges.
function _edge_tangents(::Type{T},p::DFace{D}) where {D,T}
  ed_vs = _nfaces_vertices(T,p,1)
  ts = [(t = vs[2]-vs[1])/norm(t) for vs in ed_vs ]
  ts
end

function _edge_tangents(::Type{T},p::DFace{0}) where T
  VectorValue{0,T}[]
end

function _admissible_face_vertex_permutations(p::DFace{D}) where D
  perms = Vector{Vector{Int}}[]
  _admissible_face_vertex_permutations_fill!(perms,p,Val{0}())
  perms
end

function  _admissible_face_vertex_permutations_fill!(perms,p::DFace{D},::Val{d}) where {D,d}
  faceids = p.dimranges[d+1]
  for faceid in faceids
    f = DFace{d}(p,faceid)
    f_perms = _precompute_vertex_perms_if_possible(f)
    push!(perms,f_perms)
  end
  _admissible_face_vertex_permutations_fill!(perms,p,Val{d+1}())
  nothing
end

function  _admissible_face_vertex_permutations_fill!(perms,p::DFace{D},::Val{D}) where D
  f_perms = _precompute_vertex_perms_if_possible(p)
  push!(perms,f_perms)
  nothing
end

function _precompute_vertex_perms_if_possible(p::DFace{D}) where D
  if D in (0,1)
    perms = _admissible_permutations_simplex(p)
  elseif D == 2
    perms = _admissible_permutations(p)
  else
    perms = _admissible_permutations_identity(p)
  end
  perms
end

# It generates all the admissible permutations of nodes that lead to an
# admissible polytope
function _admissible_permutations(p::DFace{D}) where D
  if D > 3
    @warn "Computing permutations for a polytope of dim > 3 is overkill"
  end
  if D in (0,1) || all( map(i->i==TET_AXIS,Tuple(p.extrusion)[2:end]) )
    perms = _admissible_permutations_simplex(p)
  elseif all( map(i->i==HEX_AXIS,Tuple(p.extrusion)[2:end]))
    perms = _admissible_permutations_n_cube(p)
  else
    @notimplemented "admissible vertex permutations only implemented for simplices and n-cubes"
  end
  perms
end

function _admissible_permutations_identity(p::DFace{D}) where D
  vs = first(_dimfrom_fs_dimto_fs(p, D, 0))
  num_vs = length(vs)
  l = [i for i = 1:num_vs]
  perms = [l,]
  perms
end

function _admissible_permutations_simplex(p::DFace{D}) where D
  vs = first(_dimfrom_fs_dimto_fs(p, D, 0))
  num_vs = length(vs)
  l = [i for i = 1:num_vs]
  perms = Combinatorics.permutations(l)
  collect(perms)
end

function _admissible_permutations_n_cube(p::DFace{D}) where D
  vs = first(_dimfrom_fs_dimto_fs(p, D, 0))
  num_vs = length(vs)
  l = [i for i = 1:num_vs]
  perms = Combinatorics.permutations(l)
  vertices = p.nfaces[first(p.dimranges)]
  permuted_vertices = similar(vertices)
  admissible_perms = Vector{Int}[]
  grads = _setup_aux_grads(vertices)
  vol = -1
  for perm in perms
    for (j,cj) in enumerate(perm)
      permuted_vertices[j] = vertices[cj]
    end
    jac = _setup_aux_jacobian(grads,permuted_vertices)
    vol_i = convert(Int,abs(det(jac)))
    if vol < 0
      vol = vol_i
    end
    if vol_i == vol
      push!(admissible_perms,perm)
    end
  end
  admissible_perms
end

function _setup_aux_grads(vertices::Vector{NFace{D}}) where D
  grads = zeros(Point{D,Int},length(vertices))
  m = zero(Mutable(Point{D,Int}))
  for (i,vertex) in enumerate(vertices)
    x = vertex.anchor
    for di in 1:D
      v = 1
      for k in 1:D
        if k == di && x[k] == 0
          v *= -1
        end
      end
      m[di] = v
    end
    grads[i] = m
  end
  grads
end

function _setup_aux_jacobian(grads,permuted_vertices::Vector{NFace{D}}) where D
  p0 = zero(Point{D,Int})
  m = zero(Mutable(outer(p0,p0)))
  for (i,pvertex) in enumerate(permuted_vertices)
    x = pvertex.anchor
    g = grads[i]
    for di in 1:D
      v = g[di]
      for dj in 1:D
        m[di,dj] += x[dj]*v
      end
    end
  end
  TensorValue(m)
end

# Some particular cases

"""
    const VERTEX = Polytope()
"""
const VERTEX = Polytope()

"""
    const SEGMENT = Polytope(HEX_AXIS)
"""
const SEGMENT = Polytope(HEX_AXIS)

# TODO use larger names
"""
    const TRI = Polytope(TET_AXIS,TET_AXIS)
"""
const TRI = Polytope(TET_AXIS,TET_AXIS)

"""
    const QUAD = Polytope(HEX_AXIS,HEX_AXIS)
"""
const QUAD = Polytope(HEX_AXIS,HEX_AXIS)

"""
    const TET = Polytope(TET_AXIS,TET_AXIS,TET_AXIS)
"""
const TET = Polytope(TET_AXIS,TET_AXIS,TET_AXIS)

"""
    const HEX = Polytope(HEX_AXIS,HEX_AXIS,HEX_AXIS)
"""
const HEX = Polytope(HEX_AXIS,HEX_AXIS,HEX_AXIS)

"""
    const WEDGE = Polytope(TET_AXIS,TET_AXIS,HEX_AXIS)
"""
const WEDGE = Polytope(TET_AXIS,TET_AXIS,HEX_AXIS)

"""
    const PYRAMID = Polytope(HEX_AXIS,HEX_AXIS,TET_AXIS)
"""
const PYRAMID = Polytope(HEX_AXIS,HEX_AXIS,TET_AXIS)


text:
using CompScienceMeshes, BEAST
using Base.Test

o, x, y, z = euclidianbasis(3)

Γ = meshsphere(1.0, 0.11)
X = lagrangecxd0(Γ)
@show numfunctions(X)

κ = 1.0; γ = im*κ
a = Helmholtz3D.singlelayer(gamma=γ)
b = Helmholtz3D.doublelayer_transposed(gamma=γ) +0.5Identity()

uⁱ = Helmholtz3D.planewave(wavenumber=κ, direction=z)
f = strace(uⁱ,Γ)
g = ∂n(uⁱ)

@hilbertspace u
@hilbertspace v

eq1 = @discretise a[v,u] == f[v] u∈X v∈X
eq2 = @discretise b[v,u] == g[v] u∈X v∈X

x1 = solve(eq1)
x2 = solve(eq2)

fcr1, geo1 = facecurrents(x1, X)
fcr2, geo2 = facecurrents(x2, X)

include(Pkg.dir("CompScienceMeshes","examples","plotlyjs_patches.jl"))
p1 = patch(geo1, real.(norm.(fcr1)))
p2 = patch(geo1, real.(norm.(fcr2)))

## test the results
Z = assemble(a,X,X);
m, n = 1, numfunctions(X)
chm, chn = chart(Γ,cells(Γ)[m]), chart(Γ,cells(Γ)[n])
ctm, ctn = center(chm), center(chn)
R = norm(cartesian(ctm)-cartesian(ctn))
G = exp(-im*κ*R)/(4π*R)
Wmn = volume(chm) * volume(chn) * G
@show abs(Wmn-Z[m,n]) / abs(Z[m,n])
@test abs(Wmn-Z[m,n]) / abs(Z[m,n]) < 2.0e-3

r = assemble(f,X)
m = 1
chm = chart(Γ,cells(Γ)[m])
ctm = center(chm)
sm = volume(chm) * f(ctm)
r[m]
@show abs(sm - r[m]) / abs(r[m])
@test abs(sm - r[m]) / abs(r[m]) < 1e-3


text:
module DemoPackageTEH

export print_greeting

print_greeting() = println("Hi body! New message.")

end


text:
function add_dual_vars_in_dual_cones(dual_model::MOI.ModelLike, primal_model::MOI.ModelLike,
                                     primal_dual_map::PrimalDualMap{T}, dual_names::DualNames,
                                     con_types::Vector{Tuple{DataType, DataType}}) where T
    dual_obj_affine_terms = Dict{VI, T}()
    for (F, S) in con_types
        for ci in MOI.get(primal_model, MOI.ListOfConstraintIndices{F,S}()) # Constraints of type {F, S}
            # Add dual variable to dual cone
            # Fill a dual objective dictionary
            # Fill the primal_con_dual_var dictionary
            ci_dual = add_dual_variable(dual_model, primal_model, dual_names,
                                        primal_dual_map.primal_con_dual_var, dual_obj_affine_terms, ci)
            push_to_primal_con_dual_con!(primal_dual_map.primal_con_dual_con, ci, ci_dual)
            push_to_primal_con_constants!(primal_model, primal_dual_map.primal_con_constants, ci)
        end
    end
    return dual_obj_affine_terms
end

# Utils for primal_con_constants dict
function push_to_primal_con_constants!(primal_model::MOI.ModelLike, primal_con_constants::Dict{CI, Vector{T}},
                                      ci::CI{F, S}) where {T, F <: MOI.AbstractScalarFunction, S <: MOI.AbstractScalarSet}
    push!(primal_con_constants, ci => get_scalar_term(primal_model, ci))
    return
end

function push_to_primal_con_constants!(primal_model::MOI.ModelLike, primal_con_constants::Dict{CI, Vector{T}},
                                       ci::CI{F, S}) where {T, F <: MOI.AbstractVectorFunction, S <: MOI.AbstractVectorSet}
    return # No constants need to be passed to the DualOptimizer in this case, don't need to push zero to the dict
end

# Utils for primal_con_dual_con dict
function push_to_primal_con_dual_con!(primal_con_dual_con::Dict{CI, CI}, ci::CI, ci_dual::CI)
    push!(primal_con_dual_con, ci => ci_dual)
    return
end

function push_to_primal_con_dual_con!(primal_con_dual_con::Dict{CI, CI}, ci::CI, ci_dual::Nothing)
    return # Don't put in the dict a nothing value
end

# Utils for dual_obj_affine_terms dict
function push_to_dual_obj_aff_terms!(primal_model::MOI.ModelLike, dual_obj_affine_terms::Dict{VI, T}, vi::VI,
                                     ci::CI{F, S}, i::Int) where {T, F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    s = get_set(primal_model, ci)
    value = set_dot(i, s, T)*get_scalar_term(primal_model, i, ci)
    if !iszero(value) # If value is different than 0 add to the dictionary
        push!(dual_obj_affine_terms, vi => value)
    end
    return
end

function push_to_dual_obj_aff_terms!(primal_model::MOI.ModelLike, dual_obj_affine_terms::Dict{VI, T}, vi::VI,
                                     ci::CI{VVF, S}, i::Int) where {T, S <: MOI.AbstractSet}
    return # It is zero so don't push to the dual_obj_affine_terms
end

function _add_dual_variable(dual_model::MOI.ModelLike, primal_model::MOI.ModelLike, dual_names::DualNames,
                            primal_con_dual_var::Dict{CI, Vector{VI}}, dual_obj_affine_terms::Dict{VI, T},
                            ci::CI{F, S}) where {T, F <: MOI.AbstractFunction,
                                                    S <: MOI.AbstractSet}

    # Change to add_constrained_varibales in MOI 0.9.0
    # because of https://github.com/guilhermebodin/Dualization.jl/issues/9
    vis = MOI.add_variables(dual_model, get_ci_row_dimension(primal_model, ci)) # Add as many variables as the dimension of the constraint
    push!(primal_con_dual_var, ci => vis) # Add the map of the added dual variable to the relationated constraint
    # Get constraint name
    ci_name = MOI.get(primal_model, MOI.ConstraintName(), ci)
    # Add each vi to the dictionary
    for (i, vi) in enumerate(vis)
        push_to_dual_obj_aff_terms!(primal_model, dual_obj_affine_terms, vi, ci, i)
        set_dual_variable_name(dual_model, vi, i, ci_name,
                               dual_names.dual_variable_name_prefix)
    end
    return vis
end

function add_dual_variable(dual_model::MOI.ModelLike, primal_model::MOI.ModelLike, dual_names::DualNames,
                           primal_con_dual_var::Dict{CI, Vector{VI}}, dual_obj_affine_terms::Dict{VI, T},
                           ci::CI{F, S}) where {T, F <: MOI.AbstractFunction,
                                                   S <: MOI.AbstractSet}

    vis = _add_dual_variable(dual_model, primal_model, dual_names,
                             primal_con_dual_var, dual_obj_affine_terms, ci)
    return add_dual_cone_constraint(dual_model, primal_model, vis, ci)
end

function set_dual_variable_name(dual_model::MOI.ModelLike, vi::VI, i::Int, ci_name::String, prefix::String)
    MOI.set(dual_model, MOI.VariableName(), vi, prefix*ci_name*"_$i")
    return
end


text:
#=

A Pythagorean triplet is a set of three natural numbers, a < b < c, for which,
a2 + b2 = c2

For example, 3^2 + 4^2 = 9 + 16 = 25 = 5^2.

There exists exactly one Pythagorean triplet for which a + b + c = 1000.
Find the product abc.

=#

p = 12

function pythagorean_triplet(perimeter)
    for b in 2:Int(ceil((perimeter-1)/2))
        for a in 1:b-1
            if 2*perimeter*(a+b) == perimeter^2 + 2*a*b
                #println(a, ", ", b, ", ", perimeter-a-b)
                return a, b, perimeter-a-b
            end
            #println(a, ", ", b, ", ", perimeter-a-b)
        end
    end
end

res = pythagorean_triplet(1000)
prod(res)

max_stars_repo_name: pedvide/project-euler

max_issues_repo_name: pedvide/project-euler

max_forks_repo_name: pedvide/project-euler


text:
using Test, Neighborhood, StaticArrays, Random, Distances
using Neighborhood: datatype, getmetric, bruteforcesearch
using Neighborhood.Testing

Random.seed!(54525)
data = [rand(SVector{3}) for i in 1:1000]
query = SVector(0.99, 0.99, 0.99)

# Theiler window and skip predicate related
data[48:52] .= [SVector(0.0, 0.0, i*0.001) for i in 1:5]
nidxs = 48:2:52
queries = [data[i] for i in nidxs]
theiler1 = Theiler(1, nidxs)
theiler2 = Theiler(2, nidxs)

r = 0.1
k = 5

@testset "Utils" begin include("util.jl") end
@testset "Neighborhood.Testing" begin include("Testing.jl") end
@testset "Brute force" begin include("bruteforce.jl") end
include("nearestneighbors.jl")


text:
"""
    AzimuthalQuadrature{N,T<:Real}

Holds information of the azimuthal quadrature, such as the azimuthal angles `ϕs`, the
azimuthal spacings `δ` as well as effective azimuthal spacings `δs` and the azimuthal
weights `ωₐ`.
"""
struct AzimuthalQuadrature{N,T<:Real,N2,N4}
    δ::T
    δs::Vector{T}
    ϕs::Vector{T}
    ωₐ::Vector{T}
end

nazim(::AzimuthalQuadrature{N}) where {N} = N
nazim2(::AzimuthalQuadrature{N,T,N2}) where {N,T,N2} = N2
nazim4(::AzimuthalQuadrature{N,T,N2,N4}) where {N,T,N2,N4} = N4

AzimuthalQuadrature(n_azim::Int, δ::Real) = AzimuthalQuadrature(Val(n_azim), δ)

function AzimuthalQuadrature(::Val{N}, δ::T) where {N,T<:Real}
    N > 0 || throw(DomainError(N, "number of azimuthal angles must be positive."))
    iszero(rem(N, 4)) || throw(DomainError(N, "number of azimuthal angles must be a " *
                                              "multiple of 4."))
    δ > 0 || throw(DomainError(δ, "azimuthal spacing must be positive."))

    N2 = div(N, 2)
    N4 = div(N, 4)

    δs, ϕs, ωₐ = ntuple(_ -> Vector{T}(undef, N2), 3)

    return AzimuthalQuadrature{N,T,N2,N4}(δ, δs, ϕs, ωₐ)
end

function init_weights!(aq::AzimuthalQuadrature)
    @unpack ϕs, ωₐ = aq
    n_azim_4 = nazim4(aq)

    for i in 1:n_azim_4
        if isone(i)
            ωₐ[i] = ϕs[i+1] - ϕs[i]
        elseif isequal(i, n_azim_4)
            ωₐ[i] = π - ϕs[i] - ϕs[i-1]
        else
            ωₐ[i] = ϕs[i+1] - ϕs[i-1]
        end
        ωₐ[i] /= 4π

        j = suplementary_idx(aq, i)
        ωₐ[j] = ωₐ[i]
    end
    return nothing
end

# TODO: use north-west y cosas asi?
right_dir(::AzimuthalQuadrature{N,T,N2,N4}) where {N,T,N2,N4} = 1:N4
left_dir(::AzimuthalQuadrature{N,T,N2,N4}) where {N,T,N2,N4} = (N4+1):N2
both_dir(::AzimuthalQuadrature{N,T,N2}) where {N,T,N2} = 1:N2

points_right(::AzimuthalQuadrature{N,T,N2,N4}, i) where {N,T,N2,N4} = i ≤ N4
points_left(aq::AzimuthalQuadrature, i) = !points_right(aq, i)

suplementary_idx(::AzimuthalQuadrature{N,T,N2}, i) where {N,T,N2} = N2 - i + 1


text:
module SpecialIterators

import Base:start,next,done,iteratorsize,iteratoreltype,eltype,length,size

include("Exterate.jl")
include("Ziperate.jl")

export ziperate,exterate

end

max_stars_repo_name: sebastianpech/SpecialIterators

max_issues_repo_name: sebastianpech/SpecialIterators

max_forks_repo_name: sebastianpech/SpecialIterators


text:
function calc_doppler_component{T}(lambda::AbstractArray{T,1}, flux::AbstractArray{T,1};
                                sigmasq_obs::AbstractArray{T,1} = 1e-16*ones(length(lambda)),	sigmasq_cor::T = 1.0, rho::T = 1.0,
                                half_chunck_size::Integer = 100)
   lambda.*calc_gp_on_segments(predict_deriv,lambda, flux, sigmasq_obs=sigmasq_obs,	sigmasq_cor=sigmasq_cor, rho=rho, half_chunck_size=half_chunck_size)
end

function calc_gp_on_segments{T}(predict_gp::Function, lambda::AbstractArray{T,1}, flux::AbstractArray{T,1};
                                sigmasq_obs::AbstractArray{T,1} = 1e-16*ones(length(lambda)),	sigmasq_cor::T = 1.0, rho::T = 1.0,
                                half_chunck_size::Integer = 100)
  @assert length(lambda) == length(flux) == length(sigmasq_obs)
  output = Array(Float64,length(lambda))
  num_seg = convert(Int64,ceil(length(lambda)/half_chunck_size)-1)
  for i in 1:num_seg
    idx_begin = 1+half_chunck_size*(i-1)
    idx_end = min(half_chunck_size*(i+1), length(lambda))
    write_idx_begin = idx_begin + div(half_chunck_size,2)
    write_idx_end = idx_end - div(half_chunck_size,2)
    if i==1 write_idx_begin=1 end
    if i==num_seg
      idx_begin = max(1,idx_end-2*half_chunck_size)
      write_idx_end=length(lambda)
    end
    #println("# i= ",i,": ", idx_begin, " - ", idx_end, " -> ", write_idx_begin, " - ", write_idx_end)
    #output[write_idx_begin:write_idx_end] = predict_gp(view(lambda,idx_begin:idx_end), view(flux,idx_begin:idx_end), view(lambda,write_idx_begin:write_idx_end), sigmasq_obs=view(sigmasq_obs,idx_begin:idx_end), sigmasq_cor=sigmasq_cor, rho=rho)
    output[write_idx_begin:write_idx_end] = predict_gp(lambda[idx_begin:idx_end], flux[idx_begin:idx_end], lambda[write_idx_begin:write_idx_end], sigmasq_obs=sigmasq_obs[idx_begin:idx_end], sigmasq_cor=sigmasq_cor, rho=rho)
  end
  output
end

function calc_doppler_component{T}(lambda::AbstractArray{T,1}, flux::AbstractArray{T,2};
                                sigmasq_obs::AbstractArray{T,1} = 1e-16*ones(length(xobs)),	sigmasq_cor::T = 1.0, rho::T = 1.0)
  doppler_basis = calc_doppler_component(lambda,vec(mean(flux,2)),sigmasq_obs=sigmasq_obs,sigmasq_cor=sigmasq_cor,rho=rho)
end

function calc_doppler_quadratic_term{T}(lambda::AbstractArray{T,1}, flux::AbstractArray{T,1};
                                sigmasq_obs::AbstractArray{T,1} = 1e-16*ones(length(lambda)),	sigmasq_cor::T = 1.0, rho::T = 1.0,
                                half_chunck_size::Integer = 100)
   0.5*lambda.^2.*calc_gp_on_segments(predict_deriv2,lambda, flux, sigmasq_obs=sigmasq_obs,	sigmasq_cor=sigmasq_cor, rho=rho, half_chunck_size=half_chunck_size)
end

function calc_doppler_quadratic_term{T}(lambda::AbstractArray{T,1}, flux::AbstractArray{T,2};
                                     sigmasq_obs::AbstractArray{T,1} = 1e-16*ones(length(xobs)),	sigmasq_cor::T = 1.0, rho::T = 1.0)
  doppler_quad_term = calc_doppler_quadratic_term(lambda,vec(mean(flux,2)),sigmasq_obs=sigmasq_obs,sigmasq_cor=sigmasq_cor,rho=rho)
end

doppler_comp = calc_doppler_component(lambda,obs, sigmasq_obs=1/150000.0*ones(length(lambda)), rho=0.1)
genpca_out = fit_gen_pca_eford(obs,doppler_comp)

doppler_quadratic_term = calc_doppler_quadratic_term(lambda,obs, sigmasq_obs=1/150000.0*ones(length(lambda)), rho=0.1)
gen2pca_out = fit_gen2_pca_eford(obs,doppler_comp,doppler_quadratic_term)

max_stars_repo_name: eford/RVSpectraKitLearn

max_issues_repo_name: eford/RVSpectraKitLearn

max_forks_repo_name: eford/RVSpectraKitLearn


text:
mutable struct BufferedIRCurve{M, C<:AbstractIRCurve{M}} <: AbstractIRCurve{M}
    dates_buffer::Dict{Date, Float64}
    yf_buffer::Dict{Rational, Float64}
    curve::C
    lock::Base.Threads.SpinLock
end

BufferedIRCurve(curve::AbstractIRCurve) = BufferedIRCurve(Dict{Date, Float64}(), Dict{Rational, Float64}(), curve, Base.Threads.SpinLock())

for fun in (:curve_get_name, :curve_get_daycount, :curve_get_compounding, :curve_get_method, :curve_get_date, :curve_get_dtm, :curve_get_zero_rates, :curve_get_model_parameters, :curve_get_spline_fit_on_rates, :curve_get_spline_fit_on_discount_factors)
    @eval begin
        ($fun)(curve::BufferedIRCurve) = ($fun)(curve.curve)
    end
end

curve_get_dict_parameter(curve::BufferedIRCurve, sym::Symbol) = curve_get_dict_parameter(curve.curve, sym)

function zero_rate(curve::BufferedIRCurve, maturity::Date)
    lock(curve.lock)
    try
        if !haskey(curve.dates_buffer, maturity)
            curve.dates_buffer[maturity] = zero_rate(curve.curve, maturity)
        end
        return curve.dates_buffer[maturity]
    finally
        unlock(curve.lock)
    end
end

function zero_rate(curve::BufferedIRCurve, maturity::YearFraction{T}) where {T<:Union{Integer, Rational}}
    yf_value = value(maturity)

    lock(curve.lock)
    try
        if !haskey(curve.yf_buffer, yf_value)
            curve.yf_buffer[yf_value] = zero_rate(curve.curve, maturity)
        end
        return curve.yf_buffer[yf_value]
    finally
        unlock(curve.lock)
    end
end

function zero_rate(curve::BufferedIRCurve, maturity_vec::Vector{T}) where {T<:Union{Date, YearFraction}}
    n = length(maturity_vec)
    result = Vector{Float64}(undef, n)
    for i in 1:n
        @inbounds result[i] = zero_rate(curve, maturity_vec[i])
    end
    return result
end


text:
using JuMP

"""
    hrep(model::JuMP.Model)

Builds an H-representation from the feasibility set of the JuMP model `model`.
Note that if non-linear constraint are present in the model, they are ignored.
"""
hrep(model::JuMP.Model) = LPHRep(model)
LPHRep(model::JuMP.Model) = LPHRep(backend(model))
polyhedron(model::JuMP.Model, args...) = polyhedron(hrep(model), args...)


text:
using FastTransforms, Compat
using Compat.Test, Compat.LinearAlgebra
import Compat.Statistics

@testset "Chebyshev--Jacobi transform" begin
    println("Testing the accuracy of asymptotics")

    Nr = 5
    v = zeros(Nr)
    Na,Nb = 5,5
    V = zeros(Na,Nb)

    for N in round.([Int],10 .^ range(1,stop=3,length=3))
        for αi=1:Na, βi=1:Nb
            α,β = -0.45+(αi-1)/Na,-0.45+(βi-1)/Nb
            p1,p2 = plan_cjt(rand(N),α,β),plan_icjt(rand(N),α,β)
            for i=1:Nr
                c = rand(N)
                v[i] = log(norm(p2*(p1*c)-c,Inf)/(20N^(1+2max(α,β))*eps()))
            end
            V[αi,βi] = Statistics.mean(v)
        end
        @test Statistics.mean(V) < 2
    end

    println("Testing the special cases length(c) = 0,1,2:")

    c = [1.5]
    @test cjt(c,0.12,0.34) == c
    @test icjt(c,0.12,0.34) == c

    pop!(c)

    @test cjt(c,0.12,0.34) == c
    @test icjt(c,0.12,0.34) == c

    @test cjt(c,0.0,0.0) == c
    @test icjt(c,0.0,0.0) == c

    c = [1.0;2.0]
    @test norm(jjt(c,0.12,0.34,0.12,0.34)-c,Inf) ≤ 2eps()

    # test same parameters

    c=[0.,1.]
    @test jjt(c,1.5,-0.5,1.5,-0.5) == c

    println("Testing the special cases (α,β) = (±0.5,±0.5)")

    # Chebyshev coefficients of exp
    c_cheb = [1.2660658777520084,1.13031820798497,0.27149533953407656,0.044336849848663804,0.005474240442093732,0.0005429263119139438,4.497732295429515e-5,3.1984364624019905e-6,1.9921248066727958e-7,1.1036771725517344e-8,5.505896079673747e-10,2.4979566169849825e-11,1.03915223067857e-12,3.9912633564144015e-14,1.4237580108256572e-15,4.740926102561496e-17,1.4801800572082976e-18,4.3499194949441696e-20,1.2074289272797528e-21,3.175356737059445e-23,7.93367197163804e-25,1.8879484042289162e-26,4.2886738765925875e-28,9.318985281777577e-30,1.9406469749017399e-31,3.8798022492260015e-33,7.45850288739152e-35,1.3807477824110646e-36,2.4648623717710964e-38,4.2485421925059915e-40,7.079001176212895e-42,1.141486778254089e-43,1.783151037543284e-45,2.7011422638738024e-47,3.9714338657415316e-49,5.672351695612112e-51,7.876788130533199e-53,1.0642416282087358e-54,1.4000817885528905e-56,1.794689085390544e-58,2.243019482662972e-60,2.7349926016641967e-62,3.2554929833033215e-64,3.78495689402033e-66,4.300544446926769e-68,4.77780574839433e-70,5.192665133046628e-72,5.523442810988436e-74,5.717985746408847e-76]

    # Jacobi coefficients of exp at (±0.5,±0.5)
    c_00 = [1.2660658777520084,2.26063641596994,0.723987572090871,0.14187791951572415,0.020020079331085648,0.0022061767595233264,0.00019937999439479753,1.5268992715616027e-5,1.014420290055232e-6,5.9506859054175615e-8,3.124851418974203e-9,1.4852160704995775e-10,6.447143371527455e-12,2.5753240175931934e-13,9.526907336854026e-15,3.281725133829881e-16,1.0576496993856426e-17,3.20238467156949e-19,9.14299140846687e-21,2.4694550834939014e-22,6.3281713790321685e-24,1.542622093731937e-25,3.5857224408995815e-27,7.964666406806272e-29,1.6939043515668796e-30,3.4556190450066486e-32,6.773312266576111e-34,1.2775611026088504e-35,2.3221236341654026e-37,4.072731101981327e-39,6.901079409575277e-41,1.1310395206213411e-42,1.794875997727078e-44,2.760732681473006e-46,4.119631242417262e-48,5.969295899580336e-50,8.405881613439669e-52,1.1512859570474932e-53,1.5347890791678689e-55,1.9929132418968355e-57,2.522290550037395e-59,3.1134869971034086e-61,3.7506700259829663e-63,4.411970068315684e-65,5.0705898146983576e-67,5.696604001638848e-69,6.259278760298182e-71,6.729592140384465e-73,7.039950057933299e-75]
    c_01 = [1.8312249817444934,1.4018135475190467,0.4211095858436539,0.07969774446521205,0.011002819207328322,0.0011944708771926122,0.00010677917238679192,8.110005868771403e-6,5.353106062476433e-7,3.1237733951100556e-8,1.6333110219427631e-9,7.735005972383578e-11,3.3473853404557083e-12,1.333595312027857e-13,4.92207038322879e-15,1.6920924742289326e-16,5.443658341283791e-18,1.6456374329092366e-19,4.691719175403533e-21,1.2655773772197326e-22,3.2393803393291905e-24,7.888321906112733e-26,1.8318188278743866e-27,4.065263937281931e-29,8.638847091039724e-31,1.7610248038105726e-32,3.449351261471638e-34,6.50183837007916e-36,1.1810743750838164e-37,2.070295858087742e-39,3.506179552173043e-41,5.743539156119835e-43,9.110325158556408e-45,1.4006615828866467e-46,2.0892357224279913e-48,3.026093616078655e-50,4.2597272086552854e-52,5.832159509538512e-54,7.772313549445791e-56,1.0089104305392274e-57,1.2765228632360985e-59,1.5752735945134057e-61,1.8971383534453687e-63,2.2310498815373167e-65,2.563461449357282e-67,2.8792582164274203e-69,3.1629293838222696e-71,3.399629156416382e-73,3.5199750289666496e-75]
    c_10 = [0.7009067737595234,0.8588228684508934,0.30287798624721707,0.062180175050512106,0.009017260123757326,0.0010117058823307142,9.260082200800561e-5,7.1589868468446235e-6,4.791096838075886e-7,2.8269125103075062e-8,1.49154039703144e-9,7.117154732612197e-11,3.099758031071747e-12,1.2417287055653363e-13,4.604836953625235e-15,1.5896326596009483e-16,5.1328386525726346e-18,1.5567472386602532e-19,4.451272233063337e-21,1.2038777062741688e-22,3.088791039702978e-24,7.537899031206638e-26,1.753903613025195e-27,3.899402469524341e-29,8.300196424629072e-31,1.694594241196076e-32,3.3239610051044727e-34,6.273772656009344e-36,1.1410492590815863e-37,2.002435243893585e-39,3.3948998574022335e-41,5.566856050093576e-43,8.838434818714371e-45,1.3600710985863591e-46,2.0303955199892708e-48,2.943202283501681e-50,4.146154404784383e-52,5.68070006093642e-54,7.575577242232898e-56,9.84002811357608e-58,1.2457676868012966e-59,1.5382134025900027e-61,1.8535316725375976e-63,2.1809201867783675e-65,2.5071283653410757e-67,2.8173457852114277e-69,3.0963493764759127e-71,3.3299629839680836e-73,3.5199750289666496e-75]
    c_11 = [1.13031820798497,0.7239875720908708,0.21281687927358628,0.040040158662171296,0.005515441898808317,0.0005981399831843923,5.34414745046561e-5,4.057681160220927e-6,2.677808657437903e-7,1.5624257094871016e-8,8.168688387747676e-10,3.8682860229164727e-11,1.6739606114355753e-12,6.668835135797819e-14,2.4612938503724105e-15,8.46119759508514e-17,2.722026970834067e-18,8.228692267620182e-20,2.3459823293192055e-21,6.328171379032168e-23,1.619753198418534e-24,3.944294684989539e-26,9.159366367827211e-28,2.0326852218802552e-29,4.31952380625831e-31,8.805305946548944e-33,1.7247074885219479e-34,3.2509730878315635e-36,5.905460097872925e-38,1.0351619114362916e-39,1.7531112569630786e-41,2.8718015963633246e-43,4.5552089244304594e-45,7.003373112109346e-47,1.044626782426559e-48,1.5130586904191403e-50,2.129879020537863e-52,2.916099250418951e-54,3.8861808216988285e-56,5.04458110007479e-58,6.382648344061989e-60,7.876407054564243e-62,9.485735646882059e-64,1.1155297592778338e-65,1.2817359024131178e-67,1.4396345613848593e-69,1.5814715428162693e-71,1.700107488097128e-73,1.7781317156635653e-75]

    @test norm(cjt(c_00,-0.5,-0.5)-c_cheb,Inf) < eps()
    @test norm(cjt(c_01,-0.5,0.5)-c_cheb,Inf) < eps()
    @test norm(cjt(c_10,0.5,-0.5)-c_cheb,Inf) < eps()
    @test norm(cjt(c_11,0.5,0.5)-c_cheb,Inf) < eps()

    @test norm(icjt(c_cheb,-0.5,-0.5)-c_00,Inf) < eps()
    @test norm(icjt(c_cheb,-0.5,0.5)-c_01,Inf) < eps()
    @test norm(icjt(c_cheb,0.5,-0.5)-c_10,Inf) < eps()
    @test norm(icjt(c_cheb,0.5,0.5)-c_11,Inf) < eps()

    c = exp.(collect(1:1000)./(-30))

    println("Testing increment/decrement operators for α,β ≤ -0.5")

    α,β = -0.9,-0.6

    @test norm(FastTransforms.fromsquare!(FastTransforms.tosquare!(copy(c),α,β),α,β)-c,Inf) < 10eps()
    @test norm(icjt(cjt(c,α,β),α,β)-c,Inf) < 20000eps()
    p1,p2 = plan_cjt(c,α,β),plan_icjt(c,α,β)
    @test norm(p2*(p1*c)-c,Inf) < 20000eps()

    println("Testing increment/decrement operators for α ≤ -0.5, β > -0.5")

    α,β = -0.9,3.1

    @test norm(FastTransforms.fromsquare!(FastTransforms.tosquare!(copy(c),α,β),α,β)-c,Inf) < 10eps()
    @test norm(icjt(cjt(c,α,β),α,β)-c,Inf) < 2000eps()
    p1,p2 = plan_cjt(c,α,β),plan_icjt(c,α,β)
    @test norm(p2*(p1*c)-c,Inf) < 2000eps()

    println("Testing increment/decrement operators for α > -0.5, β ≤ -0.5")

    α,β = 2.4,-0.75

    @test norm(FastTransforms.fromsquare!(FastTransforms.tosquare!(copy(c),α,β),α,β)-c,Inf) < 1200eps()
    @test norm(icjt(cjt(c,α,β),α,β)-c,Inf) < 1e6eps()
    p1,p2 = plan_cjt(c,α,β),plan_icjt(c,α,β)
    @test norm(p2*(p1*c)-c,Inf) < 1e6eps()

    println("Testing increment/decrement operators for α > -0.5, β > -0.5")

    α,β = 2.4,3.6

    @test norm(FastTransforms.fromsquare!(FastTransforms.tosquare!(copy(c),α,β),α,β)-c,Inf) < 300eps()
    @test norm(icjt(cjt(c,α,β),α,β)-c,Inf) < 2e5eps()
    p1,p2 = plan_cjt(c,α,β),plan_icjt(c,α,β)
    @test norm(p2*(p1*c)-c,Inf) < 2e5eps()

    println("Testing for complex coefficients")

    α,β = 0.12,0.34
    c = complex.(rand(100),rand(100))

    @test cjt(c,α,β) == complex.(cjt(real(c),α,β),cjt(imag(c),α,β))
    @test icjt(c,α,β) == complex.(icjt(real(c),α,β),icjt(imag(c),α,β))
    @test jjt(c,α,β,α,β) == complex.(jjt(real(c),α,β,α,β),jjt(imag(c),α,β,α,β))
    @test norm(jjt(c,α,β,α,β)-c,Inf) < 200eps()

    println("Testing for Vector{Float32}")

    c64 = rand(100)
    c32 = map(Float32,c64)

    cL64 = cjt(c64,0.,0.)
    cL32 = cjt(c32,0.f0,0.f0)

    @test norm(cL32-cL64,Inf) < 40eps(Float32)
    @test cL32 == cjt(c32,0.,0.)

    println("Testing for Matrix of coefficients")

    c = rand(100,100)
    @test maximum(abs,jjt(c,α,β,α,β)-c) < 10000eps()

    @testset "cjt bug (#57)" begin
        @test cjt([1.,2],.5,.5) == [1.,3]
        @test FastTransforms.incrementαβ!([1.0], -0.5, -0.5) == [1.0]
        @test FastTransforms.incrementαβ!([1.0, 2.0], -0.5, -0.5) == [1.0,2/3]
        @test FastTransforms.decrementαβ!([1.0], 0.5, 0.5) == [1.0]
        @test FastTransforms.decrementαβ!([1.0, 2.0], 0.5, 0.5) == [1.0,6]
    end
end


text:
using Base.Test
using ChaosQuantified

include("helpers.jl")
include("Rqa.jl")
include("entropy.jl")


text:
export StackOfStarsTrajectory, stackOfStarsNodes, stackOfStarsDensity

"""
    StackOfStarsTrajectory(numProfiles, numSamplingPerProfile
                  ; TE::Float64=0.0
                  , AQ::Float64=1.e-3
                  , numSlices=1
                  , angleOffset= :equispaced
                  , kargs...)

returns a 2d radial trajectory.

# Arguments
* `numProfiles::Int64`            - number of profiles
* `numSamplingPerProfile::Int64`  - number of sampling points per profile
* (`TE::Float64=0.0`)             - echo time in s
* (`AQ::Float64=1.e-3`)           - readout duration in s (per profile)
* (`numSlices=1`)                 - number of slices
* (`angleOffset= :equispaced`)    - spacing of profile angles (`:equispaced` sampling, `:golden` angle sampling or `:random` sampling)
"""
function StackOfStarsTrajectory(::Type{T}, numProfiles, numSamplingPerProfile
                  ; TE=0.0
                  , AQ=1.e-3
                  , numSlices=1
                  , angleOffset = :equispaced
                  , kargs...) where T
  nodes = cartesian3dNodes(T, numProfiles, numSamplingPerProfile; numSlices=numSlices, angleOffset = :equispaced)
  times = readoutTimes(T, numProfiles, numSamplingPerProfile, numSlices; TE=TE, AQ=AQ)
  return  Trajectory("StackOfStars", nodes, times, TE, AQ, numProfiles, numSamplingPerProfile, numSlices, false, true)
end

function stackOfStarsNodes(::Type{T}, numProfiles, numSamplingPerProfile
              ; numSlices=1, angleOffset = :equispaced, kargs...) where T

  angles = collect(pi.*(0:numProfiles-1)/numProfiles)
  nodes = zeros(T, 3, numSamplingPerProfile, numProfiles, numSlices)
  if angleOffset == :golden
    angles = [i*getGoldenAngleRad()  for i=0:numProfiles-1 ]
  elseif angleOffset == :random
    angles = sort(pi .* rand(numProfiles))
  elseif angleOffset == :equispaced
    angles = collect(pi.*(0:numProfiles-1)/numProfiles)
  end

  pos = collect((0:numSamplingPerProfile-1)/numSamplingPerProfile .- 0.5)
  posZ = collect((0:numSlices-1)/numSlices .- 0.5)
  for j = 1:numSlices
    for l = 1:numProfiles
      for k = 1:numSamplingPerProfile
        nodes[1,k,l,j] = (-1)^l * pos[k]*cos(angles[l])
        nodes[2,k,l,j] = (-1)^l * pos[k]*sin(angles[l])
        nodes[3,k,l,j] = posZ[j]
      end
    end
  end
  return reshape(nodes, 3, numSamplingPerProfile*numProfiles*numSlices)
end

function stackOfStarsDensity(::Type{T}, numSamplingPerProfile::Int64, numProfiles::Int64, numSlices::Int64) where T
  density = zeros(T, numSamplingPerProfile, numProfiles, numSlices)
  pos = collect((0:numSamplingPerProfile-1)/numSamplingPerProfile .- 0.5)
  for j= 1:numSlices
    for l = 1:numProfiles
      for k = 1:numSamplingPerProfile
        density[k,l,j] = abs(pos[k])
      end
    end
  end
  return vec(density)
end


text:
@typeparams struct Permutation <: AbstractVector{Int}
    p::{}   # Forward permutation
    ip::{}  # Backward permutation
end

Base.size(p::Permutation) = size(p.p)
Base.getindex(p::Permutation, i) = p.p[i]
Base.@propagate_inbounds function Base.setindex!(p::Permutation, pi, i)
    @unpack p,ip = p
    ip[pi] = i
    return setindex!(p,pi,i)
end

@typeparams struct PermutationTracker{K,V} <: AbstractVector{Tuple{K,V}}
    v::{<:AbstractVector{V}}     # Data
    k2i::{}                      # Key -> current index in d
    i2k::{<:AbstractVector{K}} # Current index in d -> key
end
PermutationTracker(v,k2i,i2k) = PermutationTracker{keytype(k2i), eltype(v)}(v,k2i,i2k)

Base.size(t::PermutationTracker) = size(t.v)
Base.getindex(t::PermutationTracker, i::Int) = (t.i2k[i],t.v[i])
function Base.setindex!(t::PermutationTracker{K,V}, (k,vk)::Tuple{K,V}, i::Int) where {K,V}
    @unpack v,k2i,i2k = t
    v[i] = vk
    Permutation(i2k,k2i)[i] = k
    return t
end


text:
using Test
using ComparativeAutograder

include("./stdio.jl")
include("./compare.jl")
include("./testcase.jl")
include("./testsuite.jl")
include("./runner.jl")
include("./proc.jl")
# include("./exercises/runtests.jl")


text:
# Julia wrapper for header: /opt/local/include/gdal_csv.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function GDALDefaultCSVFilename(pszBasename::Ptr{Uint8})
    ccall((:GDALDefaultCSVFilename,libgdal),Ptr{Uint8},(Ptr{Uint8},),pszBasename)
end


text:
function score(str::AbstractString)

end


text:
@info "Testing Conversions..."
@testset "Conversions Utilities" begin
    _voxels = zeros(Float32, 32, 32, 32, 2) |> gpu
    _voxels[1:15, 2:10, 18:32, :] .= 1

    _v = VoxelGrid(_voxels) |> gpu
    _m =
        load_trimesh(
            joinpath(@__DIR__, "../assets/teapot.obj"),
            joinpath(@__DIR__, "../assets/sphere.obj"),
        ) |> gpu
    _p = PointCloud(sample_points(_m, 1024)) |> gpu

    res = 28
    points = 512
    thresh = 0.9

    for algo in [:Exact, :MarchingCubes, :MarchingTetrahedra]
        m1 = TriMesh(_p, res; algo = algo)
        @test m1 isa TriMesh{Float32,UInt32,CuArray}

        m2 = TriMesh(_v, thresh = thresh, algo = algo)
        @test m2 isa TriMesh{Float32,UInt32,CuArray}

        p1 = PointCloud(_v, points, thresh = thresh, algo = algo)
        @test p1 isa PointCloud{Float32}
        @test p1.points isa CUDA.CuArray
        @test size(p1.points) == (3, points, 2)
    end

    p2 = PointCloud(_m, points)
    @test p2 isa PointCloud{Float32}
    @test size(p2.points) == (3, points, 2)
    @test p2.points isa CUDA.CuArray

    v1 = VoxelGrid(_m, res)
    @test v1 isa VoxelGrid{Float32}
    @test v1.voxels isa CUDA.CuArray
    @test size(v1.voxels) == (res, res, res, 2)

    v2 = VoxelGrid(_p, res)
    @test v2 isa VoxelGrid{Float32}
    @test v2.voxels isa CUDA.CuArray
    @test size(v2.voxels) == (res, res, res, 2)
end


text:
@reexport module Blaze
#using Reexport

import ..AbstractEchelleObservation, ..AbstractEchelleObservationSet

abstract type AbstractBlazeType end
abstract type AbstractBlaze1DType <: AbstractBlazeType end
abstract type AbstractBlaze2DType <: AbstractBlazeType end
export AbstractBlazeType, AbstractBlaze1DType, AbstractBlaze2DType

# Default functions for AbstractBlazeTypes
min_log_lambda(blaze::AbstractBlazeType) = blaze.log_lambda_lo
max_log_lambda(blaze::AbstractBlazeType) = blaze.log_lambda_hi
order(blaze::AbstractBlazeType) = blaze.order
pixel_range(blaze::AbstractBlazeType) = blaze.pixel_range
export order, min_log_lambda, max_log_lambda, pixel_range

# Constants from calling module
import ..default_blaze_degree# , ..max_blaze_degree
import ..max_num_pixels_per_order
import ..default_pixel_buffer_lo, ..default_pixel_buffer_hi

# Util functions on data from EchelleTools
import ..num_pixels_per_order, ..num_orders, ..num_obs

include("common.jl")
include("blaze_map.jl")
include("blaze_polynomial.jl")

end # module Blaze


text:
struct Pending end
struct Rescinded end
struct Waiting end

const WaiterFlags = Union{Pending,Rescinded,Waiting}
const WaiterState{T} = Union{Pending,Rescinded,Waiting,T}

struct Waiter{T,State<:Reagents.Ref{WaiterState{T}}}
    state::State
    task::Base.Task
end

Waiter{T}() where {T} =
    Waiter(Reagents.Ref{WaiterState{T}}(Pending()), current_task())::Waiter{T}

backedge!(::Waiter, ::Bag) = false

const CatalystFlags = Union{Pending,Rescinded}
const CatalystState{T} = Union{Pending,Rescinded,T}

struct Deleted end
deleteref() = Reagents.Ref{Union{ImmutableList{Bag},Deleted}}(nothing)

struct Catalyst{T,State<:Reagents.Ref{CatalystState{T}}}
    deleteref::typeof(deleteref())
    state::State
end

Catalyst{T}(deleteref = deleteref()) where {T} =
    Catalyst(deleteref, Reagents.Ref{CatalystState{T}}(Pending()))::Catalyst{T}

function backedge!(offer::Catalyst, @nospecialize(bag::Bag))
    old = offer.deleteref[]
    while true
        old isa Deleted && return true
        old, success = cas!(offer.deleteref, old, ImmutableListNode{Bag}(bag, old))
        success && break
    end
    _cleanupbags!(old)  # cleanup previous offer
    return false
end

function _delete!(offer::Catalyst)
    offer.state[] = Rescinded()
    old = offer.deleteref[]
    while true
        old isa Deleted && return
        old, success = cas!(offer.deleteref, old, Deleted())
        success && break
    end
    _cleanupbags!(old)
end

function _cleanupbags!(@nospecialize(list::ImmutableList{Bag}))
    list === nothing && return
    for bag in list
        for _ in bag  # TODO: better cleanup
        end
    end
end

const Offer{T} = Union{Waiter{T},Catalyst{T}}

offerid(offer::Waiter) = objectid(offer.state)
offerid(::Catalyst) = UInt(1)
offerid(::Nothing) = UInt(0)

"""
    Reaction

* `caslist::ImmutableList{CAS}`
* `offers::ImmutableList{Offer}`
* `postcommithooks::ImmutableList{Any}`

* `restart_on_failure::Bool` (default: `false`): A flag controlling if the
  offer should be invalidated on the failure.  Consider `(r | s) ⨟ t` where `s`
  and `t` both contain a `Swap`. When a `Swap` in `s` is triggered from the
  dual, it may indicate that `r` is now reactable (i.e., `s` is used as a
  condition variable).  In this case, it is now required to register the offer
  at `t`.  To support this behavior, `restart_on_failure` is set to `true` in
  the continuations of `Choice` (`|`). Then, in `tryreact!`, the offers of the
  dual messsage in `Swap` with `restart_on_failure` flag set are aborted.  In
  particular, if the state is in `Waiting`, the task is re-scheduled so that
  other branches of `|` can be re-evaluated.
"""
struct Reaction
    caslist::ImmutableList{CAS}
    offers::ImmutableList{Offer}
    postcommithooks::ImmutableList{Any}
    restart_on_failure::Bool
end

Reaction() = Reaction(nothing, nothing, nothing, false)

hascas(rx::Reaction) = rx.caslist !== nothing

uintptr(x) = UInt(pointer_from_objref(x.ref))

withcas(rx::Reaction, cas::CAS) = @set rx.caslist = pushsortedby(uintptr, rx.caslist, cas)
withoffer(rx::Reaction, offer::Offer) = @set rx.offers = pushfirst(rx.offers, offer)
withpostcommit(rx::Reaction, @nospecialize(f)) =
    @set rx.postcommithooks = pushfirst(rx.postcommithooks, f)

combine(rx1::Reaction, rx2::Reaction) = Reaction(
    combinesortedby(uintptr, rx1.caslist, rx2.caslist),
    combine(rx1.offers, rx2.offers),
    combine(rx1.postcommithooks, rx2.postcommithooks),
    rx1.restart_on_failure | rx2.restart_on_failure,
)

setcasing!(::Nothing) = true
function setcasing!(list::ImmutableListNode{CAS})
    (; ref, old) = list.value
    (_, success) = @atomicreplace ref.value old => CASing()
    success || return false
    if setcasing!(list.tail)
        return true
    else
        @atomic ref.value = old
        return false
    end
end

setall!(::Nothing) = true
function setall!(list::ImmutableListNode{CAS})
    (; ref, new) = list.value
    @atomic ref.value = new
    return setall!(list.tail)
end

commit!(rx::Reaction) = setcasing!(rx.caslist) && setall!(rx.caslist)

has_stale_cas(rx::Reaction) = has_stale_cas(rx.caslist)
has_stale_cas(::Nothing) = false
function has_stale_cas(list::ImmutableListNode{CAS})
    (; ref, old) = list.value
    value = @atomic ref.value
    return value !== old
end

struct Commit end

hascas(::Commit) = false

const Reactable = Union{Reactor,Commit}

maysync(actr::Reactor) = maysync(actr.reagent) || maysync(actr.continuation)
maysync(::Commit) = false

maysync(::Reagent) = false


text:
# module Ecto

module Schema

export primary_key, schema, field
import Base: getindex, ==

include("assoc.jl")

type Metadata
    state::Symbol # :built :loaded :deleted
    prefix::String # ""
    source::String # "source"
    context::Void
end

immutable t
    modul::Module
    struct::Assoc
    metadata::Metadata
    function t(modul::Module, struct::Assoc)
        metadata = Metadata(:built, "", "source", nothing)
        new(modul, struct, metadata)
    end
end
type PutMeta
    opts::Dict
end

type PrimaryKey
    name::Symbol
    typ::Symbol
    opts::Dict{Symbol,Any}
end

function ==(lhs::Schema.t, rhs::Schema.t)
    lhs.struct == lhs.struct
end

function getindex(schema::Schema.t, key::Symbol)
    if :fields == key
        get_attribute(schema.modul, :ecto_fields)
    elseif :associations == key
        get_attribute(schema.modul, :ecto_assocs)
    elseif :read_after_writes == key
        get_attribute(schema.modul, :ecto_fields)
    elseif :primary_key == key
        get_attribute(schema.modul, :ecto_primary_keys)
    elseif :types == key
        get_attribute(schema.modul, :ecto_fields)
    end
# `__schema__(:source)` - Returns the source as given to `schema/2`;
# `__schema__(:prefix)` - Returns optional prefix for source provided by
# `@schema_prefix` schema attribute;
# `__schema__(:primary_key)` - Returns a list of primary key fields (empty if there is none);
#
# `__schema__(:fields)` - Returns a list of all non-virtual field names;
# `__schema__(:type, field)` - Returns the type of the given non-virtual field;
# `__schema__(:types)` - Returns a map of all non-virtual
# field names and their type;
#
# `__schema__(:associations)` - Returns a list of all association field names;
# `__schema__(:association, assoc)` - Returns the association reflection of the given assoc;
#
# `__schema__(:embeds)` - Returns a list of all embedded field names;
# `__schema__(:embed, embed)` - Returns the embedding reflection of the given embed;
#
# `__schema__(:read_after_writes)` - Non-virtual fields that must be read back
# from the database after every write (insert or update);
#
# `__schema__(:autogenerate_id)` - Primary key that is auto generated on insert;
end

function primary_key(pk::Bool)
    __MODULE__ = current_module()
    __init_attributes__(__MODULE__)
    __attributes__[__MODULE__][:primary_key] = pk
end

function schema(block::Function, source::String)
    schema(block, source, true, :id)
end

function schema(block::Function, source::String, meta::Bool, typ::Symbol)
    __MODULE__ = current_module()
    __init_attributes__(__MODULE__)
    register_attribute(__MODULE__, :changeset_fields, accumulate= true)
    register_attribute(__MODULE__, :struct_fields, accumulate= true)
    prefix = get_attribute(__MODULE__, :schema_prefix)
    if meta
        put_struct_field(__MODULE__, :__meta__, Metadata(:built, prefix, source, nothing))
    end
    primary_key = get_attribute(__MODULE__, :primary_key)
    if primary_key == nothing
        primary_key = PrimaryKey(:id, typ, Dict(:autogenerate=>true))
    end
    if false == primary_key
        primary_key_fields = []
    else
        (name, typ, opts) = (primary_key.name, primary_key.typ, primary_key.opts)
        Schema.__field__(__MODULE__, name, typ, merge(Dict(:primary_key=>true), opts))
        primary_key_fields = [name]
    end
    block()
end

## API
function field(name::Symbol, typ::Symbol; kw...)
    __MODULE__ = current_module()
    opts = Dict(kw)
    Schema.__field__(__MODULE__, name, typ, opts)
end

function timestamps
end

function has_many
end

function has_one
end

function belongs_to
end

function many_to_many
end

## Callbacks
function __field__(modul::Module, name::Symbol, typ::Symbol, opts::Dict)
    pk = haskey(opts, :primary_key) ? opts[:primary_key] : false
    default = default_for_type(typ, opts)
    put_attribute(modul, :changeset_fields, (name, typ))
    put_struct_field(modul, name, default)
    virtual = haskey(opts, :virtual) ? opts[:virtual] : false
    if !virtual
        raw = haskey(opts, :read_after_writes) ? opts[:read_after_writes] : false
        raw && put_attribute(modul, :ecto_raw, name)
        gen = haskey(opts, :autogenerate) ? opts[:autogenerate] : false
        if isa(gen, Bool)
            gen && store_type_autogenerate!(modul, name, typ, pk)
        else
            store_mfa_autogenerate!(mod, name, typ, gen)
        end
        pk && put_attribute(modul, :ecto_primary_keys, name)
        put_attribute(modul, :ecto_fields, (name, typ))
    end
end

function store_mfa_autogenerate!(modul::Module, name::Symbol, typ::Symbol, mfa)
    put_attribute(modul, :ecto_autogenerate, (name, mfa))
end

function store_type_autogenerate!(modul::Module, name::Symbol, typ::Symbol, pk::Bool)
    id = autogenerate_id(typ)
    put_attribute(modul, :ecto_autogenerate_id, (name, id))
    put_attribute(modul, :ecto_autogenerate, (name, (:typ, :autogenerate, [])))
end

function autogenerate_id(id)
    id in [:id, :binary_id] ? id : nothing
end

function default_for_type(typ, opts)
    haskey(opts, :default) ? opts[:default] : nothing
end

## attributes
__attributes__ = Dict{Module, Any}()

function register_attribute(modul::Module, key::Symbol; kw...)
end

function put_attribute(modul::Module, key::Symbol, obj::Any)
    push!(__attributes__[modul][key], obj)
end

function get_attribute(modul::Module, key::Symbol)
    __attributes__[modul][key]
end

function put_struct_field(modul::Module, name::Symbol, assoc::Any)
    put_attribute(modul, :struct_fields, (name,assoc))
end

function __init_attributes__(modul::Module)
    if !haskey(__attributes__, modul)
        __attributes__[modul] = Dict{Symbol,Any}(
            # beforehands
            :primary_key => nothing,
            :schema_prefix => "",
            :foreign_key_type => nothing,
            :timestamps_opts => nothing,
            :derive => nothing,
            # attributes
            :changeset_fields => Vector{Tuple{Symbol,Any}}(),
            :struct_fields => Vector{Tuple{Symbol,Any}}(),
            :ecto_primary_keys => Vector(),
            :ecto_fields => Vector{Tuple{Symbol,Any}}(),
            :ecto_assocs => Vector(),
            :ecto_embeds => Vector(),
            :ecto_raw => Vector(),
            :ecto_autogenerate => Vector{Tuple{Symbol,Any}}(),
            :ecto_autoupdate => Vector(),
            :ecto_autogenerate_id => Vector{Tuple{Symbol,Any}}()
        )
    end
end

end # module Ecto.Schema

import Base: rem, |>
import .MirrorTypes: plain_to_mirrorstruct, mirrorstruct_to_schema

function rem(modul::Module; vec...)::Schema.t
    Schema.t(modul, Schema.Assoc(vec))
end

function rem(typ::Type; kw...)::Schema.t
    data = nothing
    mirror = plain_to_mirrorstruct(typ)
    modul = mirrorstruct_to_schema(mirror)
    struct = Assoc(kw)
    Schema.t(modul, struct)
end

function |>(schema::Schema.t, func::Function, args...; kw...)
    func(schema, args...; kw...)
end

function |>(schema::Schema.t, meta::Schema.PutMeta)::Schema.t
    for (k,v) in meta.opts
        setfield!(schema.metadata, k, v)
    end
    schema
end

function put_meta(; kw...)::Schema.PutMeta
    Schema.PutMeta(Dict(kw))
end

## Errors
type NoPrimaryKeyFieldError
    message
    NoPrimaryKeyFieldError(schema::Schema.t) = new("schema has no primary key")
end

type NoPrimaryKeyValueError
    message
end

type ChangeError
    message
end


text:
module UDPConnections

export UDPConnection, udp_connect, udp_send, udp_receive, udp_close

using Sockets

mutable struct UDPConnection
    sock::UDPSocket
    target_ip::IPv4
    send_port::Int
    receive_port::Int
end

"""
    UDPConnection

Create an object for communicating over UDP using `udp_send` and `udp_receive`.
Don't forget to close it with `udp_close(conn)` or `close(conn.sock)`!
"""
function UDPConnection(target_ip::IPv4 = ip"127.0.0.1", send_port::Int = 2000, receive_port::Int = 2001)
    return UDPConnection(UDPSocket(), target_ip, send_port, receive_port)
end

"""
    udp_connect
"""
function udp_connect(conn::UDPConnection)
    # TCP equivalent:
    # sock = listen(ip"127.0.0.1", 9000)
    # connect(target_ip, send_port)
    # Create a socket and bind it in order to listen. 0.0.0.0 means "any".
    if !bind(conn.sock, ip"0.0.0.0", conn.receive_port)
        warn("Could not bind. Closing socket.")
        close(conn.sock)
    end
end

"""
    udp_connect(conn) do conn ...
"""
function udp_connect(f::Function, conn::UDPConnection)
    try
        udp_connect(conn)
        f(conn)
    finally
        udp_close(conn)
    end
end

"""
    udp_connect(target_ip, send_port, receive_port) do conn ...
"""
function udp_connect(f::Function, target_ip::IPv4, send_port::Int, receive_port::Int)
    conn = UDPConnection(target_ip, send_port, receive_port)
    try
        f(conn)
    finally
        udp_close(conn)
    end
end

"""
    udp_send

Send some data over to the target.

```julia
udp_send(conn, var1, var2, varn)
```
"""
function udp_send(conn::UDPConnection, args...)
    # TCP equivalent
    # write(conn.sock, Ref(args))
    raw = vcat(map(a -> reinterpret(UInt8, isa(a, Array) ? a : [a]), args)...)
    send(conn.sock, conn.target_ip, conn.send_port, raw)
end

"""
    udp_receive

Receive some data from the target. This uses the input arguments as "examples"
to create the mapping to the output arguments.

```julia
var1, var2, varn = udp_receive(conn, var1, var2, varn)
```
"""
function udp_receive(conn::UDPConnection, args...)
    # TCP equivalent
    # return read(conn.sock, Ref(args))[]
    raw = recv(conn.sock) # A bunch of UInt8s will come in.
    outs = [args] # Create an array so we can point to the args tuple.
    n = min(sizeof(args), sizeof(raw)); # Figure out how much to read.
    unsafe_copy!(reinterpret(Ptr{UInt8}, pointer(outs)), pointer(raw), n)
    return outs[1]
end

"""
    udp_close

Close out the UDP socket. Always do this when done. Use try/catch to make sure.
"""
function udp_close(conn::UDPConnection)
    close(conn.sock)
end

end


text:
mutable struct EnergyMarket{R,T,P} # Assumes completely inelastic demand

    # Parameters
    demand::Array{Float64,3}  # (MW, r x t x p)
    pricecap::Vector{Float64} # ($/MW, r)

    # Variables
    shortfall::Array{VariableRef,3} # Load shortfall variable (MW, r x t x p)

    # Expressions
    shortfallcost::Matrix{ExpressionRef} # Shortfall costs ($, r x p)
    totalshortfallcost::Vector{ExpressionRef} # Shortfall costs ($, r x p)

    # Constraints
    minshortfall::Array{GreaterThanConstraintRef,3} # Minimum load shortfall (r x t x p)
    marketclearing::Array{EqualToConstraintRef,3} # Power balance (r x t x p)

    function EnergyMarket{}(demand::Array{Float64,3}, pricecap::Vector{Float64})

        R, T, P = size(demand)

        @assert all(x -> x >= 0, demand)
        @assert all(x -> x >= 0, pricecap)
        @assert length(pricecap) == R

        new{R,T,P}(demand, pricecap)

    end

end

function setup!(
    market::EnergyMarket{R,T,P}, m::Model,
    ops::Operations{R,G1,G2,G3,I,T,P}, periodweights::Vector{Float64},
    s::AbstractScenario
) where {R,G1,G2,G3,I,T,P}

    invprob = s.investmentproblem
    Rs = invprob.regionnames
    Ts = string.(1:T)
    Ps = invprob.periodnames

    # Variables

    market.shortfall =
        @variable(m, [r in 1:R, t in 1:T, p in 1:P])
    varnames!(market.shortfall, "shortfall_energy_$(s.name)", Rs, Ts, Ps)

    # Expressions

    market.shortfallcost =
        @expression(m, [r in 1:R, p in 1:P],
                    sum(market.shortfall[r,t,p] * market.pricecap[r]
                        for t in 1:T))

    market.totalshortfallcost =
        @expression(m, [r in 1:R],
            sum(market.shortfallcost[r,p] * periodweights[p] for p in 1:P))

    # Constraints

    market.minshortfall =
        @constraint(m, [r in 1:R, t in 1:T, p in 1:P],
                    market.shortfall[r,t,p] >= 0)

    market.marketclearing =
        @constraint(m, [r in 1:R, t in 1:T, p in 1:P],
                    energy(ops, r, t, p) + market.shortfall[r,t,p] ==
                    market.demand[r,t,p] + ops.transmission.exports[r,t,p])

    return

end

welfare(x::EnergyMarket) = -sum(x.shortfallcost)


text:
## Point3-vector identity (typealias)
@test Point3(2, 1, 3) == Vec3(2, 1, 3)

## Vector operators
@test Vec3(1, 2, 3) + Vec3(3, 4, 0) == Vec3(4, 6, 3)
@test Vec3(1, 2, 3) - Vec3(3, 4, 0) == Vec3(-2, -2, 3)
@test Vec3(1, 2, 3) * 1 == Vec3(1, 2, 3) # identity
@test Vec3(1, 2, 3) * 0 == Vec3(0, 0, 0) # null
@test Vec3(1, 2, 3) * 2 == Vec3(2, 4, 6) # scalar multiplication
@test Vec3(1, 2, 3) / 2 == Vec3(0.5, 1, 1.5)

## Euclidean norm/magnitude: norm()
@test norm(Vec3(1, 2, 0)) == sqrt(5)
@test norm(Vec3(0, 0, 0)) == 0

## Bounding boxes: BoundingCube
BCT_point_1 = Vec3(1, 1, 1)
BCT_point_2 = Vec3(2, 3, 2)
BCT_point_3 = Vec3(4, 5, 3)
BCT = BoundingCube(BCT_point_1, BCT_point_2, BCT_point_3)

### BoundingCube attributes
@test BCT == BoundingCube(1, 4, 1, 5, 1, 3)
@test height(BCT) == 4
@test width(BCT) == 3
@test depth(BCT) == 2
@test xmin(BCT) == 1
@test xmax(BCT) == 4
@test center(BCT) == Vec3(2.5, 3, 2)
@test xrange(BCT) == (1, 4)
@test yrange(BCT) == (1, 5)
@test zrange(BCT) == (1, 3)

### BoundingCube operations
BCT_1 = BoundingCube(2, 3, 4, 5, 6, 7)
BCT_2 = BoundingCube(6, 7, 8, 9, 10, 11)

#### BoundingCube (+)
@test BCT_1 + BCT_2 == BoundingCube(2, 7, 4, 9, 6, 11)

#### BoundingCube (&)
@test BCT_1 & BCT_2 == BoundingCube(NaN, NaN, NaN, NaN, NaN, NaN)

BCT_3 = BoundingCube(1, 5, 1, 5, 1, 5)
BCT_4 = BoundingCube(2, 8, 2, 8, 2, 8)

@test BCT_3 & BCT_4 == BoundingCube(2, 5, 2, 5, 2, 5)

#### deform()
@test deform(BCT_1, -1, 2, -3, 4, -5, 6) == BoundingCube(1, 5, 1, 9, 1, 13)

#### shift()
@test shift(BCT_1, -1, 2, -3) == BoundingCube(1, 2, 6, 7, 3, 4)

#### scale()
@test BCT_1 * 3 == BoundingCube(1, 4, 3, 6, 5, 8)

#### isinside()
@test is_inside(BCT_1, Point3(2.5, 4.5, 6.5))
@test is_inside(BCT_1, Point3(2, 4, 6))
@test is_inside(BCT_1, Point3(1, 3, 5)) == false


text:
# This file contains functions for the Cost Efficiency Additive DEA model
"""
    CostAdditiveDEAModel
An data structure representing a cost additive DEA model.
"""
struct CostAdditiveDEAModel <: AbstractCostDEAModel
    n::Int64
    m::Int64
    s::Int64
    weights::Symbol
    monetary::Bool
    rts::Symbol
    disposY::Symbol
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    lambda::SparseMatrixCSC{Float64, Int64}
    techeff::Vector
    alloceff::Vector
    normalization::Vector
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    deacostadditive(X, Y, W, model)
Compute cost efficiency using additive data envelopment analysis for
inputs `X`, outputs `Y` and price of inputs `W`.

Model specification:
- `:Ones`: standard additive DEA model.
- `:MIP`: Measure of Inefficiency Proportions. (Charnes et al., 1987; Cooper et al., 1999)
- `:Normalized`: Normalized weighted additive DEA model. (Lovell and Pastor, 1995)
- `:RAM`: Range Adjusted Measure. (Cooper et al., 1999)
- `:BAM`: Bounded Adjusted Measure. (Cooper et al, 2011)
- `:Custom`: User supplied weights.

# Optional Arguments
- `rts=:VRS`: chooses variable returns to scale. For constant returns to scale choose `:CRS`.
- `rhoX`: matrix of weights of inputs. Only if `model=:Custom`.
- `disposal=:Strong`: chooses strong disposal of outputs. For weak disposal choose `:Weak`.
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6 8];

julia> Y = [1; 1; 1; 1; 1; 1; 1; 1];

julia> W = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1];

julia> deacostadd(X, Y, W, :Ones)
Cost Additive DEA Model
DMUs = 8; Inputs = 2; Outputs = 1
Orientation = Input; Returns to Scale = VRS
Weights = Ones
──────────────────────────────
   Cost  Technical  Allocative
──────────────────────────────
1   0.0        0.0         0.0
2   1.0        0.0         1.0
3   1.0        0.0         1.0
4   3.0        3.0         0.0
5   6.0        6.0         0.0
6   3.0        2.0         1.0
7   3.0        3.0         0.0
8   5.6        5.2         0.4
──────────────────────────────
```
"""
function deacostadd(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector},
    W::Union{Matrix,Vector}, model::Symbol = :Default;
    rts::Symbol = :VRS,
    rhoX::Union{Matrix,Vector,Nothing} = nothing,
    dispos::Symbol = :Strong, monetary::Bool = false,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::CostAdditiveDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    nw, mw = size(W, 1), size(W, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nw != nx
        throw(DimensionMismatch("number of rows in W and X ($nw, $nx) are not equal"));
    end
    if mw != m
        throw(DimensionMismatch("number of columns in W and X ($mw, $m) are not equal"));
    end

    if dispos != :Strong && dispos != :Weak
        throw(ArgumentError("`dispos` must be :Strong or :Weak"));
    end

    # Default behaviour
    if model == :Default
        # If no weights are specified use :Ones
        if rhoX === nothing
            model = :Ones
        else
            model = :Custom
        end
    end

    # Get weights based on the selected model
    if model != :Custom
        # Display error if both model and weights are specified
        if rhoX !== nothing
            throw(ArgumentError("`rhoX` not allowed if model != :Custom"));
        end

        # Get weights for selected model
        rhoX, rhoY = deaaddweights(X, Y, model, orient = :Input)
    end

    if size(rhoX) != size(X)
        throw(DimensionMismatch("size of rhoX and X ($(size(rhoX)), $(size(X))) are not equal"));
    end

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(:LP)
    end

    # Get minimum cost targets and lambdas
    n = nx

    Xtarget, clambdaeff = deamincost(X, Y, W, rts = rts, dispos = dispos, optimizer = optimizer)
    Ytarget = Y[:,:]

    # Cost, technical and allocative efficiency
    cefficiency  = vec(sum(W .* X, dims = 2) .- sum(W .* Xtarget, dims = 2))
    normalization = vec(minimum(W ./ rhoX, dims = 2))

    if model == :Custom
        rhoY = ones(size(Y))
        techefficiency = efficiency(deaadd(X, Y, model, orient = :Input, rts = rts, rhoX = rhoX, rhoY = rhoY, disposY = dispos, optimizer = optimizer))
    else
        techefficiency = efficiency(deaadd(X, Y, model, orient = :Input, rts = rts, disposY = dispos, optimizer = optimizer))
    end

    if monetary
        techefficiency = techefficiency .* normalization
    else
        cefficiency = cefficiency ./ normalization
    end

    allocefficiency = cefficiency - techefficiency

    return CostAdditiveDEAModel(n, m, s, model, monetary, rts, dispos, names, cefficiency, clambdaeff, techefficiency, allocefficiency, normalization, Xtarget, Ytarget)

end

function Base.show(io::IO, x::CostAdditiveDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    eff = efficiency(x)
    disposY = x.disposY
    dmunames = names(x)

    techeff = efficiency(x, :Technical)
    alloceff = efficiency(x, :Allocative)

    if !compact
        print(io, "Cost Additive DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Orientation = Input")
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")
        print(io, "Weights = ", string(x.weights))
        print(io, "\n")

        if disposY == :Weak println(io, "Weak disposability of outputs") end

        show(io, CoefTable(hcat(eff, techeff, alloceff), ["Cost", "Technical", "Allocative"], dmunames))
    end

end


text:
function test_GRPFParams()
    ga = GRPFParams()
    gb = GRPFParams(100, 500000, 3, 5000, 1e-9, false)
    gc = GRPFParams(200, 10000, 3, 10000, 1e-8, true)

    # equality
    @test ga == gb
    @test ga != gc

    # constructors
    @test GRPFParams(5000, 1e-3) == GRPFParams(100, 500000, 3, 5000, 1e-3, false)
    @test GRPFParams(5000, 1e-3, true) == GRPFParams(100, 500000, 3, 5000, 1e-3, true)

    @test_logs (:warn,
        "GRPFParams `tess_sizehint` is greater than `maxnodes`") GRPFParams(100, 10, 3,
                                                                    5000, 1e-9, false)
end

function test_functionstructs()
    rmin, rmax = -1.0, 1.0
    imin, imax = -1.0, 1.0

    g2f = RP.Geometry2Function(rmin, rmax, imin, imax)
    @test eltype(g2f) == Float64

    fval = complex(-1.0, 1.0)
    f2g = RP.fcn2geom(fval, rmin, rmax, imin, imax)

    @test reim(f2g) == (RP.MINCOORD, RP.MAXCOORD)
    @test g2f(f2g) == fval

    fval = complex(0.2, -0.1)
    f2g = RP.fcn2geom(fval, rmin, rmax, imin, imax)
    @test g2f(f2g) ≈ fval   atol=1e-15  # floating point limitation

    p = RP.Point2D(reim(f2g)...)
    @test g2f(p) ≈ fval     atol=1e-15

    fval2 = complex(0.4, 0.8)
    f2g2 = RP.fcn2geom(fval2, rmin, rmax, imin, imax)
    p2 = RP.Point2D(reim(f2g2)...)
    e = RP.DelaunayEdge(p, p2)
    f1, f2 = g2f(e)
    @test f1 ≈ fval     atol=1e-15
    @test f2 ≈ fval2     atol=1e-15

    rmin, rmax = big(-1.0), big(1.0)
    imin, imax = big(-1.0), big(1.0)
    g2f = RP.Geometry2Function(rmin, rmax, imin, imax)
    @test eltype(g2f) == BigFloat
end

@testset "RootsAndPoles.jl" begin
    test_GRPFParams()
    test_functionstructs()
end


text:
#  multi-value dictionary (multidict)

struct MultiDict{K,V}
    d::Dict{K,Vector{V}}

    MultiDict{K,V}() where {K,V} = new{K,V}(Dict{K,Vector{V}}())
    MultiDict{K,V}(kvs) where {K,V} = new{K,V}(Dict{K,Vector{V}}(kvs))
    MultiDict{K,V}(ps::Pair{K,Vector{V}}...) where {K,V} = new{K,V}(Dict{K,Vector{V}}(ps...))
end

MultiDict() = MultiDict{Any,Any}()
MultiDict(kv::Tuple{}) = MultiDict()
MultiDict(kvs) = multi_dict_with_eltype(kvs, eltype(kvs))

multi_dict_with_eltype(kvs, ::Type{Tuple{K,Vector{V}}}) where {K,V} = MultiDict{K,V}(kvs)
function multi_dict_with_eltype(kvs, ::Type{Tuple{K,V}}) where {K,V}
    md = MultiDict{K,V}()
    for (k,v) in kvs
        insert!(md, k, v)
    end
    return md
end
multi_dict_with_eltype(kvs, t) = MultiDict{Any,Any}(kvs)

MultiDict(kv::AbstractArray{Pair{K,V}}) where {K,V}  = MultiDict(kv...)
function MultiDict(ps::Pair{K,V}...) where {K,V}
    md = MultiDict{K,V}()
    for (k,v) in ps
        insert!(md, k, v)
    end
    return md
end

## Functions

## Most functions are simply delegated to the wrapped Dict

@delegate MultiDict.d [ Base.haskey, Base.get, Base.get!, Base.getkey,
                        Base.getindex, Base.length, Base.isempty, Base.eltype,
                        Base.iterate, Base.keys, Base.values]

Base.sizehint!(d::MultiDict, sz::Integer) = (sizehint!(d.d, sz); d)
Base.copy(d::MultiDict) = MultiDict(d)
Base.empty(d::MultiDict{K,V}) where {K,V} = MultiDict{K,V}()
Base.:(==)(d1::MultiDict, d2::MultiDict) = d1.d == d2.d
Base.delete!(d::MultiDict, key) = (delete!(d.d, key); d)
Base.empty!(d::MultiDict) = (empty!(d.d); d)

function Base.insert!(d::MultiDict{K,V}, k, v) where {K,V}
    if !haskey(d.d, k)
        d.d[k] = V[]
    end
    push!(d.d[k], v)
    return d
end

function Base.in(pr::(Tuple{Any,Any}), d::MultiDict{K,V}) where {K,V}
    k = convert(K, pr[1])
    v = get(d,k,Base.secret_table_token)
    (v !== Base.secret_table_token) && (pr[2] in v)
end

function Base.pop!(d::MultiDict, key, default)
    vs = get(d, key, Base.secret_table_token)
    if vs === Base.secret_table_token
        if default !== Base.secret_table_token
            return default
        else
            throw(KeyError(key))
        end
    end
    v = pop!(vs)
    (length(vs) == 0) && delete!(d, key)
    return v
end
Base.pop!(d::MultiDict, key) = pop!(d, key, Base.secret_table_token)

Base.push!(d::MultiDict, kv::Pair) = insert!(d, kv[1], kv[2])
#Base.push!(d::MultiDict, kv::Pair, kv2::Pair) = (push!(d.d, kv, kv2); d)
#Base.push!(d::MultiDict, kv::Pair, kv2::Pair, kv3::Pair...) = (push!(d.d, kv, kv2, kv3...); d)

Base.push!(d::MultiDict, kv) = insert!(d, kv[1], kv[2])
#Base.push!(d::MultiDict, kv, kv2...) = (push!(d.d, kv, kv2...); d)

Base.count(d::MultiDict) = length(keys(d)) == 0 ? 0 : mapreduce(k -> length(d[k]), +, keys(d))
Base.size(d::MultiDict) = (length(keys(d)), count(d::MultiDict))

# enumerate

struct EnumerateAll
    d::MultiDict
end
enumerateall(d::MultiDict) = EnumerateAll(d)

Base.length(e::EnumerateAll) = count(e.d)

function Base.iterate(e::EnumerateAll)
    V = eltype(eltype(values(e.d)))
    vs = V[]
    dstate = iterate(e.d.d)
    vstate = iterate(vs)
    dstate === nothing || vstate === nothing && return nothing
    k = nothing
    while vstate === nothing
        ((k, vs), dst) = dstate
        dstate = iterate(e.d.d, dst)
        vstate = iterate(vs)
    end
    v, vst = vstate
    return ((k, v), (dstate, k, vs, vstate))
end

function Base.iterate(e::EnumerateAll, s)
    dstate, k, vs, vstate = s
    dstate === nothing || vstate === nothing && return nothing
    while vstate === nothing
        ((k, vs), dst) = dstate
        dstate = iterate(e.d.d, dst)
        vstate = iterate(vs)
    end
    v, vst = vstate
    return ((k, v), (dstate, k, vs, vstate))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@doc raw"""
The Cauchy or Paul wavelet.

In contrast to the `MorletWavelet`, the Cauchy or Paul wavelet is a proper wavelet with a
similar good localization in time and scale.

There are different definitions of the Cauchy wavelet around. Here one is implemented where
the center frequency is always 1 irrespective of the value of α:
```math
 g(t) = h(t) = \left(1-i\,\frac{2\pi\,t}{\alpha}\right)^{-(1+\alpha)}\,,
```
and its reproducing kernel
```math
 P_{g,h}(b, a) = \Gamma(2\alpha+1)\,a^{\alpha+1}\,\left(1+a-\frac{i\,b}{a}\right)^{-(2\alpha+1)}\,.
```
"""
struct CauchyWavelet <: GenericContinuousWavelet
    α::Float64;
    norm::Float64;
    ϵ::Float64

    @doc raw"""
        CauchyWavelet(α::Real; ϵ::Real=1e-2)

    Constructs a new Cauchy wavelet, whith the given α specifying the time-frequency resolution.
    The optional keyword arguent ϵ specifies the cutoff at which the kernel evaluation gets
    truncated. It is defined as the fraction of total power loss of the mother wavelet. Smaller
    values of ϵ will increase the precision of the wavelet transform on the cost of longer kernels
    leading to slower convolutions.
    """
    function CauchyWavelet(α::Real; ϵ::Real=1e-2)
        norm = exp(-2*log(2π) - logabsgamma(2α+1)[1]/2 + logabsgamma(α+1)[1]
          + (2α+2)*log(2)/2 + log(α));
        new(α, norm, ϵ);
    end
end

function eval_analysis(wav::CauchyWavelet, t::Float64)
    (1 - 2im*π*t/wav.α)^(-1-wav.α);
end

eval_synthesis(wav::CauchyWavelet, t::Float64) = wav.norm * eval_analysis(wav, t);

function eval_repkern(wav::CauchyWavelet, a::Float64, b::Float64)
    c = wav.α*log(a) + logabsgamma(2*wav.α-1)[1] - (1+2*wav.α)*log(2π);
    exp(c) * (1+a-2im*π*b/wav.α)^(-1-2*wav.α);
end

function cutoff_time(wav::CauchyWavelet)
    wav.α*sqrt(wav.ϵ^(-2/(wav.α+1))-1)/(2π);
end;

function cutoff_freq(wav::CauchyWavelet)
    1 + 1/( wav.α^2 * (wav.ϵ^(-2 / (wav.α + 1)) - 1) / ((2π)^2) );
end;


text:
@testset "Simulator: Hopper (3D)" begin
    # Reference trajectory
    s = get_simulation("hopper_3D", "flat_3D_lc", "flat")
    model = s.model
    env = s.env

    ref_traj = get_trajectory(s.model, s.env,
    	joinpath(module_dir(), "src/dynamics/hopper_3D/gaits/gait_in_place.jld2"),
    	load_type = :joint_traj)
    T = ref_traj.H
    h = ref_traj.h

    for t = 1:T
    	r = ContactImplicitMPC.residual(model, env, ref_traj.z[t], ref_traj.θ[t], 0.0)
    	@test norm(r, Inf) < 1.0e-5
    end

    # initial conditions
    q1 = ref_traj.q[2]
    v1 = (ref_traj.q[2] - ref_traj.q[1]) ./ h

    # simulator
    sim = ContactImplicitMPC.simulator(s, T, h=h)

    # simulate
    status = ContactImplicitMPC.simulate!(sim, q1, v1)
    @test status
end