text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Spot_julia"
version = v"2.9.7"
julia_version = v"1.6.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://www.lrde.epita.fr/dload/spot/spot-2.9.7.tar.gz","1eea67e3446cdbbbb705ee6e26fd869020cdb7d82c563fead9cb4394b9baa04c"),
    GitSource("https://github.com/MaximeBouton/spot_julia.git", "6ffcf4b64f64fc9e3363db22f4cc57a957d28128")
    ]

    # Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/spot-2.9.7/
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-python
make -j${nproc}
make install

# build with cmake
cd $WORKSPACE/srcdir/spot_julia/spot_julia

# Override compiler ID to silence the horrible "No features found" cmake error
if [[ $target == *"apple-darwin"* ]]; then
  macos_extra_flags="-DCMAKE_CXX_COMPILER_ID=AppleClang -DCMAKE_CXX_COMPILER_VERSION=10.0.0 -DCMAKE_CXX_STANDARD_COMPUTED_DEFAULT=11"
fi
Julia_PREFIX=$prefix
mkdir build
cd build
cmake -DJulia_PREFIX=$Julia_PREFIX \
    -DCMAKE_SPOT_LIB_DIR=${libdir} \
    -DCMAKE_SPOT_INCLUDE_DIR=${includedir} \
    -DCMAKE_FIND_ROOT_PATH=$prefix \
    -DJlCxx_DIR=${libdir}/cmake/JlCxx \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    ${macos_extra_flags} \
    ${windows_extra_flags} \
    -DCMAKE_BUILD_TYPE=Release \
    ..
VERBOSE=ON cmake --build . --config Release --target install -- -j${nproc}
install_license $WORKSPACE/srcdir/spot_julia/LICENSE.md
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc = "glibc"),
    Platform("x86_64", "windows"),
    Platform("x86_64", "macos")
]
platforms = expand_cxxstring_abis(platforms)

# # uncomment when pushing to yggdrasil
# # These are the platforms we will build for by default, unless further
# # platforms are passed in on the command line
# include("../../L/libjulia/common.jl")
# platforms = libjulia_platforms(julia_version)
# platforms = filter!(!Sys.iswindows, platforms) # Singular does not support Windows
# platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("genltl", :genltl),
    LibraryProduct("libspot", :libspot),
    ExecutableProduct("ltl2tgta", :ltl2tgta),
    ExecutableProduct("ltlsynt", :ltlsynt),
    ExecutableProduct("ltlcross", :ltlcross),
    LibraryProduct("libspotgen", :libspotgen),
    ExecutableProduct("autcross", :autcross),
    ExecutableProduct("genaut", :genaut),
    ExecutableProduct("ltl2tgba", :ltl2tgba),
    ExecutableProduct("randaut", :randaut),
    ExecutableProduct("autfilt", :autfilt),
    ExecutableProduct("ltlfilt", :ltlfilt),
    ExecutableProduct("ltlgrind", :ltlgrind),
    ExecutableProduct("ltldo", :ltldo),
    LibraryProduct("libbddx", :libbddx),
    LibraryProduct("libspotltsmin", :libspotltsmin),
    ExecutableProduct("randltl", :randltl),
    ExecutableProduct("dstar2tgba", :dstar2tgba),
    LibraryProduct("libspot_julia", :libspot_julia)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("libcxxwrap_julia_jll"),
    BuildDependency(PackageSpec(name="libjulia_jll", version=julia_version))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"9")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    tensordecomp(k::Int, d::Int)

Decomposition of the integers  0:(k^p - 1) in powers of k.

"""
function tensordecomp(k::Int, d::Int)

    sequences = zeros(Int, k^d, d)

    for n = 0:k^d-1
        i = d
        m = n

        while i > 0
            i = i - 1
            j = i
            f = floor(Int, m / k^i)

            while j > 0 && f == 0
                j = j - 1
                f = floor(Integer, m / k^j)
            end

            if f > 0
                sequences[n + 1, j + 1] = f
                i = j
            elseif f == 0
                sequences[n + 1, 1] = m
                i = 0
            end
            m = m - f * k^i
        end
    end

    return sequences
end

"""
    even_sampling_rules(dim::Int, split_factor::Int) -> Array{Float64, 2}

Generate rules for evenly distributed points within a simplex. To do this,
we perform a shape-preserving splitting of the simplex, given a splitting
factor.

