text:
# ---------------------------------------------------------------------------- #
#
#   export.jl
#
#   List of functions and constants that are exposed outside of the package
#
#   egads.jl
#   Spring 2018
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

# functions
export EG_alloc, EG_calloc, EG_reall, EG_free, EG_revision, EG_open,
    EG_loadModel, EG_saveModel, EG_deleteObject, EG_makeTransform,
    EG_getTransformation, EG_getContext, EG_setOutLevel, EG_getInfo,
    EG_copyObject, EG_flipObject, EG_close, EG_setUserPointer,
    EG_getUserPointer, EG_attributeAdd, EG_attributeDel, EG_attributeNum,
    EG_attributeGet, EG_attributeRet, EG_attributeDup, EG_getGeometry,
    EG_makeGeometry, EG_getRange, EG_evaluate, EG_invEvaluate,
    EG_invEvaluateGuess, EG_arcLength, EG_curvature, EG_approximate,
    EG_fitTriangles, EG_otherCurve, EG_isSame, EG_isoCline,
    EG_convertToBSpline, EG_convertToBSplineRange, EG_skinning,
    EG_tolerance, EG_getTolerance, EG_setTolerance, EG_getTopology,
    EG_makeTopology, EG_makeLoop, EG_getArea, EG_makeFace, EG_getBodyTopos,
    EG_indexBodyTopo, EG_objectBodyTopo, EG_inTopology, EG_inFace,
    EG_getEdgeUV, EG_getPCurve, EG_getBody, EG_makeSolidBody,
    EG_getBoundingBox, EG_getMassProperties, EG_isEquivalent,
    EG_sewFaces, EG_replaceFaces, EG_mapBody, EG_matchBodyFaces,
    EG_setTessParam, EG_makeTessGeom, EG_getTessGeom, EG_makeTessBody,
    EG_remakeTess, EG_mapTessBody, EG_locateTessBody, EG_getTessEdge,
    EG_getTessFace, EG_getTessLoops, EG_getTessQuads,
    EG_makeQuads, EG_getQuads, EG_getPatch, EG_quadTess, EG_insertEdgeVerts,
    EG_deleteEdgeVert, EG_moveEdgeVert, EG_openTessBody, EG_initTessBody,
    EG_statusTessBody, EG_setTessEdge, EG_setTessFace, EG_localToGlobal,
    EG_getGlobal, EG_solidBoolean, EG_intersection, EG_imprintBody, EG_filletBody,
    EG_chamferBody, EG_hollowBody, EG_extrude, EG_rotate, EG_sweep, EG_loft,
    EG_blend, EG_ruled, EG_blend_init, EG_blend_pos, EG_blend_sens, EG_blend_seval,
    EG_sens_free, EG_spline1dEval, EG_spline1dEval_psens, EG_spline1dDeriv,
    EG_spline1dFit, EG_spline1dFit_psens, EG_spline2dEval, EG_spline2dEval_psens,
    EG_spline2dDeriv, EG_spline2dAprx, eggGenerate, eggMorph, eggInfo, eggFree,
    EMP_Init, EMP_Done, EMP_ThreadCreate, EMP_ThreadExit, EMP_ThreadWait,
    EMP_ThreadSpin, EMP_ThreadID, EMP_ThreadDestroy, EMP_LockCreate, EMP_LockSet,
    EMP_LockTest, EMP_LockRelease, EMP_LockDestroy, EMP_for, EMP_sum, EMP_min,
    prm_CreateU, prm_CreateUV, prm_SmoothU, prm_SmoothUV, prm_NormalizeU, prm_NormalizeUV,
    prm_LimitGridSize, prm_FixedCfit, prm_FixedGrid, prm_BestCfit, prm_BestGrid,
    prm_EvalCfit, prm_EvalGrid, prm_FreeCfit, prm_FreeGrid

# own functions
export getBodyFromModel, cleanup, EG_join, EG_getBodyFaces

# types
export ego

# null pointer
export NULL_E

# constants
export MAGIC, MTESSPARAM, CONTXT, TRANSFORM, TESSELLATION,NIL, EMPTY
export REFERENCE, PCURVE, CURVE, SURFACE, NODE, EDGE, LOOP, FACE, SHELL
export BODY, MODEL, LINE, CIRCLE, ELLIPSE, PARABOLA, HYPERBOLA, TRIMMED
export BEZIER, BSPLINE, OFFSET, PLANE, SPHERICAL, CYLINDRICAL, REVOLUTION
export TOROIDAL, CONICAL, EXTRUSION, SREVERSE, NOMTYPE, SFORWARD, ONENODE
export TWONODE, OPEN, CLOSED, DEGENERATE, WIREBODY, FACEBODY, SHEETBODY
export SOLIDBODY, ATTRINT, ATTRREAL, ATTRSTRING, ATTRCSYS, ATTRPTR, SUBTRACTION
export INTERSECTION,FUSION, BOX, SPHERE, CONE, CYLINDER, TORUS, UISO, VISO
export NODEOFF, EDGEOFF, FACEDUP, FACECUT, FACEOFF

export EGADS_READERR, EGADS_TESSTATE, EGADS_EXISTS, EGADS_ATTRERR, EGADS_TOPOCNT
export EGADS_OCSEGFLT, EGADS_BADSCALE, EGADS_NOTORTHO, EGADS_DEGEN
export EGADS_CONSTERR, EGADS_TOPOERR, EGADS_GEOMERR, EGADS_NOTBODY
export EGADS_WRITERR, EGADS_NOTMODEL, EGADS_NOLOAD, EGADS_RANGERR, EGADS_NOTGEOM
export EGADS_NOTTESS, EGADS_EMPTY, EGADS_NOTTOPO, EGADS_REFERCE, EGADS_NOTXFORM
export EGADS_NOTCNTX, EGADS_MIXCNTX, EGADS_NODATA, EGADS_NONAME, EGADS_INDEXERR
export EGADS_MALLOC, EGADS_NOTOBJ, EGADS_NULLOBJ, EGADS_NOTFOUND, EGADS_SUCCESS
export EGADS_OUTSIDE, EGADSMAJOR, EGADSMINOR


text:
#### Export functions based on primitive data

"""
    export_data(data::Array, filename::S, headers::Vector{S})

Export `data` array to filename `filename` with headers `headers`.
"""
function export_data(data::Array, filename::S, headers::Vector{S}) where S<:AbstractString
  s = size(data)
  open(filename, "w") do f
    write(f, join(headers, ",")*"\n")
    for i in 1:s[1], j in 1:s[2]
      if j==s[2]
        @printf(f, "%18.8f\n", data[i,j])
      else
        @printf(f, "%18.8f,", data[i,j])
      end
    end
  end
end

"""
    import_data!(data::Array, filename::S, headers::Vector{S})

Import `data` array from filename `filename` with headers `headers`, exported from `export_csv`.
"""
function import_data!(data::Array{FT}, filename::S, headers::Vector{S}) where {S<:AbstractString,FT}
  data_all = readdlm(filename, ',')
  data .= FT.(data_all[2:end,:]) # Remove headers
end

ext: jl

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
using BinDeps

@BinDeps.setup

mpi = library_dependency("juliampi", aliases=["libjuliampi", "libjuliampi.so", "libjuliampi.dylib"], os = :Unix)
builddir = joinpath(BinDeps.depsdir(mpi), "build")
prefix = joinpath(BinDeps.depsdir(mpi), "usr")
src = joinpath(BinDeps.depsdir(mpi), "src")

if isdir(builddir)
    rm(builddir, recursive = true)
end
if isdir(prefix)
    rm(prefix, recursive = true)
end
if !isdir(src)
    mkdir(src)
end

# Currently these are all MPI 2.1, but we should version them
MPI_constants = [
    # types
    "MPI_BYTE",
    "MPI_INTEGER1",
    "MPI_INTEGER2",
    "MPI_INTEGER4",
    "MPI_INTEGER8",
    "MPI_REAL4",
    "MPI_REAL8",
    "MPI_COMPLEX8",
    "MPI_COMPLEX16",

    # comms
    "MPI_COMM_NULL",
    "MPI_COMM_SELF",
    "MPI_COMM_WORLD",

    # comm type
    "MPI_COMM_TYPE_SHARED",

    # reductions operations
    "MPI_OP_NULL",
    "MPI_BAND",
    "MPI_BOR",
    "MPI_BXOR",
    "MPI_LAND",
    "MPI_LOR",
    "MPI_LXOR",
    "MPI_MAX",
    "MPI_MAXLOC",
    "MPI_MIN",
    "MPI_NO_OP",
    "MPI_MINLOC",
    "MPI_PROD",
    "MPI_REPLACE",
    "MPI_SUM",

    # request
    "MPI_REQUEST_NULL",

    # info
    "MPI_INFO_NULL",

    # proc
    "MPI_PROC_NULL",

    # status
    "MPI_STATUS_SIZE", # size of Status array
    "MPI_SOURCE", # index in Status array
    "MPI_TAG",    # index in Status array
    "MPI_ERROR",  # index in Status array

    # lock_type
    "MPI_LOCK_EXCLUSIVE",
    "MPI_LOCK_SHARED",

    # misc constants
    "MPI_ANY_SOURCE",
    "MPI_ANY_TAG",
    "MPI_TAG_UB", # upper bound for tag value
    "MPI_UNDEFINED",
]

MPI_functions = [
    "MPI_ABORT",
    "MPI_ACCUMULATE",
    "MPI_ALLGATHER",
    "MPI_ALLGATHERV",
    "MPI_ALLREDUCE",
    "MPI_ALLTOALL",
    "MPI_ALLTOALLV",
    "MPI_BARRIER",
    "MPI_BCAST",
    "MPI_BSEND",
    "MPI_CANCEL",
    "MPI_CART_CREATE",
    "MPI_CART_COORDS",
    "MPI_CART_SHIFT",
    "MPI_COMM_DUP",
    "MPI_COMM_FREE",
    "MPI_COMM_GET_PARENT",
    "MPI_COMM_RANK",
    "MPI_COMM_SIZE",
    "MPI_COMM_SPLIT",
    "MPI_COMM_SPLIT_TYPE",
    "MPI_DIMS_CREATE",
    "MPI_EXSCAN",
    "MPI_FETCH_AND_OP",
    "MPI_FINALIZE",
    "MPI_FINALIZED",
    "MPI_GATHER",
    "MPI_GATHERV",
    "MPI_GET",
    "MPI_GET_ACCUMULATE",
    "MPI_GET_ADDRESS",
    "MPI_GET_COUNT",
    "MPI_GET_PROCESSOR_NAME",
    "MPI_INFO_CREATE",
    "MPI_INFO_DELETE",
    "MPI_INFO_FREE",
    "MPI_INFO_GET",
    "MPI_INFO_GET_VALUELEN",
    "MPI_INFO_SET",
    "MPI_INIT",
    "MPI_INITIALIZED",
    "MPI_INTERCOMM_MERGE",
    "MPI_IPROBE",
    "MPI_IRECV",
    "MPI_ISEND",
    "MPI_OP_CREATE",
    "MPI_OP_FREE",
    "MPI_PACK",
    "MPI_PACK_SIZE",
    "MPI_PROBE",
    "MPI_PUT",
    "MPI_RECV",
    "MPI_RECV_INIT",
    "MPI_REDUCE",
    "MPI_REQUEST_FREE",
    "MPI_RSEND",
    "MPI_SCAN",
    "MPI_SCATTER",
    "MPI_SCATTERV",
    "MPI_SEND",
    "MPI_SEND_INIT",
    "MPI_SSEND",
    "MPI_TEST",
    "MPI_TESTALL",
    "MPI_TESTANY",
    "MPI_TESTSOME",
    "MPI_TYPE_CONTIGUOUS",
    "MPI_TYPE_CREATE_STRUCT",
    "MPI_TYPE_COMMIT",
    "MPI_UNPACK",
    "MPI_WAIT",
    "MPI_WAITALL",
    "MPI_WAITANY",
    "MPI_WAITSOME",
    "MPI_WIN_ATTACH",
    "MPI_WIN_ALLOCATE_SHARED",
    "MPI_WIN_CREATE",
    "MPI_WIN_CREATE_DYNAMIC",
    "MPI_WIN_DETACH",
    "MPI_WIN_FENCE",
    "MPI_WIN_FLUSH",
    "MPI_WIN_FREE",
    "MPI_WIN_LOCK",
    "MPI_WIN_SHARED_QUERY",
    "MPI_WIN_SYNC",
    "MPI_WIN_UNLOCK",
    "MPI_WTICK",
    "MPI_WTIME",
    ]

MPI_C_constants = [
    "MPI_IN_PLACE",
]
MPI_C_types = ["Comm", "Info", "Win"]

open(joinpath(src,"gen_constants.f90"),"w") do f
    print(f,"""
program gen_constants
  implicit none
  include 'mpif.h'
""")
    for name in MPI_constants
        println(f,"  print '(\"const $name = Cint(\",i0,\")\")', $name")
    end
    print(f,"""
end program gen_constants
""")
end

open(joinpath(src,"gen_functions.c"),"w") do f
    print(f,"""
#include <stdio.h>
#include <stddef.h>
#include "jlmpi_f2c.h"
#include "mpi.h"

// check that the size of MPI_Aint is consistent
typedef int check_sizeof_MPI_Aint[sizeof(MPI_Aint) == sizeof(ptrdiff_t) ? 1 : -1];

#define STRING1(s) #s
#define STRING(s) STRING1(s)

int main(int argc, char *argv[]) {
""")
    for name in MPI_functions
        println(f,"  printf(\"const $name = (\\\"%s\\\",libmpi)\\n\", STRING(JLMPI_GLOBAL_($(lowercase(name)),$name)));")
    end

    for name in MPI_C_constants
        println(f,"  printf(\"const $name = %d\\n\", ((int)$name));")
    end

    for name in MPI_C_types
        println(f,"  printf(\"primitive type C$name %d end\\n\", (int)(sizeof(MPI_$name) * 8));")
    end
    print(f,"""
  return 0;
}
""")
end

# Add cmake flags to specify the MPI library location
configure = `cmake`
flags = ["MPI_C_COMPILER",
         "MPI_C_COMPILE_FLAGS",
         "MPI_C_INCLUDE_PATH",
         "MPI_C_LINK_FLAGS",
         "MPI_C_LIBRARIES",
         "MPI_Fortran_COMPILER",
         "MPI_Fortran_COMPILE_FLAGS",
         "MPI_Fortran_INCLUDE_PATH",
         "MPI_Fortran_LINK_FLAGS",
         "MPI_Fortran_LIBRARIES",
         "MPI_INCLUDE_PATH",
         "MPI_LIBRARIES"]
for flag in flags
    try
        val = ENV["JULIA_$flag"]
        global configure = `$configure -D$flag=$val`
    catch
        nothing
    end
end
configure = `$configure -DCMAKE_INSTALL_PREFIX=$src -DCMAKE_LIB_INSTALL_PREFIX=$prefix/lib ..`

provides(SimpleBuild,
    (@build_steps begin
        CreateDirectory(builddir)
        (@build_steps begin
            ChangeDirectory(builddir)
            (@build_steps begin
                configure
                `make`
                `make install`
            end)
        end)
    end), [mpi], os = :Unix)

@BinDeps.install Dict(:mpi => :mpi)


text:
"""
# plotcoef

Multiple regression coefficient plot for Turing models.

$(SIGNATURES)

### Required arguments
```julia
* `models`                             : Vector of `SampleModel`s to compare
* `pars`                               : Vector of parameters to include in comparison
```
### Optional arguments
```julia
* `fig=""`                             : File to store plot
* `title=""`                           : Title for plot
* `samples=NUTS(0.65)`                 : Sampler used for sampling
* `nsamples=2000`                      : No of samples (1000 warmup, 1000 used)
* `nchains=4`                          : No of chains used
* `func=nothing`                       : Funtion to apply to sample df
* `quap_samples=10000`                 : Default no of quap simulation samples
```
Currently the only function available is `quap` which will constrct
a quap model and simulate 10000 draws.

### Return values
```julia
* `(res, fig)`                         : (particles, plot)
```

"""
function plotcoef(
  models::Vector{T},
  pars::Vector{Symbol};
  fig="", title="",
  sampler=NUTS(0.65), nsamples=2000, nchains=4,
  func=nothing, quap_samples=10000) where {T <: DynamicPPL.Model}

  mnames = String.(nameof.(models))
  levels = length(models) * (length(pars) + 1)
  colors = [:blue, :red, :green, :darkred, :black]

  s = Vector{NamedTuple}(undef, length(models))
  for (mindx, mdl) in enumerate(models)
    if isnothing(func)
        chns = mapreduce(c -> sample(mdl, sampler, nsamples),
            chainscat, 1:nchains)
        df = DataFrame(Array(chns), names(chns, [:parameters]))
        m, l, u = estimparam(df)
        d = Dict{Symbol, NamedTuple}()
        for (indx, par) in enumerate(names(chns, [:parameters]))
            d[par] = (mean=m[indx], lower=l[indx], upper=u[indx])
        end
        s[mindx] =   (; d...)
    else
        quap_mdl = quap(mdl)
        post = rand(quap_mdl.distr, 10_000)
        df = DataFrame(post', [keys(quap_mdl.coef)...])
        m, l, u = estimparam(df)
        d = Dict{Symbol, NamedTuple}()
        for (indx, par) in enumerate([keys(quap_mdl.coef)...])
            d[par] = (mean=m[indx], lower=l[indx], upper=u[indx])
        end
        s[mindx] =   (; d...)
    end
  end

  xmin = 0; xmax = 0.0
  for i in 1:length(s)
    for par in pars
      syms = Symbol.(keys(s[i]))
      if Symbol(par) in syms
        mp = s[i][par].mean
        xmin = min(xmin, s[i][par].lower)
        xmax = max(xmax, s[i][par].upper)
      end
    end
  end

  ylabs = String[]
  for j in 1:length(models)
    for i in 1:length(pars)
      l = length(String(pars[i]))
      str = repeat(" ", 9-l) * String(pars[i])
      append!(ylabs, [str])
    end
    l = length(mnames[j])
    str = mnames[j] * repeat(" ", 9-l)
    append!(ylabs, [str])
  end

  ys = [string(ylabs[i]) for i = 1:length(ylabs)]
  p = plot(xlims=(xmin, xmax), leg=false, framestyle=:grid)
  title!(title)
  yran = range(1, stop=length(ylabs), length=length(ys))
  yticks!(yran, ys)

  line = 0
  for (mindx, model) in enumerate(models)
    line += 1
    hline!([line] .+ length(pars), color=:darkgrey, line=(2, :dash))
    for (pindx, par) in enumerate(pars)
      line += 1
      syms = Symbol.(keys(s[mindx]))
      if Symbol(par) in syms
        ypos = (line - 1)
        mp = s[mindx][Symbol(par)].mean
        lower = s[mindx][Symbol(par)].lower
        upper = s[mindx][Symbol(par)].upper
        plot!([lower, upper], [ypos, ypos], leg=false, color=colors[pindx])
        scatter!([mp], [ypos], color=colors[pindx])
      end
    end
  end
  if length(fig) > 0
    savefig(p, fig)
  end
  (s, p)
end

"""
# plotcoef

Multiple regression coefficient plot for a single Turing model.

$(SIGNATURES)

### Required arguments
```julia
* `model`                              : SampleModel to display
* `pars`                               : Vector of parameters to include in comparison
* `fig`                                : File to store the produced plot
```
### Optional arguments
```julia
* `fig=""`                             : File to store plot
* `title=""`                           : Title for plot
* `samples=NUTS(0.65)`                 : Sampler used for sampling
* `nsamples=2000`                      : No of samples (1000 warmup, 1000 used)
* `nchains=4`                          : No of chains used
* `func=nothing`                       : Funtion to apply to sample df
```
Currently the only function available is `quap` which will constrct
a quap model and simulate 10000 draws.

### Return values
```julia
* `(res, fig)`                         : (particles, plot)
```

"""
function plotcoef(
  mdl::DynamicPPL.Model,
  pars::Vector{Symbol};
  fig="", title="",
  sampler=NUTS(0.65), nsamples=2000, nchains=4,
  func=nothing, quap_samples=10000)

  plotcoef([mdl], pars;
    fig, title, sampler, nsamples, nchains, func, quap_samples)

 end

export
  plotcoef


text:
# Test of minimum spanning tree algorithms

using Graphs
using Base.Test

g = simple_inclist(7, is_directed=false)

wedges = [
    (1, 2, 7.),
    (2, 3, 8.),
    (1, 4, 5.),
    (2, 4, 9.),
    (2, 5, 7.),
    (3, 5, 3.),
    (4, 5, 15.),
    (4, 6, 6.),
    (5, 6, 8.),
    (5, 7, 9.),
    (6, 7, 11.) ]

m = length(wedges)
eweights = zeros(m)

for i = 1 : m
    we = wedges[i]
    add_edge!(g, we[1], we[2])
    eweights[i] = we[3]
end

@assert num_vertices(g) == 7
@assert num_edges(g) == m

re, rw = prim_minimum_spantree(g, eweights, 1)
@test length(re) == 6
@test length(rw) == 6

function verify_redge(re, rw, p, w)
    return re.source == p[1] && re.target == p[2] && rw == w
end

@test verify_redge(re[1], rw[1], (1, 4), 5.)
@test verify_redge(re[2], rw[2], (4, 6), 6.)
@test verify_redge(re[3], rw[3], (1, 2), 7.)
@test verify_redge(re[4], rw[4], (2, 5), 7.)
@test verify_redge(re[5], rw[5], (5, 3), 3.)
@test verify_redge(re[6], rw[6], (5, 7), 9.)

re, rw = kruskal_minimum_spantree(g, eweights)

@test verify_redge(re[1], rw[1], (3, 5), 3.)
@test verify_redge(re[2], rw[2], (1, 4), 5.)
@test verify_redge(re[3], rw[3], (4, 6), 6.)
@test verify_redge(re[4], rw[4], (1, 2), 7.)
@test verify_redge(re[5], rw[5], (2, 5), 7.)
@test verify_redge(re[6], rw[6], (5, 7), 9.)


text:
import JuLIP: read_dict, write_dict
export @D, @DD, @GRAD, @pot

using JuLIP: AtomicNumber

# ===========================================================================
#     implement some fun little macros for easier access
#     to the potentials
# ===========================================================================

# REMARK on @pot:
# ---------------
# Julia 0.4 version:
# call(pp::Potential, varargs...) = evaluate(pp, varargs...)
# call(pp::Potential, ::Type{Val{:D}}, varargs...) = evaluate_d(pp, varargs...)
# call(pp::Potential, ::Type{Val{:DD}}, varargs...) = evaluate_dd(pp, varargs...)
# call(pp::Potential, ::Type{Val{:GRAD}}, varargs...) = grad(pp, varargs...)
# unfortunately, in 0.5 `call` doesn't take anabstractargument anymore,
# which means that we need to specify for every potential how to
# create this syntactic sugar. This is what `@pot` is for.

"""
Annotate a type with `@pot` to setup the syntax sugar
for `evaluate, evaluate_d, evaluate_dd, grad`.

## Usage:

For example, the declaration
```julia
"documentation for `LennardJones`"
mutable struct LennardJones <: PairPotential
   r0::Float64
end
@pot LennardJones
```
creates the following aliases:
```julia
lj = LennardJones(1.0)
lj(args...) = evaluate(lj, args...)
@D lj(args...) = evaluate_d(lj, args...)
@DD lj(args...) = evaluate_dd(lj, args...)
@GRAD lj(args...) = grad(lj, args...)
```

Usage of `@pot` is not restricted to pair potentials, but can be applied to
*any* type.
"""
macro pot(fsig)
   @assert fsig isa Symbol
   sym = esc(:x)
   tsym = esc(fsig)
   quote
      @inline ($sym::$tsym)(args...) = evaluate($sym, args...)
      @inline ($sym::$tsym)(::Type{Val{:D}}, args...) = evaluate_d($sym, args...)
      @inline ($sym::$tsym)(::Type{Val{:DD}}, args...) = evaluate_dd($sym, args...)
      @inline ($sym::$tsym)(::Type{Val{:GRAD}}, args...) = grad($sym, args...)
   end
end

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

# next create macros that translate
"""
`@D`: Use to evaluate the derivative of a potential. E.g., to compute the
Lennard-Jones potential,
```julia
lj = LennardJones()
r = 1.0 + rand(10)
ϕ = lj(r)
ϕ' = @D lj(r)
```
see also `@DD`.
"""
macro D(fsig::Expr)
    @assert fsig.head == :call
    insert!(fsig.args, 2, Val{:D})
    for n = 1:length(fsig.args)
        fsig.args[n] = esc(fsig.args[n])
    end
    return fsig
end

"`@DD` : analogous to `@D`"
macro DD(fsig::Expr)
    @assert fsig.head == :call
    for n = 1:length(fsig.args)
        fsig.args[n] = esc(fsig.args[n])
    end
    insert!(fsig.args, 2, Val{:DD})
    return fsig
end

"`@GRAD` : analogous to `@D`, but escapes to `grad`"
macro GRAD(fsig::Expr)
    @assert fsig.head == :call
    for n = 1:length(fsig.args)
        fsig.args[n] = esc(fsig.args[n])
    end
    insert!(fsig.args, 2, Val{:GRAD})
    return fsig
end

# ----------------------------------------------------------------------
#    Managing a list of species
# ----------------------------------------------------------------------

abstract type AbstractZList end

"""
`ZList` and `SZList{NZ}` : simple data structures that store a list
of species and convert between atomic numbers and the index in the list.
Can be constructed via
* `ZList(zors)` : where `zors` is an Integer  or `Symbol` (single species)
* `ZList(zs1, zs2, ..., zsn)`
* `ZList([sz1, zs2, ..., zsn])`
* All of these take a kwarg `static = {true, false}`; if `true`, then `ZList`
will return a `SZList{NZ}` for (possibly) faster access.
"""
struct ZList <: AbstractZList
   list::Vector{AtomicNumber}
end

Base.length(zlist::AbstractZList) = length(zlist.list)

struct SZList{N} <: AbstractZList
   list::SVector{N, AtomicNumber}
end

function ZList(zlist::AbstractVector{<: Number};
               static = false, sorted=true)
   sortfun = sorted ? sort : identity
   return (static ? SZList(SVector( (AtomicNumber.(sortfun(zlist)))... ))
                  :  ZList( convert(Vector{AtomicNumber}, sortfun(zlist)) ))
end

ZList(s::Symbol; kwargs...) =
      ZList( [ atomic_number(s) ]; kwargs... )

ZList(S::AbstractVector{Symbol}; kwargs...) =
      ZList( atomic_number.(S); kwargs... )

ZList(args...; kwargs... ) =
      ZList( [args...]; kwargs...)

i2z(Zs::AbstractZList, i::Integer) = Zs.list[i]

function z2i(Zs::AbstractZList, z::AtomicNumber)
   if Zs.list[1] == JuLIP.Chemistry.__zAny__
      return 1
   end
   for j = 1:length(Zs.list)
      if Zs.list[j] == z
         return j
      end
   end
   error("z = $z not found in ZList $(Zs.list)")
end

zlist(V) = V.zlist
i2z(V, i::Integer) = i2z(zlist(V), i)
z2i(V, z::AtomicNumber ) = z2i(zlist(V), z)
numz(V) = length(zlist(V))

write_dict(zlist::ZList) = Dict("__id__" => "JuLIP_ZList",
                                  "list" => Int.(zlist.list))

read_dict(::Val{:JuLIP_ZList}, D::Dict) = ZList(D)
ZList(D::Dict) = ZList(D["list"])

write_dict(zlist::SZList) = Dict("__id__" => "JuLIP_SZList",
                                 "list" => Int.(zlist.list))
read_dict(::Val{:JuLIP_SZList}, D::Dict) = SZList(D)
SZList(D::Dict) = ZList([D["list"]...], static = true)


text:
struct ContinuousView{T, N, p, Tnode <: AbstractNode{<: Any, N}, Tblocks <: AbstractArray{<: AbstractLeafNode{T, N, p}, N}, Tindices <: Tuple} <: AbstractArray{T, N}
    parent::Tnode # needed for new node allocation
    blocks::Tblocks
    indices::Tindices
    blockoffset::NTuple{N, Int}
end

struct BlockLocalIndex{N}
    globalindex::NTuple{N, Int}
    blockindex::CartesianIndex{N}
    localindex::Int
end
@inline function BlockLocalIndex(x::ContinuousView, I::Int...)
    dims = TreeSize(parent(x))[end]
    blockindex, localindex = block_local_index(dims, I...)
    BlockLocalIndex(I, CartesianIndex(blockindex .- blockoffset(x)), Base._sub2ind(dims, localindex...))
end
BlockLocalIndex(x::ContinuousView, I::CartesianIndex) = BlockLocalIndex(x, Tuple(I)...)

Base.size(x::ContinuousView) = map(length, x.indices)
Base.axes(x::ContinuousView) = Base.IdentityUnitRange.(x.indices)
Base.parent(x::ContinuousView) = x.parent
rootnode(x::ContinuousView) = parent(x)

Base.propertynames(x::ContinuousView{<: Any, <: Any, <: Any, <: Any, <: AbstractArray{<: StructLeafNode}}) = (:parent, :blocks, :indices, fieldnames(eltype(x))...)
@inline function Base.getproperty(x::ContinuousView{<: Any, N, <: Any, <: Any, <: AbstractArray{<: StructLeafNode}}, name::Symbol) where {N}
    name == :parent && return getfield(x, :parent)
    name == :blocks && return getfield(x, :blocks)
    name == :indices && return getfield(x, :indices)
    name == :blockoffset && return getfield(x, :blockoffset)
    T = fieldtype(leafeltype(parent(x)), name)
    PropertyArray{T, N, name}(x)
end

blockoffset(x::ContinuousView) = x.blockoffset

for f in (:(Base.getindex), :isactive, :allocate!, :unsafe_getindex)
    @eval @inline function $f(x::ContinuousView, I::BlockLocalIndex)
        @_propagate_inbounds_meta
        @inbounds begin
            block = x.blocks[I.blockindex]
            $f(block, I.localindex)
        end
    end
end

@inline function Base.setindex!(x::ContinuousView, v, I::BlockLocalIndex)
    @_propagate_inbounds_meta
    block = x.blocks[I.blockindex]
    if isnull(block)
        leaf = _setindex!_getleaf(TreeView(parent(x)), v, I.globalindex...)
        x.blocks[I.blockindex] = leaf
    else
        block[I.localindex] = v
    end
    x
end

for f in (:(Base.getindex), :isactive, :allocate!, :unsafe_getindex)
    @eval @inline function $f(x::ContinuousView{<: Any, N}, I::Vararg{Int, N}) where {N}
        @boundscheck checkbounds(x, I...)
        @inbounds $f(x, BlockLocalIndex(x, I...))
    end
end

@inline function Base.setindex!(x::ContinuousView{<: Any, N}, v, I::Vararg{Int, N}) where {N}
    @boundscheck checkbounds(x, I...)
    @inbounds setindex!(x, v, BlockLocalIndex(x, I...))
end

function _continuousview(parentdims::Tuple, A::TreeView{<: Any, N}, I::Vararg{Union{Int, AbstractUnitRange, Colon}, N}) where {N}
    indices = to_indices(A, I)
    node = rootnode(A)
    dims = TreeSize(node)[end]
    start = block_index(dims, first.(indices)...)
    stop = block_index(dims, last.(indices)...)
    ContinuousView(node, generateblocks(CartesianIndex(start):CartesianIndex(stop), A), indices, start .- 1)
end
function _spotview(parentdims::Tuple, A::TreeView{<: Any, N}, I::Vararg{Int, N}) where {N}
    node = rootnode(A)
    dims = TreeSize(node)[end]
    start = CartesianIndex(block_index(dims, I...))
    stop = start + oneunit(start)
    ContinuousView(node,
                   generateblocks(SArray{NTuple{N, 2}}(start:stop), A),
                   UnitRange.(I, I.+dims),
                   Tuple(start) .- 1)
end

gettreeview(x::TreeView) = x
gettreeview(x::AbstractTreeArray) = x.tree
for f in (:continuousview, :spotview)
    _f = Symbol(:_, f)
    @eval begin
        @inline function $f(A::Union{AbstractTreeArray, TreeView}, I::Union{Int, AbstractUnitRange, Colon}...)
            @boundscheck checkbounds(A, I...)
            @inbounds $_f(size(A), gettreeview(A), I...)
        end
        @inline function $f(A::Union{AbstractTreeArray, TreeView}, inds::CartesianIndices)
            @boundscheck checkbounds(A, inds)
            @inbounds $_f(size(A), gettreeview(A), inds.indices...)
        end
        @inline function $f(A::Union{AbstractTreeArray, TreeView}, I::CartesianIndex)
            @boundscheck checkbounds(A, I)
            @inbounds $_f(size(A), gettreeview(A), Tuple(I)...)
        end
        @inline function $f(A::PropertyArray{<: Any, <: Any, name}, I...) where {name}
            @boundscheck checkbounds(A, I...)
            @inbounds getproperty($f(A.parent, I...), name)
        end
    end
end

dropleafindex(x::TreeIndex{depth}) where {depth} = TreeIndex(ntuple(i -> x.I[i], Val(depth-1)))
function generateblocks(blockindices, A::TreeView)
    node = rootnode(A)
    nblocks = nleafblocks(node)
    broadcast(blockindices) do i
        @_inline_meta
        blockindex = Tuple(i)
        checkbounds(Bool, CartesianIndices(nblocks), blockindex...) || return null(leaftype(node))
        I = blockindex .* TreeSize(node)[end] # global index
        treeindex = dropleafindex(TreeLinearIndex(A, I...))
        @inbounds rootnode(A[treeindex]) # should be @inbounds?
    end
end

@inline function Base.getindex(x::PropertyArray{<: Any, <: Any, name}, i::BlockLocalIndex) where {name}
    @_propagate_inbounds_meta
    getproperty(x.parent[i], name)
end

@inline function Base.setindex!(x::PropertyArray{<: Any, <: Any, name}, v, i::BlockLocalIndex) where {name}
    @_propagate_inbounds_meta
    allocated = allocate!(x.parent, i)
    set!(allocated, name, v)
    x
end

@inline function unsafe_getindex(x::PropertyArray{<: Any, <: Any, name}, i::BlockLocalIndex) where {name}
    @_propagate_inbounds_meta
    getproperty(unsafe_getindex(x.parent, i), name)
end


text:
function readintegers()
    [parse(Int, x) for x in split(readline())]
end

function solve()
    n, m = readintegers()
    a = readintegers()
    n - sum(a) < 0 ? -1 : n - sum(a)
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
export External

"""
    External(n1::Node, n2::Node, [portname::AbstractString])

`External` objects `show` a FastHenry .external command.
"""
struct External <: Element
  node1 :: Node
  node2 :: Node
  portname :: String
  External(n1::Node, n2::Node, portname::AbstractString = "") =
      new(n1,n2,portname)
end


text:
# # tests for lasso: 1/2\|Ax-b\|^2 + λ \|x\|_1

@testset "Lasso ($T)" for T in [Float32, Float64, ComplexF32, ComplexF64]
    using Test
    using LinearAlgebra
    using CIAOAlgorithms
    using ProximalOperators
    using Base.Iterators: take
    using Random

    Random.seed!(0)
    R = real(T)

    # problem definition
    N, n = 6, 3 # A in R^{N x n}
    p = 2 # nonzeros in the solution

    y_star = rand(R, N)
    y_star ./= norm(y_star) # y^star
    C = rand(R, N, n) .* 2 .- 1
    CTy = abs.(C' * y_star)
    # indices with decreasing order by abs
    perm = sortperm(CTy, rev = true)

    rho, λ = R(10), R(1)
    alpha = zeros(T, n)
    for i = 1:n
        if i <= p
            alpha[perm[i]] = λ / CTy[perm[i]]
        else
            alpha[perm[i]] = (CTy[perm[i]] < 0.1 * λ) ? λ : λ * rand() / CTy[perm[i]]
        end
    end
    A = C * diagm(0 => alpha)   # scaling the columns of Cin
    # generate the solution
    x_star = zeros(T, n)
    for i = 1:n
        if i <= p
            x_star[perm[i]] = rand() * rho / sqrt(p) * sign(dot(A[:, perm[i]], y_star))
        end
    end
    b = A * x_star + y_star

    # cost function
    cost_lasso(x) = norm(A * x - b)^2 / 2 + λ * norm(x, 1)

    f_star = cost_lasso(x_star)

    # preparations for the solver
    F = Vector{LeastSquares}(undef, 0)
    L = Vector{R}(undef, 0)
    for i = 1:N
        tempA = A[i:i, :]
        f = LeastSquares(tempA, b[i:i], R(N))
        Lf = opnorm(tempA)^2 * N
        push!(F, f)
        push!(L, Lf)
    end
    g = NormL1(λ)
    x0 = zeros(T, n)

    maxit = 1000
    tol = 1e-4

    @testset "Finito" begin
        # sweeping 1, 2, 3 for randomined, cyclic and shuffled sampling strategies, respectively.

        ## test the solver
        # basic finito
        @testset "basic Finito" for sweeping in collect(1:3)
            solver = CIAOAlgorithms.Finito{R}(maxit = maxit, sweeping = sweeping)
            x_finito, it_finito = solver(x0, F = F, g = g, L = L, N = N)
            @test cost_lasso(x_finito) - f_star < tol
            @test eltype(x_finito) == T
        end

        # limited memory finito
        @testset "LFinito" for sweeping in collect(2:3)
            # @testset "cyclical" begin
            solver =
                CIAOAlgorithms.Finito{R}(maxit = maxit, sweeping = sweeping, LFinito = true)
            x_finito, it_finito = solver(x0, F = F, g = g, L = L, N = N)
            @test cost_lasso(x_finito) - f_star < tol
            @test eltype(x_finito) == T
        end

        # # adaptive variant
        @testset "adaptive finito" for sweeping in collect(1:3)
            solver = CIAOAlgorithms.Finito{R}(
                maxit = maxit,
                tol = R(1e-5),
                sweeping = sweeping,
                adaptive = true,
            )
            x_finito, it_finito = solver(x0, F = F, g = g, L = L, N = N)
            @test cost_lasso(x_finito) - f_star < tol
            @test eltype(x_finito) == T
        end

        # basic finito with minibatch
        vec_ref = [(1, 2), (2, 2), (3, 3)] # different samplings and batch sizes
        @testset "Finito_minibatch" for (sweeping, batch) in vec_ref
            solver = CIAOAlgorithms.Finito{R}(
                maxit = maxit,
                sweeping = sweeping,
                minibatch = (true, batch),
            )
            x_finito, it_finito = solver(x0, F = F, g = g, L = L, N = N)
            @test cost_lasso(x_finito) - f_star < tol
            @test eltype(x_finito) == T
        end

        # limited memory finito with minibatch
        vec_ref = [(2, 1), (2, 2), (3, 3)] # different samplings and batch sizes
        @testset "LFinito_minibatch" for (sweeping, batch) in vec_ref
            solver = CIAOAlgorithms.Finito{R}(
                maxit = maxit,
                sweeping = sweeping,
                LFinito = true,
                minibatch = (true, batch),
            )
            x_finito, it_finito = solver(x0, F = F, g = g, L = L, N = N)
            @test cost_lasso(x_finito) - f_star < tol
            @test eltype(x_finito) == T
        end

        # test with user defined stepsizes
        @testset "γ and L as scalars" begin
            @testset "randomized" begin
                γ = N / maximum(L)
                solver = CIAOAlgorithms.Finito{R}(maxit = maxit, γ = γ)
                x_finito, it_finito = solver(x0, F = F, g = g, L = L, N = N)
                @test cost_lasso(x_finito) - f_star < tol
            end
            @testset "cyclic" begin
                solver = CIAOAlgorithms.Finito{R}(maxit = maxit)
                x_finito, it_finito = solver(x0, F = F, g = g, L = maximum(L), N = N)
                @test cost_lasso(x_finito) - f_star < tol
            end
        end

        ## test the iterator
        vec_ref = [(1, false, false), (2, false, false), (3, false, true), (3, true, false)]
        @testset "the iterator" for (sweeping, LFinito, adaptive) in vec_ref
            solver = CIAOAlgorithms.Finito{R}(
                sweeping = sweeping,
                LFinito = LFinito,
                adaptive = adaptive,
            )
            iter = CIAOAlgorithms.iterator(solver, x0, F = F, g = g, L = L, N = N)
            @test iter.x0 === x0

            for state in take(iter, 2)
                @test solution(state) === state.z
                @test eltype(solution(state)) == T
            end
        end
    end

    @testset "SVRG" begin

        ## test the solver
        γ = 1 / (7 * maximum(L))
        @testset "SVRG-Base" begin
            solver = CIAOAlgorithms.SVRG{R}(maxit = maxit, γ = γ)
            x_SVRG, it_SVRG = solver(x0, F = F, g = g, N = N)
            @test cost_lasso(x_SVRG) - f_star < tol
            @test eltype(x_SVRG) == T
        end
        @testset "SVRG++" begin
            solver = CIAOAlgorithms.SVRG{R}(maxit = 16, γ = γ, m = 1, plus = true)
            x_SVRG, it_SVRG = solver(x0, F = F, g = g, N = N)
            @test cost_lasso(x_SVRG) - f_star < tol
            @test eltype(x_SVRG) == T
        end

        # test the iterator
        @testset "the iterator" begin
            solver = CIAOAlgorithms.SVRG{R}(γ = γ)
            iter = CIAOAlgorithms.iterator(solver, x0, F = F, g = g, N = N)
            @test iter.x0 === x0

            for state in take(iter, 2)
                @test solution(state) === state.z_full
                @test eltype(solution(state)) == T
            end
            next = iterate(iter) # next = (state, state)
            # one iteration with the solver
            solver = CIAOAlgorithms.SVRG{R}(γ = γ, maxit = 1)
            x_finito, it_finito = solver(x0, F = F, g = g, L = L, N = N)
            @test solution(next[2]) == x_finito
        end
    end

    @testset "SAGA" begin

        ## test the solver
        @testset "SAGA-Base" begin
            solver = CIAOAlgorithms.SAGA{R}(maxit = maxit)
            x_SAGA, it_SAGA = solver(x0, F = F, g = g, N = N, L = L)
            @test cost_lasso(x_SAGA) - f_star < tol
            @test eltype(x_SAGA) == T
        end
        @testset "SAGA-stepsize" begin
            γ = 1 / (3 * maximum(L))
            solver = CIAOAlgorithms.SAGA{R}(maxit = maxit, γ = γ)
            x_SAGA, it_SAGA = solver(x0, F = F, g = g, N = N)
            @test cost_lasso(x_SAGA) - f_star < tol
            @test eltype(x_SAGA) == T
        end

        # test the iterator
        @testset "the iterator" begin
            γ = 1 / (3 * maximum(L))
            solver = CIAOAlgorithms.SAGA{R}(γ = γ)
            iter = CIAOAlgorithms.iterator(solver, x0, F = F, g = g, N = N)
            @test iter.x0 === x0

            for state in take(iter, 2)
                @test solution(state) === state.z
                @test eltype(solution(state)) == T
            end
            next = iterate(iter) # next = (state, state)
            # one iteration with the solver
            solver = CIAOAlgorithms.SAGA{R}(γ = γ, maxit = 1)
            x_finito, it_finito = solver(x0, F = F, g = g, L = L, N = N)
            @test solution(next[2]) == x_finito
        end
    end

    @testset "SAG" begin
        # note that proximal SAG may not be theoretically convergent
        maxit = 10000
        ## test the solver
        @testset "SAG-Base" begin
            solver = CIAOAlgorithms.SAG(R, maxit = maxit)
            x_SAG, it_SAG = solver(x0, F = F, g = g, N = N, L = L)
            @test cost_lasso(x_SAG) - f_star < tol
            @test eltype(x_SAG) == T
        end
        @testset "SAG" begin
            γ = 1 / (16 * maximum(L))
            solver = CIAOAlgorithms.SAG(R, maxit = maxit, γ = γ)
            x_SAG, it_SAG = solver(x0, F = F, g = g, N = N)
            @test cost_lasso(x_SAG) - f_star < tol
            @test eltype(x_SAG) == T
        end

        # test the iterator
        @testset "the iterator" begin
            γ = 1 / (16 * maximum(L))
            solver = CIAOAlgorithms.SAG(R, γ = γ)
            iter = CIAOAlgorithms.iterator(solver, x0, F = F, g = g, N = N)
            @test iter.x0 === x0

            for state in take(iter, 2)
                @test solution(state) === state.z
                @test eltype(solution(state)) == T
            end
            next = iterate(iter) # next = (state, state)
            # one iteration with the solver
            solver = CIAOAlgorithms.SAG(R, γ = γ, maxit = 1)
            x_finito, it_finito = solver(x0, F = F, g = g, L = L, N = N)
            @test solution(next[2]) == x_finito
        end
    end
end


text:
module Lighthouse

using Statistics, Dates, LinearAlgebra, Random, Logging
using Base.Threads
using StatsBase: StatsBase
using TensorBoardLogger
using Makie
using Printf

include("plotting.jl")

include("utilities.jl")
export majority

include("metrics.jl")
export confusion_matrix, accuracy, binary_statistics, cohens_kappa, calibration_curve

include("classifier.jl")
export AbstractClassifier

include("learn.jl")
export LearnLogger, learn!, upon, evaluate!, predict!

end # module


text:
# Part 1: digits never decreasing, has a double digit, is 6 digits long.
has_double_digit(ds) = any([ds[i] == ds[i+1] for i in 1:5])
never_decreasing(ds) = all([ds[i] >= ds[i+1] for i in 1:5])
valid_password(ds) = has_double_digit(ds) && never_decreasing(ds)
passwords = filter(x -> valid_password(digits(x)), 165432:707912)
println("Part 1: Number of possible passwords: $(length(passwords))")

# Part 2: now we need at least one group of 2 adjacent digits to not be part of larger group
has_isolated_double_digit(ds) =
    (ds[1] == ds[2] != ds[3]) ||
    any(ds[i] != ds[i+1] == ds[i+2] != ds[i+3] for i in 1:3) ||
    (ds[4] != ds[5] == ds[6])
passwords = [password for password in passwords if has_isolated_double_digit(digits(password))]
println("Part 2: Number of possible passwords: $(length(passwords))")

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
# Data topmute
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: December 2017
#

export marineTopmute2D, judiMarineTopmute2D
export model_topmute, judiTopmute, find_water_bottom, depth_scaling, judiDepthScaling, laplace

############################################ Data space preconditioners ################################################

function marineTopmute2D(Dobs::judiVector, muteStart::Integer; mute=Array{Any}(undef, 3), flipmask=false)
    # Data topmute for end-on spread marine streamer data
    Din = deepcopy(Dobs)

    for j=1:Din.nsrc

        # Design mute window
        x0 = 1f0
        xend = length(Din[j].geometry.xloc[1])
        nt = Din[j].geometry.nt[1]
        nrec = length(Din[j].geometry.xloc[1])
        drec = abs(Din[j].geometry.xloc[1][1] - Din[j].geometry.xloc[1][2])
        offsetDirectWave = 1.5f0*Din[j].geometry.t[1]
        idxOffset = Int(round(offsetDirectWave/drec))
        dx = round(idxOffset - idxOffset/10f0)

        if j==1 && ~isassigned(mute)
            z0 = muteStart - Int(round(muteStart/10))
            slope = 1.05f0*(nt - z0)/dx
            mute[1] = x0
            mute[2] = z0
            mute[3] = slope
        else#if j==1 && isassigned(mute)
            x0 = mute[1]
            z0 = mute[2]
            slope = mute[3]
        end

        mask = ones(Float32,nt,nrec)
        mask[1:z0,:] .= 0f0

        # Linear mute
        if (nrec-x0 < dx)
            x = nrec
            zIntercept = Int(round(z0+slope*(x-x0)))
            zax = z0+1:1:zIntercept
        else
            x = x0+dx
            zax = z0+1:1:nt
        end
        if length(zax) > 1
            xax = Array{Int}(round.(range(x0,stop=x,length=length(zax))))
        else
            xax = Int(round(x0))
        end
        for k=1:length(zax)
            mask[zax[k],xax[k]:end] .= 0f0
        end
        flipmask == true && (mask = reverse(mask, dims=2))
        Din.data[j] = Din.data[j].*mask
    end
    return Din
end

function judiMarineTopmute2D(muteStart,geometry;params=Array{Any}(undef, 3),flipmask=false)
# JOLI wrapper for the linear depth scaling function
    nsrc = length(geometry.xloc)
    N = 0
    for j=1:nsrc
        N += geometry.nt[j]*length(geometry.xloc[j])
    end
    D = joLinearFunctionFwd_T(N,N,
                             v -> marineTopmute2D(v,muteStart;mute=params, flipmask=flipmask),
                             w -> marineTopmute2D(w,muteStart;mute=params, flipmask=flipmask),
                             Float32,Float32,name="Data mute")
    return D
end

############################################ Model space preconditioners ###############################################

function model_topmute(n::Tuple{Int64,Int64}, mute_end::Array{Integer,1}, length::Int64, x_orig)
# Model domain topmute for a velocity model of dimensions n = [nx, nz].
    x = copy(reshape(x_orig,n))
    for j=1:n[1]
        mute_start = mute_end[j] - length
        filter = zeros(Float32, n[2])
        filter[1:mute_start-1] .= 0f0
        filter[mute_end[j]+1:end] .= 1f0
        taper_length = mute_end[j] - mute_start + 1
        taper = (1f0 .+ sin.((Float32(pi)*(0:taper_length-1))/(taper_length - 1) .- Float32(pi)/2f0))/2f0
        filter[mute_start:mute_end[j]] = taper
        M = diagm(0=>filter)
        global x[j,:] = x[j,:].*filter
    end
    return vec(x)
end

function model_topmute(n::Tuple{Int64,Int64}, mute_end::Int64, length::Int64, x_orig)
# Model domain topmute for a velocity model of dimensions n = [nx, nz].
    x = copy(reshape(x_orig,n))
    mute_start = mute_end - length
    filter = zeros(Float32, n[2])
    filter[1:mute_start-1] .= 0f0
    filter[mute_end+1:end] .= 1f0
    taper_length = mute_end - mute_start + 1
    taper = (1f0 .+ sin.((Float32(pi)*(0:taper_length-1))/(taper_length - 1) .- Float32(pi)/2f0))/2f0
    filter[mute_start:mute_end] = taper
    M = diagm(0=>filter)
    for j=1:n[1]
        global x[j,:] = x[j,:].*filter
    end
    return vec(x)
end

model_topmute(n::Tuple{Int64,Int64}, mute_end::Array{Float32, 2}, length, x) = vec(mute_end) .* vec(x)

function judiTopmute(n, mute_start, length)
    # JOLI wrapper for model domain topmute
    N = prod(n)
    T = joLinearFunctionFwd_T(N,N,
                             v -> model_topmute(n, mute_start, length, v),
                             w -> model_topmute(n, mute_start, length, w),
                             Float32,Float32,name="Model topmute")
    return T
end

function find_water_bottom(m)
    #return the indices of the water bottom of a seismic image
    n = size(m)
    idx = zeros(Integer,n[1])
    eps = 1e-4
    for j=1:n[1]
        k=1
        while true
            if abs(m[j,k]) > eps
                idx[j] = k
                break
            end
            k += 1
        end
    end
    return idx
end

function depth_scaling(m,model)
# Linear depth scaling function for seismic images
    m = reshape(m,model.n)
    filter = sqrt.(0f0:model.d[2]:model.d[2]*(model.n[2]-1))
    F = diagm(0=>filter)
    for j=1:model.n[1]
        m[j,:] = F*m[j,:]
    end
    return vec(m)
end

function judiDepthScaling(model)
# JOLI wrapper for the linear depth scaling function
    N = prod(model.n)
    D = joLinearFunctionFwd_T(N, N,
                             v -> depth_scaling(v,model),
                             w -> depth_scaling(w,model),
                             Float32,Float32,name="Depth scaling")
end

# function laplace(model::TimeModeling.Model)
# # 2D Laplace operator
#
#     # 2nd derivative in x direction
#     d1 = ones(Float32,model.n[1]-1)*1f0/(model.d[1]^2)
#     d2 = ones(Float32,model.n[1])*-2f0/(model.d[1]^2)
#     d3 = d1
#     d = (d1, d2, d3)
#     position = (-1, 0, 1)
#
#     Dx = spdiagm(d,position,model.n[1],model.n[1])
#     Ix = speye(Float32,model.n[1])
#
#     # 2nd derivative in z direction
#     d1 = ones(Float32,model.n[2]-1)*1f0/(model.d[2]^2)
#     d2 = ones(Float32,model.n[2])*-2f0/(model.d[2]^2)
#     d3 = d1
#     d = (d1, d2, d3)
#     position = (-1, 0, 1)
#
#     Dz = spdiagm(d,position,model.n[2],model.n[2])
#     Iz = speye(Float32,model.n[2])
#
#     # 2D Laplace operator
#     D = kron(Dz,Ix) + kron(Iz,Dx)
# end


text:
using Flux
using Flux: onehot, onehotbatch, logitcrossentropy, reset!, throttle
using Statistics: mean
using Random
using Unicode
using Parameters: @with_kw

@with_kw mutable struct Args
    lr::Float64 = 1e-3     # learning rate
    N::Int = 15            # Number of perceptrons in hidden layer
    test_len::Int = 100    # length of test data
    langs_len::Int = 0     # Number of different languages in Corpora
    alphabet_len::Int = 0  # Total number of characters possible, in corpora
    throttle::Int = 10     # throttle timeout
end

function get_processed_data(args)
    corpora = Dict()

    for file in readdir("corpus")
        lang = Symbol(match(r"(.*)\.txt", file).captures[1])
        corpus = split(String(read("corpus/$file")), ".")
        corpus = strip.(Unicode.normalize.(corpus, casefold=true, stripmark=true))
        corpus = filter(!isempty, corpus)
        corpora[lang] = corpus
    end

    langs = collect(keys(corpora))
    args.langs_len = length(langs)
    alphabet = ['a':'z'; '0':'9'; ' '; '\n'; '_']
    args.alphabet_len = length(alphabet)

    # See which chars will be represented as "unknown"
    unique(filter(x -> x ∉ alphabet, join(vcat(values(corpora)...))))

    dataset = [(onehotbatch(s, alphabet, '_'), onehot(l, langs)) for l in langs for s in corpora[l]] |> shuffle

    train, test = dataset[1:end-args.test_len], dataset[end-args.test_len+1:end]
    return train, test
end

function build_model(args)
    scanner = Chain(Dense(args.alphabet_len, args.N, σ), LSTM(args.N, args.N))
    encoder = Dense(args.N, args.langs_len)
    return scanner, encoder
end

function model(x, scanner, encoder)
    state = scanner.(x.data)[end]
    reset!(scanner)
    encoder(state)
end

function train(; kws...)
    # Initialize Hyperparameters
    args = Args(; kws...)
    # Load Data
    train_data, test_data = get_processed_data(args)

    @info("Constructing Model...")
    scanner, encoder = build_model(args)

    loss(x, y) = logitcrossentropy(model(x, scanner, encoder), y)
    testloss() = mean(loss(t...) for t in test_data)

    opt = ADAM(args.lr)
    ps = params(scanner, encoder)
    evalcb = () -> @show testloss()
    @info("Training...")
    Flux.train!(loss, ps, train_data, opt, cb = throttle(evalcb, args.throttle))
end

cd(@__DIR__)
train()

max_stars_repo_name: shawn-kern/model-zoo

max_issues_repo_name: shawn-kern/model-zoo

max_forks_repo_name: shawn-kern/model-zoo


text:
const stft_docs = """
The short-time fourier transform can be thought of as a filterbank, with each
filter centered at each fft bin frequency. The window function determines the
frequency response of each bandpass filter (each filter is an FIR filter of the
window function modulated by the center frequency).

This implementation has a few differences from the `stft` implementation in
DSP.jl that I think make it more consistent with the filterbank interpretation,
and make the output easier to reason about:

1. The first window is centered at the first sample of `x`, and in general the
   nth window is centered at `(n-1)*hop`. This also ensures that no input samples
   are multiplied by zero, which would make them unrecoverable.
2. `x` is zero-padded to ensure that the center of the last window is no
   earlier than the last sample of `x`. Similarly this ensures that no information
   from `x` is lost, and also makes the width of the result easier to reason about.
   (it is only a function of the hop size, not the window size).
3. Each window is centered at index 1 before taking the FFT, avoiding a linear
   phase shift in each frame. One benefit of this is that when the same frequency
   component present in `x` is picked up by multiple filters, the phase will be
   coherent across the bands (assuming it is not demodulated).
4. It provides the `demod` option to work with demodulated filter bank outputs.
   This is useful because the phase progresses much more slowly (proportional to
   the difference between the signal frequency and the center frequency, rather
   than at the signal frequency itself).
5. An `istft2` implementation is provided. Note that `istft2` is currently not
   fully compatible with DSP.jl's `stft` function. In order to make it compatible
   we'd need to match the window offset and also perform an `fftshift`.

`stft2` and `istft2` should satisfy the property that for any real-valued signal
`x` and `n`:

    istft2(stft2(x, n), n)[1:length(x)] ≈ x

Note that the round-trip may introduce some zero-padding at the end.
"""

# TODO: support windows longer than nfft - alternatively support windows with
# padding in the middle where the window length is still nfft. They should be
# equivilent ways to do the same thing, though maybe `stft` and `istft` will need
# to handle them specially...
# TODO: add options for boundary conditions so you don't get big bursts
# in the STFT at the beginning and end for truncated signals
"""
    stft2(x, nfft, hop=nfft÷2;
          window=rect,
          demod=false)

Perform a short-time fourier transform on a real-valued time signal `x`. The
result is a complex-valued matrix of dimension:

    (nfft÷2+1, cld(length(x)-1, hop)+1)

- `nfft` is the FFT size.
- `hop` is the hop size from frame to frame. You can think of this as the
   sampling period of the continuous filterbank.
- `window` adds an analysis window. It can be an `AbstractArray`
  of length `nfft`, or a function that can be called as
  `f(nfft; zerophase=true)` to generate such a window. Note that generally the
  window function should be periodic and centered around index `1`, which is
  what the `zerophase` option in DSP.jl's window functions does. See the DSP.jl
  documentation for more details.
- `demod` should be set to `true` if the STFT phase should be demodulated (i.e
  a "sliding window" stft). This is equivilent to de-modulating the output of
  each filter by its center frequency.

$stft_docs
"""
function stft2(x, nfft, hop=nfft÷2;
               window=rect,
               demod=false)
    # we could handle this, but it's simpler not to, and this is the common case
    @assert iseven(nfft)
    @assert hop > 0

    L = length(x)
    N = nfft÷2+1
    # the first window will be centered at the first sample. The last window
    # should be centered at the last sample or later, so that COLA will apply,
    # and we don't need to worry about numerical issues because some samples
    # got scaled way down by the window function.
    M = cld(L-1, hop)+1

    X = zeros(Complex{eltype(x)}, N, M)

    win = if window isa Function
        window(nfft; zerophase=true)
    elseif window isa AbstractArray
        @assert length(window) == nfft
        window
    else
        throw(ArgumentError("window must be an Array or Function"))
    end

    xbuf = zeros(eltype(x), nfft)
    Xbuf = similar(X, N)

    p = plan_rfft(xbuf)

    # each time we copy into the buffer for the FFT, we have to copy the
    # positive part (including the window center) and the negative part
    # separately. We want the window center to be at the beginning of the fft
    # buffer, or else we add a half-buffer circular delay to each frame

    # we could simplify this (at the expense of more copying) by just copying
    # the data in one block and calling `ifftshift`.

    xidx = 1
    nfft2 = nfft÷2

    for c in 1:M
        fill!(xbuf, zero(eltype(xbuf)))
        # copy the positive half to the beginning of the fft buffer
        # note the center of the window could be past the last sample
        # of `x`, in which case `npos` should be 0.
        npos = max(0, min(nfft2, L-xidx+1))
        copyto!(xbuf, 1, x, xidx, npos)
        # copy the negative half to the end of the fft buffer
        nneg = min(nfft2, xidx-1)
        # note that if xidx is past the end of `x`, we don't actually have
        # nneg samples available to copy. it's still useful to compute the
        # indices for the copying though
        copystart = xidx-nneg
        ncopy = max(0, min(nneg, L-copystart+1))
        copyto!(xbuf, nfft-nneg+1, x, copystart, ncopy)
        # zero-out any parts of the buffer that weren't filled. This should
        # only happen at the beginning and end of the process
        # xbuf[npos+1:end-nneg] .= zero(eltype(xbuf))
        if win !== nothing
            xbuf .*= win
        end
        mul!(Xbuf, p, xbuf)
        copyto!(X, (c-1)*N+1, Xbuf, 1, N)
        xidx += hop
    end

    if demod
        # we want to convert to a sliding-window STFT  by multiplying by
        # exp(-jnω). This is equivilent to demodulating each band of the STFT
        # filter bank
        n = (0:M-1)*hop
        ω = range(0, π, length=N)
        # TODO we could do this without allocating a whole copy of X by doing
        # the multipliation frame-wise
        X .*= exp.(-im .* n' .* ω)
    end

    X
end

"""
    istft2(X, nfft, hop=nfft÷2;
           out_nfft=nfft,
           out_hop=Int(out_nfft * hop / nfft),
           window=rect,
           demod=false)
Perform an inverse short-time fourier transform on the complex-valued matrix
`X`.

- `nfft` is the FFT size that was used to compute the STFT.
- `hop` is the hop size of the frames in `X`
- `out_nfft` is the desired iFFT size to use for reconstruction.
- `out_hop` is the hop size for the resynthesized windows. You can set this
  to something different from `hop` for basic time stretching (though you will
  likely also need to manually update the phase values).
- `window` adds a resynthesis window. It can be an `AbstractArray`
  of length `nfft`, or a function that can be called as
  `f(nfft; zerophase=true)` to generate such a window. We currently only handle
  the painless case where the round-trip window (elementwise product of the
  analysis and synthesis windows) has the constant-overlap-add property. Note
  that generally the window function should be periodic and centered around
  index `1`, which is what the `zerophase` option in DSP.jl's window functions
  does. See the DSP.jl documentation for more details.
- `demod` should be set to `true` if the STFT phase has been demodulated (i.e
  a "sliding window" stft).

$stft_docs
"""
function istft2(X, nfft, hop=nfft÷2;
                out_nfft=nfft,
                # out_hop will throw an error if the given sizes aren't
                # compatible. TODO: give a better error message or handle this
                # more gracefully
                out_hop=Int(out_nfft * hop / nfft),
                window=rect,
                demod=false)
    # we could handle this case, but it's simpler not to
    @assert iseven(nfft)
    @assert iseven(out_nfft)
    @assert hop > 0
    @assert out_hop > 0

    N = size(X, 1)
    M = size(X, 2)

    # the output signal starts in the center of the first window, and
    # ends no later than the center of the last window
    L = (M-1) * out_hop + 1
    out = zeros(real(eltype(X)), L)

    # TODO: if the element type is e.g. Float32, we should use the same for
    # the window
    # TODO: scaling should probably be handled more automatically, and
    # also we should handle the non-COLA case. Also we're not scaling
    # the windows above because we're assuming they're COLA and add to 1
    # ALSO there are still some artifacts at the very beginning and end,
    # where the samples don't get the right number of windows adding up,
    # so the scale factor is off.
    win = if window isa Function
        window(out_nfft, zerophase=true)
    elseif window isa AbstractArray
        @assert length(window) == out_nfft
        window
    else

    end

    if demod
        # we want to convert to a sliding-signal STFT (which is what most stft
        # implementations do natively) by multiplying by exp(jnω). This is
        # equivilent to modulating the STFT filter bank by the band's center
        # frequency.
        n = (0:M-1)*hop
        ω = range(0, π, length=N)
        # TODO we could do this without allocating a whole copy of X by doing
        # the multipliation frame-wise
        X .*= exp.(im .* n' .* ω)
    end

    # this is probably pretty close to supporting nfft != N, but not really
    # tested in that case
    Xbuf = similar(X, out_nfft÷2+1)

    # fill with zeros, so that if we partially fill it later, it ends up being
    # zero-padded
    fill!(Xbuf, zero(eltype(X)))
    outbuf = similar(out, out_nfft)
    p = plan_irfft(Xbuf, out_nfft)
    copylen = min(length(Xbuf), N)

    nfft2 = out_nfft÷2

    outoffset = 0
    for c in 1:M
        copyto!(Xbuf, 1, X, (c-1)*N+1, copylen)
        mul!(outbuf, p, Xbuf)
        outbuf .*= win
        npos = min(nfft2, L-outoffset)
        @simd for i in 1:npos
            @inbounds out[outoffset+i] += outbuf[i]
        end
        nneg = min(nfft2, outoffset)
        @simd for i in 1:nneg
            @inbounds out[outoffset-nneg+i] += outbuf[end-nneg+i]
        end
        outoffset += out_hop
    end

    out
end

max_stars_repo_name: ssfrr/SFRDSP

max_issues_repo_name: ssfrr/SFRDSP

max_forks_repo_name: ssfrr/SFRDSP


text:
using Documenter, QuHamiltonian

makedocs(
     modules = [QuHamiltonian],
     format  = :html,
     sitename = "QuHamiltonian",
     pages = Any[
         "Introduction"   => "index.md"
     ],
     # Use clean URLs, unless built as a "local" build
     html_prettyurls = !("local" in ARGS),
     # html_canonical = "https://juliadocs.github.io/Documenter.jl/latest/",
)

deploydocs(
     repo = "git@github.com:Roger-luo/QuHamiltonian.jl.git",
     target = "build",
     julia = "1.0",
     osname = "linux",
     deps = nothing,
     make = nothing,
)

max_stars_repo_name: Roger-luo/QuHamiltonian

max_issues_repo_name: Roger-luo/QuHamiltonian

max_forks_repo_name: Roger-luo/QuHamiltonian


text:
module ThreeBodyUnitTests

using Test
using UnitfulAstrodynamics
using UnitfulAstrodynamics.TwoBody.Systems: Earth, Sun

@testset "ThreeBody" begin

    # Hardcode Gravity parameters for the Sun,
    # and the Earth-Moon System
    μₛ = 1.32712440018e20u"m^3/s^2"
    μₑ = 4.035032351966808e14u"m^3/s^2"

    # Dimensional initial conditions for spacecraft
    r = [2e9, 7000, 2000]u"km"
    v = [0.001, 0.08, 0.02]u"km/s"
    t = 500u"d"
    a = 1.0u"AU"

    # Construct nondimensional state
    sys = ThreeBodyState(μₛ, μₑ, a, r, v, t);

    @test true

    # This should run!
    μ = nondimensionalize(mass_parameter(Sun), mass_parameter(Earth)) |> upreferred
    r, v, T = halo(μ; Az = 1e-3, L = 2)
    sys = NondimensionalThreeBodyState(r, v, μ, T, 1.0u"AU", 500u"d")
    sys = redimensionalize(sys, mass_parameter(Sun), mass_parameter(Earth))
    sys = nondimensionalize(sys)

    @test true

end

end


text:
"""
filter_matrix_meshwidth(f, domain, N)

Assemble discrete filtering matrix from a continuous filter `f` width constant width
\$h(x) = \\Delta x / 2\$.
"""
function filter_matrix_meshwidth(f::TopHatFilter, domain::PeriodicIntervalDomain, N)
    x = discretize(domain, N)
    Δx = 1 // N * (domain.right - domain.left)
    h = f.width
    h₀ = h(x[1])

    all(≈(h(x[1])), h.(x)) || error("Filter width must be constant")
    if Δx .≈ 2h₀
        # Three point stencil
        inds = [-1, 0, 1]
        stencil = [1 // 24, 11 // 12, 1 // 24]
    elseif Δx .≈ h₀
        # Three point stencil
        inds = [-1, 0, 1]
        stencil = [1 // 6, 4 // 6, 1 // 6]
    else
        error("Filter width must be equal to mesh width")
    end

    # Five point stencil
    # inds = [-2, -1, 0, 1, 2]
    # stencil = [-6 // 2033, 77 // 1440, 863 // 960, 77 // 1440, -6 // 2033]

    # Construct banded matrix
    diags = [i => fill(s, N - abs(i)) for (i, s) ∈ zip(inds, stencil)]
    W = spdiagm(diags...)

    # Periodic extension of three point stencil
    W[1, end] = stencil[1]
    W[end, 1] = stencil[end]

    # Periodic extension of five point stencil
    # W[1, [end - 1, end]] = stencil[[1, 2]]
    # W[2, end] = stencil[1]
    # W[end - 1, 1] = stencil[end]
    # W[end, [1, 2]] = stencil[[end-1, end]]

    W
end

function filter_matrix_meshwidth(f::TopHatFilter, domain::ClosedIntervalDomain, N)
    x = discretize(domain, N)
    Δx = 1 // N * (domain.right - domain.left)
    h = f.width
    h₀ = h(x[1])

    all(≈(h(x[1])), h.(x)) || error("Filter width must be constant")
    Δx ≈ 2h₀ || error("Filter width must be equal to mesh width")

    # Three point stencil
    inds = [-1, 0, 1]
    stencil = [1 / 24, 11 / 12, 1 / 24]

    # Five point stencil
    # inds = [-2, -1, 0, 1, 2]
    # stencil = [-6 / 2033, 77 / 1440, 863 / 960, 77 / 1440, -6 / 2033]

    # Construct banded matrix
    diags = [i => fill(s, N + 1 - abs(i)) for (i, s) ∈ zip(inds, stencil)]
    W = spdiagm(diags...)

    # Boundary weights for three point stencil
    W[1, [1, 2]] = [3 / 4, 1 / 4]
    W[end, [end, end - 1]] = [3 / 4, 1 / 4]

    # Boundary weights for five point stencil
    # W[1, [1, 2, 3]] = [2 / 3, 5 / 12, -1 / 12]
    # W[2, [1, 2, 3, 4]] = [1/24, 11/12, 1/24, 0]
    # W[end - 1, [end, end - 1, end - 2, end - 3]] = [1/24, 11/12, 1/24, 0]
    # W[end, [end, end - 1, end - 2]] = [2 / 3, 5 / 12, -1 / 12]

    W
end


text:
function load_distributions_samples()
return quote

using .Distributions: Uniform, Normal

"""
    _sample_unit_nsphere_muller!(D::Vector{Vector{N}}, n::Int, p::Int;
                                 rng::AbstractRNG=GLOBAL_RNG,
                                 seed::Union{Int, Nothing}=nothing) where {N}

Draw samples from a uniform distribution on an ``n``-dimensional unit
sphere using Muller's method.

### Input

- `D`    -- output, vector of points
- `n`    -- dimension of the sphere
- `p`    -- number of random samples
- `rng`  -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed` -- (optional, default: `nothing`) seed for reseeding

### Output

A vector of `nsamples` vectors.

### Algorithm

This function implements Muller's method of normalised Gaussians [1] to uniformly
sample over the ``n``-dimensional sphere ``S^n`` (which is the bounding surface
of the ``n``-dimensional unit ball).

Given ``n`` canonical Gaussian random variables ``Z₁, Z₂, …, Z_n``, the
distribution of the vectors

```math
\\dfrac{1}{α}\\left(z₁, z₂, …, z_n\\right)^T,
```
where ``α := \\sqrt{z₁² + z₂² + … + z_n²}``, is uniform over ``S^n``.

[1] Muller, Mervin E. *A note on a method for generating points uniformly on
    n-dimensional spheres.* Communications of the ACM 2.4 (1959): 19-20.
"""
function _sample_unit_nsphere_muller!(D::Vector{Vector{N}}, n::Int, p::Int;
                                      rng::AbstractRNG=GLOBAL_RNG,
                                      seed::Union{Int, Nothing}=nothing) where {N}
    rng = reseed(rng, seed)
    Zdims = [Normal() for _ in 1:n] # normal distributions, one for each dimension
    v = Vector{N}(undef, n) # sample direction
    @inbounds for j in 1:p
        α = zero(N)
        for i in 1:n
            v[i] = rand(rng, Zdims[i])
            α += v[i]^2
        end
        D[j] = v ./ sqrt(α)
    end
    return D
end

"""
    _sample_unit_nball_muller!(D::Vector{Vector{N}}, n::Int, p::Int;
                               rng::AbstractRNG=GLOBAL_RNG,
                               seed::Union{Int, Nothing}=nothing) where {N}

Draw samples from a uniform distribution on an ``n``-dimensional unit ball
using Muller's method.

### Input

- `D`    -- output, vector of points
- `n`    -- dimension of the ball
- `p`    -- number of random samples
- `rng`  -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed` -- (optional, default: `nothing`) seed for reseeding

### Output

A vector of `nsamples` vectors.

### Algorithm

This function implements Muller's method of normalised Gaussians [1] to
sample from the interior of the ball.

Given ``n`` Gaussian random variables ``Z₁, Z₂, …, Z_n`` and a uniformly
distributed random variable ``r`` with support in ``[0, 1]``, the distribution
of the vectors

```math
\\dfrac{r^{1/n}}{α} \\left(z₁, z₂, …, z_n\\right)^T,
```
where ``α := \\sqrt{z₁² + z₂² + … + z_n²}``, is uniform over the
``n``-dimensional unit ball.

[1] Muller, Mervin E. *A note on a method for generating points uniformly on
    n-dimensional spheres.* Communications of the ACM 2.4 (1959): 19-20.
"""
function _sample_unit_nball_muller!(D::Vector{Vector{N}}, n::Int, p::Int;
                                    rng::AbstractRNG=GLOBAL_RNG,
                                    seed::Union{Int, Nothing}=nothing) where {N}

    rng = reseed(rng, seed)
    Zdims = [Normal() for _ in 1:n] # normal distributions, one for each dimension
    Zrad = Uniform() # distribution to pick random radius
    one_over_n = one(N)/n
    v = Vector{N}(undef, n) # sample direction
    @inbounds for j in 1:p
        α = zero(N)
        for i in 1:n
            v[i] = rand(rng, Zdims[i])
            α += v[i]^2
        end
        r = rand(rng, Zrad)
        β = r^one_over_n / sqrt(α)
        D[j] = v .* β
    end
    return D
end

end # quote
end # function load_distributions_samples()


text:
module Issue760

using Gridap
using Test

cells = (2,2)
domain = (0,1,0,1)
model = CartesianDiscreteModel(domain,cells)

Ω = Interior(model)
Ω1 = Interior(model,[1,2])
Ω2 = Interior(model,[2,4])

f = CellField(1,Ω1)
x = get_cell_points(Ω2)
arr = f(x)

#print_op_tree(arr)

@test arr.maps.value.values_pos.values[1] == [1,1,1,1]
@test arr.maps.value.values_neg.values[1] == [0,0,0,0]

end # module


text:
using LinearAlgebra, SparseArrays, InfiniteArrays, FillArrays, LazyArrays, Statistics, DSP, Test
    import InfiniteArrays: OrientedInfinity, OneToInf, InfUnitRange, InfStepRange

@testset "∞" begin
    @test ∞ ≠ 1
    @test ∞ == ∞
    @test ∞ == Inf

    @test isless(1, ∞)
    @test !isless(Inf, ∞)
    @test !isless(∞, Inf)
    @test !isless(∞, 1)

    @test !isless(∞, ∞)
    @test !(∞ < ∞)
    @test ∞ ≤ ∞
    @test !(∞ > ∞)
    @test ∞ ≥ ∞

    @test ∞ + ∞ ≡ ∞
    @test ∞ + 1 ≡ ∞
    @test *(∞) ≡ ∞
    @test ∞*∞ ≡ ∞

    # oriented infinity
    @test OrientedInfinity(∞) ≡ convert(OrientedInfinity, ∞) ≡ OrientedInfinity() ≡
        OrientedInfinity(false)

    @test -∞ ≡ OrientedInfinity(true)
    @test +∞ ≡ ∞

    @test ∞ == +∞ == OrientedInfinity(∞)
    @test ∞ ≠  -∞
    @test 1-∞ == -∞

    @test (-∞)*(-∞) ≡ ∞*OrientedInfinity(∞) ≡ OrientedInfinity(∞)*∞

    @test  isless(-∞, 1)
    @test !isless(-∞, -Inf)
    @test !isless(-Inf, -∞)
    @test !isless(1, -∞)

    @test max(1,∞) == max(∞,1) == ∞
    @test min(1,∞) == min(∞,1) == 1
    @test maximum([1,∞]) == ∞
    @test minimum([1,∞]) == 1

    @test OrientedInfinity(true) + OrientedInfinity(true) == OrientedInfinity(true)
    @test OrientedInfinity(false) + OrientedInfinity(false) == OrientedInfinity(false)
    @test OrientedInfinity(true)+1 == OrientedInfinity(true)
    @test OrientedInfinity(false)+1 == OrientedInfinity(false)

    @test exp(im*π/4)*∞ == Inf+im*Inf
    @test exp(im*π/4)+∞ == ∞
end

@testset "ranges" begin
    @test size(10:1:∞) == (∞,)
    @testset "colon" begin
        @test @inferred(10:1:∞) === @inferred(range(10; step=1, length=∞))
        @inferred(1:0.2:∞)  === @inferred(range(1; step=0.2, length=∞))
        @inferred(1.0:0.2:∞)  === @inferred(range(1.0; step=0.2, length=∞))
        @inferred(1:∞) === @inferred(range(1; length=∞))
    end
    @test_throws ArgumentError 2:-.2:∞
    @test_throws ArgumentError 0.0:-∞

    @testset "indexing" begin
        L32 = @inferred(Int32(1):∞)
        L64 = @inferred(Int64(1):∞)
        @test @inferred(L32[1]) === Int32(1) && @inferred(L64[1]) === Int64(1)
        @test L32[2] == 2 && L64[2] == 2
        @test L32[3] == 3 && L64[3] == 3
        @test L32[4] == 4 && L64[4] == 4

        @test @inferred((1.0:∞)[1]) === 1.0
        @test @inferred((1.0f0:∞)[1]) === 1.0f0
        @test @inferred((Float16(1.0):∞)[1]) === Float16(1.0)

        @test @inferred((0.1:0.1:∞)[2]) === 0.2
        @test @inferred((0.1f0:0.1f0:∞)[2]) === 0.2f0
        @test @inferred((1:∞)[1:4]) === 1:4
        @test @inferred((1.0:∞)[1:4]) === 1.0:4
        @test (2:∞)[1:4] == 2:5
        @test (1:∞)[2:5] === 2:5
        @test (1:∞)[2:2:5] === 2:2:4
        @test (1:2:∞)[2:6] === 3:2:11
        @test (1:2:∞)[2:3:7] === 3:6:13

        @test isempty((1:∞)[5:4])
        @test_throws BoundsError (1:∞)[8:-1:-2]
    end
    @testset "length" begin
        @test length(.1:.1:∞) == ∞
        @test length(1.1:1.1:∞) == ∞
        @test length(1.1:1.3:∞) == ∞
        @test length(1:1:∞) == ∞
        @test length(1:.2:∞) == ∞
        @test length(1.:.2:∞) == ∞
        @test length(2:-.2:-∞) == ∞
        @test length(2.:-.2:-∞) == ∞
    end

    @test_throws ArgumentError (2:.2:-∞)
    @test_throws ArgumentError (2.:.2:-∞)
    @test_throws ArgumentError (1:-∞)

    @testset "intersect" begin
        @test intersect(1:∞, 2:3) == 2:3
        @test intersect(2:3, 1:∞) == 2:3
        @test intersect(1:∞, 2:∞) == 2:∞
        @test intersect(2:∞, 1:∞) == 2:∞
        @test intersect(-3:∞, 2:8) == 2:8
        @test intersect(1:∞, -2:3:15) == 1:3:15
        @test intersect(1:∞, -2:3:∞) == 1:3:∞
        @test intersect(1:11, -2:2:∞) == intersect(-2:2:∞,1:11) == 2:2:10
        @test intersect(1:∞, -2:1:15) == 1:15
        @test intersect(1:∞, 15:-1:-2) == 1:15
        @test intersect(1:∞, 15:-4:-2) == 3:4:15
        @test intersect(-20:∞, -10:3:-2) == -10:3:-2
        @test isempty(intersect(-5:5, -6:13:∞))
        @test isempty(intersect(1:∞, 15:4:-2))

        @test @inferred(intersect(0:3:∞, 0:4:∞)) == intersect(0:4:∞, 0:3:∞) == 0:12:∞

        @test intersect(24:-3:0, 0:4:∞) == 0:12:24
        @test_throws ArgumentError intersect(1:6:∞, 0:4:∞) # supporting empty would break type inferrence

        @test intersect(1:∞,3) == 3:3
        @test intersect(1:∞, 2:∞, UnitRange(3,7), UnitRange(4,6)) == UnitRange(4,6)

        @test intersect(1:∞, 2) === intersect(2, 1:∞) === 2:2
        @test_skip intersect(1.0:∞, 2) == intersect(2, 1.0:∞) == [2.0] # gives infinite loop
    end

    @testset "sort/sort!/partialsort" begin
        @test sort(1:∞) == 1:∞
        @test sort!(1:∞) == 1:∞
    end
    @testset "in" begin
        @test 0 in UInt(0):100:∞
        @test 0 in 0:-100:-∞
        @test ∞ ∉ UInt(0):100:∞

        @test !(3.5 in 1:∞)
        @test (3 in 1:∞)
        @test (3 in 5:-1:-∞)

        let r = 0.0:0.01:∞
            @test (r[30] in r)
        end
        let r = (-4*Int64(maxintfloat(Int === Int32 ? Float32 : Float64))):∞
            @test (3 in r)
            @test (3.0 in r)
        end
    end
    @testset "in() works across types, including non-numeric types (#21728)" begin
        @test 1//1 in 1:∞
        @test 1//1 in 1.0:∞
        @test !(5//1 in 6:∞)
        @test !(5//1 in 6.0:∞)
        @test Complex(1, 0) in 1:∞
        @test Complex(1, 0) in 1.0:∞
        @test Complex(1.0, 0.0) in 1:∞
        @test Complex(1.0, 0.0) in 1.0:∞
        @test_skip !(Complex(1, 1) in 1:∞)  # this is an infinite-loop at the moment
        @test_skip !(Complex(1, 1) in 1.0:∞)
        @test_skip !(Complex(1.0, 1.0) in 1:∞)
        @test_skip !(Complex(1.0, 1.0) in 1.0:∞)
        @test !(π in 1:∞)
    end

    @testset "indexing range with empty range (#4309)" begin
        @test (3:∞)[5:4] == 7:6
        @test (0:2:∞)[7:6] == 12:2:10
    end
    # indexing with negative ranges (#8351)
    for a=[3:∞, 0:2:∞], b=[0:1, 2:-1:0]
        @test_throws BoundsError a[b]
    end

    @testset "sums of ranges" begin
        @test sum(1:∞) ≡ mean(1:∞) ≡ median(1:∞) ≡ ∞
        @test sum(0:∞) ≡ mean(1:∞) ≡ median(1:∞) ≡ ∞
        @test sum(0:2:∞) ≡ mean(0:2:∞) ≡ median(0:2:∞) ≡ OrientedInfinity(∞)
        @test sum(0:-2:-∞) ≡ mean(0:-2:-∞) ≡ median(0:-2:-∞) ≡ -∞
    end

    @testset "broadcasted operations with scalars" begin
        @test broadcast(-, 1:∞, 2) ≡ -1:∞
        @test broadcast(-, 1:∞, 0.25) ≡ 1-0.25:∞
        @test broadcast(+, 1:∞, 2) ≡ 3:∞
        @test broadcast(+, 1:∞, 0.25) ≡ 1+0.25:∞
        @test broadcast(+, 1:2:∞, 1) ≡ 2:2:∞
        @test broadcast(+, 1:2:∞, 0.3) ≡ 1+0.3:2:∞
        @test broadcast(-, 1:2:∞, 1) ≡ 0:2:∞
        @test broadcast(-, 1:2:∞, 0.3) ≡ 1-0.3:2:∞
        @test broadcast(-, 2, 1:∞) ≡ 1:-1:-∞
    end

    @testset "near-equal ranges" begin
        @test 0.0:0.1:∞ != 0.0f0:0.1f0:∞
    end

    @testset "comparing InfiniteUnitRanges and OneToInf" begin
        @test 1:2:∞ == 1:2:∞ != 1:3:∞ != 2:3:∞ == 2:3:∞ != 2:∞
        @test 1:1:∞ == 1:∞ == 1:∞ == OneToInf() == OneToInf()
    end

    @testset "Base.OneTo (misleading) overrides" begin
        @test Base.OneTo{BigInt}(∞) isa OneToInf{BigInt}
        @test Base.OneTo(∞) isa OneToInf{Int}
    end

    @testset "issue #6973" begin
        r1 = 1.0:0.1:∞
        r2 = 1.0f0:0.2f0:∞
        r3 = 1:2:∞
        @test r1 + r1 == 2*r1
        @test_broken r1 + r2 == 2.0:0.3:∞
        @test (2r1)-3r1 == -1:(2step(r1)-3step(r1)):-∞
        @test_broken (r1 + r2) - r2 == r1
        @test r1 + r3 == 2.0:2.1:∞
        @test r3 + r3 == 2 * r3
    end

    @testset "Preservation of high precision upon addition" begin
        let r = (-0.1:0.1:∞) + broadcast(+, -0.3:0.1:∞, 1e-12)
            @test_broken r[3] == 1e-12
        end
    end

    @testset "issue #8584" begin
        @test (0:1//2:∞)[1:2:3] == 0:1//1:1
    end

    @testset "issue #9962" begin
        @test eltype(0:1//3:∞) <: Rational
        @test (0:1//3:∞)[1] == 0
        @test (0:1//3:∞)[2] == 1//3
    end
    @testset "converting ranges (issue #10965)" begin
        @test promote(0:∞, UInt8(2):∞) === (0:∞, 2:∞)
        @test convert(InfUnitRange{Int}, 0:∞) === 0:∞
        @test convert(InfUnitRange{Int128}, 0:∞) === Int128(0):∞

        @test InfUnitRange{Int16}(1:∞) ≡ AbstractVector{Int16}(1:∞) ≡
                AbstractArray{Int16}(1:∞) ≡ Int16(1):∞

        @test OneToInf{Int16}(OneToInf()) ≡ AbstractVector{Int16}(OneToInf()) ≡
                AbstractArray{Int16}(OneToInf()) ≡ OneToInf{Int16}()

        @test promote(0:1:∞, UInt8(2):UInt8(1):∞) === (0:1:∞, 2:1:∞)
        @test convert(InfStepRange{Int,Int}, 0:1:∞) === 0:1:∞
        @test convert(InfStepRange{Int128,Int128}, 0:1:∞) === Int128(0):Int128(1):∞

        @test promote(0:1:∞, 2:∞) === (0:1:∞, 2:1:∞)
        @test convert(InfStepRange{Int128,Int128}, 0:∞) === Int128(0):Int128(1):∞
        @test convert(InfStepRange, 0:∞) === 0:1:∞
        @test convert(InfStepRange{Int128,Int128}, 0.:∞) === Int128(0):Int128(1):∞

        @test_broken promote(0f0:inv(3f0):∞, 0.:2.:∞) === (0:1/3:∞, 0.:2.:∞)

        @test promote(0:1/3:∞, 0:∞) === (0:1/3:∞, 0.:1.:∞)
    end

    @testset "inf-range[inf-range]" begin
        @test (1:∞)[1:∞] == 1:∞
        @test (1:∞)[OneToInf()] == 1:∞
        @test (1:∞)[2:∞] == 2:∞
        @test_throws BoundsError (1:∞)[-1:∞]
        @test (1:-1:-∞)[1:∞] == 1:-1:-∞
    end

    @testset "OneToInf" begin
        let r = OneToInf()
            @test !isempty(r)
            @test length(r) == ∞
            @test size(r) == (∞,)
            @test step(r) == 1
            @test first(r) == 1
            @test last(r) == ∞
            @test minimum(r) == 1
            @test maximum(r) == ∞
            @test r[2] == 2
            @test r[2:3] === 2:3
            @test_throws BoundsError r[0]
            @test broadcast(+, r, 1) === 2:∞
            @test 2*r === 2:2:∞
            @test r + r === 2:2:∞

            @test r - r === Zeros{Int}(∞)

            @test intersect(r, Base.OneTo(2)) == Base.OneTo(2)
            @test intersect(r, 0:5) == 1:5
            @test intersect(r, 2) === intersect(2, r) === 2:2
        end
    end

    @testset "show" begin
        @test summary(1:∞) == "InfUnitRange{Int64} with indices OneToInf()"
        @test Base.inds2string(axes(1:∞)) == "OneToInf()"
    end
end

@testset "fill" begin
    for A in (Zeros(∞), Fill(1,∞), Ones(∞))
        @test length(A) == ∞
    end
    @test size(Zeros(∞,5)) === (∞,5)
    @test size(Zeros(5,∞)) === (5,∞)
end

@testset "diagonal" begin
    D = Diagonal(1:∞)
    @test D[1:10,1:10] == Diagonal(1:10)
end

@testset "concat" begin
    A = Vcat(1:10, 1:∞)
    @test @inferred(length(A)) == ∞
    @test @inferred(A[5]) == A[15] == 5

    A = Vcat(Ones(1,∞), Zeros(2,∞))
    @test @inferred(size(A)) == (3,∞)
    @test @inferred(A[1,5]) == 1
    @test @inferred(A[3,5]) == 0
    @test_throws BoundsError A[4,1]

    A = Vcat(Ones{Int}(1,∞), Diagonal(1:∞))
    @test @inferred(size(A)) ≡ (∞,∞)
    @test @inferred(A[1,5]) ≡ 1
    @test @inferred(A[5,5]) ≡ 0
    @test @inferred(A[6,5]) ≡ 5
    @test_throws BoundsError A[-1,1]

    A = Vcat(Ones{Float64}(1,∞), Diagonal(1:∞))
    @test @inferred(size(A)) ≡ (∞,∞)
    @test @inferred(A[1,5]) ≡ 1.0
    @test @inferred(A[5,5]) ≡ 0.0
    @test @inferred(A[6,5]) ≡ 5.0
    @test_throws BoundsError A[-1,1]

    A = Vcat(1, Zeros(∞))
    @test @inferred(A[1]) ≡ 1.0
    @test @inferred(A[2]) ≡ 0.0

    A = Hcat(Ones(∞), Zeros(∞,2))
    @test @inferred(size(A)) == (∞,3)
    @test @inferred(A[5,1]) == 1
    @test @inferred(A[5,3]) == 0
    @test_throws BoundsError A[1,4]

    A = Hcat(Ones{Int}(∞), Diagonal(1:∞))
    @test @inferred(size(A)) ≡ (∞,∞)
    @test @inferred(A[5,1]) ≡ 1
    @test @inferred(A[5,5]) ≡ 0
    @test @inferred(A[5,6]) ≡ 5
    @test_throws BoundsError A[-1,1]
end

# This should be generalized, but it at the moment
# it is restricted to a single Number. Support smart
# addition for any number of Number/SVector's would be better
# allowibng for the tail to be variable lenth
@testset "Vcat special case" begin
    @test Vcat(1,Zeros{Int}(∞)) + Vcat(3,Zeros{Int}(∞)) ≡
          Vcat(1,Zeros{Int}(∞)) .+ Vcat(3,Zeros{Int}(∞)) ≡
          Vcat(4,Zeros{Int}(∞))

          size(Vcat(1:∞)) ≡ (∞,)
end

@testset "Fill indexing" begin
    B = Ones(∞,∞)
    @test IndexStyle(B) == IndexCartesian()
    V = view(B,:,1)
    @test_broken size(V) == (∞,1)
    V = view(B,1,:)
    @test size(V) == (∞,)
    V = view(B,1:1,:)
    @test size(V) == (1,∞)
end

@testset "∞ BroadcastArray" begin
    A = 1:∞
    B = BroadcastArray(exp, A)
    @test length(B) == ∞
    @test B[6] == exp(6)
    @test exp.(A) ≡ B
    B = Diagonal(1:∞) .+ 1
    @test B isa BroadcastArray{Int}
    @test B[1,5] ≡ 1
    @test B[6,6] == 6+1
    B = Diagonal(1:∞) - Ones{Int}(∞,∞) # lowers to broadcast
    @test B isa BroadcastArray{Int}
    @test B[1,5] ≡ -1
    @test B[6,6] == 6-1
end

@testset "Taylor ODE" begin
    e₁ = Vcat(1, Zeros(∞))
    D = Hcat(Zeros(∞), Diagonal(1:∞))
    @test typeof(Eye(∞)) == Eye{Float64,OneToInf{Int}}
    @test Base.BroadcastStyle(typeof(Eye(∞))) == LazyArrays.LazyArrayStyle{2}()
    L = Vcat(e₁', Eye(∞) + D)
    @test L[1:3,1:3] == [1.0 0.0 0.0;
                         1.0 1.0 0.0;
                         0.0 1.0 2.0]
end

@testset "Cumsum" begin
    @test cumsum(Ones(∞)) ≡ 1.0:1.0:∞
    @test cumsum(Fill(2,∞)) ≡ 2:2:∞
    @test cumsum(Ones{Int}(∞)) ≡ Base.OneTo(∞)
    @test cumsum(Ones{BigInt}(∞)) ≡ Base.OneTo{BigInt}(∞)

    x = Vcat([3,4], Ones{Int}(5), 3, Fill(2,∞))
    y = @inferred(cumsum(x))
    @test y isa Vcat
    @test y[1:12] == cumsum(x[1:12])

    @test cumsum(x).arrays[2] ≡ 8:12
    @test last(y.arrays) == sum(x[1:9]):2:∞
    r = (3:4:∞)
    @test cumsum(r)[1:20] == cumsum(r[1:20])
end

@testset "Sub-array" begin
    @test Ones(∞)[3:∞] ≡ Ones(∞)
    @test Ones{Int}(∞)[4:6] ≡ Ones{Int}(3)
    @test (1:∞)[3:∞] ≡ 3:∞
end

@testset "Broadcast Fill Lowers" begin
    @test broadcast(+, Zeros{Int}(∞) , Fill(1,∞)) isa Fill
    @test broadcast(+, Zeros{Int}(∞) , Zeros(∞)) isa Zeros
    @test broadcast(*, Ones(∞), Ones(∞)) ≡ Ones(∞)
end

@testset "maximum/minimum Vcat" begin
    x = Vcat(1:2, [1,1,1,1,1], 3, Fill(4,∞))
    @test maximum(x) == 4
    @test minimum(x) == 1

    @test_throws ArgumentError maximum(exp.(1:∞))
end

@testset "conv" begin
    @test conv(1:∞, [2]) ≡ conv([2], 1:∞) ≡ 2:2:∞
    @test conv(1:2:∞, [2]) ≡ conv([2], 1:2:∞) ≡ 2:4:∞
    @test conv(1:∞, Ones(∞))[1:5] == conv(Ones(∞),1:∞)[1:5] == [1,3,6,10,15]
    @test conv(Ones(∞), Ones(∞)) ≡ 1.0:1.0:∞
    @test conv(Ones{Int}(∞), Ones{Int}(∞)) ≡ Base.OneTo(∞)
    @test conv(Ones{Bool}(∞), Ones{Bool}(∞)) ≡ Base.OneTo(∞)
end

@testset "show" begin
    @test repr(Vcat(1:∞)) == "[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, …]"
    @test repr(Vcat(2,1:∞)) == "[2, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, …]"
end

@testset "Array constructor errors" begin
    @test_throws ArgumentError Array{Float64}(undef, ∞)
    @test_throws ArgumentError Array{Float64}(undef, ∞, ∞)
    @test_throws ArgumentError Array{Float64}(undef, 1, ∞)
    @test_throws ArgumentError Array{Float64}(undef, ∞, 1)

    @test_throws ArgumentError Vector{Float64}(undef, ∞)
    @test_throws ArgumentError Matrix{Float64}(undef, ∞, ∞)
    @test_throws ArgumentError Matrix{Float64}(undef, 1, ∞)
    @test_throws ArgumentError Matrix{Float64}(undef, ∞, 1)
end


text:
# This is the main output file for PS3 Question 1
# We write down matrices T, Z, Q, H of each process
# α_{t} = T * α_{t-1} + η_{t} ~ N(0, Q)
# y_{t} = Z * α_{t} + ϵ_{t} ~ N(0, H)

## Packages
using Random
using Distributions
using Plots
using Optim
using ForwardDiff
using LinearAlgebra

## Files
include("KalmanFilter.jl")

rng = Random.seed!(1234);
n = 100; # number of periods

############################## Part 1 (AR(1)) #############################
## Question 1.1: AR(1) process (x_{t} = ρ * x_{t-1} + ϵ_{t})
# Parameters
ρ_AR1 = 0.6;
σ_AR1 = 0.5;

# Simulate AR(1) process
y_AR1 = repeat([reshape([0.0], 1, 1)], n); # [0.0] is a Vector{Float64}, 'reshape' changes it to Matrix{Float64}
[y_AR1[i] = [ρ_AR1] * y_AR1[i - 1] + σ_AR1 * randn(rng, 1) for i in 2:n];

## Getting Matrices to apply kalman filter
T_AR1 = reshape([ρ_AR1], 1, 1);
Z_AR1 = reshape([1.0], 1, 1);
Q_AR1 = reshape([σ_AR1^2], 1, 1);
H_AR1 = reshape([0.0], 1, 1);
α0_AR1 = reshape(mean(y_AR1), 1, 1);

## Applying kalman_filter
αhat_AR1, v_AR1, yhat_AR1, P_AR1, F_AR1, G_AR1, K_AR1  = KalmanFilter(T_AR1, Z_AR1, Q_AR1, H_AR1, n, α0_AR1, y_AR1)

## Plots
vec_y_AR1 = [y_AR1[i][1, 1] for i in 1:n]
vec_yhat_AR1 = [yhat_AR1[i][1, 1] for i in 2:n]

plot(vec_y_AR1, xlabel = " Time ", label =" Simulated AR(1) ", legend=:bottomright)
plot!(vec_yhat_AR1, label = "Estimate from Kalman Filter")
savefig("figs/AR1.png")

## Writing the likelhood function
# For the AR1 process we need to estimate two Parameters ρ and σ:
function log_L_AR1(θ::Vector)  # θ = [ρ, log_σ]
    ρ1 = θ[1]
    σ1 = exp(θ[2])
    T1 = reshape([ρ1], 1, 1)
    Q1 = reshape([σ1^2], 1, 1)

    v1 = [KalmanFilter(T1, Z_AR1, Q1, H_AR1, n, α0_AR1, y_AR1)[2][i][1] for i in 1:n]
    F1 = [KalmanFilter(T1, Z_AR1, Q1, H_AR1, n, α0_AR1, y_AR1)[5][i][1] for i in 1:n]
    Ft = abs.(det.(vec(F1)))

    return sum([-0.5 * log(Ft[i]) - 0.5 * v1[i]' * inv(Ft[i]) * v1[i] for i in 1:n])
end

## Optimization
res_AR1 = optimize(θ->-log_L_AR1(θ), [-0.6, -0.69])
ρ_hat_AR1, σ_hat_AR1 = Optim.minimizer(res_AR1)[1], exp(Optim.minimizer(res_AR1)[2]) # (0.61, 0.48)

############################## Part 2 (AR(2)) #############################
## Question 1.2: AR(2) process (x_{t} = ρ1 * x_{t-1} + ρ2 * x_{t-2} + ϵ_{t})
# Parameters
ρ1_AR2 = 0.3
ρ2_AR2 = 0.4
σ_AR2 = 0.5

# Simulate AR(2) process
y_AR2 = repeat([reshape([0.0], 1, 1)], n);
[y_AR2[i] = [ρ1_AR2] * y_AR2[i - 1] + [ρ2_AR2] * y_AR2[i - 2] + σ_AR2 * randn(rng, 1) for i in 3:n];

## Getting Matrices to apply kalman filter
T_AR2 = [ρ1_AR2 ρ2_AR2; 1.0 0.0];
Z_AR2 = [1.0 0.0];
Q_AR2 = [σ_AR2^2 0.0; 0.0 0.0];
H_AR2 = reshape([0.0], 1, 1);
α0_AR2 = [mean(y_AR2); mean(y_AR2)];

αhat_AR2, v_AR2, yhat_AR2, P_AR2, F_AR2, G_AR2, K_AR2 = KalmanFilter(T_AR2, Z_AR2, Q_AR2, H_AR2, n, α0_AR2, y_AR2)

vec_y_AR2 = [y_AR2[i][1, 1] for i in 3:n]
vec_y_hat_AR2 = [yhat_AR2[i][1, 1] for i in 3:n]

plot(vec_y_AR2, xlabel = " Time ", label =" Simulated AR(2) ", legend=:bottomright)
plot!(vec_y_hat_AR2, label = "Estimate from Kalman Filter")
savefig("figs/AR2.png")

## Writing the likelhood function
# For the AR2 process we need to estimate three Parameters ρ1, ρ2, and σ:
function log_L_AR2(θ::Vector)  # θ = [ρ1, ρ2, log_σ]
    ρ1 = θ[1]
    ρ2 = θ[2]
    σ1 = exp(θ[3])
    T1 = [ρ1 ρ2; 1.0 0.0]
    Q1 = [σ1^2 0.0; 0.0 0.0]

    v1 = [KalmanFilter(T1, Z_AR2, Q1, H_AR2, n, α0_AR2, y_AR2)[2][i][1] for i in 1:n]
    F1 = [KalmanFilter(T1, Z_AR2, Q1, H_AR2, n, α0_AR2, y_AR2)[5][i][1] for i in 1:n]
    Ft = abs.(det.(vec(F1)))

    return sum([-0.5 * log(Ft[i]) - 0.5 * v1[i]' * inv(Ft[i]) * v1[i] for i in 3:n])
end

## Optimization
res_AR2 = optimize(θ->-log_L_AR2(θ), [0.3, 0.7, -0.6])
ρ1_hat_AR2,ρ2_hat_AR2,σ_hat_AR2 = Optim.minimizer(res_AR2)[1],Optim.minimizer(res_AR2)[2],exp(Optim.minimizer(res_AR2)[3])
# (0.26, 0.48, 0.54)

############################## Part 3 (MA(1)) #############################
## Question 1.3: MA(1) process (x_{t} = ϵ_{t} + ρ * ϵ_{t-1})
# Parameters
ρ_MA1 = 0.7
σ_MA1 = 0.8

# Simulate MA(1) process
ϵ_MA1 = rand(Normal(0, σ_MA1), n)
y_MA1 = repeat([reshape([0.0], 1, 1)], n);
[y_MA1[i] = reshape([ϵ_MA1[i]], 1,1) + [ρ_MA1]*ϵ_MA1[i-1] for i in 2:n]

## Getting Matrices to apply kalman filter
T_MA1 = [0.0 0.0; 1.0 0.0]
Z_MA1 = [1.0 ρ_MA1]
Q_MA1 = [σ_MA1^2 0.0; 0.0 0.0]
H_MA1 = reshape([0.0], 1, 1)

α0_MA1 = reshape([0.0; 0.0],2,1);
αhat_MA1, v_MA1, yhat_MA1, P_MA1, F_MA1, G_MA1, K_MA1 = KalmanFilter(T_MA1, Z_MA1, Q_MA1, H_MA1, n, α0_MA1, y_MA1)

vec_y_MA1 = [y_MA1[i][1, 1] for i in 2:n]
vec_y_hat_MA1= [yhat_MA1[i][1, 1] for i in 2:n]

plot(vec_y_MA1, xlabel = " Time ", label =" Simulated MA(1) ", legend=:bottomright)
plot!(vec_y_hat_MA1, label = "Estimate from Kalman Filter")
savefig("figs/MA1.png")

## Writing the likelhood function
# For the MA1 process we need to estimate two Parameters ρ and σ:
function log_L_MA1(θ::Vector)  # θ = [ρ, log_σ]
    ρ1 = θ[1]
    σ1 = exp(θ[2])
    Z1 = [1.0 ρ1]
    Q1 = [σ1^2 0.0; 0.0 0.0]

    v1 = [KalmanFilter(T_MA1, Z1, Q1, H_MA1, n, α0_MA1, y_MA1)[2][i][1] for i in 1:n]
    F1 = [KalmanFilter(T_MA1, Z1, Q1, H_MA1, n, α0_MA1, y_MA1)[5][i][1] for i in 1:n]
    Ft = abs.(det.(vec(F1)))

    return sum([-0.5 * log(Ft[i]) - 0.5 * v1[i]' * inv(Ft[i]) * v1[i] for i in 2:n])
end

## Optimization
res_MA1 = optimize(θ->-log_L_MA1(θ), [0.6, -0.6])
ρ_hat_MA1,σ_hat_MA1 = Optim.minimizer(res_MA1)[1],exp(Optim.minimizer(res_MA1)[2]) # (0.61, 0.73)

############################## Part 4 (Random Walk) #############################
## Question 1.4: Random Walk process
# Parameters
σ_ϵ = 0.8
σ_η = 0.6

# Simulate MA(1) process
ϵ_RW = rand(Normal(0, σ_ϵ), n)
η_RW = rand(Normal(0, σ_η), n)

μ_RW = repeat([reshape([0.0], 1, 1)], n)
[μ_RW[i] = μ_RW[i-1] + reshape([η_RW[i]], 1,1) for i in 2:n]

y_RW = repeat([reshape([0.0], 1, 1)], n)
[y_RW[i] = μ_RW[i] + reshape([ϵ_RW[i]], 1,1) for i in 1:n]

## Getting Matrices to apply kalman filter
T_RW = reshape([1.0], 1, 1)
Z_RW = copy(T_RW)
Q_RW = reshape([σ_η^2], 1, 1)
H_RW = reshape([σ_ϵ^2], 1, 1)

α0_RW = reshape([0.0],1,1);
αhat_RW, v_RW, yhat_RW, P_RW, F_RW, G_RW, K_RW = KalmanFilter(T_RW, Z_RW, Q_RW, H_RW, n, α0_RW, y_RW)

vec_y_RW = [y_RW[i][1, 1] for i in 2:n]
vec_y_hat_RW= [yhat_RW[i][1, 1] for i in 2:n]

plot(vec_y_RW, xlabel = " Time ", label =" Simulated Random Walk ", legend=:bottomright)
plot!(vec_y_hat_RW, label = "Estimate from Kalman Filter")
savefig("figs/RW.png")

## Writing the likelhood function
# For the Random Walk process we need to estimate two Parameters σ_ϵ and σ_η:
function log_L_RW(θ::Vector)  # θ = [log_σ_ϵ, log_σ_η]
    σ_ϵ = exp(θ[1])
    σ_η = exp(θ[2])
    Q1 = reshape([σ_η^2], 1, 1)
    H1 = reshape([σ_ϵ^2], 1, 1)

    v1 = [KalmanFilter(T_RW, Z_RW, Q1, H1, n, α0_RW, y_RW)[2][i][1] for i in 1:n]
    F1 = [KalmanFilter(T_RW, Z_RW, Q1, H1, n, α0_RW, y_RW)[5][i][1] for i in 1:n]
    Ft = abs.(det.(vec(F1)))

    return sum([-0.5 * log(Ft[i]) - 0.5 * v1[i]' * inv(Ft[i]) * v1[i] for i in 2:n])
end

## Optimization
res_RW = optimize(θ->-log_L_RW(θ), [-0.3, -0.6])
σ_ϵ_hat, σ_η_hat= exp(Optim.minimizer(res_RW)[1]), exp(Optim.minimizer(res_RW)[2]) # (0.70, 0.65)


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module Products
using Nemo

export ModuleProducts
export ∏, Product, F!, RisingFactorial, ↑, FallingFactorial, ↓, MultiFactorial
export V000407, V124320, V265609, V000142, V081125, V001147
export V000165, V032031, V007559, V008544, V007696, V001813, V008545, V047053

"""

* ∏, Product, F!, RisingFactorial, ↑, FallingFactorial, ↓, MultiFactorial
* V000407, V124320, V265609, V000142, V081125, V001147, V000165, V032031, V007559, V008544, V007696, V001813, V008545, V047053
"""
const ModuleProducts = ""

"""

If ``a ≤ b`` then return the product of ``i`` ∈ ``a:b`` else return ``1``.
"""
function ∏(a, b)
    n = b - a
    if n < 24
        p = fmpz(1)
        for k ∈ a:b
            p *= k
        end
        return fmpz(p)
    end
    m = div(a + b, 2)
    ∏(a, m) * ∏(m + 1, b)
end

"""

If ``a ≤ b`` then return the product of ``i`` ∈ ``a:b`` else return ``1``.
"""
Product(a, b) = ∏(a, b)

"""

Return the accumulated product of an array.
"""
function ∏(A)
    function prod(a, b)
        n = b - a
        if n < 24
            p = fmpz(1)
            for k ∈ a:b
                p *= A[k]
            end
            return fmpz(p)
        end
        m = div(a + b, 2)
        prod(a, m) * prod(m + 1, b)
    end
    A == [] && return fmpz(1)
    prod(1, length(A))
end

"""

Return the accumulated product of an array.
"""
Product(A) = ∏(A)

"""

Return ``\frac{n!}{⌊n/2⌋!}``.
"""
V081125(n::Int) = ∏(div(n, 2) + 1, n)

"""

Return the rising factorial which is the product of ``i`` ∈ ``n:(n + k - 1)``.
"""
RisingFactorial(n::Int, k::Int) = ∏(n, n + k - 1)

"""

Return the rising factorial which is the product of ``i`` ∈ ``n:(n + k - 1)``. A convenient infix syntax for the rising factorial is n ↑ k.
"""
↑(n, k) = RisingFactorial(n, k)

"""

Return the rising factorial i.e. the product of ``i`` ∈ ``n:(n + k - 1)``.
"""
V265609(n::Int, k::Int) = RisingFactorial(n, k)

# *** deprecated, use n ↑ k instead ***
#"""
#Return 'Pochhammer(n, k)', which is ambiguous in the literature, as the
#RisingFactorial(n,k).
#"""
#Pochhammer(n::Int, k::Int) = ∏(n, n + k - 1)
# *** deprecated ***

"""

Return the falling factorial which is the product of ``i`` ∈ ``(n - k + 1):n``.
"""
FallingFactorial(n::Int, k::Int) = ∏(n - k + 1, n)

"""

Return the falling factorial which is the product of ``i`` ∈ ``(n - k + 1):n``. A convenient infix syntax for the falling factorial is n ↓ k.
"""
↓(n, k) = FallingFactorial(n, k)

"""

Return the number of permutations of n letters, ``n! = ∏(1, n)``, the factorial of ``n``. (The notation is a shortcut breaking Julia conventions.)
"""
F!(n) = Nemo.factorial(ZZ(n))

"""

Return the factorial numbers.
"""
V000142(n::Int) = F!(n)

"""

Return the central rising factorial ``(n+1) ↑ (n+1) = (2n+1)! / n!``.
"""
V000407(n::Int) = (n + 1)↑(n + 1)
# V000407(n::Int) = ∏(n + 1, 2n + 1)

"""

Return the restricted rising factorial which is zero for ``n < 0`` or ``k > n``.
"""
V124320(n::Int, k::Int) = (n < 0 || k > n) ? 0 : ∏(n, n + k - 1)

"""

Return the multi-factorial which is the function ``n → ∏(a + b, a(n-1) + b)``
"""
MultiFactorial(a::Int, b::Int) = n -> ∏([a * k + b for k ∈ 0:(n-1)])

"""

Return the double factorial of odd numbers, ``1×3×5×...×(2n-1) = (2n-1)!!``.
"""
V001147(n::Int) = MultiFactorial(2, 1)(n)

"""

Return the double factorial of even numbers: ``2^n n! = (2n)!!``.
"""
V000165(n::Int) = MultiFactorial(2, 2)(n)

"""

Return the triple factorial numbers with shift 1, ``3^n n! = (3n)!!!``.
"""
V007559(n::Int) = MultiFactorial(3, 1)(n)

"""

Return the triple factorial numbers with shift 2.
"""
V008544(n::Int) = MultiFactorial(3, 2)(n)

"""

Return the triple factorial numbers with shift 3.
"""
V032031(n::Int) = MultiFactorial(3, 3)(n)

"""

Return the quadruple factorial numbers with shift 1.
"""
V007696(n::Int) = MultiFactorial(4, 1)(n)

"""

Return the quadruple factorial numbers with shift 2, ``(2n)!/n!``.
"""
V001813(n::Int) = MultiFactorial(4, 2)(n) # = ∏(n + 1, 2n)

"""

Return the quadruple factorial numbers with shift 3.
"""
V008545(n::Int) = MultiFactorial(4, 3)(n)

"""

Return the quadruple factorial numbers ``4^n n!``.
"""
V047053(n::Int) = MultiFactorial(4, 4)(n)

# RisingFactorial(n, k)
# n\k [0  1   2    3     4      5        6         7          8]
# --------------------------------------------------------------
# [0] [1, 0,  0,   0,    0,     0,       0,        0,         0]
# [1] [1, 1,  2,   6,   24,   120,     720,     5040,     40320]
# [2] [1, 2,  6,  24,  120,   720,    5040,    40320,    362880]
# [3] [1, 3, 12,  60,  360,  2520,   20160,   181440,   1814400]
# [4] [1, 4, 20, 120,  840,  6720,   60480,   604800,   6652800]
# [5] [1, 5, 30, 210, 1680, 15120,  151200,  1663200,  19958400]
# [6] [1, 6, 42, 336, 3024, 30240,  332640,  3991680,  51891840]
# [7] [1, 7, 56, 504, 5040, 55440,  665280,  8648640, 121080960]
# [8] [1, 8, 72, 720, 7920, 95040, 1235520, 17297280, 259459200]

#START-TEST-########################################################

using Test, SeqTests

function test()
    @testset "FallingFact" begin
        @test FallingFactorial(100, 100) == factorial(BigInt(100))
        @test (100↓100) == factorial(BigInt(100))
        @test FallingFactorial(333, 333) == factorial(BigInt(333))
        @test FallingFactorial(111, 0) == 1
    end
    @testset "RisingFact" begin
        @test RisingFactorial(11, 11) == 14079294028800
        @test (11↑11) == 14079294028800
        @test RisingFactorial(
            33,
            33
        ) == 31344295059422473624824839739793024055460338073600000000
        @test RisingFactorial(111, 0) == 1
    end
    @testset "MultiFact" begin
        a = [MultiFactorial(2, 1)(n) for n ∈ 0:6]
        b = [1, 1, 3, 15, 105, 945, 10395]
        @test all(a .== b)
    end

    if data_installed()
        V = [
            V000142,
            V000165,
            V007696,
            V001813,
            V047053,
            V001147,
            V008545,
            V081125,
            V000407,
            V032031,
            V007559,
            V008544
        ]

        @testset "Products" begin
            for v ∈ V
                SeqTest(v, 'V')
            end
        end
    end
end

function demo()
    for n ∈ 0:6
        println([RisingFactorial(n, k) for k ∈ 0:7])
    end

    for n ∈ 0:6
        println([FallingFactorial(n, k) for k ∈ 0:7])
    end

    println(∏([FallingFactorial(n, n) for n ∈ 0:12]))
    println(Product([FallingFactorial(n, n) for n ∈ 0:12]))

    for n ∈ 0:4, k ∈ 0:4
        println(n, k, " → ", n↑k, " ", n↓k)
    end
    println()
end

"""

for n ∈ 1:10000 F!(n) end
    1.480517 seconds (10.00 k allocations: 156.250 KiB)
for n ∈ 1:1000 A000407(n) end
    0.177942 seconds (997.77 k allocations: 15.225 MiB)
for n ∈ 1:200, k ∈ 1:200 Pochhammer(n, k) end
    1.784024 seconds (4.43 M allocations: 67.596 MiB, 33.35% gc time)
"""
function perf()
    GC.gc()
    @time (for n ∈ 1:10000
        F!(n)
    end)
    GC.gc()
    @time (for n ∈ 1:1000
        V000407(n)
    end)
    GC.gc()
    @time (for n ∈ 1:200, k ∈ 1:200
        RisingFactorial(n, k)
    end)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

#=

Test Summary: | Pass  Total
FallingFact   |    4      4
Test Summary: | Pass  Total
RisingFact    |    4      4
Test Summary: | Pass  Total
MultiFact     |    1      1

fmpz[1, 0, 0, 0, 0, 0, 0, 0]
fmpz[1, 1, 2, 6, 24, 120, 720, 5040]
fmpz[1, 2, 6, 24, 120, 720, 5040, 40320]
fmpz[1, 3, 12, 60, 360, 2520, 20160, 181440]
fmpz[1, 4, 20, 120, 840, 6720, 60480, 604800]
fmpz[1, 5, 30, 210, 1680, 15120, 151200, 1663200]
fmpz[1, 6, 42, 336, 3024, 30240, 332640, 3991680]
fmpz[1, 0, 0, 0, 0, 0, 0, 0]
fmpz[1, 1, 0, 0, 0, 0, 0, 0]
fmpz[1, 2, 2, 0, 0, 0, 0, 0]
fmpz[1, 3, 6, 6, 0, 0, 0, 0]
fmpz[1, 4, 12, 24, 24, 0, 0, 0]
fmpz[1, 5, 20, 60, 120, 120, 0, 0]
fmpz[1, 6, 30, 120, 360, 720, 720, 0]

127313963299399416749559771247411200000000000
127313963299399416749559771247411200000000000

00 → 1 1
01 → 0 0
02 → 0 0
03 → 0 0
04 → 0 0ü
10 → 1 1
11 → 1 1
12 → 2 0
13 → 6 0
14 → 24 0
20 → 1 1
21 → 2 2
22 → 6 2
23 → 24 0
24 → 120 0
30 → 1 1
31 → 3 3
32 → 12 6
33 → 60 6
34 → 360 0
40 → 1 1
41 → 4 4
42 → 20 12
43 → 120 24
44 → 840 24

  1.469480 seconds (10.00 k allocations: 156.250 KiB)
  0.183021 seconds (997.77 k allocations: 15.225 MiB)
  1.953392 seconds (4.43 M allocations: 67.596 MiB, 31.24% gc time)

=#

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
# Performs cross validation using CatBoost's built-in CV system via `cv`
module CrossValidation

using CatBoost
using DataFrames

cv_data = [["France", 1924, 44], ["USA", 1932, 37], ["Switzerland", 1928, 25],
           ["Norway", 1952, 30], ["Japan", 1972, 35], ["Mexico", 1968, 112]]

labels = [1, 1, 0, 0, 0, 1]

cat_features = [0]

cv_dataset = Pool(; data=cv_data, label=labels, cat_features=cat_features)

params = Dict("iterations" => 100, "depth" => 2, "loss_function" => "Logloss",
              "verbose" => false)

scores = cv(cv_dataset; fold_count=2, params)

display(scores)

end # module