text:
@testset "OneHotArray" begin
  using Transformers.Basic: OneHot, onehotsize
  idx = rand(5:10)
  ot = OneHot(10)
  o = OneHot{10}(idx)

  @test onehotsize(o) == 10
  @test o == OneHot(10, idx)
  @test o == ot(idx)

  @test size(o) == (10,)
  for i = 1:10
    @test o[i] == (i == idx)
  end
  @test_throws BoundsError o[11]

  @test UInt32(o) == UInt32(idx)
  @test UInt64(o) == UInt64(idx)
  @test Int32(o) == Int32(idx)
  @test Int64(o) == Int64(idx)

  @test convert(UInt32, o) == UInt32(idx)
  @test convert(UInt64, o) == UInt64(idx)
  @test convert(Int32, o) == Int32(idx)
  @test convert(Int64, o) == Int64(idx)
  @test convert(Int8, o) == Int8(idx)
  @test_throws InexactError convert(Bool, o)
  @test convert(ot, idx) == o
  @test convert(ot, UInt8(idx)) == o
  @test convert(ot, UInt16(idx)) == o
  @test convert(ot, UInt32(idx)) == o
  @test convert(ot, UInt64(idx)) == o
  @test convert(ot, UInt128(idx)) == o
  @test convert(ot, Int8(idx)) == o
  @test convert(ot, Int16(idx)) == o
  @test convert(ot, Int32(idx)) == o
  @test convert(ot, Int64(idx)) == o
  @test convert(ot, Int128(idx)) == o

  @test convert(ot, true) == OneHot(10, 1)
  @test convert(ot, false) == OneHot(10, 0)

  @test convert(OneHot(15), o) == OneHot(15, idx)
  @test convert(OneHot(idx+1), o) == OneHot(idx+1, idx)
  @test_throws Basic.OneHotEncodeError convert(OneHot(idx-1), o)

  @test zero(o) == OneHot(10, 0)
  @test !iszero(o)
  @test iszero(zero(o))

  @test_throws Basic.OneHotEncodeError convert(ot, -1)
  @test_throws Basic.OneHotEncodeError convert(ot, 12)
  @test_throws Basic.OneHotEncodeError convert(ot, typemax(Int64))
  @test_throws Basic.OneHotEncodeError convert(ot, typemax(UInt64))

  ind = rand(1:10, 5, 5)
  oa = OneHotArray{10}(ind)

  @test onehotsize(oa) == 10
  @test OneHotArray(15, oa) |> onehotsize == 15

  @test size(oa) == (10, 5, 5)
  @test oa[:, 5, 5] == oa.onehots[5, 5]
  @test oa[:, 5, :] == OneHotArray(oa.onehots[5, :])
  @test oa[ind[1], 1, 1]

  @test vcat(o, o) == vcat(collect(o), collect(o))
  @test hcat(o, o).onehots == [o,o]
  @test cat(o, o; dims=1) == vcat(collect(o), collect(o))
  @test cat(o, o; dims=2).onehots == [o,o]

  @test vcat(oa, oa) == vcat(collect(oa), collect(oa))
  @test hcat(oa, oa).onehots == vcat(oa.onehots, oa.onehots)
  @test cat(oa, oa; dims=1) == vcat(collect(oa), collect(oa))
  @test cat(oa, oa; dims=2).onehots == vcat(oa.onehots, oa.onehots)
  @test cat(oa, oa; dims=3).onehots == cat(oa.onehots, oa.onehots; dims=2)

  @test reshape(oa, 10, 25) isa OneHotArray
  @test reshape(oa, 10, :) isa OneHotArray
  @test reshape(oa, :, 25) isa OneHotArray
  @test reshape(oa, 50, :) isa Base.ReshapedArray{Bool, 2}
  @test reshape(oa, 5, 10, 5) isa Base.ReshapedArray{Bool, 3}

end


text:
### SymataSyntax

global symatasyntax_inited = false

function _init_symatasyntax()
    if ! symatasyntax_inited
        try
            @eval using SymataSyntax
            global symatasyntax_inited = true
            return true
        catch e
            @warn("Unable to load 'SymataSyntax'. If the module is not installed, " *
                  "read https://github.com/jlapeyre/SymataSyntax.jl`.")
            rethrow(e)
        end
    else
        nothing
    end
    return true
end

### MmaSyntax

@sjdoc MmaSyntax """
    MmaSyntax()

enters a Mathematica-syntax mode. Type `ctrl-d` to return to the usual Symata mode in the terminal.
In Jupyter, type `MmaSyntax[False]` to return to Symata syntax.

In a terminal, the Mathematica-syntax mode is a modified version of the mathics REPL (terminal interface).
Input expressions must be valid Mathematica expressions. Output is printed as Mathematica expressions.
In Jupyter, Mathematica syntax is expected as input.

This mode uses the mathics parser to read input and the mathics formatter to write output. But, all evaluation
is done by Symata. Symata often calls the python package sympy, and builtin Julia functions.

This mode requires the Julia package `SymataSyntax`, which requires the python package `mathics`.

The only relation to Mathematica is the syntax. This mode does not interact with Mathematica in any way and is independent
of Mathematica and Wolfram language and all other Wolfram company products.
"""
@mkapprule MmaSyntax nargs => 0:1

@doap function MmaSyntax()
    (! _init_symatasyntax()) && return
    if isymata_inited()
        isymata_mma_mode(true)
    else
        Base.invokelatest(SymataSyntax.mmasyntax_REPL)
    end
end

@doap function MmaSyntax(b::Bool)
    (! _init_symatasyntax()) && return
    if isymata_inited()
        isymata_mma_mode(b)
    else
        if b
            SymataSyntax.mmasyntax_REPL()
        else
            symprintln("Use ctrl-d to exit the Mma syntax REPL")
        end
    end
end

@mkapprule SymataSyntax

@doap function SymataSyntax()
    (! _init_symatasyntax()) && return
    if isymata_inited()
        isymata_mma_mode(false)
    else
        symprintln("Use ctrl-d to exit the Mma syntax REPL")
    end
end

### Mathics

@sjdoc Mathics """
    Mathics()

enters the mathics shell. Type `ctrl-d` to return to the Symata.

The mathics process is separate from the Symata process. In particular, setting a variable in mathics
does nothing in Symata and vice-versa.
"""

@mkapprule Mathics nargs => 0

@doap function Mathics()
    (! _init_symatasyntax()) && return
    SymataSyntax.mathics_REPL()
end

### GetMma

@mkapprule GetMma

@sjdoc GetMma """
    GetMma(filename)

read and evaluate Mathematica-syntax expressions from file `filename`.
"""

@doap function GetMma(fname::String)
    (! _init_symatasyntax()) && return
    SymataSyntax.read_file(fname)
end

ext: jl


text:
#=--- TerminalGraphics / integrations / Cairo.jl --------------------------=#

#=-------------------------------------------------------------------------=#
#
#   Integrations for 'Cairo.jl'
#
#   When Cairo is loaded, automatically allow displaying a CairoSurface on
#   the console by drawing to a temporary png file.
#
#=-------------------------------------------------------------------------=#

import Base: display
display(d::TerminalGraphics.SixelDisplay, surface::Cairo.CairoSurface) = begin
    # define a temporary filename and plot to it
    fn = tempname() * ".png"
    isfile(fn) &&
        error("The intended temporary file '$fn' already exists. Try again.")
    Cairo.write_to_png(surface, fn)
    # draw the 'png'
    plottype = string(typeof(plt))
    title = string(typeof(surface))
    size = "size ($(surface.width), $(surface.height))"
    println(d.io, "Drawing $title of $size, temporarily at $fn")
    TerminalGraphics.draw(d.io, fn)
    # delete the 'png'
    isfile(fn) && Base.Filesystem.rm(fn)
    isfile(fn) &&
        warn("The temporary file '$fn' couldn't be deleted. Please inspect.")
end

#=--- end of file ---------------------------------------------------------=#


text:
using Turing, ReverseDiff, Memoization
using Statistics
using Random
Turing.setadbackend(:reversediff)
Turing.setrdcache(true)
Turing.turnprogress(false)
using Distributions

@model g(rn,T) = begin
    sigma ~ Uniform(0.0,20.0)

    rn ~ MvNormal(zeros(T),sigma)
end

for i in 1:10
# create a sample from a normal
# sample from a Gaussian
    N = 1000
    d = Normal(0.0, i)
    data=rand(d,N)
    chn = sample(g(data,N),NUTS(0.65),1000)
    Turing.emptyrdcache()

    println("Std of chain: ",mean(Array(chn[:sigma])),"Std of data: ",std(data))
end


text:
"Provide routines for making fake Gaussian models which have analytic Fourier transforms, for the purposes of testing."
module gauss

export imageGauss, FTGauss

using LinearAlgebra
using ..constants

"
    imageGauss(ll::AbstractVector{Float64}, mm::AbstractVector{Float64}, p::Vector{Float64}, k::Real; theta=0)

Given two arrays of ``l`` and ``m`` coordinates, corresponding to ``x`` and ``y``, fill an array of the Gaussian image
following the MATLAB convention for images (each row corresponds to a different y value).

`p0` is a vector of `[mu_RA, mu_DEC, sigma_x, sigma_y, rho]` in units of arcseconds.
``\\rho`` is the correlation of the Gaussian, ranging from 0 to 1.

```math
\\rho = \\frac{\\sigma_{xy}}{\\sigma_x \\sigma_y}
```

``\\mu_\\alpha`` and ``\\mu_\\delta`` are the locations of the centroid emission relative to the
image origin (RA=0, DEC=0).

``k`` is a scaling pre-factor to adjust the amplitude of the Gaussian.

The image intensity as a function of ``l`` and ``m`` (sky plane coordinates) is

```math
    I(l,m) = \\frac{k}{2 \\pi \\sqrt{|\\boldsymbol{\\Sigma}|}} \\exp \\left \\{ -\\frac{1}{2} \\boldsymbol{R}^\\mathrm{T} \\boldsymbol{\\Sigma}^{-1} \\boldsymbol{R} \\right \\}
```

Where ``\\boldsymbol{\\Sigma}`` is
```math
\\boldsymbol{\\Sigma} = \\left [ \\begin{array}{cc}
\\sigma_x & \\sigma_{xy} \\\\
\\sigma_{xy} & \\sigma_y \\\\
\\end{array} \\right ]
```