Returns the convex expansion coefficients of the points of the resulting
subsimplices in terms of the vertices of the original simplex.
"""
function even_sampling_rules(dim::Int, split_factor::Int)

    sequences::Array{Int, 2} = tensordecomp(split_factor, dim)
    n_seq = size(sequences, 1)

    χ1 = sequences .* (dim + 1)
    χ2 = repeat(transpose(collect(1:dim)), n_seq, 1)
    χ::Array{Int, 2} = χ1 .+ χ2
    χ = sort(χ, dims=2)

    # Define multiplicity matrix M
    M = zeros(Float64, size(χ, 1), size(χ, 2) + 1)
    M[:, 1] = χ[:, 1]
    M[:, 2:(end - 1)] = χ[:, 2:end] - χ[:, 1:(end - 1)]
    M[:, end] = (dim+1)*split_factor * ones(size(χ, 1)) - χ[:, end]

    M = M ./ (split_factor * (dim + 1))

    return copy(transpose(M))
end

"""
Evenly sample points within a simplex by performing a shape-preserving
subdivision of the simplex with a given `split_factor`. If the simplex
lives in a space of dimension `dim`, the resulting number of points is
`split_factor`^(dim).
"""
function evenly_sample(simplex::AbstractArray{Float64, 2}, split_factor::Int)
    dim = size(simplex, 2)
    centroids_exp_coeffs = copy(transpose(even_sampling_rules(dim, split_factor)))
    centroids_exp_coeffs * simplex
end

ext: jl


text:
import POSIXClock

function main()
  interval = 100000 # 100us

  t = POSIXClock.timespec(0,0)
  POSIXClock.gettime!(t,POSIXClock.CLOCK_MONOTONIC)
  @show t

  n = 0
  # Lock future memory allocations, disable GC
  ccall(:mlockall, Cint, (Cint,), 2) # CAUTION: will crash Julia on future memory allocations
  gc_enable(false)

  @time while n <= div(5000000000,interval)
    n+=1
    POSIXClock.nanosleep!(t,interval)
  end # time should be around 5 seconds

  # Unlock memory allocations, enable GC
  ccall(:munlockall, Cint, ())
  gc_enable(true)

  return nothing
end

main()


text:
module FastFilter

@static isdefined(Core.Compiler, :AbstractInterpreter) || error("FastFilter.jl only works with Julia versions 1.6 and higher")

# imports
# =======

import Base:
    uniontypes,
    get_world_counter,
    _methods_by_ftype

import Core:
    Const,
    MethodMatch,
    MethodInstance

import Core.Compiler:
    NativeInterpreter,
    specialize_method,
    InferenceResult,
    InferenceState,
    typeinf,
    widenconst

# utils
# =====

function infer(@nospecialize(tt),
               world = get_world_counter(),
               interp = NativeInterpreter(world),
               )
    mms = _methods_by_ftype(tt, -1, world)
    length(mms) === 1 || return nothing

    mm = first(mms)::MethodMatch

    linfo = specialize_method(mm.method, mm.spec_types, mm.sparams)
    result = InferenceResult(linfo)
    frame = InferenceState(result, #=cached=# true, interp)

    typeinf(interp, frame)

    return result
end

# filter
# ======

function fastfilter(@nospecialize(f), a::Array{T,N}) where {T,N}
    if @generated
        ft = f
        et = T
        isa(et, Union) || return :(filter(f, a))
        ets = uniontypes(et)

        filtered = []
        for et in ets
            tt = Tuple{ft,et}
            result = infer(tt)
            if isa(result, InferenceResult)
                result = result.result
                isa(result, Const) && result.val === false && continue
            end
            push!(filtered, et)
        end
        et′ = Union{filtered...}

        return quote
            j = 1
            b = Vector{$(et′)}(undef, length(a))
            isempty(b) && return b
            fallback = first(b)
            for ai in a
                c = f(ai)
                @inbounds b[j] = (c ? ai : fallback)::$(et′)
                j = ifelse(c, j+1, j)
            end
            resize!(b, j-1)
            sizehint!(b, length(b))
            b
        end
    else
        # fallback when there're not enough type information for code generation
        return filter(f, a)
    end
end

function fastfilter(@nospecialize(pred), s::AbstractSet)
    if @generated
        ft = pred
        et = eltype(s)
        isa(et, Union) || return :(filter(pred, s))
        ets = uniontypes(et)

        filtered = []
        for et in ets
            tt = Tuple{ft,et}
            result = infer(tt)
            if isa(result, InferenceResult)
                result = result.result
                isa(result, Const) && result.val === false && continue
            end
            push!(filtered, et)
        end
        et′ = Union{filtered...}

        return :(Base.mapfilter(pred, push!, s, Set{$(et′)}()))
    else
        return filter(pred, s)
    end
end

function fastreplace(@nospecialize(new::Base.Callable), A; count = nothing)
    if @generated
        # when all elements in `A` aren't guaranteed to be replaced, just fallback
        cnt = count
        count === Nothing || return :(replace(new, A; count))
        fallback = :(replace(new, A; count = typemax(Int)))

        ft = new
        et = eltype(A)
        isa(et, Union) || return fallback
        ets = uniontypes(et)

        replaced = []
        for et in ets
            tt = Tuple{ft,et}
            result = infer(tt)
            isa(result, InferenceResult) || return fallback # unsuccessful inference
            push!(replaced, widenconst(result.result))
        end
        et′ = Union{replaced...}

        return :(Base._replace!(new, Base._similar_or_copy(A, $(et′)), A, Base.check_count(typemax(Int))))
    else
        if isnothing(count)
            count = typemax(Int)
        end
        return replace(new, A; count)
    end
end

export
    fastfilter,
    fastreplace

end


text:
@snn_kw struct NoisyIFParameter{FT=Float32} <: AbstractIFParameter
    σ::FT = 0
end

@snn_kw mutable struct NoisyIF{VFT=Vector{Float32}} <: AbstractIF
    param::NoisyIFParameter = NoisyIFParameter()
    N::Int32 = 100
    randncache::VFT = randn(N)
end

"""
Noisy Integrate-And-Fire Neuron
"""
NoisyIF

function integrate!(p::NoisyIF, param::NoisyIFParameter, dt::Float32)
    @unpack N, randncache = p
    @unpack σ = param
    randn!(randncache)
    @inbounds for i = 1:N
        v[i] += dt * (ge[i] + gi[i] - (v[i] - El) + I[i] + σ / √dt * randncache[i]) / τm
        ge[i] += dt * -ge[i] / τe
        gi[i] += dt * -gi[i] / τi
    end
    @inbounds for i = 1:N
        fire[i] = v[i] > Vt
        v[i] = ifelse(fire[i], Vr, v[i])
    end
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tuple{typeof(get_rebugger_modeswitch_dict)})
    precompile(Tuple{typeof(rebugrepl_init)})
end


text:
using GenericLinearAlgebra, ControlSystems, Optim, LineSearches, ForwardDiff, Plots, Printf
include("expm.jl") # Native matrix exponential

struct IdConditions
    θ0::Array{Float64,1}    # [Rp,C,Rc,L]
    x00::Array{Float64,1}
    θid::Array{Int,1}       # e.g. [1,1,1,0]
    x0id::Bool              # identify x0
    x0solve::Bool # solve for x0 using its preseriodicity at every optimization step
    G_::Function # discrete time transfer function, of form `G_(θ,h;discretize=true)`
    IdConditions(θ0,x00,θid,x0id,x0solve,G_) = (length(θ0) == length(θid)) ? new(θ0,x00,θid,x0id,x0solve,G_) : error("Length of `θ0` != `θid`.")
    # FIXME add check on G_ parameter expectation
end

IdConditions(θ0, x00, θid, x0id, x0solve) = IdConditions(θ0,x00,θid,x0id,x0solve,G) # default to 4-element windkessel model, G

"""
Identifies subset of {system parameters θ, initial state x0}.
"""
function identify(θ0,x00,u,y,h;θid=[1,1,1,1],x0id=true,x0solve=false,G_=G,θpos=true,v=false)
    id = IdConditions(θ0,x00,θid,x0id,x0solve,G_)
    p0 = p(θ0,x00,id)
    f(p) = cost(p,id,u,y,h)
    g(f) = (grad,p) -> ForwardDiff.gradient!(grad,f,p)
    g! = g(f)

    function gethess(f)
        return (hess,p) -> begin
            ForwardDiff.hessian!(hess,f,p)
            for c in 1:length(p)
               for r in 1:(c-1)
                    hess[r,c] = (hess[r,c]+hess[c,r])/2  # Enforces hess >= 0
                    hess[c,r] = hess[r,c]
               end
            end
        end
    end
    h! = gethess(f)

    df = TwiceDifferentiable(f,g!,h!,p0)
    opts = Optim.Options(store_trace = true, extended_trace = true, show_trace = v)
    local res
    if θpos
        up = Float64[]
        lp = vcat(zeros(+(id.θid...)),fill(-Inf,length(x00)))
        dfc = TwiceDifferentiableConstraints(lp, up)
        res = optimize(df,dfc,p0,IPNewton(),opts)
    else
        res = optimize(df,p0,Newton(alphaguess = LineSearches.InitialPrevious(), linesearch = LineSearches.BackTracking()),opts)
    end

    phat = Optim.minimizer(res)
    return phat,res
end

"""
Identifies subset of {system parameters θ, initial state x0} using NelderMead.
"""
function identifyglobal(θ0,x00,u,y,h;θid=[1,1,1,1],x0id=true,x0solve=false,v=false)
    id = IdConditions(θ0,x00,θid,x0id,x0solve)
    p0 = p(θ0,x00,id)
    f(p) = cost(p,id,u,y,h)
    opts = Optim.Options(iterations=Int(1e5),store_trace = true, extended_trace = true, show_trace = v)
    res = optimize(f, p0, NelderMead(),opts) # FIXME: constrain to θ > 0

    phat = Optim.minimizer(res)
    return phat,res
end

"""
Returns the cost (to be minimized).
"""
function cost(p,id::IdConditions,u,y,h)
    θ,x0 = deal(p,id)
    if id.x0solve x0 = solvex0(id.G_(θ,h),u) end
    ym, = lsim(id.G_(θ,h),u,h.*(1:length(u)),x0=x0)
    return dot(y-ym,y-ym)
end

"""
Simulate discrete-time system G(θ,h) with input vector u, starting in initial condition x0.
"""
function simulate(θ,x0,u,h; xout=false)
    y, _, x = lsim(G(θ,h),u,h.*(1:length(u)),x0=x0)
    return xout ? (y, x) : y
end

"""
Generates parallel 4-element Windkessel instance with parameter values of θ=[Rp,C,Rc,L], ZOH discretized with period h.
"""
function G(θ,h;discretize=true)
    Rp,C,Rc,L = tuple(θ...)
    # Parallel form as used by Burattini, Westerhof, Stergiopulos
    # s = tf("s")
    # Gc = ss((s*L*Rc)/(s*L+Rc) + Rp/(1+s*C*Rp))

    # With x1 equal to the charge (volume) on C.
    A = [-1/(C*Rp) 0; 0 -Rc/L]
    B = [1; Rc]
    C = [1/C -Rc/L]
    D = [Rc]
    Gc = ss(A,B,C,D)
    G = c2d(Gc,h,:zoh)
    return discretize ? G : Gc
end

"""
Generates 3-element Windkessel instance with parameter values of θ=[Rp,C,Rc], ZOH discretized with period h. Corresponds to 4 element with L = Inf.
"""
function G3(θ,h;discretize=true)
    Rp,C,Rc = tuple(θ...)
    # With x1 equal to the charge (volume) on C.
    A = [-1/(C*Rp)]
    B = [1]
    C = [1/C]
    D = [Rc]
    Gc = ss(A,B,C,D)
    G = c2d(Gc,h,:zoh)
    return discretize ? G : Gc
end

"""
Generates 2-element Windkessel instance with parameter values of θ=[Rp,C], ZOH discretized with period h. Corresponds to element with L = 0.
"""
function G2(θ,h;discretize=true)
    Rp,C = tuple(θ...)
    # With x1 equal to the charge (volume) on C.
    A = [-1/(C*Rp)]
    B = [1]
    C = [1/C]
    D = [0]
    Gc = ss(A,B,C,D)
    G = c2d(Gc,h,:zoh)
    return discretize ? G : Gc
end

"""
Build parameter vector for optimization
"""
function p(θ,x0,id::IdConditions)
    p = θ[findall(Bool.(id.θid))]
    if id.x0id append!(p,x0) end
    return Float64.(p)
end

"""
Retrieve original parametrization.
"""
function deal(p,id::IdConditions)
    pidx = cumsum(id.θid)
    θ = [Bool(id_) ? p_ : θ0_ for (p_,θ0_,id_) in zip(p[pidx],id.θ0,id.θid)]
    x0 = id.x0id ? p[sum(id.θid).+(1:length(id.x00))] : id.x00
    return θ,x0
end

"""
Parameters from stergiopulos1999total.
"""
function params_stergiopulos1999total()
    Rp = .79                # peripheral (systemic) resistance [mmHg*ml^-1*s]
    C = 1.22                # total arterial compliance [ml/mmHg]
    Rc = .056               # characeristic aoritic resistance
    L = .0051               # total arterial impedance [mmHg*s^2*ml^-1]
    θ = [Rp,C,Rc,L]
    x0 = [8.066, -8.296]    # taken from steady-state sim results
    return θ, x0
end

"""
Solve for x0 by leveraging the periodicity of the input signal.
"""
function solvex0(G, u::Array{Float64,1})
    n = length(u)
    un = u[end]
    _, _, x = lsim(G,u)
    xn0 = x[end,:]

    x0 = (I-G.A^n)\(G.A^(n-1)*G.B*un+xn0)
    return x0[:]
end

max_stars_repo_name: hpigot/windkessel-id

max_issues_repo_name: hpigot/windkessel-id

max_forks_repo_name: hpigot/windkessel-id


text:
using GeophysicalModelGenerator
using Test

@testset "Data import.jl" begin
    include("test_data_import.jl")
end
@testset "Data types.jl" begin
    include("test_data_types.jl")
end
@testset "Paraview" begin
    include("test_paraview.jl")
end
@testset "Gravity model" begin
    include("test_voxel_gravity.jl")
end
@testset "Utils" begin
    include("test_utils.jl")
end
@testset "Transformations" begin
    include("test_transformation.jl")
end

@testset "LaMEM" begin
    include("test_lamem.jl")
end

@testset "STL" begin
    include("test_stl.jl")
end

# Cleanup
foreach(rm, filter(endswith(".vts"), readdir()))
foreach(rm, filter(endswith(".vtu"), readdir()))
rm("./markers/",recursive=true)


text:
using JuliaConPackage
using Test

@testset "JuliaConPackage.jl" begin
    # Write your tests here.
end


text:
"""
file tfw
"""
function save_tfw(output::String, GSD::Float64, lx::Float64, uy::Float64)
    fname = splitext(output)[1]
    io = open(fname * ".tfw", "w")
    write(io, "$(Float64(GSD))\n")
    write(io, "0.000000000000000\n")
    write(io, "0.000000000000000\n")
    write(io, "-$(Float64(GSD))\n")
    L = @sprintf("%f", lx)
    U = @sprintf("%f", uy)
    write(io, "$L\n")
    write(io, "$U\n")
    close(io)
end

"""
Save point cloud extracted file .las.
"""
function savepointcloud(params::ParametersOrthophoto, temp::String)

    flushprintln("Point cloud: saving ...")

    # update number of points in header
    params.mainHeader.records_count = params.numPointsProcessed

    # update header bounding box
    flushprintln("Point cloud: update bbox ...")
    params.mainHeader.x_min = params.tightBB.x_min
    params.mainHeader.y_min = params.tightBB.y_min
    params.mainHeader.z_min = params.tightBB.z_min
    params.mainHeader.x_max = params.tightBB.x_max
    params.mainHeader.y_max = params.tightBB.y_max
    params.mainHeader.z_max = params.tightBB.z_max

    # write las file
    pointtype = LasIO.pointformat(params.mainHeader) # las point format

    if params.numPointsProcessed != 0 # if n == 0 nothing to save
        # in temp : list of las point records
        open(temp) do s
            # write las
            open(params.outputfile, "w") do t
                write(t, LasIO.magic(LasIO.format"LAS"))
                write(t, params.mainHeader)

                # LasIO.skiplasf(s)
                for i = 1:params.numPointsProcessed
                    p = read(s, pointtype)
                    write(t, p)
                end
            end
        end
    end

    rm(temp) # remove temp
    flushprintln("Point cloud: done ...")
end

"""
Save orthoprojection image.
"""
function saveimage(params::ParametersOrthophoto)
    flushprintln("Image: saving ...")
    # save tfw
    save_tfw(params.outputimage, params.GSD, params.refX, params.refY)

    # save image
    save(params.outputimage, Images.colorview(RGB, params.RGBtensor))

    flushprintln("Image: done ...")
end

function newPointRecord(
    point::Point,
    type::LasIO.DataType,
    header::LasIO.LasHeader;
    affineMatrix = Matrix(Common.I, 4, 4),
)

    position = Common.apply_matrix(affineMatrix, point.position)
    x = LasIO.xcoord(position[1], header)
    y = LasIO.ycoord(position[2], header)
    z = LasIO.zcoord(position[3], header)
    intensity = UInt16(0)
    flag_byte = UInt8(0)
    raw_classification = UInt8(0)
    scan_angle = Int8(0)
    user_data = UInt8(0)
    pt_src_id = UInt16(0)

    if type == LasIO.LasPoint0
        return type(
            x,
            y,
            z,
            intensity,
            flag_byte,
            raw_classification,
            scan_angle,
            user_data,
            pt_src_id,
        )

    elseif type == LasIO.LasPoint1
        gps_time = Float64(0)
        return type(
            x,
            y,
            z,
            intensity,
            flag_byte,
            raw_classification,
            scan_angle,
            user_data,
            pt_src_id,
            gps_time,
        )

    elseif type == LasIO.LasPoint2
        red = point.color[1]
        green = point.color[2]
        blue = point.color[3]
        return type(
            x,
            y,
            z,
            intensity,
            flag_byte,
            raw_classification,
            scan_angle,
            user_data,
            pt_src_id,
            red,
            green,
            blue,
        )

    elseif type == LasIO.LasPoint3
        gps_time = Float64(0)
        red = rgb[1]
        green = rgb[2]
        blue = rgb[3]
        return type(
            x,
            y,
            z,
            intensity,
            flag_byte,
            raw_classification,
            scan_angle,
            user_data,
            pt_src_id,
            gps_time,
            red,
            green,
            blue,
        )

    end

end


text:
function neighquadlets(x, q, quadletlist, quadletnum, atomtype, ilist, alist)

    dim = size(x,1);
    ncq = size(q,1);
    inum = length(ilist);
    ijklnum = sum(quadletnum(ilist+1)-quadletnum(ilist));

    ai = zeros(Int64,ijklnum);
    aj = zeros(Int64,ijklnum);
    ak = zeros(Int64,ijklnum);
    al = zeros(Int64,ijklnum);
    ti = zeros(Int64,ijklnum);
    tj = zeros(Int64,ijklnum);
    tk = zeros(Int64,ijklnum);
    tl = zeros(Int64,ijklnum);
    xij = zeros(dim,ijklnum);
    xik = zeros(dim,ijklnum);
    xil = zeros(dim,ijklnum);
    qi = zeros(ncq,ijklnum);
    qj = zeros(ncq,ijklnum);
    qk = zeros(ncq,ijklnum);
    ql = zeros(ncq,ijklnum);

    count = 0;
    for ii=1:inum
      i = ilist[ii]; # atom i
      n1 = quadletnum[i];
      n2 = quadletnum[i+1];
      n = n2-n1;
      gj = quadletlist((n1+1):n2,1); # a list of neighbors around atom i
      j = alist[gj];
      gk = quadletlist((n1+1):n2,2); # a list of neighbors around atom i
      k = alist[gk];
      gl = quadletlist((n1+1):n2,3); # a list of neighbors around atom i
      l = alist[gl];

      ind = (count+1):(count+n);
      count = count + n;

      ai[ind] = i;
      aj[ind] = j;
      ak[ind] = k;
      al[ind] = l;
      ti[ind] = atomtype[i];
      tj[ind] = atomtype[j];
      tk[ind] = atomtype[k];
      tl[ind] = atomtype[l];
      xij[:,ind] = x[:,gj] - x[:,i];
      xik[:,ind] = x[:,gk] - x[:,i];
      xil[:,ind] = x[:,gl] - x[:,i];
      if ncq>0
          qi[ind] = q[i];
          qj[ind] = q[j];
          qk[ind] = q[k];
          ql[ind] = q[l];
      end
    end

    return xij, xik, xil, ai, aj, ak, al, ti, tj, tk, tl, qi, qj, qk, ql

    end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using LinearAlgebra
using Plots
using Images
using ImageView

n = 1000
X = zeros(n, n)
X[250:750, 250:750] .= 1
U, S, V = svd(X)
imshow(X)
p1 = plot(S, yaxis = :log)

nAngles = 12 # sweep through 12 angles, from 0:4:44

for j = 2:nAngles
    Y = imrotate(X, -(j-1)*4*pi/180, axes(X)) # rotate by (j-1)*4
    Y[isnan.(Y)] .= 0
    U2, S2, V2 = svd(Y)
    plot!(S2)
end


text:
function RandomGraph_duplication_noise_experiments(testing=true)

    if testing
        default_n = [20]
        default_p = [0.5]
        default_sp = [.25]

        alphas = [.5]
        betas = [1.0]

        trial_count = 1
        parameter_sets = [
            ([50],default_p,default_sp),#size exp
            (default_n,[.9],default_sp),
            (default_n,default_p,[.25]),
        ]

        result_path = all_results_path*"/RG_DupNoise/test/"
    else
        #  --  Parameter Sets  --  #
        default_n = [250]
        default_p = [0.5]
        default_sp = [.25]

        seeded=true
        trial_count = 20
        parameter_sets = [
            (default_n,collect(0.0:.1:1.0),default_sp),
            (default_n,default_p,[0.05,0.1,.2,.25,.3,.4,.5]),
            ([100,250,500,1000,1250,1500],default_p,default_sp)#,#size exp
        ]
        alphas = [.5,1.0]
        betas = [0.0,1.0,1e1,1e2]
        result_path = all_results_path*"/RG_DupNoise/"
    end
    seeded=true
    if !isdir(result_path)
        mkpath(result_path)
    end

    test_set = 0
    for (n_sizes,edge_inclusion_p,step_percentage) in parameter_sets

        #Lambda TAME results
        #    Gram Matching KlauAlgo Post Processing
        X = distributed_random_trials(trial_count,DuplicationNoise(),seeded;postProcessing=KlauAlgo(rtype=1,maxiter=500),method=ΛTAME_M(),
                                    graph=RandomGeometric(),degreedist=LogNormal(log(5),1),matchingMethod=ΛTAME_GramMatching(),
                                    n_sizes, edge_inclusion_p,step_percentage,profile=true,alphas,betas)
        filename=filter(x-> !isspace(x),"LambdaTAME_graphType:RG_alphas:$(alphas)_betas:$(betas)_degdist:LogNormal-log5_n:$(n_sizes)_noiseModel:Duplication_p:$(edge_inclusion_p)_sp:$(step_percentage)_postProcess:KlauAlgo_trialcount:$(trial_count).json") #Bug: trialcount shortened because filename is too long
        save_results(result_path*filename,X)
        #    Gram Matching LocalSearch Post Processing
        X = distributed_random_trials(trial_count,DuplicationNoise(),seeded;postProcessing=LocalSearch(),method=ΛTAME_M(),
                                    graph=RandomGeometric(),degreedist=LogNormal(log(5),1),matchingMethod=ΛTAME_GramMatching(),
                                    n_sizes, edge_inclusion_p,step_percentage,profile=true,alphas,betas)
        filename=filter(x-> !isspace(x),"LambdaTAME_graphType:RG_alphas:$(alphas)_betas:$(betas)_degdist:LogNormal-log5_n:$(n_sizes)_noiseModel:Duplication_p:$(edge_inclusion_p)_sp:$(step_percentage)_postProcess:LocalSearch_trialcount:$(trial_count).json") #Bug: trialcount shortened because filename is too long
        save_results(result_path*filename,X)

        #=#   LambdaTAME - ROM  - KlauAlgo Post Processing
        X = distributed_random_trials(trial_count,DuplicationNoise(),seeded;postProcessing=KlauAlgo(),method=ΛTAME_M(),alphas,betas,
                                            graph=RandomGeometric(),degreedist=LogNormal(log(5),1),matchingMethod=ΛTAME_rankOneMatching(),
                                            n_sizes, edge_inclusion_p,step_percentage)
        filename=filter(x-> !isspace(x),"LambdaTAME_ROM_graphType:RG_alphas:$(alphas)_betas:$(betas)_degdist:LogNormal-log5_n:$(n_sizes)_noiseModel:Duplication_p:$(edge_inclusion_p)_sp:$(step_percentage)_postProcess:KlauAlgo_trialcount:$(trial_count).json") #Bug: trialcount shortened because filename is too long
        save_results(result_path*filename,X)
        =#

        #Low Rank TAME results
        #    Local Post Processing
        X = distributed_random_trials(trial_count,DuplicationNoise(),seeded;postProcessing=LocalSearch(),method=LowRankTAME_M(),
                                    graph=RandomGeometric(),degreedist=LogNormal(log(5),1), n_sizes, edge_inclusion_p,step_percentage,profile=true,alphas,betas)
        filename=filter(x-> !isspace(x),"LowRankTAME_graphType:RG_alphas:$(alphas)_betas:$(betas)_degdist:LogNormal-log5_n:$(n_sizes)_noiseModel:Duplication_p:$(edge_inclusion_p)_sp:$(step_percentage)_postProcess:LocalSearch_trialcount:$(trial_count).json") #Bug: trialcount shortened because filename is too long
        save_results(result_path*filename,X)
        #    KlauAlgo Post Processing
        X = distributed_random_trials(trial_count,DuplicationNoise(),seeded;postProcessing=KlauAlgo(),method=LowRankTAME_M(),
                                    graph=RandomGeometric(),degreedist=LogNormal(log(5),1), n_sizes, edge_inclusion_p,step_percentage,profile=true,alphas,betas)
        filename=filter(x-> !isspace(x),"LowRankTAME_graphType:RG_alphas:$(alphas)_betas:$(betas)_degdist:LogNormal-log5_n:$(n_sizes)_noiseModel:Duplication_p:$(edge_inclusion_p)_sp:$(step_percentage)_postProcess:KlauAlgo_trialcount:$(trial_count).json") #Bug: trialcount shortened because filename is too long
        save_results(result_path*filename,X)

        #= # low rank matching - KlauAlgo Post Processing
        X = distributed_random_trials(trial_count,DuplicationNoise(),seeded;postProcessing=KlauAlgo(),method=LowRankTAME_M(),
                                    graph=RandomGeometric(),degreedist=LogNormal(log(5),1), n_sizes, edge_inclusion_p,step_percentage,profile=true,low_rank_matching=true,alphas,betas)
        filename=filter(x-> !isspace(x),"LowRankTAME-lrm_graphType:RG_alphas:$(alphas)_betas:$(betas)_degdist:LogNormal-log5_n:$(n_sizes)_noiseModel:Duplication_p:$(edge_inclusion_p)_sp:$(step_percentage)_postProcess:KlauAlgo_trialcount:$(trial_count).json") #Bug: trialcount shortened because filename is too long
        save_results(result_path*filename,X)
        #    low rank matching - LocalSearch Post Processing
        X = distributed_random_trials(trial_count,DuplicationNoise(),seeded;postProcessing=LocalSearch(),method=LowRankTAME_M(),
                                    graph=RandomGeometric(),degreedist=LogNormal(log(5),1), n_sizes, edge_inclusion_p,step_percentage,profile=true,low_rank_matching=true,alphas,betas)
        filename=filter(x-> !isspace(x),"LowRankTAME-lrm_graphType:RG_alphas:$(alphas)_betas:$(betas)_degdist:LogNormal-log5_n:$(n_sizes)_noiseModel:Duplication_p:$(edge_inclusion_p)_sp:$(step_percentage)_postProcess:LocalSearch_trialcount:$(trial_count).json") #Bug: trialcount shortened because filename is too long
        save_results(result_path*filename,X)
        =#

        X = distributed_random_trials(trial_count,DuplicationNoise(),seeded;postProcessing=LocalSearch(),method=LowRankEigenAlign_M(),
                                    graph=RandomGeometric(),degreedist=LogNormal(log(5),1), n_sizes, edge_inclusion_p,step_percentage,profile=true)
        filename=filter(x-> !isspace(x),"LowRankEigenAlign_graphType:RG_degdist:LogNormal-log5_n:$(n_sizes)_noiseModel:Duplication_p:$(edge_inclusion_p)_sp:$(step_percentage)_postProcess:LocalSearch_trialcount:$(trial_count).json")
        save_results(result_path*filename,X)

        X = distributed_random_trials(trial_count,DuplicationNoise(),seeded;postProcessing=KlauAlgo(),method=LowRankEigenAlign_M(),
                                    graph=RandomGeometric(),degreedist=LogNormal(log(5),1), n_sizes, edge_inclusion_p,step_percentage,profile=true)
        filename=filter(x-> !isspace(x),"LowRankEigenAlign_graphType:RG_degdist:LogNormal-log5_n:$(n_sizes)_noiseModel:Duplication_p:$(edge_inclusion_p)_sp:$(step_percentage)_postProcess:KlauAlgo_trialcount:$(trial_count).json")
        save_results(result_path*filename,X)

        println("finished parameter set $(test_set)")
        test_set += 1
    end
end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Given a start and end date, how many steps/periods are in between
guess(a::YearDate, b::YearDate, c::Year) = Int64(div(value(b - a), value(c)))
guess(a::YearDate, b::YearDate, c::Semester) = Int64(div(2 * value(b - a), value(c)))
guess(a::YearDate, b::YearDate, c::Quarter) = Int64(div(4 * value(b - a), value(c)))
guess(a::YearDate, b::YearDate, c::Month) = Int64(div(12 * value(b - a), value(c)))
guess(a::SemesterDate, b::SemesterDate, c::Year) = Int64(div(value(b - a), semesters(c)))
guess(a::SemesterDate, b::SemesterDate, c::Semester) = Int64(div(value(b - a), value(c)))
guess(a::SemesterDate, b::SemesterDate, c::Quarter) = Int64(div(2 * value(b - a), value(c)))
guess(a::SemesterDate, b::SemesterDate, c::Month) = Int64(div(6 * value(b - a), value(c)))
guess(a::QuarterDate, b::QuarterDate, c::Year) = Int64(div(value(b - a), quarters(c)))
guess(a::QuarterDate, b::QuarterDate, c::Semester) = Int64(div(value(b - a), quarters(c)))
guess(a::QuarterDate, b::QuarterDate, c::Quarter) = Int64(div(value(b - a), value(c)))
guess(a::QuarterDate, b::QuarterDate, c::Month) = Int64(div(3 * value(b - a), value(c)))
guess(a::MonthDate, b::MonthDate, c::Year) = Int64(div(value(b - a), months(c)))
guess(a::MonthDate, b::MonthDate, c::Semester) = Int64(div(value(b - a), months(c)))
guess(a::MonthDate, b::MonthDate, c::Quarter) = Int64(div(value(b - a), months(c)))
guess(a::MonthDate, b::MonthDate, c::Month) = Int64(div(value(b - a), value(c)))
guess(a::WeekDate, b::WeekDate, c::Year) = Int64(div(value(b - a), weeks(c)))
guess(a::WeekDate, b::WeekDate, c::Semester) = Int64(div(value(b - a), weeks(c)))
guess(a::WeekDate, b::WeekDate, c::Quarter) = Int64(div(value(b - a), weeks(c)))
guess(a::WeekDate, b::WeekDate, c::Month) = Int64(div(value(b - a), weeks(c)))
guess(a::WeekDate, b::WeekDate, c::Week) = Int64(div(value(b - a), value(c)))
guess(a::WeekDate, b::WeekDate, c::Day) = Int64(7 * div(value(b - a), value(c)))
guess(a::DayDate, b::DayDate, c::Year) = Int64(div(value(b - a), days(c)))
guess(a::DayDate, b::DayDate, c::Semester) = Int64(div(value(b - a), days(c)))
guess(a::DayDate, b::DayDate, c::Quarter) = Int64(div(value(b - a), days(c)))
guess(a::DayDate, b::DayDate, c::Month) = Int64(div(value(b - a), days(c)))
guess(a::DayDate, b::DayDate, c::Week) = Int64(div(value(b - a), days(c)))
guess(a::DayDate, b::DayDate, c::Day) = Int64(div(value(b - a), value(c)))
guess(a::UndatedDate, b::UndatedDate, c) = Int64(div(value(b - a), value(c)))
function len(a, b, c)
    lo, hi, st = min(a, b), max(a, b), abs(c)
    i = guess(a, b, c)
    v = lo + st * i
    prev = v  # Ensure `v` does not overflow
    while v <= hi && prev <= v
        prev = v
        v += st
        i += 1
    end
    return i - 1
end
Base.length(r::StepRange{<:SimpleDate}) =
    isempty(r) ? Int64(0) : len(r.start, r.stop, r.step) + 1
# Period ranges hook into Int64 overflow detection
#Base.length(r::StepRange{<:Dates.Period}) =
#    length(StepRange(value(r.start), value(r.step), value(r.stop)))

# Overload Base.steprange_last because `rem` is not overloaded for `TimeType`s
function Base.steprange_last(start::T, step, stop) where {T<:SimpleDate}
    if isa(step, AbstractFloat)
        throw(ArgumentError("StepRange should not be used with floating point"))
    end
    z = zero(step)
    step == z && throw(ArgumentError("step cannot be zero"))

    if stop == start
        last = stop
    else
        if (step > z) != (stop > start)
            last = Base.steprange_last_empty(start, step, stop)
        else
            diff = stop - start
            if (diff > zero(diff)) != (stop > start)
                throw(OverflowError())
            end
            remain = stop - (start + step * len(start, stop, step))
            last = stop - remain
        end
    end
    last
end

import Base.in
function in(x::T, r::StepRange{T}) where {T<:SimpleDate}
    n = len(first(r), x, step(r)) + 1
    n >= 1 && n <= length(r) && r[n] == x
end

Base.iterate(r::StepRange{<:SimpleDate}) =
    length(r) <= 0 ? nothing : (r.start, (length(r), 1))
Base.iterate(r::StepRange{<:SimpleDate}, (l, i)) =
    l <= i ? nothing : (r.start + r.step * i, (l, i + 1))

for p in (:Semester, :Undated)
    @eval begin
        +(x::$p, r::AbstractRange{<:SimpleDate}) = (x+first(r)):step(r):(x+last(r))
        +(r::AbstractRange{<:SimpleDate}, x::$p) = x + r
        -(r::AbstractRange{<:SimpleDate}, x::$p) = (first(r)-x):step(r):(last(r)-x)
        *(x::$p, r::AbstractRange{<:Real}) = (x*first(r)):(x*step(r)):(x*last(r))
        *(r::AbstractRange{<:Real}, x::$p) = x * r
        /(r::AbstractRange{<:P}, x::P) where {P<:$p} = (first(r)/x):(step(r)/x):(last(r)/x)
    end
end

# Combinations of types and periods for which the range step is regular
Base.RangeStepStyle(::Type{<:OrdinalRange{<:SimpleDate,<:Dates.FixedPeriod}}) =
    Base.RangeStepRegular()


text:
# Autogenerated wrapper script for ONNXRuntime_jll for aarch64-apple-darwin
export libonnxruntime

JLLWrappers.@generate_wrapper_header("ONNXRuntime")
JLLWrappers.@declare_library_product(libonnxruntime, "@rpath/libonnxruntime.1.10.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libonnxruntime,
        "lib/libonnxruntime.1.10.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using CategoricalArrays

# CategoricalArrays
v = CategoricalArray(repeat(["b", "a"], inner = 5))
@test isequal(rcopy(CategoricalArray,RObject(v)), v)
v = CategoricalArray(repeat(["b", "a"], inner = 5), ordered=true)
@test isequal(rcopy(CategoricalArray,RObject(v)), v)
@test CategoricalArrays.levels(rcopy(CategoricalArray,R"factor(c('c','a','a'))")) == ["a","c"]
@test CategoricalArrays.isordered(rcopy(CategoricalArray,R"ordered(c('c','a','a'))"))

a = Array{Union{String, Missing}}(repeat(["b", "a"], inner = 5))
a[repeat([true, false], outer = 5)] .= missing
v = CategoricalArray(a)
@test isequal(rcopy(CategoricalArray,RObject(v)), v)
v = CategoricalArray(a, ordered=true)
@test isequal(rcopy(CategoricalArray,RObject(v)), v)
@test CategoricalArrays.levels(rcopy(CategoricalArray,R"factor(c('c',NA,'a'))")) == ["a","c"]
@test CategoricalArrays.isordered(rcopy(CategoricalArray,R"ordered(c('c',NA,'a'))"))


text:
# countrystructure.jl - Definición de tipos contenedores para bases de datos de
# variaciones intermensuales del IPC VarCPIBase definidas en cpibase.jl

import Base: show, summary, convert, getindex, eltype

"""
    CountryStructure{N, T <: AbstractFloat}

Tipo abstracto que representa el conjunto de bases del IPC de un país.
"""
abstract type CountryStructure{N, T <: AbstractFloat} end

"""
    UniformCountryStructure{N, T, B} <: CountryStructure{N, T}

Estructura que representa el conjunto de bases del IPC de un país,
posee el campo `base`, que es una tupla de la estructura `VarCPIBase`. Todas
las bases deben tener el mismo tipo de índice base.
"""
struct UniformCountryStructure{N, T, B} <: CountryStructure{N, T}
    base::NTuple{N, VarCPIBase{T, B}}
end

# Este tipo se puede utilizar con datos cuyos primeros índices no sean todos 100
"""
    MixedCountryStructure{N, T} <: CountryStructure{N, T}

Estructura que representa el conjunto de bases del IPC de un país,
posee el campo `base`, que es una tupla de la estructura `VarCPIBase`, cada una
con su propio tipo de índices base B. Este tipo es una colección de un tipo abstracto.
"""
struct MixedCountryStructure{N, T} <: CountryStructure{N, T}
    base::NTuple{N, VarCPIBase{T, B} where B}
end

# Util methods for constructing CountryStructures
UniformCountryStructure(bases::Vararg{VarCPIBase{T, B}, N}) where {N, T, B} = UniformCountryStructure{N, T, B}(bases)
MixedCountryStructure(bases::Vararg{VarCPIBase}) = MixedCountryStructure(bases)

# Util methods from FullCPIBase -> CountryStructure
function UniformCountryStructure(fullbases::Vararg{FullCPIBase{T, B}, N}) where {N, T, B}
    varbases = map(VarCPIBase, fullbases)
    UniformCountryStructure{N, T, B}(varbases)
end
MixedCountryStructure(fullbases::Vararg{FullCPIBase}) = MixedCountryStructure(map(VarCPIBase, fullbases))

# Resumen y método para mostrar

function summary(io::IO, cst::CountryStructure)
    datestart, dateend = _formatdate.((first(cst.base).dates[begin], last(cst.base).dates[end]))
    print(io, typeof(cst), ": ", datestart, "-", dateend)
end

function show(io::IO, cst::CountryStructure)
    l = length(cst.base)
    suffix = l == 1 ? " base" : " bases"
    println(io, typeof(cst), " con ", l, suffix)
    for base in cst.base
        println(io, "└─→ ", sprint(summary, base))
    end
end
# possible chars '├', '└', '│', '─', '⋮'

# Conversión entre tipos de datos flotantes

# Este método crea una copia a través de los métodos de conversión de bases
function convert(::Type{T}, cst::CountryStructure) where {T <: AbstractFloat}
    # Convert each base to type T
    conv_b = convert.(T, cst.base)
    getunionalltype(cst)(conv_b)
end

# Funciones de acceso a bases del IPC

"""
    getindex(cst::CountryStructure, i::Int)

Devuelve la base número `i` de un contenedor `CountryStructure`.
"""
getindex(cst::CountryStructure, i::Int) = cst.base[i]

# Función de ayuda para obtener los índices que corresponden a una fecha
# específica de una base
function _base_index(cst, date, retfirst=true)
    for (b, base) in enumerate(cst.base)
        dates = base.dates
        dateindex = findfirst(dates .== date)
        if !isnothing(dateindex)
            return b, dateindex
        end
    end
    # return first or last
    if retfirst
        1, 1
    else
        length(cst.base), size(cst.base[end].v, 1)
    end
end

"""
    getindex(cst::CountryStructure, startdate::Date, finaldate::Date)

Devuelve una copia del `CountryStructure` con las bases modificadas para tener
observaciones entre las fechas indicada por `startdate` y `finaldate`.
"""
function getindex(cst::CountryStructure, startdate::Date, finaldate::Date)

    # Obtener base y fila de inicio
    start_base, start_index = _base_index(cst, startdate, true)
    final_base, final_index = _base_index(cst, finaldate, false)

    bases = deepcopy(cst.base[start_base:final_base])
    if start_base == final_base
        # copy same base and slice
        @debug "Fechas en la misma base"
        # @info bases[1]
        onlybase = bases[1]
        newbase = VarCPIBase(
            onlybase.v[start_index:final_index, :],
            copy(onlybase.w), onlybase.dates[start_index:final_index], copy(onlybase.baseindex))

        return getunionalltype(cst)(newbase)
    else
        # different bases
        @debug "Fechas en diferentes bases"
        firstbase = first(bases)
        lastbase = last(bases)
        newstart = VarCPIBase(
            firstbase.v[start_index:end, :],
            copy(firstbase.w), firstbase.dates[start_index:end], copy(firstbase.baseindex))
        newfinal = VarCPIBase(
            lastbase.v[begin:final_index, :],
            copy(lastbase.w), lastbase.dates[begin:final_index], copy(lastbase.baseindex))

        if final_base - start_base > 1
            # more than one base
            @debug "Más de dos bases"
            newbases = (newstart, bases[start_base+1:final_base-1], newfinal)
        else
            # only two bases
            @debug "Dos bases"
            newbases = (newstart, newfinal)
            return getunionalltype(cst)(newbases)
        end
    end

end

"""
    getindex(cst::CountryStructure, finaldate::Date)

Devuelve una copia del `CountryStructure` hasta la fecha indicada por `finaldate`.
"""
function getindex(cst::CountryStructure, finaldate::Date)
    startdate = cst.base[1].dates[1]
    getindex(cst, startdate, finaldate)
end

# Métodos getunionalltype: estos sirven para obtener el tipo concreto de un CountryStructure y poder así construir nuevos objetos de forma genérica.

"""
    getunionalltype(::UniformCountryStructure)

Devuelve el tipo `UniformCountryStructure`. Utilizado al llamar
`getunionalltype` sobre un `CountryStructure` para obtener el tipo concreto
`UnionAll`.
"""
getunionalltype(::UniformCountryStructure) = UniformCountryStructure

"""
    getunionalltype(::MixedCountryStructure)

Devuelve el tipo `MixedCountryStructure`. Utilizado al llamar `getunionalltype`
sobre un `CountryStructure` para obtener el tipo concreto `UnionAll`.
"""
getunionalltype(::MixedCountryStructure) = MixedCountryStructure

## Utilidades

"""
    eltype(::CountryStructure{N, T})

Tipo de dato de punto flotante del contenedor de la estructura de país
`CountryStructure`.
"""
eltype(::CountryStructure{N, T}) where {N,T} = T

"""
    periods(cst::CountryStructure)

Computa el número de períodos (meses) en las bases de variaciones intermensuales
de la estructura de país.
"""
periods(cst::CountryStructure) = sum(size(b.v, 1) for b in cst.base)

"""
    infl_periods(cst::CountryStructure)

Computa el número de períodos de inflación de la estructura de país. Corresponde
al número de observaciones intermensuales menos las primeras 11 observaciones de
la primera base del IPC.
"""
infl_periods(cst::CountryStructure) = periods(cst) - 11

"""
    infl_periods(cst::CountryStructure)

Fechas correspondientes a la trayectorias de inflación computadas a partir un
`CountryStructure`.
"""
infl_dates(cst::CountryStructure) =
    first(cst.base).dates[12]:Month(1):last(cst.base).dates[end]

# TO-DO: Resolver problema de infl_dates con bases disconexas en el CountryStructure


text:
module UIVariables

#------------------------------------------------------------------------------
# EXports
#------------------------------------------------------------------------------
export  BasicValidation,
        Variable,
        on

using ..Twinkle
using Parameters
import UUIDs

function typedict(x::Union{AbstractUIVariable, AbstractUIValidation})
    res = Dict()

    for fn in filter(x->!startswith(string(x), "_"), fieldnames(typeof(x)))
        key = fn
        val = getfield(x, fn)
        if (val === nothing)
            continue
        end
        if (isstructtype(typeof(val)) && !(val isa String))
            @show key
            res[key] = typedict(val)
        else
            res[key] = val;
        end
    end
    # if (isstructtype(typeof(x)))
    # Dict(fn=>getfield(x, fn) for fn ∈ filter(x->!startswith(string(x), "_"), fieldnames(typeof(x))))

    return res;
end

# ----------------------------------------------------------------------
# Validation
# ----------------------------------------------------------------------
@with_kw mutable struct BasicValidation <: AbstractUIValidation
    type::String = "any"
    min::Union{Nothing, Float64} = nothing
    max::Union{Nothing, Float64} = nothing
    valid_values::Union{Nothing, Vector{Any}} = nothing
    invalid_values::Union{Nothing, Vector{Any}} = nothing
end

function BasicValidation(min::Any, max::Any)
    # _, min, max = promote(0.0, min, max)
    return BasicValidation(type="number", min=min, max=max)
end

# ----------------------------------------------------------------------
# Variables
# ----------------------------------------------------------------------
@with_kw mutable struct Variable <: AbstractUIVariable
    name::String = "defaultVariableName"
    type::String = "any"
    value::Union{Nothing, Any} = nothing;
    validation::Union{Nothing, BasicValidation} = nothing

    _app::Union{Nothing, Twinkle.AbstractUIApp} = nothing
    _on_change_func::Union{Nothing, Function} = nothing
end

function on(func::Function, var::AbstractUIVariable)
    var._on_change_func = func
end

function Base.getindex(var::Variable)
    if (var.type == "integer" || var.type == "int")
        if (isa(var.value, String))
            return Int(floor(parse(Float64, var.value)))
        elseif (isa(var.value, Int64))
            return var.value
        elseif (isa(var.value, Float64))
            return Int(floor(var.value))
        else
            return var.value
        end
    elseif (var.type == "float64")
        if (isa(var.value, String))
            return parse(Float64, var.value)
        elseif (isa(var.value, Int))
            return convert(Float64, var.value)
        elseif (isa(var.value, Float64))
            return var.value
        else
            return var.value
        end
    else
        return var.value
    end
end

function Base.setindex!(var::Variable, val::Any)
    # @info "setting variable $(var.name) to $(val) - app: $(var._app)"
    var.value = val;
    Twinkle.updateVariable!(var._app, var)
end

function Base.setindex!(var::Variable, val::Any, type::String)
    if (type == "png")
        try
            val = Twinkle.base64png(val)
        catch e
            @error "Can't convert to SVG [$e]"
        end
    elseif (type == "svg")
        try
            val = Twinkle.base64svg(val)
        catch e
            @error "Can't convert to SVG [$e]"
        end
    end
    var.value = val;
    Twinkle.updateVariable!(var._app, var)
end

end # module Variables


text:
simplestep(s) = d -> d>0 ? s : -s
sinstep(s,w) = d -> (abs(d) < w ? sin(π/2*d/w) : sign(d))*s

function faultline!(hf::Heightfield, iters, step::Function)
    nx,nz = size(hf.y)
    x = linspace(-1,1,nx)*hf.w/2
    z = linspace(-1,1,nz)*hf.h/2
    d = sqrt(hf.w^2 + hf.h^2)
    for i = 1:iters
        v = 2π*rand()
        a,b = sin(v),cos(v)
        c = (2rand()-1)*d/2

        for xi = eachindex(x)
            for zi = eachindex(z)
                hf.y[xi,zi] += step(a*x[xi] + b*z[zi] - c)
            end
        end
    end
end

export faultline!, simplestep, sinstep


text:
using Documenter

using Fridge

makedocs(sitename="Fridge.jl",
    format = Documenter.HTML(),
    modules=[Fridge], # add your module
    pages=[
        "index.md",
        "main function" => "mainFunction.md",
        "search algorithms" => "searchAlgorithms.md",
        "recipe webscraper"=> "recipeWebscraper.md",
    ])

deploydocs(
            repo = "github.com/wardvanbelle/Fridge.jl",
        )


text:
using MCMCChain
using Test

# Tests for missing values.

@testset "utils" begin
    x = [1, missing, 3, 2]
    @test MCMCChain.cummean(x) == [1., 1., 2., 2.]
end

@testset "stats" begin
    chn = Chains(randn(1000, 2, 2))

    # Call describe without missing values.
    describe(devnull, chn)
    s1, s2 = MCMCChain.summarystats(chn), MCMCChain.quantile(chn)

    # Add missing values.
    chn = Chains(cat(chn.value, ones(1, 2, 2) .* missing, dims = 1))

    # Call describe with missing values.
    describe(devnull, chn)
    m1, m2 = MCMCChain.summarystats(chn), MCMCChain.quantile(chn)

    @test all(s1.value[:,1:4,:] .== m1.value[:,1:4,:])
    @test all(s1.value[:,5,:] .+ 1 .== m1.value[:,5,:])
end

@testset "diagnostic functions" begin
    chn = Chains(randn(5000, 2, 2))

    # Add missing values.
    chn_m = Chains(cat(chn.value, ones(1, 2, 2) .* missing, dims = 1))

    # Currently we only check if diag. functions throw an assertion if a value is missing.
    @test_throws AssertionError gelmandiag(chn_m)

    @test all(gewekediag(chn).value .== gewekediag(chn_m).value)
    @test all(heideldiag(chn).value .== heideldiag(chn_m).value)
    @test all(rafterydiag(chn).value .== rafterydiag(chn_m).value)

    @test_throws AssertionError discretediag(chn_m)
end


text:
function gnorm!(c::AbstractVector, γ)
    T = eltype(c)
    if γ == zero(T)
        c[1] = exp(c[1])
        return c
    end

    gain = one(T) + γ*c[1]
    for i=2:length(c)
        @inbounds c[i] /= gain
    end
    c[1] = gain^(one(T)/γ)

    c
end

gnorm(c::AbstractVector, γ=0.0) = gnorm!(copy(c), γ)

function gnorm!{T<:GeneralizedLogCepstrum}(state::SpectralParamState{T})
    assert_not_ready_to_filt(state)
    assert_gain_unnormalized(state)

    γ = glog_gamma(paramdef(state))
    state.data = gnorm!(rawdata(state), γ)
    state.gain_normalized = true
    state
end

function gnorm!{T<:Union{AllPoleLogCepstrum,StandardLogCepstrum}}(state::SpectralParamState{T})
    assert_not_ready_to_filt(state)
    assert_gain_unnormalized(state)
    state
end

function gnorm{T<:MelGeneralizedCepstrum}(state::SpectralParamState{T})
    gnorm!(copy(state))
end


text:
export WeightedSelector

using StatsBase: sample, Weights

"""
    WeightedSelector(is_normalized::Bool)

`is_normalized` is used to indicating if the feeded action values
are alrady normalized to have a sum of `1.0`.
"""
struct WeightedSelector <: AbstractDiscreteActionSelector
    is_normalized::Bool
end

"""
    (p::WeightedSelector)(values::AbstractVector)

!!! note
    Action `values` are normalized to have a sum of 1.0
    and then used as the probability to sample a random action.
"""
(s::WeightedSelector)(values::AbstractVector{T}; kw...) where {T} =
    sample(Weights(values, s.is_normalized ? one(T) : sum(wsum)))


text:
module QuantumFluidSpectra

using Tullio
using FFTW
using SpecialFunctions
using PaddedViews
using UnPack
using TensorCast

abstract type Field end
struct Psi{D} <: Field
    ψ::Array{Complex{Float64},D}
    X::NTuple{D}
    K::NTuple{D}
end

include("arrays.jl")
include("analysis.jl")

export Psi, xvecs, kvecs
export auto_correlate, cross_correlate
export bessel_reduce, sinc_reduce, gv
export log10range, convolve

export xk_arrays, fft_differentials
export gradient, velocity, current
export energydecomp, helmholtz, kinetic_density
export incompressible_spectrum, compressible_spectrum, qpressure_spectrum
export incompressible_density, compressible_density, qpressure_density
export ic_density, iq_density, cq_density
export density_spectrum, trap_spectrum

end


text:
################################################################################
#
#       GrpAb/Elem.jl : Elements in finitely generated abelian groups
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016, 2017: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016, 2017 Tommy Hofmann, Claus Fieker
#
################################################################################

export GrpAbFinGen, GrpAbFinGenElem, parent, isfinite, is_infinite, rank,
       getindex, show, +, *, ngens, snf_with_transform, nrels,
       -, ==, order, exponent,
       quo, sub, rels, has_image, haspreimage, is_snf, is_cyclic, hom, kernel,
       psylow_subgroup

import Base.+, Nemo.snf, Nemo.parent, Base.rand, Nemo.is_snf

function Base.deepcopy_internal(x::GrpAbFinGenElem, dict::IdDict)
  return GrpAbFinGenElem(parent(x), Base.deepcopy_internal(x.coeff, dict))
end

################################################################################
#
#  Constructors
#
################################################################################

# This destroy's the input. If you don't want this, use A(::fmpz_mat)

function GrpAbFinGenElem(A::GrpAbFinGen, a::fmpz_mat)
  if is_snf(A)
    return elem_snf(A, a)
  else
    return elem_gen(A, a)
  end
end

function elem_gen(A::GrpAbFinGen, a::fmpz_mat)
  assure_has_hnf(A)
  reduce_mod_hnf_ur!(a, A.hnf)
  z = GrpAbFinGenElem()
  z.parent = A
  z.coeff = a
  return z
end

function reduce_mod_snf!(a::fmpz_mat, v::Vector{fmpz})
  GC.@preserve a begin
    for i = 1:length(v)
      d = v[i]
      if !iszero(d)
        for j = 1:nrows(a)
          t = ccall((:fmpz_mat_entry, libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), a, j - 1, i - 1)
          ccall((:fmpz_mod, libflint), Nothing, (Ptr{fmpz}, Ptr{fmpz}, Ref{fmpz}), t, t, d)
        end
        #a[1, i] = mod(a[1, i], A.snf[i])
      end
    end
  end
end

function elem_snf(A::GrpAbFinGen, a::fmpz_mat)
  reduce_mod_snf!(a, A.snf)
  z = GrpAbFinGenElem()
  z.parent = A
  z.coeff = a
  return z
end

################################################################################
#
#  Generators
#
################################################################################

@doc Markdown.doc"""
    gens(G::GrpAbFinGen) -> Vector{GrpAbFinGenElem}

The sequence of generators of $G$.
"""
gens(G::GrpAbFinGen) = GrpAbFinGenElem[G[i] for i = 1:ngens(G)]

@doc Markdown.doc"""
    gen(G::GrpAbFinGen, i::Int) -> Vector{GrpAbFinGenElem}

The $i$-th generator of $G$.
"""
gen(G::GrpAbFinGen, i::Int) = G[i]

################################################################################
#
#  Parent
#
################################################################################

@doc Markdown.doc"""
    parent(x::GrpAbFinGenElem) -> GrpAbFinGen

Returns the parent of $x$.
"""
function parent(x::GrpAbFinGenElem)
  return x.parent
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, a::GrpAbFinGenElem)
  @show_special_elem(io, a)

  if get(io, :compact, false)
    print(io, a.coeff)
  else
    s = get_attribute(parent(a), :name)
    s === nothing
    if s === nothing
      print(io, "Element of\n")
      print(io, parent(a))
      print(io, "\nwith components ", a.coeff)
    else
      print(io, "Element of ", s, " with components ", a.coeff)
    end
  end
end

################################################################################
#
#  Hashing
#
################################################################################

function Base.hash(a::GrpAbFinGenElem, s::UInt)
  return hash(a.coeff, s)
end

################################################################################
#
#  Indexing
#
################################################################################

@doc Markdown.doc"""
    getindex(x::GrpAbFinGenElem, i::Int) -> fmpz

Returns the $i$-th component of the element $x$.
"""
function getindex(x::GrpAbFinGenElem, i::Int)
  return x.coeff[1, i]
end

################################################################################
#
#  Comparison
#
################################################################################

function ==(a::GrpAbFinGenElem, b::GrpAbFinGenElem)
  a.parent == b.parent || error("Elements must belong to the same group")
  return a.coeff == b.coeff
end

################################################################################
#
#  Arithmetic
#
################################################################################

function +(x::GrpAbFinGenElem, y::GrpAbFinGenElem, L::GrpAbLattice = GroupLattice)
  if x.parent === y.parent
    n = GrpAbFinGenElem(x.parent, x.coeff + y.coeff)
    return n
  end

  b, m = can_map_into(L, x.parent, y.parent)
  if b
    return GrpAbFinGenElem(y.parent, x.coeff*m) + y
  end

  b, m = can_map_into(L, y.parent, x.parent)
  if b
    return x + GrpAbFinGenElem(x.parent, y.coeff*m)
  end

  b, G, m1, m2 = can_map_into_overstructure(L, x.parent, y.parent)
  if b
    return GrpAbFinGenElem(G, x.coeff * m1 + y.coeff * m2)
  end

  error("Cannot coerce elements into common structure")
end

op(x::GrpAbFinGenElem, y::GrpAbFinGenElem, L::GrpAbLattice = GroupLattice) = +(x, y, L)

function -(x::GrpAbFinGenElem, y::GrpAbFinGenElem)
  x.parent == y.parent || error("Elements must belong to the same group")
  n = GrpAbFinGenElem(x.parent, x.coeff - y.coeff)
  return n
end

function -(x::GrpAbFinGenElem)
  n = GrpAbFinGenElem(x.parent, -x.coeff)
  return n
end

function *(x::fmpz, y::GrpAbFinGenElem)
  n = x*y.coeff
  return GrpAbFinGenElem(y.parent, n)
end

function *(x::Integer, y::GrpAbFinGenElem)
  n = x*y.coeff
  return GrpAbFinGenElem(y.parent, n)
end

*(x::GrpAbFinGenElem, y::fmpz) = y*x

*(x::GrpAbFinGenElem, y::Integer) = y*x

################################################################################
#
#  Neutral element
#
################################################################################

iszero(a::GrpAbFinGenElem) = iszero(a.coeff)

isone(a::GrpAbFinGenElem) = iszero(a.coeff)

is_identity(a::GrpAbFinGenElem) = iszero(a.coeff)

################################################################################
#
#  Parent object overloading
#
################################################################################

@doc Markdown.doc"""
    (A::GrpAbFinGen)(x::Vector{fmpz}) -> GrpAbFinGenElem

Given an array `x` of elements of type `fmpz` of the same length as ngens($A$),
this function returns the element of $A$ with components `x`.
"""
function (A::GrpAbFinGen)(x::Vector{fmpz})
  ngens(A) != length(x) && error("Lengths do not coincide")
  y = matrix(FlintZZ, 1, ngens(A), x)
  z = GrpAbFinGenElem(A, y)
  return z
end

@doc Markdown.doc"""
    (A::GrpAbFinGen)(x::Vector{Integer}) -> GrpAbFinGenElem

Given an array `x` of elements of type `Integer` of the same length as
ngens($A$), this function returns the element of $A$ with components `x`.
"""
function (A::GrpAbFinGen)(x::AbstractVector{T}) where T <: Integer
  ngens(A) != length(x) && error("Lengths do not coincide")
  z = A(map(fmpz, x))
  return z
end

@doc Markdown.doc"""
    (A::GrpAbFinGen)(x::fmpz_mat) -> GrpAbFinGenElem

Given a matrix over the integers with $1$ row and `ngens(A)` columns,
this function returns the element of $A$ with components `x`.
"""
function (A::GrpAbFinGen)(x::fmpz_mat)
  ngens(A) != ncols(x) && error("Lengths do not coincide")
  nrows(x) != 1 && error("Matrix should have only one row")
  z = GrpAbFinGenElem(A, Base.deepcopy(x))
  return z
end

function (A::GrpAbFinGen)()
  y = zero_matrix(FlintZZ, 1, ngens(A))
  z = GrpAbFinGenElem(A, y)
  return z
end

@doc Markdown.doc"""
    getindex(A::GrpAbFinGen, i::Int) -> GrpAbFinGenElem

Returns the element of $A$ with components $(0,\dotsc,0,1,0,\dotsc,0)$,
where the $1$ is at the $i$-th position.
"""
function getindex(A::GrpAbFinGen, i::Int)
  (i < 0 || i > ngens(A)) && error("Index ($i) out of range (1:$(ngens(A)))")
  if i==0
    return GrpAbFinGenElem(A, zero_matrix(FlintZZ, 1, ngens(A)))
  end
  z = zero_matrix(FlintZZ, 1, ngens(A))
  for j in 1:ngens(A)
    z[1, j] = fmpz()
  end
  z[1, i] = fmpz(1)
  return GrpAbFinGenElem(A, z)
end

################################################################################
#
#  Order
#
################################################################################

@doc Markdown.doc"""
    order(A::GrpAbFinGenElem) -> fmpz

Returns the order of $A$. It is assumed that the order is finite.
"""
function order(a::GrpAbFinGenElem)
  G, m = snf(a.parent)
  b = m\a
  o = fmpz(1)
  for i=1:ngens(G)
    if iszero(G.snf[i])
      if !iszero(b[i])
        error("Element has infinite order")
      else
        continue
      end
    end
    o = lcm(o, divexact(G.snf[i], gcd(G.snf[i], b[i])))
  end
  return o
end

##############################################################################
#
#  Random elements
#
##############################################################################

#this allows some more complicated rand(G, (2,2)) and similar.
#TODO: figure out how this SHOULD be done

rand(rng::AbstractRNG, a::Random.SamplerTrivial{GrpAbFinGen, GrpAbFinGenElem}) = rand(a.self)

@doc Markdown.doc"""
    rand(G::GrpAbFinGen) -> GrpAbFinGenElem

Returns an element of $G$ chosen uniformly at random.
"""
rand(A::GrpAbFinGen) = is_snf(A) ? rand_snf(A) : rand_gen(A)

function rand_snf(G::GrpAbFinGen)
  if !isfinite(G)
    error("Group is not finite")
  end
  return G([rand(1:G.snf[i]) for i in 1:ngens(G)])
end

function rand_gen(G::GrpAbFinGen)
  S, mS = snf(G)
  return image(mS, rand(S))
end

@doc Markdown.doc"""
    rand(G::GrpAbFinGen, B::fmpz) -> GrpAbFinGenElem

For a (potentially infinite) abelian group $G$, return an element
chosen uniformly at random with coefficients bounded by $B$.
"""
rand(G::GrpAbFinGen, B::fmpz) = is_snf(G) ? rand_snf(G, B) : rand_gen(G, B)

@doc Markdown.doc"""
    rand(G::GrpAbFinGen, B::Integer) -> GrpAbFinGenElem

For a (potentially infinite) abelian group $G$, return an element
chosen uniformly at random with coefficients bounded by $B$.
"""
rand(G::GrpAbFinGen, B::Integer) = is_snf(G) ? rand_snf(G, B) : rand_gen(G, B)

function rand_snf(G::GrpAbFinGen, B::fmpz)
  z = G([rand(1:(iszero(G.snf[i]) ? B : min(B, G.snf[i]))) for i in 1:ngens(G)])
  return z
end

function rand_snf(G::GrpAbFinGen, B::Integer)
  return rand(G, fmpz(B))
end

function rand_gen(G::GrpAbFinGen, B::fmpz)
  S, mS = snf(G)
  return image(mS, rand(S, fmpz(B)))
end

function rand_gen(G::GrpAbFinGen, B::Integer)
  return rand(G, fmpz(B))
end

################################################################################
#
#  Iterator
#
################################################################################

function Base.iterate(G::GrpAbFinGen)
  if order(G) > typemax(UInt)
    error("Group too large for iterator")
  end

  return _elem_from_enum(G, UInt(0)), UInt(1)
end

function Base.iterate(G::GrpAbFinGen, st::UInt)
  if st >= order(G)
    return nothing
  end

  a = _elem_from_enum(G, st)
  return a, st + 1
end

function _elem_from_enum(G::GrpAbFinGen, st::UInt)
  if G.is_snf
    el = fmpz[]
    s = fmpz(st)
    for i = 1:ngens(G)
      push!(el, s % G.snf[i])
      s = div(s - el[end], G.snf[i])
    end
    return G(el)
  end

  S, mS = snf(G)
  return image(mS, _elem_from_enum(S, st))
end

Base.IteratorSize(::Type{GrpAbFinGen}) = Base.HasLength()

Base.length(G::GrpAbFinGen) = Int(order(G))

Base.eltype(::Type{GrpAbFinGen}) = GrpAbFinGenElem


text:
using OVERT: get_symbols, find_variables, substitute!

""" this function turns :(w == -(x+1)) to :(w == -1(x+1)) """
function fix_negation!(expr::Expr)
    if length(expr.args[3].args) == 2
        if expr.args[3].args[1] == :-
            arg = expr.args[3].args[2]
            expr.args[3].args = [:*, -1, arg]
        end
    end
    return expr
end

""" find if this expr constains max """
is_max_expr(expr::Union{Symbol, Expr}) = :max in get_symbols(expr)
is_min_expr(expr::Union{Symbol, Expr}) = :min in get_symbols(expr)

"given expr = :(2x+3y-3z+1) return [:x, :y, :c], [2, 3, -3], 1"
get_linear_coeffs(expr::Symbol) = [expr], [1.0], 0
function get_linear_coeffs(expr::Expr)
    vars = find_variables(expr)
    coeffs = zeros(length(vars))
    scalar = deepcopy(expr)
    for v in vars
        substitute!(scalar, v, 0)
    end
    scalar = eval(scalar)

    for i = 1:length(vars)
        expr_copy = deepcopy(:($expr - $scalar))
        for v in vars
            if vars[i] == v
                substitute!(expr_copy, v, 1)
            else
                substitute!(expr_copy, v, 0)
            end
        end
        coeffs[i] = eval(:($expr_copy))
    end
    return vars, coeffs, scalar
end


text:
using MyExample2
using Test

@testset "MyExample2.jl" begin
    # Write your tests here.
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

using MathOptInterface, Test

const MOI = MathOptInterface
const MOIU = MOI.Utilities

const NUM_BENCHMARKS = length(MOI.Benchmarks.BENCHMARKS)

@testset "suite" begin
    suite = MOI.Benchmarks.suite() do
        return MOIU.MockOptimizer(MOIU.Model{Float64}())
    end
    @test length(suite.data) == NUM_BENCHMARKS

    suite = MOI.Benchmarks.suite(exclude = [r"delete_"]) do
        return MOIU.MockOptimizer(MOIU.Model{Float64}())
    end
    # Note: update this value whenever more benchmarks are added to
    # `src/Benchmarks/Benchmarks.jl`.
    @test 6 <= length(suite.data) <= NUM_BENCHMARKS - 3
end

@testset "Perform benchmark" begin
    params = joinpath(@__DIR__, "baseline_params.json")
    baseline = joinpath(@__DIR__, "baseline_baseline.json")
    @test !isfile(params)
    @test !isfile(baseline)
    @testset "create_baseline" begin
        suite = MOI.Benchmarks.suite() do
            return MOIU.MockOptimizer(MOIU.Model{Float64}())
        end
        MOI.Benchmarks.create_baseline(
            suite,
            "baseline";
            directory = @__DIR__,
            seconds = 2,
            verbose = true,
        )
    end
    @test isfile(params)
    @test isfile(baseline)
    @testset "compare_against_baseline" begin
        suite = MOI.Benchmarks.suite() do
            return MOIU.MockOptimizer(MOIU.Model{Float64}())
        end
        MOI.Benchmarks.compare_against_baseline(
            suite,
            "baseline";
            directory = @__DIR__,
            seconds = 2,
            verbose = true,
        )
    end
    rm(params)
    rm(baseline)
    @testset "Report" begin
        report = read(joinpath(@__DIR__, "report.txt"), String)
        @test occursin("=> invariant", report)
    end
    rm(joinpath(@__DIR__, "report.txt"))
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
using DelimitedFiles
using StatsBase

"returns the index of the winner, or -1 in case there is no winner."
function winner(marked)
    for k in 1:1:size(marked,3)
        board = view(marked,:,:,k)
        for col in eachcol(board)
            if sum(col) == 5
                return k
            end
        end
        for row in eachrow(board)
            if sum(row) == 5
                return k
            end
        end
    end
    return -1
end
function answer1()
    input = "input"
    cardboards = permutedims((reshape(readdlm(input,Int,skipstart=2)', 5,5,:)),[2,1,3])
    ballots = parse.(Int,split(readline(input),","))
    marked = zeros(size(cardboards))

    for ballot in ballots
        marked .+= (cardboards .== ballot)
        first_winner = winner(marked)
        if first_winner > 0
            # @show first_winner, ballot
            # sum(-1 .* (marked[:,:,first_winner] .- 1) .* cardboards[:,:,first_winner])*ballot
            return Int(sum(-1 .* (marked[:,:,first_winner] .- 1) .* cardboards[:,:,first_winner])*ballot)
        end
    end
end
@show answer1()

"Returns a list of winners, [-1] if no new winners."
function winner(marked, candidates)
    # @show candidates
    winners = []
    for k in candidates
        board = view(marked,:,:,k)
        for col in eachcol(board)
            if sum(col) == 5
                append!(winners,k)
            end
        end
        for row in eachrow(board)
            if sum(row) == 5
                append!(winners,k)
            end
        end
    end
    if size(winners,1)>0
        return unique(winners)
    end
    return [-1]
end
function answer2()
    input = "input"
    cardboards = permutedims((reshape(readdlm(input,Int,skipstart=2)', 5,5,:)),[2,1,3])
    ballots = parse.(Int,split(readline(input),","))
    marked = zeros(size(cardboards))
    winners = ones(size(cardboards,3))
    candidates = collect(1:1:size(winners,1))[winners .== 1]
    last_winner = 0
    points_until_won = 0
    for ballot in ballots
        marked .+= (cardboards .== ballot)
        first_winner = winner(marked,candidates)
        for first_winner in winner(marked,candidates)
            if first_winner != -1
                winners[first_winner] = 0
                candidates = collect(1:1:size(winners,1))[winners .== 1]
            end
            if first_winner > 0 && first_winner != last_winner
                last_winner = first_winner
                points_until_won = Int(sum(-1 .* (marked[:,:,last_winner] .- 1) .* cardboards[:,:,last_winner])*ballot)
            end
        end
    end
    return points_until_won
end
@show answer2()


text:
@reexport using Distributions
import Distributions: mean, var, logpdf, pdf, cdf, maximum, minimum, insupport, quantile, rand, rand!, params
export DiscreteUnivariateCopula, pdf_constants
export pmf_copula!, reorder_pmf

struct DiscreteUnivariateCopula{
    DistT <: DiscreteUnivariateDistribution,
    T     <: Real
    } <: DiscreteUnivariateDistribution
    d  :: DistT
    μ :: T
    σ2 :: T
    c0 :: T
    c1 :: T
    c2 :: T
    c  :: T # normalizing constant
    maxvalue :: T
    y_sample :: Vector{T}
    pmf_vec :: Vector{T}
end

"""
    DiscreteUnivariateCopula(d, c0, c1, c2)
The distribution with density `c * P(x = x) * (c0 + c1 * x + c2 * x^2)`, where `f`
is the density of the base distribution `d` and `c` is the normalizing constant.
"""
function DiscreteUnivariateCopula(
    d  :: DistT,
    c0 :: T,
    c1 :: T,
    c2 :: T) where {DistT <: DiscreteUnivariateDistribution, T <: Real}
    μ = mean(d)
    σ2 = var(d)
    c  = inv(c0 + c1 * μ + c2 * (σ2 + abs2(μ)))
    Tc = typeof(c)
    maxvalue = quantile(Base.typename(typeof(d)).wrapper(Distributions.params(d)...), 0.999999999999)
    y_sample = collect(0:maxvalue)
    pmf_vec = zeros(length(y_sample)) # marginal pmf
    DiscreteUnivariateCopula(d, Tc(μ), Tc(σ2), Tc(c0), Tc(c1), Tc(c2), c, Tc(maxvalue), Tc.(y_sample), Tc.(pmf_vec))
end

"""
    pdf_constants(Γ::Matrix{<:Real}, res::Vector{<:Real}, i::Int64, dist::DiscreteUnivariateDistribution)
This function will fill out the appropriate constants, c0, c1, c2 for each conditional distribution to form the `DiscreteUnivariateCopula` structure.
"""
function pdf_constants(Γ::Matrix{T}, res::Vector{T}, i::Int64, dist::DiscreteUnivariateDistribution) where T <: Real
    μ = Distributions.mean(dist)
    σ2 = Distributions.var(dist)
    c_0 = μ^2 * inv(σ2)
    c__0 = μ * inv(sqrt(σ2)) * sum(crossterm_res(res, i, Γ))
    c_1 = -2μ * inv(σ2)
    c__1 = inv(sqrt(σ2)) * sum(crossterm_res(res, i, Γ))
    c_2 = inv(σ2)
    storage = zeros(i-1)
    # first multiply Γ[1:i-1, 1:i-1] *  res[1:i-1]
    mul!(storage, Γ[1:i-1, 1:i-1], res[1:i-1])
    # then multiply to get transpose(res[1:i-1]) *  Γ[1:i-1, 1:i-1] *  res[1:i-1]
    c0 = 1 + 0.5 * dot(storage, res[1:i-1]) +  0.5 * tr(Γ[i+1:end, i+1:end]) + 0.5 * Γ[i, i] * c_0  - c__0
    c1 = 0.5 * Γ[i, i] * c_1  + c__1
    c2 = 0.5 * Γ[i, i] * c_2
    DiscreteUnivariateCopula(dist, c0, c1, c2)
end

"""
    pdf_constants(γ::T, dist::DiscreteUnivariateDistribution)
This function will fill out the appropriate constants, c0, c1, c2 for the univariate marginal distribution to form the `DiscreteUnivariateCopula` structure.
When the number of observations in the cluster is 1
"""
function pdf_constants(γ::T, dist::DiscreteUnivariateDistribution) where T <: Real
    μ = Distributions.mean(dist)
    σ2 = Distributions.var(dist)
    c_0 = μ^2 * inv(σ2)
    c_1 = -2μ * inv(σ2)
    c_2 = inv(σ2)
    c0 = 1 + 0.5 * γ * c_0
    c1 = 0.5 * γ * c_1
    c2 = 0.5 * γ * c_2
    DiscreteUnivariateCopula(dist, c0, c1, c2)
end

"""
    pmf_copula!(dist::DiscreteUnivariateCopula)
This function will get the appropriate probability mass function, using our copula density.
For discrete distributions with countably infinite values in the range, we want to find some approximate maximum value that is large enough so that the probability mass vector sums to about 1.
"""
function pmf_copula!(dist::DiscreteUnivariateCopula) where T<: Real
    for k in 1:length(dist.y_sample)
        dist.pmf_vec[k] = pdf(dist, dist.y_sample[k]) # get the pmf given maximum and our pdf implementation
    end
    nothing
end

"""
    reorder_pmf(pmf::Vector{<:Real}, μ)
This function will re-order the probability mass function, by sorting the vector of probabilities in decreasing order (starting with mean μ).
"""
function reorder_pmf(pmf::Vector{T}, μ) where T <: Real
    listofj = zeros(Int64, length(pmf))
    k = Integer(floor.(μ))
    reordered_pmf = zeros(length(pmf))
    i = 1
    j = k[1]
    while(i < length(pmf) && j > 0 && j < length(pmf))
        listofj[i] = j
        reordered_pmf[i] = pmf[j + 1]
        if i%2 == 1
            j = j + i
            elseif i%2 == 0
            j = j - i
        end
        i = i + 1
    end
    if j == 0
        listofj[i] = 0
        reordered_pmf[i] = pmf[1]
        for s in i+1:length(pmf)
            listofj[s] = s - 1
            reordered_pmf[s] = pmf[s]
            end
        end
    listofj, reordered_pmf
end

"""
    rand(dist::DiscreteUnivariateCopula)
This function will simulate the discrete random variable under our copula model.
"""
function rand(dist::DiscreteUnivariateCopula) where T <: Real
    pmf_copula!(dist) # get pmf under our copula density
    listofj, reordered_pmf = reorder_pmf(dist.pmf_vec, dist.μ) # re-order the pmf
    # reorder_pmf!(dist)
    sample = rand() # generate x from uniform(0, 1)
    (random_deviate, s) = listofj[1], reordered_pmf[1] # if the cumulative probability mass is less than the P(X = listofj[1]) then leave it as the mean
    for i in 2:length(reordered_pmf)
        if sample < s
            random_deviate = listofj[i - 1]
            break
        else
            s += reordered_pmf[i]
        end
    end
    random_deviate
end

"""
    rand(dist::DiscreteUnivariateCopula, n_reps::Int64)
This function will simulate the discrete random variable under our copula model, n_reps times.
"""
function rand(dist::DiscreteUnivariateCopula, n_reps::Int64) where T <: Real
    random_deviate = zeros(Int64, n_reps)
    for l in 1:n_reps
        random_deviate[l] = rand(dist)
    end
    random_deviate
end

"""
    rand!(dist::DiscreteUnivariateCopula, sample::Vector{T})
This function will write over each entry in the specified sample vector, the simulated univariate discrete values under our copula model.
"""
function rand!(dist::DiscreteUnivariateCopula, sample::Vector{T}) where T <: Real
    for i in 1:length(sample)
        sample[i] = rand(dist)
    end
    sample
end

"""
    rand(dist::DiscreteUnivariateCopula)
This function will simulate the discrete random variable under our copula model.
"""
function rand(dist::DiscreteUnivariateCopula{Bernoulli{T}, T}) where T <: Real
    random_deviate = Float64(rand(Bernoulli(QuasiCopula.mean(dist))))
    random_deviate
end

"""
    rand(dist::DiscreteUnivariateCopula, n_reps::Int64)
This function will simulate the discrete random variable under our copula model, n_reps times.
"""
function rand(dist::DiscreteUnivariateCopula{Bernoulli{T}, T}, n_reps::Int64) where T <: Real
    random_deviate = zeros(Int64, n_reps)
    for l in 1:n_reps
        random_deviate[l] = rand(dist)
    end
    random_deviate
end

"""
    rand!(dist::DiscreteUnivariateCopula, sample::Vector{T})
This function will write over each entry in the specified sample vector, the simulated univariate discrete values under our copula model.
"""
function rand!(dist::DiscreteUnivariateCopula{Bernoulli{T}, T}, sample::Vector{T}) where T <: Real
    for i in 1:length(sample)
        sample[i] = rand(dist)
    end
    sample
end


text:
"""
    $(TYPEDEF)

Make a two-dimensional, one-legged robot walk forward as fast as possible. The robot model is
based on the work T. Erez, Y. Tassa, and E. Todorov: ["Infinite Horizon Model Predictive
Control for Nonlinear Periodic Tasks"](https://homes.cs.washington.edu/~todorov/papers/ErezRSS11.pdf), 2011.

# Spaces

* **State: (20, )**
* **Action: (3, )**
* **Observation: (11, )**
"""
mutable struct HopperV2{Sim, SSpace, OSpace} <: AbstractMuJoCoEnvironment
    sim::Sim
    statespace::SSpace
    observationspace::OSpace
    last_torso_x::Float64
    randreset_distribution::Uniform{Float64}
    function HopperV2(sim::MJSim)
        sspace = MultiShape(
            simstate=statespace(sim),
            last_torso_x=ScalarShape(Float64)
        )
        ospace = MultiShape(
            cropped_qpos = VectorShape(Float64, sim.m.nq - 1),
            qvel = VectorShape(Float64, sim.m.nv)
        )
        env = new{typeof(sim), typeof(sspace), typeof(ospace)}(
            sim, sspace, ospace, 0, Uniform(-0.005, 0.005)
        )
        reset!(env)
    end
end

HopperV2() = first(tconstruct(HopperV2, 1))

function LyceumBase.tconstruct(::Type{HopperV2}, n::Integer)
    modelpath = joinpath(@__DIR__, "hopper.xml")
    [HopperV2(s) for s in tconstruct(MJSim, n, modelpath, skip=4)]
end

@inline getsim(env::HopperV2) = env.sim

@inline LyceumBase.statespace(env::HopperV2) = env.statespace

function LyceumBase.getstate!(state, env::HopperV2)
    checkaxes(statespace(env), state)
    @uviews state begin
        shaped = statespace(env)(state)
        getstate!(shaped.simstate, env.sim)
        shaped.last_torso_x = env.last_torso_x
    end
    state
end

function LyceumBase.setstate!(env::HopperV2, state)
    checkaxes(statespace(env), state)
    @uviews state begin
        shaped = statespace(env)(state)
        setstate!(env.sim, shaped.simstate)
        env.last_torso_x = shaped.last_torso_x
    end
    env
end

@inline LyceumBase.observationspace(env::HopperV2) = env.observationspace

function LyceumBase.getobservation!(obs, env::HopperV2)
    checkaxes(observationspace(env), obs)
    qpos = env.sim.d.qpos
    @views @uviews qpos obs begin
        shaped = observationspace(env)(obs)
        copyto!(shaped.cropped_qpos, qpos[2:end])
        copyto!(shaped.qvel, env.sim.d.qvel)
        clamp!(shaped.qvel, -10, 10)
    end
    obs
end

function LyceumBase.getreward(state, action, ::Any, env::HopperV2)
    checkaxes(statespace(env), state)
    checkaxes(actionspace(env), action)
    @uviews state begin
        shapedstate = statespace(env)(state)
        alive_bonus = 1.0
        reward = (_torso_x(shapedstate, env) - shapedstate.last_torso_x) / timestep(env)
        reward += alive_bonus
        reward -= 1e-3 * sum(x->x^2, action)
        reward
    end
end

function LyceumBase.reset!(env::HopperV2)
    reset!(env.sim)
    env.last_torso_x = _torso_x(env)
    env
end

function LyceumBase.randreset!(rng::AbstractRNG, env::HopperV2)
    reset_nofwd!(env.sim)
    perturb!(rng, env.randreset_distribution, env.sim.d.qpos)
    perturb!(rng, env.randreset_distribution, env.sim.d.qvel)
    forward!(env.sim)
    env.last_torso_x = _torso_x(env)
    env
end

function LyceumBase.step!(env::HopperV2)
    env.last_torso_x = _torso_x(env)
    step!(env.sim)
    env
end

function LyceumBase.isdone(state, ::Any, env::HopperV2)
    checkaxes(statespace(env), state)
    @uviews state begin
        shapedstate = statespace(env)(state)
        torso_x = _torso_x(shapedstate, env)
        height = _torso_height(shapedstate, env)
        torso_ang = _torso_ang(shapedstate, env)
        qpos = shapedstate.simstate.qpos
        qvel = shapedstate.simstate.qvel

        done = !(
            all(isfinite, state)
            && all(x->abs(x) < 100, uview(qpos, 3:length(qpos)))
            && all(x->abs(x) < 100, uview(qvel))
            && height > 0.7
            && abs(torso_ang) < 0.2
        )
        done
    end
end

@inline _torso_x(shapedstate::ShapedView, ::HopperV2) = shapedstate.simstate.qpos[1]
@inline _torso_x(env::HopperV2) = env.sim.d.qpos[1]
@inline _torso_height(shapedstate::ShapedView, ::HopperV2) = shapedstate.simstate.qpos[2]
@inline _torso_ang(shapedstate::ShapedView, ::HopperV2) = shapedstate.simstate.qpos[3]


text:
import Base: reducedim, sum, maximum
export argmax, argmin

function reducedim(f::CuFunction, x::CuArray, dim::Int)
    x3d = reshape3d(x, dim)
    by = 1024
    while true
        gy = Int(ceil(size(x3d,2)/by))
        y = similar(x, ntuple(i -> i==dim ? gy : size(x,i), ndims(x)))
        y3d = reshape3d(y, dim)
        f(x3d, y3d, dx=size(x3d,1), dy=size(x3d,2), dz=size(x3d,3), bx=1, by=by, bz=1)
        gy == 1 && return y
        x3d = y3d
    end
end

@generated function maximum{T}(x::CuArray{T}, dim::Int)
    f = reducedim_f(T, "-1e10", "a > b ? a : b")
    quote
        reducedim($f, x, dim)
    end
end

@generated function sum{T}(x::CuArray{T}, dim::Int)
    f = reducedim_f(T, "0", "a + b")
    quote
        reducedim($f, x, dim)
    end
end

function reducedim_find(f::CuFunction, x::CuArray, dim::Int)
    x3d = reshape3d(x, dim)
    by = 1024
    while true
        gy = Int(ceil(size(x3d,2)/by))
        y = CuArray{Cint}(ntuple(i -> i==dim ? gy : size(x,i), ndims(x)))
        y3d = reshape3d(y, dim)
        f(x3d, y3d, dx=size(x3d,1), dy=size(x3d,2), dz=size(x3d,3), bx=1, by=by, bz=1)
        gy == 1 && return y
        x3d = y3d
    end
end

@generated function argmax{T}(x::CuArray{T}, dim::Int)
    f = find_f(T, "-1e10", "a < b")
    quote
        reducedim_find($f, x, dim)
    end
end

function reducedim_f{T}(::Type{T}, v0::String, op::String)
    CuFunction("""
    template<typename T>
    __inline__ __device__ void warpReduce(T &a) {
        for (int delta = warpSize/2; delta > 0; delta /= 2) {
            T b = __shfl_down(a, delta);
            a = $op;
        }
    }
    template<typename T>
    __inline__ __device__ void blockReduce(T &value) {
        static __shared__ T values[32]; // 32*32=1024
        int laneId = threadIdx.y % warpSize;
        int warpId = threadIdx.y / warpSize;

        warpReduce<T>(value);
        if (laneId == 0) values[warpId] = value;
        __syncthreads();

        value = (laneId < 32) ? values[laneId] : $v0;
        warpReduce<T>(value);
    }
    __global__ void reduce(Array<$T,3> x, Array<$T,3> y) {
        int idx_x = blockIdx.x * blockDim.x + threadIdx.x;
        int idx_y = blockIdx.y * blockDim.y + threadIdx.y;
        int idx_z = blockIdx.z * blockDim.z + threadIdx.z;

        $T a = (idx_y < x.dims[1]) ? x(idx_x,idx_y,idx_z) : $v0;
        for (int i = idx_y+blockDim.y*gridDim.y; i < x.dims[1]; i += blockDim.y*gridDim.y) {
            $T b = x(idx_x, i, idx_z);
            a = $op;
        }
        blockReduce<$T>(a);
        if (threadIdx.y == 0) y(blockIdx.x, blockIdx.y, blockIdx.z) = a;
    }
    """)
end

function find_f{T}(::Type{T}, v0::String, op::String)
    CuFunction("""
    template<typename T>
    __inline__ __device__ void warpReduce(T &a, int &i) {
        for (int delta = warpSize/2; delta > 0; delta /= 2) {
            T b = __shfl_down(a, delta);
            int j = __shfl_down(i, delta);
            if ($op) {
                a = b;
                i = j;
            }
        }
    }
    template<typename T>
    __inline__ __device__ void blockReduce(T &value, int &index) {
        static __shared__ T values[32]; // 32*32=1024
        static __shared__ int indexs[32];
        int laneId = threadIdx.y % warpSize;
        int warpId = threadIdx.y / warpSize;

        warpReduce<T>(value, index);
        if (laneId == 0) {
            values[warpId] = value;
            indexs[warpId] = index;
        }
        __syncthreads();

        value = (laneId < 32) ? values[laneId] : $v0;
        index = (laneId < 32) ? indexs[laneId] : -1;
        warpReduce<T>(value, index);
    }
    __global__ void f(Array<$T,3> x, Array<int,3> y) {
        int idx_x = blockIdx.x * blockDim.x + threadIdx.x;
        int idx_y = blockIdx.y * blockDim.y + threadIdx.y;
        int idx_z = blockIdx.z * blockDim.z + threadIdx.z;

        $T value = (idx_y < x.dims[1]) ? x(idx_x,idx_y,idx_z) : $v0;
        int index = (idx_y < x.dims[1]) ? blockIdx.y * blockDim.y + threadIdx.y : -1;
        //for (int i = idx_y+blockDim.y*gridDim.y; i < x.dims[1]; i += blockDim.y*gridDim.y) {
        //    $T b = x(idx_x, i, idx_z);
        //    a = $op;
        //}
        blockReduce<$T>(value, index);
        if (threadIdx.y == 0) y(blockIdx.x, blockIdx.y, blockIdx.z) = index + 1;
    }
    """)
end

#=
function reducedim_f{T}(::Type{T}, op::String)
    blocksize = 1024
    CuFunction("""
    __global__ void reduce(Array<$T,3> x, Array<$T,3> y) {
        static __shared__ $T temp[1024];

        int idx_x = blockIdx.x * blockDim.x + threadIdx.x;
        //int idx_y = blockIdx.y * blockDim.y + threadIdx.y;
        int idx_z = blockIdx.z * blockDim.z + threadIdx.z;
        //if (idx_y >= x.dims[1]) return;

        int tid = threadIdx.y;
        int idx_y = blockIdx.y * $blocksize * 2 + threadIdx.y;
        $T a = (idx_y < x.dims[1]) ? x(idx_x, idx_y, idx_z) : 0;
        if (idx_y+$blocksize < x.dims[1]) {
            a += x[idx_y+$blocksize];
        }

        temp[tid] = a;
        $T b;
        __syncthreads();

        if (($blocksize >= 512) && (tid < 256)) {
            b = temp[tid+256];
            a = $op;
            temp[tid] = a;
            __syncthreads();
        }
        if (($blocksize >= 256) &&(tid < 128)) {
            b = temp[tid+128];
            a = $op;
            temp[tid] = a;
            __syncthreads();
        }
        if (($blocksize >= 128) && (tid < 64)) {
            b = temp[tid+64];
            a = $op;
            temp[tid] = a;
            __syncthreads();
        }

        if (tid < 32) {
            if ($blocksize >= 64) a += temp[tid+32];
            for (int offset = warpSize/2; offset > 0; offset /= 2) {
                a += __shfl_down(a, offset);
            }
        }
        if (threadIdx.y == 0) y(blockIdx.x, blockIdx.y, blockIdx.z) = a;
    }""")
end
=#


text:
using Distributions, SpecialFunctions
import Distributions: mean
"""
GMPN: Our variation of Fearnhead & Liu 2007 (MPN algo)
            for the Gaussian task (section 2.3.1)
Update functions can be found in gbayesfilter.jl
The type of approximation can be found in approximations.jl
"""
struct GMPN{D, A} <: BayesFilter
    N::Int
    posteriors::Array{D, 1}
    prior::D
    probabilities::Array{Float64, 1}
    sigma::Float64
    changeprobability::Float64
    approximation::A
    mu_n::Array{Float64, 1}
    seed::Any
    rng::MersenneTwister
    lastprob::Array{Float64, 1}
end
function GMPN(;sigma = .1, mu_0 = 0., sigma_0 = 1.,
                    changeprobability = 0.01, N = 20, seed = 3)

    prior = GaussKnownVariance(mu_0, sigma_0, sigma)
    approximation = DropArgMinApprox(N)
    rng = MersenneTwister(seed)
    GMPN(N, [deepcopy(prior)], prior, [1.], sigma, changeprobability, approximation,
        [0.], seed, rng, [1.])
end