``\\theta`` is a desired angle to rotate the Gaussian about the origin. To make sense of the order of operations, the rotation must occur about the origin. A positive value means a counter-clockwise rotation of the Gaussian from North towards East.
"
function imageGauss(ll::AbstractVector{Float64}, mm::AbstractVector{Float64}, p::Vector{Float64}, k::Real; theta=0)

    # Both ll and mm increase with array index
    nx = length(ll)
    ny = length(mm)

    img = Array{Float64}(undef, ny, nx)
    mu = p[1:2] * arcsec #ll and mm shifts
    sigma_x = p[3] * arcsec
    sigma_y = p[4] * arcsec
    rho = p[5]
    sigma_xy = rho * sigma_x * sigma_y
    Sigma = Float64[[sigma_x^2, sigma_xy] [sigma_xy, sigma_y^2]]
    rot = [[cosd(theta), +sind(theta)] [-sind(theta), cosd(theta)]]
    # Sigma = Diagonal((p[3:4] * arcsec).^2) #Convert from arcsec to radians
    pre = 1.0 / (2pi * sqrt(det(Sigma))) * k
    for j=1:ny
        for i=1:nx
            # Rotation of Gaussian occurs about origin, so subtract off mean first.
            R = rot * (Float64[ll[i] , mm[j]] - mu)
            img[j, i] = pre * exp(-0.5 * (R' * (Sigma\R))[1]) # backslash solves for Sigma^{-1}R
        end
    end
    return img
end

"
    FTGauss(uu::Float64, vv::Float64, p::Vector{Float64}, k::Real)

Given u and v coordinates in [kλ], evaluate the analytic FT of the
aforementioned Gaussian.

`p` is a length 5 vector of `[mu_RA, mu_DEC, sigma_x, sigma_y, rho]` in units of arcseconds, corresponding to the **image** plane.

N.B. Here Sigma refers to the (same) covariance matrix in the *image* domain.

This function always returns a complex value."
function FTGauss(uu::Float64, vv::Float64, p::Vector{Float64}, k::Real; theta=0)
    uu = uu .* 1e3 #[λ]
    vv = vv .* 1e3 #[λ]

    rot = Float64[[cosd(theta), -sind(theta)] [sind(theta), cosd(theta)]]

    mu = p[1:2] * arcsec #ll and mm shifts
    sigma_x = p[3] * arcsec
    sigma_y = p[4] * arcsec
    rho = p[5]
    sigma_xy = rho * sigma_x * sigma_y
    Sigma = Float64[[sigma_x^2, sigma_xy] [sigma_xy, sigma_y^2]]

    # mu_RA, mu_DEC = p[1:2]
    # mu = Float64[mu_RA, mu_DEC] * arcsec #ll and mm shifts

    R0 = Float64[uu, vv]

    # Position to evaluate the rotated Gaussian centered at the origin. Evaluating the Gaussian as if we were
    # in uu', vv'.
    R = rot * R0
    model = k * exp(-2 * (pi^2) * (R' * Sigma * R)[1])
    # in this case, Sigma serves as an inverse matrix, since the covariance matrix corresponds to the *image* plane.
    # Sigma = Diagonal((p[3:4] * arcsec).^2) #Convert from arcsec to radians

    # Phase shift should be evaluated using the original coordinates and mu, since this is the direction that we are actually moving things. And, the interferometer samples uu,vv, not uu', vv'.
    phase_shift = exp(-2pi * 1.0im * (R0' * mu)[1]) # Not actually in polar phase form

    return model * phase_shift
end

"
    FTGauss(uu::AbstractVector{Float64}, vv::AbstractVector{Float64}, p::Vector{Float64}, k::Int)

`p` is a length 5 vector of `mu_RA`, `mu_DEC`, `sigma_RA`, `sigma_DEC`, `rho`.

Given two arrays of u and v coordinates in [kλ], fill an array with the
analytic FT of `imageGauss` evaluated at every pairwise (u,v) pair."
function FTGauss(uu::AbstractVector{Float64}, vv::AbstractVector{Float64}, p::Vector{Float64}, k::Real; theta=0)
    nu = length(uu)
    nv = length(vv)
    # Both uu and vv increase with array index
    img = Array{ComplexF64}(undef, nv, nu)
    for j=1:nv
        for i=1:nu
            img[j, i] = FTGauss(uu[i], vv[j], p, k, theta=theta)
        end
    end
    return img
end

end # module


text:
using Transformations
using PlotRecipes

export
    ChainPlot

# -------------------------------------------------------------------------
# Chain

@recipe function f(chain::Chain)
    tmap = Dict()
    for (i,t) in enumerate(chain.ts)
        tmap[input_node(t)] = i
    end
    # connect input (1) --> ts[1], and ts[end] --> output (-1)
    n = length(chain.ts)
    source = Int[n+1, n]
    destiny = Int[1, n+2]
    colors = Symbol[]
    shapes = Symbol[]
    for (i,t) in enumerate(chain.ts)
        c,s = if isa(t, Activation)
            :orangered, :circle
        elseif isa(t, Linear) || isa(t, Affine)
            :steelblue, :hexagon
        else
            :lightgreen, :rect
        end
        push!(colors, c)
        push!(shapes, s)
        for tonode in output_node(t).tonodes
            push!(source, i)
            push!(destiny, tmap[tonode])
        end
    end
    append!(colors, [:cyan,:cyan])
    append!(shapes, [:circle,:circle])
    method --> :tree
    names --> vcat(map(string, chain.ts), "Input", "Output")
    markercolor --> colors
    markershape --> shapes
    markeralpha --> 0.6
    PlotRecipes.GraphPlot((source, destiny))
end

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

"""
A ChainPlot is a bunch of TracePlots that are updated with the param/input/output
values and gradients from a Transformations.Chain.  An example:

```julia
using MLPlots, Transformations

# build a neural net
nn = nnet(2,3,[2,10],:softplus)

# build the ChainPlot.  attributes pass through to the plot
cp = ChainPlot(nn, size=(1000,1000))

for i=1:20
    # make some random updates to the Chain
    transform!(nn, rand(2))
    grad!(nn, rand(3))

    # update the ChainPlot and display it
    update!(cp)
    gui()
end
```
"""
type ChainPlot
    chain::Chain
    plt::Plot
    nt::Int # number of transformations
    np::Int # number of params
    lastx::Float64
    θtps::Vector{TracePlot}
    ∇tps::Vector{TracePlot}
    ytps::Vector{TracePlot}
    dytps::Vector{TracePlot}
end

function add_traceplots(plt, valtps, gradtps, n, idxbase, i, idxdiff, title, add_yguide, kw) #, xguide="")
    push!(valtps, TracePlot(n;
        sp=plt[idxbase+i],
        title = title,
        yguide = (add_yguide ? "Value" : ""),
        # xguide = xguide,
        kw...
    ))
    push!(gradtps, TracePlot(n;
        sp = plt[idxbase+i+idxdiff],
        yguide = (add_yguide ? "Grad" : ""),
        # xguide = xguide,
        kw...
    ))
end

function ChainPlot(chain::Chain; kw...)
    nparams = count(t->length(params(t))>0, chain.ts)
    n = length(chain.ts)

    plt = plot(layout = @layout([grid(2,nparams)
                                 grid(2,n+1)]),
                leg=false, titlefont=font(10))

    θtps = TracePlot[]
    ∇tps = TracePlot[]
    ytps = TracePlot[]
    dytps = TracePlot[]
    pidx = 1
    for (i,t) in enumerate(chain.ts)
        np = length(params(t))
        if np > 0
            add_traceplots(plt, θtps, ∇tps, np, 0, pidx, nparams, string(t), i==1, kw) #, "params $i $t")
            pidx += 1
        end

        if i == 1
            nin = input_length(t)
            add_traceplots(plt, ytps, dytps, nin, 2nparams, i, n+1, "Input", true, kw) #, "input $i $t")
        end

        nout = output_length(t)
        add_traceplots(plt, ytps, dytps, nout, 2nparams+1, i, n+1, string(t), false, kw) #, "output $i $t")
    end

    ChainPlot(chain, plt, n, nparams, 0.0, θtps, ∇tps, ytps, dytps)
end

function LearnBase.update!(cp::ChainPlot, x = cp.lastx+1)
    pidx = 1
    for (i,t) in enumerate(cp.chain.ts)
        np = length(params(t))

        if np > 0
            push!(cp.θtps[pidx], x, params(t))
            push!(cp.∇tps[pidx], x, grad(t))
            pidx += 1
        end
        idxbase = 2 * cp.np + 1
        if i == 1
            push!(cp.ytps[1], x, input_value(t))
            push!(cp.dytps[1], x, input_grad(t))
        end
        push!(cp.ytps[i+1], x, output_value(t))
        push!(cp.dytps[i+1], x, output_grad(t))
    end
    cp.lastx = x
    return
end


text:
using DiffEqBase
using Base.Test

@time @testset "Number of Parameters Calculation" begin include("numargs_test.jl") end
@time @testset "Solution Interface" begin include("solution_get_tests.jl") end
@time @testset "Extended Functions" begin include("extended_function_tests.jl") end
@time @testset "Callbacks" begin include("callbacks.jl") end
@time @testset "Constructed Parameterized Functions" begin include("constructed_pf_test.jl") end
@time @testset "Plot Variables" begin include("plot_vars.jl") end


text:
using RosDataProcess

ts = range(-2pi, stop=2pi, length=500)

a = TimeSeries(sin, ts)
b = differentiate(a)

#plot(a)
#plot(b)

max_stars_repo_name: mcorah/RosDataProcess

max_issues_repo_name: mcorah/RosDataProcess

max_forks_repo_name: mcorah/RosDataProcess


text:
# --------------------------------------------------------------------------
# ACE.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

import ACE.SphericalHarmonics: SHBasis, index_y

@doc raw"""
`struct RnYlmBasisFcn` : 1-particle basis function specification
for bases of the type ``P \otimes Y`` with `P::ScalarBasis` and `Y::SHBasis`
"""
struct RnYlmBasisFcn <: OnepBasisFcn
   n::Int
   l::Int
   m::Int
end

function RnYlmBasisFcn(t::VecOrTup)
   @assert length(t) == 3
   return RnYlmBasisFcn(t[1], t[2], t[3])
end

Base.show(io::IO, b::RnYlmBasisFcn) = print(io, "nlm[$(b.n),$(b.l),$(b.m)]")

write_dict(b::RnYlmBasisFcn) =
   Dict("__id__" => "ACE_RnYlmBasisFcn",
        "nlm" => [ b.n, b.l, b.m ] )

read_dict(::Val{:SHIPs_RnYlmBasisFcn}, D::Dict) =
   read_dict(Val{:ACE_RnYlmBasisFcn}(), D)

read_dict(::Val{:ACE_RnYlmBasisFcn}, D::Dict) =
   RnYlmBasisFcn(D["nlmz"]...)

scaling(b::RnYlmBasisFcn, p) = b.n^p + b.l^p + abs(b.m)^p

@doc raw"""
`struct RnYlm1pBasis <: OneParticleBasis`

One-particle basis of the form
```math
\phi_{nlm}({\bm r}) = J_{n}(r) Y_l^m(\hat{\br r})
```
where ``J_{n}`` denotes a radial basis.
"""
struct RnYlm1pBasis{T, TJ <: ScalarBasis{T}} <: OneParticleBasis{T}
   J::TJ           # J = R
   SH::SHBasis{T}  # SH = Ylm
   spec::Vector{RnYlmBasisFcn}
end

function RnYlm1pBasis(J::ScalarBasis{T}, zlist::SZList,
                      spec::Vector{RnYlmBasisFcn}) where {T}
   # now get the maximum L-degree to generate the SH basis
   maxL = maximum(b.l for b in spec)
   SH = SHBasis(maxL, T)
   # construct a preliminary Aindices array to get an incorrect basis
   NZ = length(zlist)
   P = RnYlm1pBasis(J, SH, zlist, spec,
                       Matrix{UnitRange{Int}}(undef, NZ, NZ))
   return P
end

function RnYlm1pBasis(J::ScalarBasis;
                         species = :X,
                         D::AbstractDegree = SparsePSHDegree())
   # get a generic basis spec
   spec = _build_PSH_1p_spec(length(J), D, species)
   # construct the basis
   zlist = ZList(species; static=true)
   return RnYlm1pBasis(J, zlist, spec)
end

cutoff(basis::RnYlm1pBasis) = cutoff(basis.J)

Base.length(basis::RnYlm1pBasis, z0::AtomicNumber) =
      numz(basis) *  length(basis.spec)

Base.length(basis::RnYlm1pBasis, iz0::Integer) =
      numz(basis) *  length(basis.spec)

Base.length(basis::RnYlm1pBasis, iz::Integer, iz0::Integer) =
      length(basis.spec)

Base.rand(basis::RnYlm1pBasis) =
      AtomState(rand(basis.zlist.list), ACE.Random.rand_vec(basis.J))

function get_basis_spec(basis::RnYlm1pBasis, z0::AtomicNumber)
   iz0 = z2i(basis, z0)
   len_iz0 = sum(length(basis, iz, iz0) for iz = 1:numz(basis))
   spec = Vector{RnYlmBasisFcn}(undef, len_iz0)
   for iz = 1:numz(basis)
      z = i2z(basis, iz)
      spec[basis.Aindices[iz, iz0]] =
               [ RnYlmBasisFcn(b.n, b.l, b.m, z)   for b in basis.spec ]
   end
   return spec
end

get_basis_spec(basis::RnYlm1pBasis, iz0::Integer, i::Integer) =
   get_basis_spec(basis, i2z(basis, iz0), i)

function get_basis_spec(basis::RnYlm1pBasis, z0::AtomicNumber, i::Integer)
   iz = ((i-1) ÷ length(basis.spec)) + 1
   inew = mod1(i, length(basis.spec))
   b = basis.spec[inew]
   return RnYlmBasisFcn(b.n, b.l, b.m, i2z(basis, iz))
end

_build_specnl(maxdeg, D, z, z0) =
         gensparse(2, maxdeg;               #     n     l    m  z
                    tup2b = t -> RnYlmBasisFcn(t[1]+1, t[2], 0, z),
                    degfun = b -> D(b, z0),
                    ordered = false)

@doc raw"""
`function _build_PSH_1p_spec`

Construct the specification for a ``P \otimes Y`` type 1-particle basis.
These must be treated differently because of the requirements that complete
``l``-blocks are represented in the basis.

See also: `_get_1p_spec`.
"""
function _build_PSH_1p_spec(maxn::Integer, D::AbstractDegree, species)
   # find out what the largest degree is that we can allow:
   specnl = []
   for s1 in species, s2 in species
      z, z0 = AtomicNumber(s1), AtomicNumber(s2)
      maxdeg = maximum( D(RnYlmBasisFcn(n, 0, 0, z), z0) for n = 1:maxn )
      # generate the `spec::Vector{RnYlmBasisFcn}` using maxn
      specnl_zz0 = _build_specnl(maxdeg, D, z, z0)
      append!(specnl, specnl_zz0)
   end
   specnl = unique([ RnYlmBasisFcn(b.n, b.l, 0, 0) for b in specnl ])
   # add the m-parameters
   return [ RnYlmBasisFcn(b.n, b.l, m, 0)
              for b in specnl for m = -b.l:b.l ]
end

_build_PSH_1p_spec(maxn::Integer, D::AbstractDegree, species::Symbol) =
         _build_PSH_1p_spec(maxn::Integer, D, (species,))

# ------------------------------------------------------
#  FIO code

==(P1::RnYlm1pBasis, P2::RnYlm1pBasis) =  ACE._allfieldsequal(P1, P2)

write_dict(basis::RnYlm1pBasis{T}) where {T} = Dict(
      "__id__" => "ACE_RnYlm1pBasis",
           "J" => write_dict(basis.J),
          "SH" => write_dict(basis.SH),
        "spec" => write_dict.(basis.spec),
       "zlist" => write_dict(basis.zlist),
   )

read_dict(::Val{:SHIPs_RnYlm1pBasis}, D::Dict) =
   read_dict(Val{:ACE_RnYlm1pBasis}(), D::Dict)

function read_dict(::Val{:ACE_RnYlm1pBasis}, D::Dict)
   J = read_dict(D["J"])
   SH = read_dict(D["SH"])
   zlist = read_dict(D["zlist"])
   spec = read_dict.(D["spec"])
   P = RnYlm1pBasis(J, SH, zlist, spec,
                 Matrix{UnitRange{Int}}(undef, length(zlist), length(zlist)))
   set_Aindices!(P)
   return P
end

# ------------------------------------------------------
#  Evaluation code

fltype(basis::RnYlm1pBasis{T}) where T = Complex{T}
rfltype(basis::RnYlm1pBasis{T}) where T = T

alloc_temp(basis::RnYlm1pBasis, args...) =
      (
        BJ = alloc_B(basis.J, args...),
        tmpJ = alloc_temp(basis.J, args...),
        BY = alloc_B(basis.SH, args...),
        tmpY = alloc_temp(basis.SH, args...),
       )

function evaluate!(B, tmp, basis::RnYlm1pBasis,
                   X::AbstractState, X0::AbstractState)
   # evaluate the r-basis and the R̂-basis for the current neighbour at R
   evaluate!(tmp.BJ, tmp.tmpJ, basis.J, norm(X.rr))
   evaluate!(tmp.BY, tmp.tmpY, basis.SH, X.rr)
   # add the contributions to the A_zklm
   @inbounds for (i, nlm) in enumerate(basis.spec)
      B[i] = tmp.BJ[nlm.n] * tmp.BY[index_y(nlm.l, nlm.m)]
   end
   return B
end

# function add_into_A!(A, tmp, basis::RnYlm1pBasis,
#                      R, iz::Integer, iz0::Integer)
#    # evaluate the r-basis and the R̂-basis for the current neighbour at R
#    evaluate!(tmp.BJ, tmp.tmpJ, basis.J, norm(R))
#    evaluate!(tmp.BY, tmp.tmpY, basis.SH, R)
#    # add the contributions to the A_zklm
#    @inbounds for (i, nlm) in enumerate(basis.spec)
#       A[i] += tmp.BJ[nlm.n] * tmp.BY[index_y(nlm.l, nlm.m)]
#    end
#    return nothing
# end

# function add_into_A!(A, inds, tmp, basis::RnYlm1pBasis,
#                      R, iz::Integer, iz0::Integer)
#    # evaluate the r-basis and the R̂-basis for the current neighbour at R
#    evaluate!(tmp.BJ, tmp.tmpJ, basis.J, norm(R))
#    evaluate!(tmp.BY, tmp.tmpY, basis.SH, R)
#    # add the contributions to the A_zklm
#    @inbounds for (i, nlm) in zip(inds, basis.spec)
#       A[i] += tmp.BJ[nlm.n] * tmp.BY[index_y(nlm.l, nlm.m)]
#    end
#    return nothing
# end

alloc_temp_d(basis::RnYlm1pBasis, args...) =
      (
        BJ = alloc_B(basis.J, args...),
        tmpJ = alloc_temp(basis.J, args...),
        BY = alloc_B(basis.SH, args...),
        tmpY = alloc_temp(basis.SH, args...),
        dBJ = alloc_dB(basis.J, args...),
        tmpdJ = alloc_temp_d(basis.J, args...),
        dBY = alloc_dB(basis.SH, args...),
        tmpdY = alloc_temp_d(basis.SH, args...),
       )

function add_into_A_dA!(A, dA, tmpd, basis::RnYlm1pBasis, R, iz::Integer, iz0::Integer)
   r = norm(R)
   R̂ = R / r
   # evaluate the r-basis and the R̂-basis for the current neighbour at R
   evaluate_d!(tmpd.BJ, tmpd.dBJ, tmpd.tmpdJ, basis.J, r)
   evaluate_d!(tmpd.BY, tmpd.dBY, tmpd.tmpdY, basis.SH, R)
   # add the contributions to the A_zklm, ∇A
   @inbounds for (i, nlm) in enumerate(basis.spec)
      iY = index_y(nlm.l, nlm.m)
      A[i] += tmpd.BJ[nlm.n] * tmpd.BY[iY]
      dA[i] = (tmpd.dBJ[nlm.n] * tmpd.BY[iY]) * R̂ + tmpd.BJ[nlm.n] * tmpd.dBY[iY]
   end
   return nothing
end


text:
export StateSpace, BoundedStateSpace
export AbstractConfig, Config
export vecs2mat, mat2vecs

# StateSpace
abstract type StateSpace{S<:State} end
## BoundedStateSpace
struct BoundedStateSpace{S,B} <: StateSpace{S}
    bounded_set::B
end
BoundedStateSpace{S}(bounded_set::B) where {S,B} = BoundedStateSpace{S,B}(bounded_set)
BoundedStateSpace(bounded_set::B) where {B} = BoundedStateSpace{typeof(rand(bounded_set)),B}(bounded_set)
BoundedStateSpace(lo::S, hi::S) where {S<:StaticVector} = BoundedStateSpace{S}(AxisAlignedBox(lo, hi))
BoundedStateSpace(lo, hi) = BoundedStateSpace(AxisAlignedBox(lo, hi))
Base.in(x, SS::BoundedStateSpace) = x in SS.bounded_set
Base.rand(SS::BoundedStateSpace{S}) where {S} = convert(S, rand(SS.bounded_set))
@recipe function f(SS::BoundedStateSpace; dims=(1, 2)) #, statespace_color=:black, statespace_linecolor=statespace_color)
    dims      --> dims
    linecolor :=  :match
    fillalpha :=  0
    label     --> ""
    SS.bounded_set
end

# Config
abstract type AbstractConfig end
const Config = Union{AbstractConfig, AbstractVector{<:Number}}

# Matrix <--> Vector of (Column) Vectors Conversion
vecs2mat(V::AbstractVector{S}) where {S<:StaticArray}   = reshape(reinterpret(eltype(S), V), (length(S), length(V)))
vecs2mat(V::AbstractVector{S}) where {S<:AbstractArray} = reduce(hcat, V)
mat2vecs(::Type{S}, M::Matrix) where {S<:StaticArray}   = reshape(reinterpret(S, M), size(M, 2))
mat2vecs(::Type{S}, M::Matrix) where {S<:AbstractArray} = [convert(S, M[:,i]) for i in 1:size(M,2)]

max_stars_repo_name: pbpf/clean

max_issues_repo_name: pbpf/clean

max_forks_repo_name: pbpf/clean


text:
import Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid

function ImmersedBoundaryGrid(grid::ConformalCubedSphereGrid, immersed_boundary)
    faces = Tuple(ImmersedBoundaryGrid(get_face(grid, i), immersed_boundary) for i = 1:6)
    FT = eltype(grid)
    face_connectivity = grid.face_connectivity

    cubed_sphere_immersed_grid = ConformalCubedSphereGrid{FT, typeof(faces), typeof(face_connectivity)}(faces, face_connectivity)

    return cubed_sphere_immersed_grid
end

import Oceananigans.Operators: Γᶠᶠᵃ

@inline function Γᶠᶠᵃ(i, j, k, ibg::ImmersedBoundaryGrid{F, TX, TY, TZ, G, I}, u, v) where {F,TX,TY,TZ,G<:ConformalCubedSphereFaceGrid,I}
    Γᶠᶠᵃ(i, j, k, ibg.grid, u, v)
end


text:
# tests for sum, mean, and var, which implemented by calling other atoms

x = create_var(5, 3)
s = sum(x)
m = mean(x)
test_expr(s, (x,))
test_expr(m, (x,))

t = sum(x, 1)
n = mean(x, 1)
test_expr(t, (x,))
test_expr(n, (x,))

u = sum(x, 2)
p = mean(x, 2)
test_expr(u, (x,))
test_expr(p, (x,))

v = var(x)
test_expr(v)

info("All sum/mean/var tests passed")


text:
include(joinpath(pwd(), "models/box.jl"))
include(joinpath(pwd(), "src/constraints/contact.jl"))

# Initial state
mrp = MRP(UnitQuaternion(RotY(pi / 10.0) * RotX(pi / 15.0)))
q1 = [0.0, 0.0, 2.5, mrp.x, mrp.y, mrp.z]
v1 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
v2 = v1 - G_func(model,q1) * h
q2 = q1 + 0.5 * h * (v1 + v2)

# Simulate contact model one step
function step_contact(model, x1, u1, h)
    # Horizon
    T = 2

    # Objective
    obj = PenaltyObjective(1.0e5, model.m)

    # Constraints
    con_contact = contact_constraints(model, T)

    # Bounds
    _uu = Inf * ones(model.m)
    _uu[model.idx_u] .= u1
    _ul = zeros(model.m)
    _ul[model.idx_u] .= u1
    ul, uu = control_bounds(model, T, _ul, _uu)

    xl, xu = state_bounds(model, T, x1 = x1)

    # Problem
    prob = trajectory_optimization_problem(model,
                   obj,
                   T,
                   h = h,
                   xl = xl,
                   xu = xu,
                   ul = ul,
                   uu = uu,
                   con = con_contact)

    # Trajectory initialization
    x0 = [x1 for t = 1:T]
    u0 = [1.0e-5 * rand(model.m) for t = 1:T-1] # random controls

    # Pack trajectories into vector
    z0 = pack(x0, u0, prob)

    @time z , info = solve(prob, copy(z0), tol = 1.0e-5, c_tol = 1.0e-5)

    @assert check_slack(z, prob) < 1.0e-4
    x, u = unpack(z, prob)

    return x[end]
end

x = [x1]
for t = 1:200
    push!(x, step_contact(model, x[end], zeros(model.nu), h))
    println("step $t")
end

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model, state_to_configuration(x), Δt = h)


text:
function solve(case, N, K, times)
    deployments = 0
    endtime = 0

    for (S, E) in sort(times)
        if endtime < E
            start = max(S, endtime)
            todeploy = ceil(Int, (E - start) / K)
            deployments += todeploy
            endtime = start + todeploy * K
        end
    end

    println("Case #$case: $deployments")
end

function main()
    T = parse(Int, readline(stdin))

    for case in 1:T
        N, K = parse.(Int, stdin |> readline |> split)

        times = []
        for i in 1:N
            S, E = parse.(Int, stdin |> readline |> split)
            push!(times, (S, E))
        end

        solve(case, N, K, times)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
@testset "chi-square-analysis" begin

    msr_path = joinpath(mktempdir(),"temp.csv")

    data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
    if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
    if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end
    data["settings"]["sbase_default"] = 100.0
    # insert the load profiles
    _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

    # transform data model
    data = _PMD.transform_data_model(data);
    _PMDSE.reduce_single_phase_loadbuses!(data)

    # solve the power flow
    pf_result = _PMD.solve_mc_pf(data, _PMD.ACPUPowerModel, ipopt_solver)

    key = "1" # or "sourcebus"
    v_pu = data["settings"]["vbases_default"][key]* data["settings"]["voltage_scale_factor"] # divider [V] to get the voltage in per units.
    v_max_err = 1.15 # maximum error of voltage measurement = 0.5% or 1.15 V
    σ_v = 1/3*v_max_err/v_pu

    p_pu = data["settings"]["sbase"] # divider [kW] to get the power in per units.
    p_max_err = 0.01 # maximum error of power measurement = 10W, or 0.01 kW
    σ_p = 1/3*p_max_err/p_pu

    # sigma_dict
    σ_dict = Dict("load" => Dict("load" => σ_p,
                                "bus"   => σ_v),
                "gen"  => Dict("gen" => σ_p/100,
                                "bus" => σ_v/100)
                )

    # write measurements based on power flow
    _PMDSE.write_measurements!(_PMD.ACPUPowerModel, data, pf_result, msr_path, exclude = ["vi","vr"], σ = σ_dict)

    # read-in measurement data and set initial values
    _PMDSE.add_measurements!(data, msr_path, actual_meas = true)

    data["se_settings"] = Dict{String,Any}("criterion" => "rwlav", "rescaler" => 1)
    se_result = _PMDSE.solve_acp_red_mc_se(data, ipopt_solver)

    @test _PMDSE.get_degrees_of_freedom(data) == 34

    chi_result = _PMDSE.exceeds_chi_squares_threshold(se_result, data)
    @test chi_result[1] == false
    @test isapprox(chi_result[2], 0.0, atol = 1e-8)
    @test isapprox(chi_result[3], 48.60, atol = 1e-2)

    _PMDSE.add_measurements!(data, msr_path, actual_meas = false)
    se_result = _PMDSE.solve_acp_red_mc_se(data, ipopt_solver)
    chi_result = _PMDSE.exceeds_chi_squares_threshold(se_result, data)
    @test chi_result[1] == true #TODO: there's no real bad data, check whether there is a problem with write_meas
    @test isapprox(chi_result[2], 963.32, atol = 1e-2)
    @test isapprox(chi_result[3], 48.60, atol = 1e-2)

end

@testset "h_functions" begin

    # NB: the length (in terms of lines of code) of this sub-test could/should be significantly result but have no time now

    msr_path = joinpath(mktempdir(),"temp.csv")
    data = _PMD.parse_file(joinpath(BASE_DIR, "test/data/extra/networks/case3_unbalanced.dss"); data_model=MATHEMATICAL)
    #reduce grid
    [delete!(data["load"], l) for (l, load) in data["load"] if l!="1"]
    _PMDSE.reduce_single_phase_loadbuses!(data)
    pf_result = _PMD.solve_mc_pf(data, _PMD.ACPUPowerModel, ipopt_solver)
    _PMDSE.write_measurements!(_PMD.ACPUPowerModel, data, pf_result, msr_path, exclude = ["vr","vi"])
    _PMDSE.add_measurements!(data, msr_path, actual_meas = true)

    # state is the result of the power flow
    vv = pf_result["solution"]["bus"]
    state = [ vv["4"]["vm"][1], vv["4"]["vm"][2], vv["4"]["vm"][3], vv["1"]["vm"][1], vv["1"]["vm"][2], vv["1"]["vm"][3], vv["2"]["vm"][1], vv["2"]["vm"][2], vv["2"]["vm"][3], vv["3"]["vm"][1],
             vv["1"]["va"][1], vv["1"]["va"][2], vv["1"]["va"][3], vv["2"]["va"][1], vv["2"]["va"][2], vv["2"]["va"][3], vv["3"]["va"][1] ]

    variable_dict = _PMDSE.build_variable_dictionary(data)
    state_array = _PMDSE.build_state_array(pf_result, variable_dict)
    @test state_array == state

    # push h functions
    functions = []
    ref_bus = 4
    _PMDSE.add_h_function!(:pd, "4", data, ref_bus, variable_dict, functions)
    _PMDSE.add_h_function!(:qd, "5", data, ref_bus, variable_dict, functions)
    _PMDSE.add_h_function!(:pg, "1", data, ref_bus, variable_dict, functions)
    _PMDSE.add_h_function!(:qg, "2", data, ref_bus, variable_dict, functions)
    _PMDSE.add_h_function!(:vm, "3", data, ref_bus, variable_dict, functions)
    _PMDSE.add_h_function!(:vm, "6", data, ref_bus, variable_dict, functions)

    @test isapprox( functions[1](state), -_DST.mean(data["meas"]["4"]["dst"][1]) ) #pd
    @test isapprox( functions[2](state), -_DST.mean(data["meas"]["5"]["dst"][1]) ) # qd
    @test isapprox( functions[3](state), pf_result["solution"]["gen"]["1"]["pg"][1], atol=1e-8 )
    @test isapprox( functions[4](state), pf_result["solution"]["gen"]["1"]["pg"][2], atol=1e-8 )
    @test isapprox( functions[5](state), pf_result["solution"]["gen"]["1"]["pg"][3], atol=1e-8 )
    @test isapprox( functions[6](state), pf_result["solution"]["gen"]["1"]["qg"][1], atol=1e-8  )
    @test isapprox( functions[7](state), pf_result["solution"]["gen"]["1"]["qg"][2], atol=1e-8  )
    @test isapprox( functions[8](state), pf_result["solution"]["gen"]["1"]["qg"][3], atol=1e-8  )
    @test isapprox( functions[9](state), pf_result["solution"]["bus"]["4"]["vm"][1] )
    @test isapprox( functions[10](state), pf_result["solution"]["bus"]["4"]["vm"][2] )
    @test isapprox( functions[11](state), pf_result["solution"]["bus"]["4"]["vm"][3] )
    @test isapprox( functions[12](state), pf_result["solution"]["bus"]["3"]["vm"][1] )

    _PMDSE.add_measurement!(data, :p, :branch, 1, pf_result["solution"]["branch"]["1"]["pt"], [0.0003])
    _PMDSE.add_measurement!(data, :q, :branch, 3, pf_result["solution"]["branch"]["3"]["qf"], [0.0003, 0.0003, 0.0003])
    _PMDSE.add_measurement!(data, :q, :branch, 1, pf_result["solution"]["branch"]["1"]["qf"], [0.0003])

    _PMDSE.add_h_function!(:p, "7", data, ref_bus, variable_dict, functions)
    _PMDSE.add_h_function!(:q, "8", data, ref_bus, variable_dict, functions)
    _PMDSE.add_h_function!(:q, "9", data, ref_bus, variable_dict, functions)

    cm3 = sqrt.(pf_result["solution"]["branch"]["3"]["pf"].^2+pf_result["solution"]["branch"]["3"]["qf"].^2)./(pf_result["solution"]["bus"]["4"]["vm"])
    cm1 = sqrt(pf_result["solution"]["branch"]["1"]["pf"][1]^2+pf_result["solution"]["branch"]["1"]["qf"][1]^2)/(pf_result["solution"]["bus"]["1"]["vm"][2])

    _PMDSE.add_measurement!(data, :cm, :branch, 3, cm3, [0.0003, 0.0003, 0.0003])
    _PMDSE.add_measurement!(data, :cm, :branch, 1, [cm1], [0.0003])
    _PMDSE.add_measurement!(data, :va, :bus, 2, pf_result["solution"]["bus"]["2"]["va"], [0.0003, 0.0003, 0.0003])
    _PMDSE.add_h_function!(:cm, "10", data, ref_bus, variable_dict, functions)
    _PMDSE.add_h_function!(:cm, "11", data, ref_bus, variable_dict, functions)
    _PMDSE.add_h_function!(:va, "12", data, ref_bus, variable_dict, functions)

    @test isapprox( functions[13](state), pf_result["solution"]["branch"]["1"]["pf"][1] )
    @test isapprox( functions[14](state), pf_result["solution"]["branch"]["3"]["qf"][1], atol=1e-8 )
    @test isapprox( functions[15](state), pf_result["solution"]["branch"]["3"]["qf"][2], atol=1e-8 )
    @test isapprox( functions[16](state), pf_result["solution"]["branch"]["3"]["qf"][3], atol=1e-8 )
    @test isapprox( functions[17](state), pf_result["solution"]["branch"]["1"]["qf"][1] )
    @test isapprox( functions[18](state), cm3[1] , atol=1e-8)
    @test isapprox( functions[19](state), cm3[2] , atol=1e-8)
    @test isapprox( functions[20](state), cm3[3], atol=1e-8 )
    @test isapprox( functions[21](state), cm1 )
    @test isapprox( functions[22](state), pf_result["solution"]["bus"]["2"]["va"][1])
    @test isapprox( functions[23](state), pf_result["solution"]["bus"]["2"]["va"][2])
    @test isapprox( functions[24](state), pf_result["solution"]["bus"]["2"]["va"][3])

    pf_result_ivr = _PMD.solve_mc_pf(data, _PMD.IVRUPowerModel, ipopt_solver)
    _PMDSE.add_measurement!(data, :cr, :branch, 3, pf_result_ivr["solution"]["branch"]["3"]["cr_fr"], [0.0003, 0.0003, 0.0003])
    _PMDSE.add_measurement!(data, :ci, :branch, 3, pf_result_ivr["solution"]["branch"]["3"]["ci_fr"], [0.0003, 0.0003, 0.0003])
    _PMDSE.add_measurement!(data, :crd, :branch, 1, pf_result_ivr["solution"]["load"]["1"]["crd_bus"], [0.0003])
    _PMDSE.add_measurement!(data, :cid, :branch, 1, pf_result_ivr["solution"]["load"]["1"]["cid_bus"], [0.0003])

    _PMDSE.add_h_function!(:cr, "13", data, ref_bus, variable_dict, functions)
    _PMDSE.add_h_function!(:ci, "14", data, ref_bus, variable_dict, functions)
    _PMDSE.add_h_function!(:crd, "15", data, ref_bus, variable_dict, functions)
    _PMDSE.add_h_function!(:cid, "16", data, ref_bus, variable_dict, functions)

    @test isapprox( functions[25](state), pf_result_ivr["solution"]["branch"]["3"]["cr_fr"][1], atol=1e-8)
    @test isapprox( functions[26](state), pf_result_ivr["solution"]["branch"]["3"]["cr_fr"][2], atol=4e-5)
    @test isapprox( functions[27](state), pf_result_ivr["solution"]["branch"]["3"]["cr_fr"][3], atol=1e-8)
    @test isapprox( functions[28](state), pf_result_ivr["solution"]["branch"]["3"]["ci_fr"][1], atol=1e-5)
    @test isapprox( functions[29](state), pf_result_ivr["solution"]["branch"]["3"]["ci_fr"][2], atol=2e-5)
    @test isapprox( functions[30](state), pf_result_ivr["solution"]["branch"]["3"]["ci_fr"][3], atol=2e-5)
    @test isapprox( functions[31](state), -pf_result_ivr["solution"]["load"]["1"]["crd_bus"][1], atol=2e-4)
    @test isapprox( functions[32](state), -pf_result_ivr["solution"]["load"]["1"]["cid_bus"][1], atol=2e-4)

    _PMDSE.add_measurement!(data, :p, :branch, 2, pf_result["solution"]["branch"]["2"]["pf"], [0.0003, 0.0003, 0.0003])
    _PMDSE.add_measurement!(data, :q, :branch, 2, pf_result["solution"]["branch"]["2"]["qf"], [0.0003, 0.0003, 0.0003])
    _PMDSE.add_measurement!(data, :cr, :branch, 2, pf_result_ivr["solution"]["branch"]["2"]["cr_fr"], [0.0003, 0.0003, 0.0003])
    _PMDSE.add_measurement!(data, :ci, :branch, 2, pf_result_ivr["solution"]["branch"]["2"]["ci_fr"], [0.0003, 0.0003, 0.0003])
    cmfr = sqrt.(pf_result_ivr["solution"]["branch"]["2"]["cr_fr"].^2 + pf_result_ivr["solution"]["branch"]["2"]["ci_fr"].^2)
    _PMDSE.add_measurement!(data, :cm, :branch, 2, cmfr, [0.0003, 0.0003, 0.0003])

    _PMDSE.add_h_function!(:p, "17", data, ref_bus, variable_dict, functions)
    _PMDSE.add_h_function!(:q, "18", data, ref_bus, variable_dict, functions)
    _PMDSE.add_h_function!(:cr, "19", data, ref_bus, variable_dict, functions)
    _PMDSE.add_h_function!(:ci, "20", data, ref_bus, variable_dict, functions)
    _PMDSE.add_h_function!(:cm, "21", data, ref_bus, variable_dict, functions)

    @test isapprox( functions[33](state), pf_result["solution"]["branch"]["2"]["pf"][1], atol=1e-8)
    @test isapprox( functions[34](state), pf_result["solution"]["branch"]["2"]["pf"][2], atol=1e-8)
    @test isapprox( functions[35](state), pf_result["solution"]["branch"]["2"]["pf"][3], atol=1e-8)
    @test isapprox( functions[36](state), pf_result["solution"]["branch"]["2"]["qf"][1], atol=1e-8)
    @test isapprox( functions[37](state), pf_result["solution"]["branch"]["2"]["qf"][2], atol=1e-8)
    @test isapprox( functions[38](state), pf_result["solution"]["branch"]["2"]["qf"][3], atol=1e-8)
    @test isapprox( functions[39](state), pf_result_ivr["solution"]["branch"]["2"]["cr_fr"][1], atol=1e-4)
    @test isapprox( functions[40](state), pf_result_ivr["solution"]["branch"]["2"]["cr_fr"][2], atol=1e-4)
    @test isapprox( functions[41](state), pf_result_ivr["solution"]["branch"]["2"]["cr_fr"][3], atol=1e-4)
    @test isapprox( functions[42](state), pf_result_ivr["solution"]["branch"]["2"]["ci_fr"][1], atol=1e-4)
    @test isapprox( functions[43](state), pf_result_ivr["solution"]["branch"]["2"]["ci_fr"][2], atol=1e-4)
    @test isapprox( functions[44](state), pf_result_ivr["solution"]["branch"]["2"]["ci_fr"][3], atol=1e-4)
    @test isapprox( functions[45](state), cmfr[1], atol=1e-6)
    @test isapprox( functions[46](state), cmfr[2], atol=1e-6)
    @test isapprox( functions[47](state), cmfr[3], atol=1e-6)

end

@testset "BadData_matrices_and_LNR" begin

    msr_path = joinpath(mktempdir(),"temp.csv")
    data = _PMD.parse_file(joinpath(BASE_DIR, "test/data/extra/networks/case3_unbalanced.dss"); data_model=MATHEMATICAL)
    #reduce grid
    [delete!(data["load"], l) for (l, load) in data["load"] if l!="1"]
    _PMDSE.reduce_single_phase_loadbuses!(data)
    pf_result = _PMD.solve_mc_pf(data, _PMD.ACPUPowerModel, ipopt_solver)
    _PMDSE.write_measurements!(_PMD.ACPUPowerModel, data, pf_result, msr_path, exclude = ["vr","vi"])
    _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
    _PMDSE.assign_start_to_variables!(data)

    data["se_settings"] = Dict{String,Any}("criterion" => "wls", "rescaler" => 1)
    se_result = _PMDSE.solve_acp_red_mc_se(data, ipopt_solver)

    _PMDSE.add_zib_virtual_meas!(data, 0.00000001, exclude = [2])
    _PMDSE.add_zib_virtual_residuals!(se_result, data)

    variable_dict = _PMDSE.build_variable_dictionary(data)
    h_array = _PMDSE.build_measurement_function_array(data, variable_dict)
    state_array = _PMDSE.build_state_array(pf_result, variable_dict)

    stored_H_matrix = h5open(joinpath(BASE_DIR, "test/data/H_matrix.h5"), "r") do file
        read(file, "H")
    end

    stored_G_matrix = h5open(joinpath(BASE_DIR, "test/data/G_matrix.h5"), "r") do file
        read(file, "G")
    end

    stored_R_matrix = h5open(joinpath(BASE_DIR, "test/data/R_matrix.h5"), "r") do file
        read(file, "R")
    end

    stored_Ω_matrix = h5open(joinpath(BASE_DIR, "test/data/Ω_matrix.h5"), "r") do file
        read(file, "Ω")
    end

    H = _PMDSE.build_H_matrix(h_array, state_array)
    R = _PMDSE.build_R_matrix(data)
    G = _PMDSE.build_G_matrix(stored_H_matrix, R)
    Ω = _PMDSE.build_omega_matrix(R, stored_H_matrix, G)

    @test all(isapprox.(H, stored_H_matrix, atol=1))
    @test all(isapprox.(R, stored_R_matrix, atol=1))
    #@test all(isapprox.(G, stored_G_matrix, atol=1))
    @test all(isapprox.(Ω, stored_Ω_matrix, atol=1))

    id_val, exc = _PMDSE.normalized_residuals(data, se_result, Ω)
    #@test !exc
    #@test id_val[1] == "3"
    #@test isapprox(id_val[2], 0.11035175, atol=1e-8)

    _PMDSE.simple_normalized_residuals(data, se_result, R)
    @test haskey(se_result["solution"]["meas"]["5"], "nr")
end


text:
# These are the official errors as specified in etcd
# https://github.com/coreos/etcd

const etcd_error_codes = Int[(const ECODE_KEY_NOT_FOUND    = 100),
                             (const ECODE_TEST_FAILED      = 101),
                             (const ECODE_NOT_FILE         = 102),
                             (const ECODE_NO_MORE_PEER     = 103),
                             (const ECODE_NOT_DIR          = 104),
                             (const ECODE_NODE_EXIST       = 105),
                             (const ECODE_KEY_IS_PRESERVED = 106),
                             (const ECODE_ROOT_R_ONLY      = 107),
                             (const ECODE_DIR_NOT_EMPTY    = 108),

                             (const ECODE_VALUE_REQUIRED          = 200),
                             (const ECODE_PREV_VALUE_REQUIRED     = 201),
                             (const ECODE_TTL_NA_N                = 202),
                             (const ECODE_INDEX_NA_N              = 203),
                             (const ECODE_VALUE_OR_TTL_REQUIRED   = 204),
                             (const ECODE_TIMEOUT_NA_N            = 205),
                             (const ECODE_NAME_REQUIRED           = 206),
                             (const ECODE_INDEX_OR_VALUE_REQUIRED = 207),
                             (const ECODE_INDEX_VALUE_MUTEX       = 208),
                             (const ECODE_INVALID_FIELD           = 209),

                             (const ECODE_RAFT_INTERNAL = 300),
                             (const ECODE_LEADER_ELECT  = 301),

                             (const ECODE_WATCHER_CLEARED     = 400),
                             (const ECODE_EVENT_INDEX_CLEARED = 401)]

const etcd_errors = {
    # command related errors
    ECODE_KEY_NOT_FOUND => "Key not found",
    ECODE_TEST_FAILED => "Compare failed", # test and set
    ECODE_NOT_FILE => "Not a file",
    ECODE_NO_MORE_PEER => "Reached the max number of peers in the cluster",
    ECODE_NOT_DIR => "Not a directory",
    ECODE_NODE_EXIST => "Key already exists", # create
    ECODE_ROOT_R_ONLY => "Root is read only",
    ECODE_KEY_IS_PRESERVED => "The prefix of given key is a keyword in etcd",
    ECODE_DIR_NOT_EMPTY => "Directory not empty",

    # Post form related errors
    ECODE_VALUE_REQUIRED => "Value is Required in POST form",
    ECODE_PREV_VALUE_REQUIRED => "PrevValue is Required in POST form",
    ECODE_TTL_NA_N => "The given TTL in POST form is not a number",
    ECODE_INDEX_NA_N => "The given index in POST form is not a number",
    ECODE_VALUE_OR_TTL_REQUIRED => "Value or TTL is required in POST form",
    ECODE_TIMEOUT_NA_N => "The given timeout in POST form is not a number",
    ECODE_NAME_REQUIRED => "Name is required in POST form",
    ECODE_INDEX_OR_VALUE_REQUIRED => "Index or value is required",
    ECODE_INDEX_VALUE_MUTEX => "Index and value cannot both be specified",
    ECODE_INVALID_FIELD => "Invalid field",

    # raft related errors
    ECODE_RAFT_INTERNAL => "Raft Internal Error",
    ECODE_LEADER_ELECT => "During Leader Election",

    # etcd related errors
    ECODE_WATCHER_CLEARED => "watcher is cleared due to etcd recovery",
    ECODE_EVENT_INDEX_CLEARED => "The event in requested index is outdated and cleared"
}


text:
@testset "Test Type Constructor" begin
    lin = Linear(2, 3)
    @test lin.ninputs == 2
    @test lin.noutputs == 3
    @test lin.nparams == 6
end

@testset "Test dx Buffer" begin
    lin = Linear(2, 3)
    dx = dx_buffer(lin)
    @test typeof(dx) == Matrix{Float64}
    @test size(dx) == (lin.noutputs, lin.ninputs)
end

@testset "Test dΘ Buffer" begin
    @testset "Test Single Output" begin
        lin = Linear(3, 1)
        dΘ = dΘ_buffer(lin)
        @test typeof(dΘ) == Matrix{Float64}
        @test size(dΘ) == (lin.noutputs, lin.ninputs)
    end
    @testset "Test Multiple Output" begin
        lin = Linear(2, 3)
        dΘ = dΘ_buffer(lin)
        dΘ_expected = hcat((speye(lin.noutputs) for i=1:lin.ninputs)...)
        @test typeof(dΘ) == typeof(dΘ_expected)
        @test dΘ.m  == dΘ_expected.m
        @test dΘ.n  == dΘ_expected.n
        @test dΘ.colptr == dΘ.colptr
        @test dΘ.rowval == dΘ.rowval
    end
end

@testset "Test Linear Evaluation" begin
    lin = Linear(3, 2)
    A = [1 2 3;
         4 5 6]
    z = z_buffer(lin)
    evaluate!(lin, [2, 4, 6], vec(A), z)
    @test z == [28; 64]
    evaluate!(lin, [1, 1, 1], vec(A), z)
    @test z == [6; 15]
    evaluate!(lin, [3, 3, 3], vec(A), z)
    @test z == [18; 45]
end

@testset "Test Jacobian Matrix Computation" begin
    @testset "Test Dense" begin
        lin = Linear(3, 2)
        A = [1 2 3;
             4 5 6]
        z = z_buffer(lin)
        dx = dx_buffer(lin)
        dΘ = dΘ_buffer(lin)
        evaluate!(lin, [2, 3, 4], vec(A), z,
                  dx=dx, dΘ=dΘ)
        @test dx == A
        @test dΘ == [2 0 3 0 4 0;
                     0 2 0 3 0 4]
    end

    @testset "Test Finite Diferences" begin
        lin = Linear(4, 5)
        A = [1.0 2.0 3.1 1.2;
             4.0 5.0 6.1 1.2;
             7.0 8.0 9.1 1.2;
             1.0 2.0 3.1 1.2;
             4.0 5.0 6.1 1.2]
        x = [1.2, 2.3, 4.3, 5]
        Θ = vec(A)

        function cal_x(y)
            z = z_buffer(lin)
            evaluate!(lin, y, Θ, z)
            return z
        end
        function cal_Θ(theta)
            z = z_buffer(lin)
            evaluate!(lin, x, theta, z)
            return z
        end
        dx_numeric = Calculus.finite_difference_jacobian(cal_x, x)
        dΘ_numeric = Calculus.finite_difference_jacobian(cal_Θ, Θ)

        z = z_buffer(lin)
        dx = dx_buffer(lin)
        dΘ = dΘ_buffer(lin)
        evaluate!(lin, x, Θ, z, dx=dx, dΘ=dΘ)
        @test dx ≈ dx_numeric
        @test dΘ ≈ dΘ_numeric
    end

    @testset "Test Single Output" begin
        lin = Linear(6, 1)
        Θ = [1, 2, 3, 4, 5, 6]
        x = [10, 20, 30, 40, 50, 60]

        z = z_buffer(lin)
        dx = dx_buffer(lin)
        dΘ = dΘ_buffer(lin)
        evaluate!(lin, x, Θ, z, dx=dx, dΘ=dΘ)
        @test dx ≈ reshape(Θ, 1, 6)
        @test dΘ ≈ reshape(x, 1, 6)
    end
end

@testset "Test Initial Guess" begin
    N = 1000
    srand(0)
    lin = Linear(5, 2)
    init = ScaledGaussianInitializer()

    mean_estimated = 0
    var_estimated = 0
    for i=1:N
        params = initial_guess(lin, init)
        mean_estimated += sum(params)/(N*lin.nparams)
        var_estimated += sum(params.^2)/(N*lin.nparams)
    end

    std_estimated = sqrt(var_estimated)

    @test std_estimated ≈ 1/sqrt(5) atol=1e-1
    @test mean_estimated ≈ 0.0 atol=1e-1
end


text:
using Test, LinearAlgebra, Random
using GenericLinearAlgebra: cholUnblocked!, cholBlocked!, cholRecursive!

@testset "Cholesky" begin

    Random.seed!(123)

    n = 50

    @testset "element type: $T" for T in (Float32, Float64, BigFloat)
        @testset "is complex: $cplx" for cplx in (true, false)
            if cplx
                A = complex.(T.(rand(n,n)), T.(rand(n,n)))
            else
                A = T.(rand(n,n))
            end
            AcA = A'A
            @test LowerTriangular(cholUnblocked!(copy(AcA), Val{:L}))   ≈ cholesky(Hermitian(AcA)).L
            @test LowerTriangular(cholBlocked!(copy(AcA), Val{:L}, 5))  ≈ cholesky(Hermitian(AcA)).L
            @test LowerTriangular(cholBlocked!(copy(AcA), Val{:L}, 10)) ≈ cholesky(Hermitian(AcA)).L
            @test LowerTriangular(cholRecursive!(copy(AcA), Val{:L}, 1)) ≈ cholesky(Hermitian(AcA)).L
            @test LowerTriangular(cholRecursive!(copy(AcA), Val{:L}, 4)) ≈ cholesky(Hermitian(AcA)).L
        end
    end
end


text:
using JuMP
model = Model();

@variable(model, free_x)

@variable(model, keyword_x, lower_bound = 1, upper_bound = 2)

has_upper_bound(keyword_x)

upper_bound(keyword_x)

@variable(model, a[1:2, 1:2])

n = 10
l = [1; 2; 3; 4; 5; 6; 7; 8; 9; 10]
u = [10; 11; 12; 13; 14; 15; 16; 17; 18; 19]

@variable(model, l[i] <= x[i = 1:n] <= u[i])

@variable(model, y[i = 1:2, j = 1:2] >= 2i + j)

@variable(model, z[i = 2:3, j = 1:2:3] >= 0)

@variable(model, w[1:5,["red", "blue"]] <= 1)

@variable(model, u[i = 1:3, j = i:5])

@variable(model, v[i = 1:9; mod(i, 3) == 0])

@variable(model, integer_z, integer = true)

@variable(model, binary_z, binary = true)

@variable(model, psd_x[1:2, 1:2], PSD)

@variable(model, sym_x[1:2, 1:2], Symmetric)

model = Model()
@variable(model, x)
@variable(model, y)
@variable(model, z[1:10]);

@constraint(model, con, x <= 4)

@constraint(model, [i = 1:3], i * x <= i + 1)

@constraint(model, [i = 1:2, j = 2:3], i * x <= j + 1)

@constraint(model, [i = 1:2, j = 1:2; i != j], i * x <= j + 1)

for i in 1:3
    @constraint(model, 6x + 4y >= 5i)
end

@constraint(model, [i in 1:3], 6x + 4y >= 5i)

@constraint(model, sum(z[i] for i in 1:10) <= 1)

using GLPK

model = Model(GLPK.Optimizer)
@variable(model, x >= 0)
@variable(model, y >= 0)
set_objective_sense(model, MOI.MIN_SENSE)
set_objective_function(model, x + y)

optimize!(model)

@show objective_value(model);

objective_sense(model)

objective_function(model)

objective_function_type(model)

vector_model = Model(GLPK.Optimizer)

A= [ 1 1 9  5;
     3 5 0  8;
     2 0 6 13]

b = [7; 3; 5]

c = [1; 3; 5; 2]

@variable(vector_model, x[1:4] >= 0)
@constraint(vector_model, A * x .== b)
@objective(vector_model, Min, c' * x)

optimize!(vector_model)

@show objective_value(vector_model);


text:
using StorageMirrorServer: RegistryMeta, query_latest_hash, download_and_verify, query_artifacts, read_packages, read_records
using StorageMirrorServer: mirror_tarball
using StorageMirrorServer: decompress
using StorageMirrorServer: timeout_call
using Tar

@testset "mirror_tarball" begin
    tmp_testdir = mktempdir()
    # tmp_testdir = abspath("..", "julia")
    # rm(tmp_testdir; force=true, recursive=true)

    registry = RegistryMeta("General", "23338594-aafe-5451-b93e-139f81909106", "https://github.com/JuliaRegistries/General")
    upstreams = ["https://mirrors.bfsu.edu.cn/julia", "https://mirrors.sjtug.sjtu.edu.cn/julia", "https://us-east.storage.juliahub.com", "https://kr.storage.juliahub.com"]
    registry_hash = query_latest_hash(registry, upstreams)

    tarball = joinpath(tmp_testdir, "registry", registry.uuid, registry_hash)
    resource = "/registry/$(registry.uuid)/$(registry_hash)"
    download_and_verify(upstreams, resource, tarball)
    @test isfile(tarball)

    packages = mktempdir() do tmpdir
        rst = []
        open(tarball, "r") do io
            Tar.extract(decompress(io), tmpdir)
        end
        # only returns packages that are not stored in static_dir
        append!(rst, read_packages(tmpdir; fetch_full_registry=false, static_dir=tmp_testdir, latest_versions_num=1) do pkg
            occursin("MbedTLS", pkg.name) # platform-dependent artifacts
        end)
        append!(rst, read_packages(tmpdir; fetch_full_registry=false, static_dir=tmp_testdir, latest_versions_num=1, recursive=false) do pkg
            occursin("TestImages", pkg.name) || # platform-independent artifacts
            occursin("StanMCMCChain", pkg.name) # vanished resources
        end)
        return rst
    end

    # @time rst_hash = timeout_call(2400) do
        rst_hash = mirror_tarball(registry, upstreams, tmp_testdir; packages=packages, registry_hash=registry_hash)
    # end
    @test rst_hash == registry_hash

    # test if all packages are downloaded
    for pkg in packages
        for (ver, hash_info) in pkg.versions
            tree_hash = hash_info["git-tree-sha1"]
            tarball = joinpath(tmp_testdir, "package", pkg.uuid, tree_hash)
            if pkg.name == "StanMCMCChain"
                # instead this one should be listed in "/failed_resources.txt"
                @test !isfile(tarball)
            else
                @test isfile(tarball)
            end
        end
    end

    # test if all artifacts are downloaded
    artifacts = query_artifacts(tmp_testdir; fetch_full=false)
    @test isempty(artifacts)

    registries_file = joinpath(tmp_testdir, "registries")
    @test isfile(registries_file) && occursin(resource, readline(registries_file))

    # test if vanished resources are listed here
    failed_resources_log = joinpath(tmp_testdir, "failed_resources.txt")
    @test isfile(failed_resources_log)
    open(failed_resources_log, "a") do io
        println(io, "/artifact/uuid/arti") # invalid line
    end
    failed_records = read_records(failed_resources_log)
    for pkg_uuid in [
        "8f1571ae-b3a1-52af-8ab1-32258739efdb", # StanMCMCChain
    ]
        @test any(x->occursin(pkg_uuid, x), failed_records)
    end
    @test all(x->!occursin("/artifact/uuid/arti", x), failed_records)

    # an immediately incremental build should does nothing
    tarball = joinpath(tmp_testdir, "registry", registry.uuid, registry_hash)
    packages = mktempdir() do tmpdir
        rst = []
        open(tarball, "r") do io
            Tar.extract(decompress(io), tmpdir)
        end
        # only returns packages that are not stored in static_dir
        append!(rst, read_packages(tmpdir; fetch_full_registry=false, static_dir=tmp_testdir, latest_versions_num=1) do pkg
            occursin("MbedTLS", pkg.name) # platform-dependent artifacts
        end)
        append!(rst, read_packages(tmpdir; fetch_full_registry=false, static_dir=tmp_testdir, latest_versions_num=1, recursive=false) do pkg
            occursin("TestImages", pkg.name) || # platform-independent artifacts
            occursin("StanMCMCChain", pkg.name) # vanished resources
        end)
        return rst
    end
    @time err_msg = timeout_call(60) do
        @capture_err mirror_tarball(registry, upstreams, tmp_testdir; packages=packages, registry_hash=registry_hash)
    end
    @test occursin("$(length(packages)) previously failed resources are skipped during this build", err_msg)
end


text:
using BinaryBuilder, Pkg

name = "FastME"
version = v"2.1.6"
sources = [
    GitSource("https://gite.lirmm.fr/atgc/FastME.git", "4d0ecf8e1f06bd22e25e57dc362594548c3e10d4"),
]

script = raw"""
cd ${WORKSPACE}/srcdir/FastME

install_license /usr/share/licenses/GPL-3.0+

update_configure_scripts --reconf
autoupdate

# Checks from macros `AC_FUNC_MALLOC` and `AC_FUNC_REALLOC` may fail when cross-compiling,
# which can cause configure to remap `malloc` and `realloc` to replacement functions
# `rpl_malloc` and `rpl_realloc`, which will cause a linking error.  For more information,
# see https://stackoverflow.com/q/70725646/2442087
FLAGS=(ac_cv_func_malloc_0_nonnull=yes ac_cv_func_realloc_0_nonnull=yes)

# optimized multi-threading
./configure --prefix=$prefix --build=${MACHTYPE} --host=${target} "${FLAGS[@]}"
make -j${nproc}
make install
install -Dvm 755 $bindir/fastme${exeext} $bindir/fastmeMP${exeext}

# optimized single-threading
./configure --disable-OpenMP --prefix=$prefix --build=${MACHTYPE} --host=${target} "${FLAGS[@]}"
make clean
make -j${nproc}
make install
"""

platforms = supported_platforms()

products = [
    ExecutableProduct("fastme", :fastme),
    ExecutableProduct("fastmeMP", :fastmeMP),
]

dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"); platforms=filter(!Sys.isbsd, platforms)),
    Dependency(PackageSpec(name="LLVMOpenMP_jll", uuid="1d63c593-3942-5779-bab2-d838dc0a180e"); platforms=filter(Sys.isbsd, platforms)),
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies, julia_compat="1.6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
#!/usr/bin/env julia

# MIT License

# Copyright (c) 2016-2017 Andrew Clemons

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

type Elf
    id::UInt32
    presents::UInt32
end

function initCircle(count)
    elves = Array{Elf, 1}()

    for i in 1:count
        elf = Elf(i, 1)
        push!(elves, elf)
    end

    elves
end

function winningElf(count)
    elves = initCircle(count)

    while length(elves) > 1
        first = shift!(elves)
        second = shift!(elves)

        first.presents = first.presents + second.presents

        push!(elves, first)
    end

    elf = pop!(elves)

    @printf "Elf %d now gets all the presents\n" elf.id
end

function winningElfV2(count)
    elves = initCircle(count)

    len = length(elves)

    offset = 1
    while len > 1
        first = elves[offset]

        opposite = offset + convert(UInt32, round((len / 2), RoundDown))

        if opposite > len
            opposite = opposite - len
        end

        second = elves[opposite]
        deleteat!(elves, opposite)

        first.presents = first.presents + second.presents

        len = length(elves)

        if len > 1
            if opposite > offset
                offset += 1
            else
                if offset > len
                    offset = 1
                end
            end
        end
    end

    elf = pop!(elves)

    @printf "Elf %d now gets all the presents\n" elf.id
end

count = parse(UInt32, chomp(readline()))

winningElf(count)
winningElfV2(count)


text:
module Dec3Part2

slopes = [1 1; 3 1; 5 1; 7 1; 1 2]

function get_level(pattern::Array{String,1}, level::Int64)
    return pattern[level]
end

function get_x_value(level::String, index::Int64)
    n_l = length(level)
    if index <= n_l
        return index
    else
        return index - n_l
    end
end
function get_next_point(current_point::Array{Int64,1}, slope::Array{Int64,1}, pattern::Array{String,1})
    x_ref, y_ref = current_point + slope
    n_level = y_ref
    level = get_level(pattern, n_level)
    x = get_x_value(level, x_ref)
    y = y_ref
    return [x,y]
end

function get_marker_point(pattern::Array{String,1}, current_point::Array{Int64,1})
    row, col = current_point
    return pattern[col][row]
end

function get_n_trees(pattern, slope)
    initial_point = [1,1]
    n_max = length(pattern)
    current_point = initial_point
    counter = 0
    while current_point[2] < n_max
        current_point = get_next_point(current_point, slope, pattern)
        marker = get_marker_point(pattern, current_point)
        if marker == '#'
            counter += 1
        end
    end
    return counter
end

function n_trees_multiple_slopes(pattern::Array{String,1}, slopes::Array{Int64,2})::Int64
    count = 1
    for slope in eachrow(slopes)
        slope = Array(slope)
        count *= get_n_trees(pattern, slope)
    end
    return count
end

function main(file_path::String)
    lines = readlines(file_path)
    return n_trees_multiple_slopes(lines, slopes)
end
end


text:
using Compat
using Base.Test
using PropertyGraph
using UUID

function getnetworksize(graph::Graph, name::String, searchsteps::Integer)
	# start with a specific person...
	# the repeatedly follow the Follows links for searchsteps + 1 traversals (or until there are no more traversals to follow)
	network = query(graph,
					(vertices, v->v.typelabel == "Person" && get(v,"Name") == name),
					newonly, # only process the newly included results (relevant when looping)
					(outgoing, e->e.typelabel == "Follows"), # follow the outgoing Follows links
					head, # the head of these edges is someone being followed
					(mergedistinct,q -> 4), # merge these new results into the results we had at the last loop 4 steps back
					(loop,q-> if q.repeatcount < searchsteps && resultcount(q) > previousresultcount(q, 5) 4 else 0 end) # if we have found new people, loop back and repeat
					)
	# check the expected network size
	networksize = query(network, count)

	return networksize
end

function socialgraphtest()
	graph = buildsocialtestgraph()

	# The social graph in this test is made up of:
	#		- 30 Vertices of type Country, with the following properties:
	#			- Name
	#		- 200 Vertices of type Person, with the following properties:
	#			- Name, Age
	#		- 200 Edges between a Person and a Country, with a type label: "LivesIn" and properties:
	#			- ForYears (indicating how many years the person has lived in the country)
	#		- 1783 Edges between a Person and a Person, with a type label: "Follows" (indicating that one person follows another)

	# count the countries
	@test query(graph, (vertices, v->v.typelabel == "Country"), count) == 30

	# count the people
	@test query(graph, (vertices, v->v.typelabel == "Person"), count) == 200

	# count the LivesIn edges
	@test query(graph, (edges, e->e.typelabel == "LivesIn"), count) == 200

	# count the Follows edges
	@test query(graph, (edges, e->e.typelabel == "Follows"), count) == 1783

	# count the vertices where the condition is never met
	@test query(graph, (vertices, v-> 1 == 2), count) == 0

	# start with the vertices for Germany and France
	# follow the edges of type "LivesIn" to find all the people within the graph who live in Germany or France
	#     only edges that have a ForYears property of >= 5 (indicating the person has lived in the associated country for 5 or more years) will be followed
	# the tail of each edge is a person, filter the people selected to those with an age property >= 30
	part1 = query(graph,
				(vertices, v->v.typelabel == "Country" && (get(v,"Name") == "Germany" || get(v,"Name") == "France")),
			 	(incoming, e->e.typelabel == "LivesIn" && get(e,"ForYears",0) >= 5),
			  	(tail, v->get(v,"Age",0) >= 30))
	# get a count of these results
	@test query(part1, count) == 10

	# start with the results of the first query
	# follow edges of type "Follows" (in an incoming direction)
	# the tail of each edge is a person who follows a person selected in the first query
	# select a distinct set of vertices (as there may have been multiple edges to each person..as a person could follow more than one person selected by the first query)
	part2 = query(part1,
				(incoming, e->e.typelabel == "Follows"),
			  	tail,
				distinct
			 )
	# get a count of the results
	testcount = query(part2, count)
	@test testcount == 78

	# check the expected network size based on Follows links of a person
	# test first just the link of this person
	searchsteps = 0
	networksize = getnetworksize(graph,"Ronny Lamonica", searchsteps)
	@test networksize == 10

	# now test incuding people followed by people Ronny follows
	searchsteps = 1
	networksize = getnetworksize(graph,"Ronny Lamonica", searchsteps)
	@test networksize == 71

	# keep going
	for searchsteps in 2:5
		networksize = getnetworksize(graph,"Ronny Lamonica", searchsteps)
		println("Network size: ", networksize, " Search Steps: ", searchsteps)
	end
end

@time socialgraphtest()


text:
"""
    relative_energy_regain(ground_truth, rec)

Calculates the relative energy regain between the `ground_truth`
and the reconstruction.
Assumes that both arrays are 2 dimensional

# Reference
* [Rainer Heintzmann, \"Estimating missing information by maximum likelihood deconvolution\"](https://www.sciencedirect.com/science/article/abs/pii/S0968432806001272)
"""
function relative_energy_regain(ground_truth, rec)
    T = eltype(ground_truth)
    # go to fourier space
    ground_truth_fft = fft(ground_truth)
    rec_fft = fft(rec)

    # a dict to store the values for certain frequencies
    # we store a list since some (rounded) frequencies occur more than once
    ΔE_R_dict = Dict{T, Vector{T}}()
    E_R_dict = Dict{T, Vector{T}}()

    # round the frequencies to 4 digits, alternative would be to bin
    round4(x) = T(round(x, digits=3))

    # iterate over the frequencies and calculate the relative energy regain
    for (i₂, f₂) in enumerate(fftfreq(size(rec_fft, 2)))
        for (i₁, f₁) in enumerate(fftfreq(size(rec_fft, 1)))
            f_res = round4(√(f₁^2 + f₂^2))
            Δ_E_R = abs2(ground_truth_fft[i₁, i₂] - rec_fft[i₁, i₂])
            E_R = abs2(ground_truth_fft[i₁, i₂])

            update_dict_list!(ΔE_R_dict, f_res, Δ_E_R)
            update_dict_list!(E_R_dict, f_res, E_R)
        end
    end

    # finally transform everything into a list of frequencies and
    # a list of relative energy regains
    freqs = T[]
    G_R_list = T[]
    for f in sort(T.(keys(ΔE_R_dict)))
        push!(freqs, f)
        mean_ΔE_r = mean(ΔE_R_dict[f])
        mean_E_r = mean(E_R_dict[f])
        push!(G_R_list, (mean_E_r - mean_ΔE_r) / mean_E_r)
    end

    return freqs, G_R_list
end

"""
    update_dict_list!(d, k, v)

Updates the dict `d` which stores a list.
If `k` is in the keys of `d` we simply push `v` to the list
otherwise create a new list `[v]`
"""
function update_dict_list!(d, k, v)
    if haskey(d, k)
        push!(d[k], v)
    else
        d[k] = [v]
    end
    return d
end

"""
    normalized_cross_correlation(ground_truth, measured)

Calculates the normalized cross correlation.

External links:
* [Wikipedia](https://en.wikipedia.org/wiki/Sombrero_function)
* [StatsBase.jl](https://juliastats.org/StatsBase.jl/stable/signalcorr/#StatsBase.crosscor)
"""
function normalized_cross_correlation(ground_truth, measured)
    fl(x) = collect(Iterators.flatten(x))
    ground_truth = fl(ground_truth)
    measured = fl(measured)

    ncc = crosscor(ground_truth, measured, [0], demean=true)[begin]
    return ncc
end

"""
    normalized_variance(a, b)

Calculates the mean variance between two array, but normalizing arra a to the same mean as array b.
"""
function normalized_variance(a, b)
    factor = sum(b)/sum(a)
    sum(abs2.(a.*factor .-b))./prod(size(a))
end

function reset_summary!(summary)
    summary["losses"] = []
    summary["best_ncc"] = -Inf
    summary["best_ncc_idx"] = 0
    summary["best_ncc_img"] = []
    summary["nccs"] = []
    summary["best_nvar"] = Inf
    summary["best_nvar_idx"] = 0
    summary["best_nvar_img"] = []
    summary["times"] = []
    summary["step_sizes"] = []
    summary["nvars"] = []
end

"""
    options_trace_deconv(ground_truth, iterations, mapping, every=1; more_options...)

    A useful routine to simplify performance checks of deconvolution on simulated data.
    Returns an Options structure to be used with the deconvolution routine as an argument to `opt_options` and
    a summary dictionary with all the performance metrics calculated, which is resetted and updated during deconvolution.
    This can then be plotted or visualized.
    The summary dictionary has the following content:

    "best_nvar"     => the lowest normalized variance compared to the `ground_truth` that was achieved.
    "best_nvar_img" => the reconstruction result corresponding to this lowest normalized variance
    "best_nvar_idx" => the corresponding index where this was achieved. `(best_nvar_idx-1)*every+1` approximated the iteration number.
    "best_ncc"      => the highest normalized crosscorrelation compared to the `ground_truth` that was achieved.
    "best_ncc_img"  => the reconstruction result corresponding to this highest normalized crosscorrelation
    "losses"        => the vector of losses evaluated at each of `every` iterations.
    "nccs"          => the vector of normalized cross correlations calculated at each of `every` iterations.
    "best_ncc_idx"  => the corresponding index where this was achieved. `(best_ncc_idx-1)*every+1` approximated the iteration number.
    "nvars"         => the vector of normalized variances calculated at each of `every` iterations.

    For an example of how to plot the results, see the file `` in the `examples` folder.
# Arguments
- `ground_truth`: The underlying ground truth data. Note that this scaling is unimportant due to the normalized norms used for comparison,
                  whereas the relative offset matters.
- `iterations`: The maximal number of iterations to performance. If covergence is reached, the result may have less iterations
- `mapping`: If mappings such as the positivity contraints (e.g. `NonNegative()`) are used in the deconvolution routing, they also
             need to be provided here. Otherwises select `nothing`.
- `every`: This option allows to select every how many iterations the evaluation is performed. Note that the results will not keep track
        of this iteration number.

the argument list can be followed by a semicolon and any number of named arguments which will be passed to the option structure.
E.g. `opt_noreg, show_noreg = options_trace_deconv(ground_truth, iterations, mapping; x_tol=0.001, f_tol=0.001, f_calls_limit=100);`

# Example
```julia-repl
julia> using DeconvOptim, TestImages, Noise, Plots;

julia> obj = Float32.(testimage("resolution_test_512"));

julia> psf = Float32.(generate_psf(size(obj), 30));

julia> img_b = conv(obj, psf);

julia> img_n = poisson(img_b, 300);

julia> iterations = 100;

julia> opt_noreg, show_noreg = options_trace_deconv(obj, iterations, Non_negative());

julia> res_noreg, o = deconvolution(img_n, psf, regularizer = nothing, opt_options=opt_noreg);

julia> opt_GR, show_GR = options_trace_deconv(obj, iterations, Non_negative());

julia> res_GR, o = deconvolution(img_n, psf, λ=1e-2, regularizer=DeconvOptim.GR(), opt_options=opt_GR);

julia> opt_TV, show_TV = options_trace_deconv(obj, iterations, Non_negative());

julia> res_TV, o = deconvolution(img_n, psf, λ=1e-3, regularizer=DeconvOptim.TV(), opt_options=opt_TV);

julia> plot()

julia> show_noreg(false,"NoReg")

julia> show_GR(false,"GR")

julia> show_TV(false,"TV")

julia> using View5D

julia> @vt (ground_truth, best_ncc_img, best_nvar_img) = show_noreg(true)
```
"""
function options_trace_deconv(ground_truth, iterations, mapping, every=1; more_options...)
    @show more_options
    summary = Dict()
    reset_summary!(summary)
    summary["ground_truth"] = ground_truth # needs to be accesible
    idx = 1
    cb = tr -> begin
        img = (mapping === nothing) ? tr[end].metadata["x"] : mapping[1](tr[end].metadata["x"])
        img *= mean(summary["ground_truth"])
        record_progress!(summary, img, idx, tr[end].value,
                        tr[end].metadata["time"], tr[end].metadata["Current step size"])
        idx += 1
        false
    end

    opt_options = Optim.Options(callback = cb, iterations=iterations, show_every=every, store_trace=true, extended_trace=true; more_options...)
    return (opt_options, summary)
end

"""
    record_progress!(summary, img, idx, loss, mytime, stepsize)
    helper function for recording the iteration progress in a summary dictionary.
"""
function record_progress!(summary, img, idx, loss, mytime, stepsize)
    # iteration always starts with index 0 (before 1st iteration)
    if idx == 0
        reset_summary!(summary)
    end
    push!(summary["losses"], loss)
    push!(summary["times"], mytime)
    push!(summary["step_sizes"], stepsize)
    # current image:
    # the line below is needed, since in the iterations, the measurement is rescaled to a mean of one.
    # see deconvolution.jl.  This rescaling is only an estimate and does not affect the norms.

    ground_truth = summary["ground_truth"]
    ncc = DeconvOptim.normalized_cross_correlation(ground_truth, img)
    push!(summary["nccs"], ncc)
    summary["best_ncc"], summary["best_ncc_img"], summary["best_ncc_idx"] = let
        if ncc > summary["best_ncc"]
            (ncc, img, idx)
        else
            (summary["best_ncc"], summary["best_ncc_img"], summary["best_ncc_idx"])
        end
    end
    nvar = normalized_variance(img, ground_truth)
    push!(summary["nvars"], nvar)
    summary["best_nvar"], summary["best_nvar_img"], summary["best_nvar_idx"] = let
        if nvar < summary["best_nvar"]
                (nvar, img, idx)
        else
                (summary["best_nvar"], summary["best_nvar_img"], summary["best_nvar_idx"])
        end
    end
end


text:
module TestQuaternions

using Test
using LinearAlgebra
using StaticArrays
using Flight.Quaternions

export test_quaternions

function test_quaternions()
    @testset verbose = true "Quaternions" begin
        @testset verbose = true "FreeQuat" begin test_Quat() end
        @testset verbose = true "UnitQuat" begin test_UnitQuat() end
    end
end

function test_Quat()

    @testset "Basics" begin

        #use inner constructor only, do not rely on outer constructors
        data = SVector{4,Float64}([-3, 2, 4, 1])
        q = FreeQuat(data)

        #getindex
        @test q[:] == data[:]

        #getproperty
        @test q.real == data[1]
        @test q.imag == data[2:end]

        #copy
        q_copy = copy(q)
        @test isa(q_copy, FreeQuat)

        #normalization
        q_copy = copy(q)
        @test norm(normalize(q)) ≈ 1 #returns a normalized copy of q

    end

    @testset "Constructors" begin

        v = [2, 3, 1, -4]
        q = FreeQuat(v)

        @test q[:] == v[:] #Vector
        @test FreeQuat(q)[:] == q[:]
        @test FreeQuat(2.0)[:] == [2.0, 0, 0, 0] #scalar
        @test FreeQuat()[:] == zeros(4) #zero arguments
        @test FreeQuat(real = 2, imag = [4,3,-2])[:] == [2, 4, 3, -2] #keyword
        @test FreeQuat(real = 2)[:] == [2, 0, 0, 0] #keyword
        @test FreeQuat(imag = [3,4,5])[:] == [0, 3, 4, 5] #keyword
        @test_throws DimensionMismatch FreeQuat([1, 3, 2, 5, 0])
        @test_throws DimensionMismatch FreeQuat([1, 5, 0])

    end

    @testset "Operators" begin

        v1 = rand(4); v2 = rand(4); v3 = rand(4)
        q1 = FreeQuat(v1); q2 = FreeQuat(v2); q3 = FreeQuat(v3)

        #equality (inherited from AbstractVector)
        # @test q1 != v1 #not equal to an AbstractVector with the same components
        @test q1 != FreeQuat(-v1)
        @test q1 == q1

        #approximate equality
        @test q1 ≈ FreeQuat(v1 + [1e-9, 0, 0, 0])
        @test !(q1 ≈ q2)

        #promotion and conversion from real scalar and AbstractVector
        @test promote(q1, 3.213)[2] == FreeQuat(3.213)
        @test convert(FreeQuat, [1, 2, 3, 4]) == FreeQuat([1,2,3,4])

        #conjugation
        @test q1'.imag ≈ -q1.imag

        #unary minus
        @test isa(-q1, FreeQuat)
        @test (-q1)[:] == -(q1[:])

        #addition and substraction
        @test isa(q1 + q2, FreeQuat)
        @test (q1 + q2)[:] == v1 + v2
        @test isa(q1 - q2, FreeQuat)
        @test (q1 - q2)[:] == v1 - v2

        #multiplication, inverse and division by quaternion
        @test isa(q1 * q2, FreeQuat)
        @test q1 * FreeQuat(1.0) ≈ q1
        @test (q1 * q2) * q3 ≈ (q1 * (q2 * q3))
        @test (q1 + q2) * q3 ≈ (q1 * q3 + q2 * q3)
        @test q1 * inv(q1) ≈ FreeQuat(1.0)
        @test q1 / q2 ≈ q1 * inv(q2)
        @test q1 \ q2 ≈ inv(q1) * q2
        @test q1 \ q2 != q2 / q1 #inv(q1)*q2 != q2*inv(q1)

        #and by real
        @test 5q1 ≈ FreeQuat(5*q1[:])
        @test q1/5 ≈ FreeQuat(q1[:]/5)
        @test 5/q1 ≈ FreeQuat(5.) * inv(q1)
        @test 5\q1 ≈ inv(FreeQuat(5.)) * q1
    end

end

function test_UnitQuat()

    @testset "Basics" begin

        q = FreeQuat([-3, 2, 4, 1])
        u_badnorm = UnitQuat(q, normalization = false)
        u = UnitQuat(q)

        #built-in normalization
        @test abs(norm(u_badnorm) - 1) > 0.5 #check normalization has been bypassed
        @test norm(u) ≈ 1

        #getindex
        @test u_badnorm[:] == q[:]

        #getproperty
        @test u_badnorm.real == q[1]
        @test u_badnorm.imag == q[2:end]

        #setindex & setproperty disallowed
        @test_throws ErrorException u[4] = -5
        @test_throws ErrorException u.real = 0
        @test_throws ErrorException u.imag = [1, 2, 3]

        #copy & normalization
        u_badnorm_copy = copy(u_badnorm)
        @test isa(u_badnorm_copy, UnitQuat)
        @test norm(normalize(u_badnorm)) ≈ 1 #returns a normalized copy of u

    end

    @testset "Constructors" begin

        v = [2.0, 3, 1, -4]
        u = UnitQuat(FreeQuat(v)) #inner constructor
        q = FreeQuat(v)

        #outer constructors
        @test u[:] == normalize(v)[:] #Vector
        @test UnitQuat(q)[:] == normalize(q)[:] #AbstractQuat
        @test UnitQuat(2.0)[:] == [1, 0, 0, 0] #scalar
        @test UnitQuat(real = v[1], imag = v[2:4])[:] == normalize(v)[:]

    end

    @testset "Operators" begin

        v1 = normalize(rand(4)); v2 = normalize(rand(4))
        u1 = UnitQuat(v1); u2 = UnitQuat(v2)

        #equality (inherited from AbstractVector)
        # @test u1 != v1 #not equal to an AbstractVector with the same components
        @test u1 != UnitQuat(-v1)
        @test u1 == UnitQuat(v1)
        @test u1 == FreeQuat(u1) #equal to a FreeQuat with the same components

        #approximate equality
        @test u1 ≈ UnitQuat(FreeQuat(v1 + [1e-9, 0, 0, 0]), normalization = false)
        @test !(u1 ≈ u2)

        #promotion and conversion to FreeQuat
        @test promote(u1, FreeQuat())[1] == FreeQuat(u1)

        #conjugation
        @test u1'.imag ≈ -u1.imag

        #unary minus
        @test isa(-u1, UnitQuat)
        @test (-u1)[:] == -(u1[:])

        #multiplication, inverse and division by UnitQuat
        @test u1 * UnitQuat(1) ≈ u1
        @test u1' == inv(u1)
        @test u1 * u1' ≈ UnitQuat(1)
        @test u1 / u2 ≈ u1 * u2'
        @test u1 \ u2 ≈ u1' * u2
        @test u1 \ u2 != u2 / u1 #inv(u1)*u2 != u2*inv(u1)

        #operations with Reals and FreeQuats
        q = FreeQuat(rand(4))
        @test u1 * FreeQuat(1.0) ≈ u1
        @test q * u1' ≈ q / FreeQuat(u1)
        @test u1 / q ≈ FreeQuat(u1) / q
        @test u1 \ q ≈ FreeQuat(u1) \ q

        @test u1 * q isa FreeQuat
        @test q * u1' isa FreeQuat
        @test u1 / q isa FreeQuat
        @test u1 \ q isa FreeQuat

        @test 5u1 isa FreeQuat
        @test 5+u1 isa FreeQuat
        @test 2.0-u1 isa FreeQuat
        @test 2/u1 isa FreeQuat
    end

end

end #module


text:
# using Stars
# using Test

@testset "st_as_sf" begin
    file_vi = "data/MOD13A2_Henan_2015_2020_10km.tif"
    ga = rast(file_vi)
    mask = ga[1] != 0
    # mask = get_mask(ga, 0)
    # mask_shrink = ga.A[:, :, 1] .!= 0

    @test size(st_as_sf(ga, ga[1] != 0), 1) == 670
    @test size(st_as_sf(ga, !isnan(ga[1])), 1) == 1073

    EVI = st_as_sf(file_vi, mask)

    ## post processing
    # r_mask = rast(mask, b)
    lon, lat = st_dim(mask)
    d_coord = maskCoords(mask)
    @test length(d_coord[:, 1]) == 670
end

@testset "st_as_stars" begin
    function test_rastDiff(r1, r2)
        r_diff = (r1 - r2) |> maximum
        delta = filter(!isnan, r_diff.A) |> maximum
        @test delta <= 1e-9
    end

    f = "data/MOD13A2_Henan_2015_2020_10km.tif"
    r = rast(f)
    r_mask = get_mask(r, 0)
    EVI = st_as_sf(r, r_mask)

    # r1 = st_as_stars(EVI, mask, b)           # bbox works
    # r2 = st_as_stars(EVI, mask, bbox2vec(b)) # bbox in vector format also works
    r3 = st_as_stars(EVI, r_mask)

    # test_rastDiff(r1, r2)
    test_rastDiff(r, r3)
end

# @test r1.A == r3.A
# @test r2.A == r3.A
# @test r2.A == r1.A
# @time r = st_as_stars(EVI, mask, b; outfile = "MOD13A2_Henan_2015_2020.tif");
# @time plot(r[:,:, 1])

# heatmap(mat')
# @time d2 = unstack(d_coord, :lon, :lat, :id); "ok"
# @time Matrix(d2)


text:
# Test framework to see if config files and fixed/free parameters work

using DiskJockey.model

import YAML
config = YAML.load(open("config.yaml"))

# Read parameters from a config file
pars = config["parameters"]
model = config["model"]

# Test the dictionary conversion routine
pars_type = convert_dict(pars, model)
println(pars_type)
println(pars_type.Sigma_c)

# Test the vector conversion routine

# Convert pars dictionary to dictionary of symbols
xargs = convert(Dict{Symbol}{Float64}, pars)

# According to model.jl, these are the total parameters. Let's fix some
# "standard", ["M_star", "r_c", "T_10", "q", "gamma", "Sigma_c", "ksi", "dpc", "incl", "PA", "vel", "mu_RA", "mu_DEC"]

# Separate out the fixed values
fixed = ["gamma"]

# Separate out the free values
free = ["M_star", "r_c", "T_10", "q", "Sigma_c", "ksi", "dpc", "incl", "PA", "vel", "mu_RA", "mu_DEC"]

p = [1.0, 50.0, 40.0, 0.6, -1.0, 0.4, 145.0, 130.0, 130.0, 0.0, 0.0, 0.0]

# Convert a proposal (vector) of free parameters with to a parameter type using fixed parameters to supply missing fields.
# Convert these to a parameter type
pars_type = convert_vector(p, model, fixed; xargs...)

println(pars_type)
println(pars_type.Sigma_c)

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
function simulate!(X, mcProcess::HestonProcess, rfCurve::ZeroRateCurve, mcBaseData::MonteCarloConfiguration{type1, type2, type3, SerialMode, type5}, T::numb) where {numb <: Number, type1 <: Integer, type2 <: Integer, type3 <: StandardMC, type5 <: Random.AbstractRNG}
    r = rfCurve.r
    S0 = mcProcess.underlying.S0
    d = dividend(mcProcess)
    Nsim = mcBaseData.Nsim
    Nstep = mcBaseData.Nstep
    σ = mcProcess.σ
    σ_zero = mcProcess.σ_zero
    λ1 = mcProcess.λ
    κ = mcProcess.κ
    ρ = mcProcess.ρ
    θ = mcProcess.θ
    @assert T > 0.0

    ####Simulation
    ## Simulate
    κ_s = κ + λ1
    θ_s = κ * θ / (κ + λ1)
    r_d = r - d
    dt = T / Nstep
    zero_drift = r_d(dt * 0.0, dt)
    isDualZero = S0 * T * σ_zero * κ * θ * λ1 * σ * ρ * 0.0 * zero_drift
    #X=Matrix{typeof(isDualZero)}(undef,Nsim,Nstep+1);
    X[:, 1] .= isDualZero
    v_m = [σ_zero^2 + isDualZero for _ = 1:Nsim]
    for j = 1:Nstep
        e1 = randn(mcBaseData.rng, Nsim)
        e2 = randn(mcBaseData.rng, Nsim)
        e2 = e1 * ρ + e2 * sqrt(1 - ρ * ρ)
        tmp_ = r_d((j - 1) * dt, dt)
        X[:, j+1] = X[:, j] + (tmp_ .- 0.5 * max.(v_m, 0) * dt) + sqrt.(max.(v_m, 0)) * sqrt(dt) .* e1
        v_m += κ_s .* (θ_s .- max.(v_m, 0)) .* dt + σ .* sqrt.(max.(v_m, 0)) .* sqrt(dt) .* e2
    end
    ## Conclude
    X .= S0 .* exp.(X)
    return
end

function simulate!(X, mcProcess::HestonProcess, rfCurve::ZeroRateCurve, mcBaseData::MonteCarloConfiguration{type1, type2, type3, SerialMode, type5}, T::numb) where {numb <: Number, type1 <: Integer, type2 <: Integer, type3 <: AntitheticMC, type5 <: Random.AbstractRNG}
    r = rfCurve.r
    S0 = mcProcess.underlying.S0
    d = dividend(mcProcess)
    Nsim = mcBaseData.Nsim
    Nstep = mcBaseData.Nstep
    σ = mcProcess.σ
    σ_zero = mcProcess.σ_zero
    λ1 = mcProcess.λ
    κ = mcProcess.κ
    ρ = mcProcess.ρ
    θ = mcProcess.θ
    @assert T > 0.0

    ####Simulation
    ## Simulate
    κ_s = κ + λ1
    θ_s = κ * θ / (κ + λ1)

    dt = T / Nstep
    isDualZero = S0 * T * σ_zero * κ * θ * λ1 * σ * ρ * 0.0
    #X=Matrix{typeof(isDualZero)}(undef,Nsim,Nstep+1);
    X[:, 1] .= isDualZero
    v_m = [σ_zero^2 + isDualZero for _ = 1:Nsim]
    Nsim_2 = div(Nsim, 2)
    r_d = r - d

    for j = 1:Nstep
        e1 = randn(mcBaseData.rng, Nsim_2)
        e2 = randn(mcBaseData.rng, Nsim_2)
        e1 = [e1; -e1]
        e2 = [e2; -e2]
        e2 = e1 .* ρ .+ e2 .* sqrt(1 - ρ * ρ)
        tmp_ = r_d((j - 1) * dt, dt)
        X[:, j+1] = X[:, j] + (tmp_ .- 0.5 .* max.(v_m, 0) .* dt) + sqrt.(max.(v_m, 0)) .* sqrt(dt) .* e1
        v_m += κ_s .* (θ_s .- max.(v_m, 0)) .* dt + (σ * sqrt(dt)) .* sqrt.(max.(v_m, 0)) .* e2
    end

    ## Conclude
    X .= S0 .* exp.(X)
    return
end


text:
"""
A test for smaller_vector().

Arguments
-----------
input:Array
    input vector

i:Integer
    index of the item to be removed

solution:Array
    answer key

msg="":AbstractString
    (optional) test message
"""
function test_smaller_vector(input::Array, i::Integer, solution::Array, msg::AbstractString="")
    println("-------------------------------------")
    print_with_color(:blue, "Test smaller_vector()\n\n")
    if msg != ""
        print_with_color(:blue, "$(msg)\n\n")
    end
    answer = smaller_vector(input, i)
    @test answer == solution
    print_with_color(:green, "VERIFIED!\n")
    println("-------------------------------------")
    return true
end
test_smaller_vector([1,2,3],1,[2,3],"test 1")
test_smaller_vector([1,2,3],2,[1,3],"test 2")
test_smaller_vector([1,2,3],3,[1,2],"test 3")
test_smaller_vector([1,2,3],4,[1,2,3],"test 4")
test_smaller_vector([1,2,3],0,[1,2,3],"test 5")


text:
type DArray{T,N,A} <: AbstractArray{T,N}
    dims::NTuple{N,Int}

    chunks::Array{RemoteRef,N}

    # pmap[i]==p ⇒ processor p has piece i
    pmap::Vector{Int}

    # indexes held by piece i
    indexes::Array{NTuple{N,UnitRange{Int}},N}
    # cuts[d][i] = first index of chunk i in dimension d
    cuts::Vector{Vector{Int}}

    function DArray(dims, chunks, pmap, indexes, cuts)
        # check invariants
        assert(size(chunks) == size(indexes))
        assert(length(chunks) == length(pmap))
        assert(dims == map(last,last(indexes)))
        new(dims, chunks, pmap, indexes, cuts)
    end
end

typealias SubDArray{T,N,D<:DArray} SubArray{T,N,D}
typealias SubOrDArray{T,N}         Union(DArray{T,N}, SubDArray{T,N})

## core constructors ##

# dist == size(chunks)
function DArray(init, dims, procs, dist)
    np = prod(dist)
    procs = procs[1:np]
    idxs, cuts = chunk_idxs([dims...], dist)
    chunks = Array(RemoteRef, dist...)
    for i = 1:np
        chunks[i] = remotecall(procs[i], init, idxs[i])
    end
    p = max(1, localpartindex(procs))
    A = remotecall_fetch(procs[p], r->typeof(fetch(r)), chunks[p])
    DArray{eltype(A),length(dims),A}(dims, chunks, procs, idxs, cuts)
end

function DArray(init, dims, procs)
    if isempty(procs)
        error("no processors")
    end
    DArray(init, dims, procs, defaultdist(dims,procs))
end
DArray(init, dims) = DArray(init, dims, workers()[1:min(nworkers(),maximum(dims))])

# new DArray similar to an existing one
DArray(init, d::DArray) = DArray(init, size(d), procs(d), [size(d.chunks)...])

similar(d::DArray, T, dims::Dims)= DArray(I->Array(T, map(length,I)), dims, procs(d))
similar(d::DArray, T)= similar(d, T, size(d))
similar{T}(d::DArray{T}, dims::Dims)= similar(d, T, dims)
similar{T}(d::DArray{T})= similar(d, T, size(d))

size(d::DArray) = d.dims
procs(d::DArray) = d.pmap

chunktype{T,N,A}(d::DArray{T,N,A}) = A

## chunk index utilities ##

# decide how to divide each dimension
# returns size of chunks array
function defaultdist(dims, procs)
    dims = [dims...]
    chunks = ones(Int, length(dims))
    np = length(procs)
    f = sort!(collect(keys(factor(np))), rev=true)
    k = 1
    while np > 1
        # repeatedly allocate largest factor to largest dim
        if np%f[k] != 0
            k += 1
            if k > length(f)
                break
            end
        end
        fac = f[k]
        (d, dno) = findmax(dims)
        # resolve ties to highest dim
        dno = last(find(dims .== d))
        if dims[dno] >= fac
            dims[dno] = div(dims[dno], fac)
            chunks[dno] *= fac
        end
        np = div(np,fac)
    end
    chunks
end

# get array of start indexes for dividing sz into nc chunks
function defaultdist(sz::Int, nc::Int)
    if sz >= nc
        iround(linspace(1, sz+1, nc+1))
    else
        [[1:(sz+1)], zeros(Int, nc-sz)]
    end
end

# compute indexes array for dividing dims into chunks
function chunk_idxs(dims, chunks)
    cuts = map(defaultdist, dims, chunks)
    n = length(dims)
    idxs = Array(NTuple{n,UnitRange{Int}},chunks...)
    cartesianmap(tuple(chunks...)) do cidx...
        idxs[cidx...] = ntuple(n, i->(cuts[i][cidx[i]]:cuts[i][cidx[i]+1]-1))
    end
    idxs, cuts
end

function localpartindex(pmap::Vector{Int})
    mi = myid()
    for i = 1:length(pmap)
        if pmap[i] == mi
            return i
        end
    end
    return 0
end

localpartindex(d::DArray) = localpartindex(d.pmap)

function localpart{T,N,A}(d::DArray{T,N,A})
    lpidx = localpartindex(d)
    if lpidx == 0
        convert(A, Array(T, ntuple(N,i->0)))::A
    else
        fetch(d.chunks[lpidx])::A
    end
end
function localindexes(d::DArray)
    lpidx = localpartindex(d)
    if lpidx == 0
        ntuple(ndims(d), i->1:0)
    else
        d.indexes[lpidx]
    end
end

# find which piece holds index (I...)
function locate(d::DArray, I::Int...)
    ntuple(ndims(d), i->searchsortedlast(d.cuts[i], I[i]))
end

chunk{T,N,A}(d::DArray{T,N,A}, i...) = fetch(d.chunks[i...])::A

## convenience constructors ##

dzeros(args...) = DArray(I->zeros(map(length,I)), args...)
dzeros(d::Int...) = dzeros(d)
dones(args...) = DArray(I->ones(map(length,I)), args...)
dones(d::Int...) = dones(d)
dfill(v, args...) = DArray(I->fill(v, map(length,I)), args...)
dfill(v, d::Int...) = dfill(v, d)
drand(args...)  = DArray(I->rand(map(length,I)), args...)
drand(d::Int...) = drand(d)
drandn(args...) = DArray(I->randn(map(length,I)), args...)
drandn(d::Int...) = drandn(d)

## conversions ##

function distribute(a::AbstractArray)
    owner = myid()
    rr = RemoteRef()
    put!(rr, a)
    DArray(size(a)) do I
        remotecall_fetch(owner, ()->fetch(rr)[I...])
    end
end

function convert{S,T,N}(::Type{Array{S,N}}, d::DArray{T,N})
    a = Array(S, size(d))
    @sync begin
        for i = 1:length(d.chunks)
            @async a[d.indexes[i]...] = chunk(d, i)
        end
    end
    a
end

function convert{S,T,N}(::Type{Array{S,N}}, s::SubDArray{T,N})
    I = s.indexes
    d = s.parent
    if isa(I,(UnitRange{Int}...)) && S<:T && T<:S
        l = locate(d, map(first, I)...)
        if isequal(d.indexes[l...], I)
            # SubDArray corresponds to a chunk
            return chunk(d, l...)
        end
    end
    a = Array(S, size(s))
    a[[1:size(a,i) for i=1:N]...] = s
    a
end

function reshape{T,S<:Array}(A::DArray{T,1,S}, d::Dims)
    if prod(d) != length(A)
        error("dimensions must be consistent with array size")
    end
    DArray(d) do I
        sz = map(length,I)
        d1offs = first(I[1])
        nd = length(I)

        B = Array(T,sz)
        nr = size(B,1)
        sztail = size(B)[2:end]

        for i=1:div(length(B),nr)
            i2 = ind2sub(sztail, i)
            globalidx = [ I[j][i2[j-1]] for j=2:nd ]

            a = sub2ind(d, d1offs, globalidx...)

            B[:,i] = A[a:(a+nr-1)]
        end
        B
    end
end

## indexing ##

function getindex(r::RemoteRef, args...)
    if r.where==myid()
        getindex(fetch(r), args...)
    else
        remotecall_fetch(r.where, getindex, r, args...)
    end
end

getindex(d::DArray, i::Int) = getindex(d, ind2sub(size(d), i))
getindex(d::DArray, i::Int...) = getindex(d, sub2ind(size(d), i...))

function getindex{T}(d::DArray{T}, I::(Int...))
    chidx = locate(d, I...)
    chunk = d.chunks[chidx...]
    idxs = d.indexes[chidx...]
    localidx = ntuple(ndims(d), i->(I[i]-first(idxs[i])+1))
    chunk[localidx...]::T
end

getindex(d::DArray) = d[1]
getindex(d::DArray, I::Union(Int,UnitRange{Int})...) = sub(d,I)

copy(d::SubOrDArray) = d

# local copies are obtained by convert(Array, ) or assigning from
# a SubDArray to a local Array.

function setindex!(a::Array, d::DArray, I::UnitRange{Int}...)
    n = length(I)
    @sync begin
        for i = 1:length(d.chunks)
            K = d.indexes[i]
            @async a[[I[j][K[j]] for j=1:n]...] = chunk(d, i)
        end
    end
    a
end

function setindex!(a::Array, s::SubDArray, I::UnitRange{Int}...)
    n = length(I)
    d = s.parent
    J = s.indexes
    if length(J) < n
        a[I...] = convert(Array,s)
        return a
    end
    offs = [isa(J[i],Int) ? J[i]-1 : first(J[i])-1 for i=1:n]
    @sync begin
        for i = 1:length(d.chunks)
            K_c = {d.indexes[i]...}
            K = [ intersect(J[j],K_c[j]) for j=1:n ]
            if !any(isempty, K)
                idxs = [ I[j][K[j]-offs[j]] for j=1:n ]
                if isequal(K, K_c)
                    # whole chunk
                    @async a[idxs...] = chunk(d, i)
                else
                    # partial chunk
                    ch = d.chunks[i]
                    @async a[idxs...] = remotecall_fetch(ch.where, ()->sub(fetch(ch), [K[j]-first(K_c[j])+1 for j=1:n]...))
                end
            end
        end
    end
    a
end

# to disambiguate
setindex!(a::Array{Any}, d::SubOrDArray, i::Int) = arrayset(a, d, i)

setindex!(a::Array, d::SubOrDArray, I::Union(Int,UnitRange{Int})...) =
    setindex!(a, d, [isa(i,Int) ? (i:i) : i for i in I ]...)

## higher-order functions ##

map(f::Callable, d::DArray) = DArray(I->map(f, localpart(d)), d)

reduce(f::Function, d::DArray) =
    mapreduce(fetch, f,
              { @spawnat p reduce(f, localpart(d)) for p in procs(d) })

function map!(f::Callable, d::DArray)
    @sync begin
        for p in procs(d)
            @spawnat p map!(f, localpart(d))
        end
    end
end

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia