text:
function try_move!(r, side)
    n = 0;
    ort_side = next_side(side)
    while isborder(r, side)
        if !isborder(r, ort_side)
            move!(r, ort_side)
            n+=1
        else
            return false
        endt
    end

    move!(r, side)
    while isborder(r, opposite_side(side))
        move!(r, side)
    end
    move_for!(r, opposite_side(side), n)
    return true
end

abstract type AbstractRobot
end

import HorizonSideRobots

move!(robot::AbstractRobot, side::HorizonSide) = HorizonSideRobots.move!(get(robot), side)
isborder(robot::AbstractRobot,  side::HorizonSide) = HorizonSideRobots.isborder(get(robot), side)
putmarker!(robot::AbstractRobot) = HorizonSideRobots.putmarker!(get(robot))
ismarker(robot::AbstractRobot) = HorizonSideRobots.ismarker(get(robot))
temperature(robot::AbstractRobot) = HorizonSideRobots.temperature(get(robot))

get(robot::AbstractRobot) = robot

abstract type AbstractBorderRobot<:AbstractRobot
end

#movements! идет до упора и возвращает число шагов или делает заданное число шагов в заданном направлении
movements!(robot::AbstractBorderRobot, side::HorizonSideRobots)
movements!(robot::AbstractBorderRobot, side::HorizonSideRobots, steps::Integer)

function movements!(robot::AbstractBorderRobot,side::HorizonSide)::Integer
    n = 0
    while try_move!(get(robot),side)
        n += 1
    end
    return n
end

function snake(robot::AbstractRobot)::Nothing    # function snake(robot::Union{AbstractRobot,Robot}::Nothing
    side = Ost
    movements!(robot,side)
    while (!isborder(robot,Nord))
        move!(robot, Nord)
        side = anti_side(side)
        movements!(robot,side)
    end
end

struct PutMerkersBorder<:AbstractBorderRobot
    robot::Robot
end

function try_move!(robot::AbstractBorderRobot,side::HorizonSide)
    try_move!(get(robot),side)
    putmarker!(r)
end
function PutMarkersBorder()
    putmarker!(robot)
    new(robot)
end

function mark_field!(robot::Robot)
    robot = PutMarkersBorder
    snake(robot)
end


text:
using LyapunovExponents
using Plots
demo = solve!(LyapunovExponents.lozi_map(); progress=1)
report(demo)
plt = plot(demo)

ext: jl


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# Initialize data structures in element container
function init_elements!(elements, mesh::P4estMesh{2}, basis::LobattoLegendreBasis)
  @unpack node_coordinates, jacobian_matrix,
          contravariant_vectors, inverse_jacobian = elements

  calc_node_coordinates!(node_coordinates, mesh, basis)

  for element in 1:ncells(mesh)
    calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis)

    calc_contravariant_vectors!(contravariant_vectors, element, jacobian_matrix)

    calc_inverse_jacobian!(inverse_jacobian, element, jacobian_matrix)
  end

  return nothing
end

# Interpolate tree_node_coordinates to each quadrant at the nodes of the specified basis
function calc_node_coordinates!(node_coordinates,
                                mesh::P4estMesh{2},
                                basis::LobattoLegendreBasis)
  # Hanging nodes will cause holes in the mesh if its polydeg is higher
  # than the polydeg of the solver.
  @assert length(basis.nodes) >= length(mesh.nodes) "The solver can't have a lower polydeg than the mesh"

  calc_node_coordinates!(node_coordinates, mesh, basis.nodes)
end

# Interpolate tree_node_coordinates to each quadrant at the specified nodes
function calc_node_coordinates!(node_coordinates,
                                mesh::P4estMesh{2},
                                nodes::AbstractVector)
  # We use `StrideArray`s here since these buffers are used in performance-critical
  # places and the additional information passed to the compiler makes them faster
  # than native `Array`s.
  tmp1    = StrideArray(undef, real(mesh),
                        StaticInt(2), static_length(nodes), static_length(mesh.nodes))
  matrix1 = StrideArray(undef, real(mesh),
                        static_length(nodes), static_length(mesh.nodes))
  matrix2 = similar(matrix1)
  baryweights_in = barycentric_weights(mesh.nodes)

  # Macros from `p4est`
  p4est_root_len = 1 << P4EST_MAXLEVEL
  p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l)

  trees = unsafe_wrap_sc(p4est_tree_t, mesh.p4est.trees)

  for tree in eachindex(trees)
    offset = trees[tree].quadrants_offset
    quadrants = unsafe_wrap_sc(p4est_quadrant_t, trees[tree].quadrants)

    for i in eachindex(quadrants)
      element = offset + i
      quad = quadrants[i]

      quad_length = p4est_quadrant_len(quad.level) / p4est_root_len

      nodes_out_x = 2 * (quad_length * 1/2 * (nodes .+ 1) .+ quad.x / p4est_root_len) .- 1
      nodes_out_y = 2 * (quad_length * 1/2 * (nodes .+ 1) .+ quad.y / p4est_root_len) .- 1
      polynomial_interpolation_matrix!(matrix1, mesh.nodes, nodes_out_x, baryweights_in)
      polynomial_interpolation_matrix!(matrix2, mesh.nodes, nodes_out_y, baryweights_in)

      multiply_dimensionwise!(
        view(node_coordinates, :, :, :, element),
        matrix1, matrix2,
        view(mesh.tree_node_coordinates, :, :, :, tree),
        tmp1
      )
    end
  end

  return node_coordinates
end

# Initialize node_indices of interface container
@inline function init_interface_node_indices!(interfaces::P4estInterfaceContainer{2},
                                              faces, orientation, interface_id)
  # Iterate over primary and secondary element
  for side in 1:2
    # Align interface in positive coordinate direction of primary element.
    # For orientation == 1, the secondary element needs to be indexed backwards
    # relative to the interface.
    if side == 1 || orientation == 0
      # Forward indexing
      i = :i_forward
    else
      # Backward indexing
      i = :i_backward
    end

    if faces[side] == 0
      # Index face in negative x-direction
      interfaces.node_indices[side, interface_id] = (:begin, i)
    elseif faces[side] == 1
      # Index face in positive x-direction
      interfaces.node_indices[side, interface_id] = (:end, i)
    elseif faces[side] == 2
      # Index face in negative y-direction
      interfaces.node_indices[side, interface_id] = (i, :begin)
    else # faces[side] == 3
      # Index face in positive y-direction
      interfaces.node_indices[side, interface_id] = (i, :end)
    end
  end

  return interfaces
end

# Initialize node_indices of boundary container
@inline function init_boundary_node_indices!(boundaries::P4estBoundaryContainer{2},
                                             face, boundary_id)
  if face == 0
    # Index face in negative x-direction
    boundaries.node_indices[boundary_id] = (:begin, :i_forward)
  elseif face == 1
    # Index face in positive x-direction
    boundaries.node_indices[boundary_id] = (:end, :i_forward)
  elseif face == 2
    # Index face in negative y-direction
    boundaries.node_indices[boundary_id] = (:i_forward, :begin)
  else # face == 3
    # Index face in positive y-direction
    boundaries.node_indices[boundary_id] = (:i_forward, :end)
  end

  return boundaries
end

# Initialize node_indices of mortar container
# faces[1] is expected to be the face of the small side.
@inline function init_mortar_node_indices!(mortars::P4estMortarContainer{2},
                                           faces, orientation, mortar_id)
  for side in 1:2
    # Align mortar in positive coordinate direction of small side.
    # For orientation == 1, the large side needs to be indexed backwards
    # relative to the mortar.
    if side == 1 || orientation == 0
      # Forward indexing for small side or orientation == 0
      i = :i_forward
    else
      # Backward indexing for large side with reversed orientation
      i = :i_backward
    end

    if faces[side] == 0
      # Index face in negative x-direction
      mortars.node_indices[side, mortar_id] = (:begin, i)
    elseif faces[side] == 1
      # Index face in positive x-direction
      mortars.node_indices[side, mortar_id] = (:end, i)
    elseif faces[side] == 2
      # Index face in negative y-direction
      mortars.node_indices[side, mortar_id] = (i, :begin)
    else # faces[side] == 3
      # Index face in positive y-direction
      mortars.node_indices[side, mortar_id] = (i, :end)
    end
  end

  return mortars
end

end # @muladd


text:
using PiGPIOMEM
using LLVM
#using BaremetalPi

nop() = LLVM.Interop.@asmcall("nop")
spin_ns(n) = for i in 1:n nop() end

function test(pin)

    for i in 1:10_000_000
        clear(pin)
        spin_ns(80)
        set(pin)
        spin_ns(80)
    end
end

#function test2()
#    init_gpio()
#    gpio_set_mode(18, :out)
#
#    while true
#        gpio_clear(18)
#        spin_ns(100)
#        gpio_set(18)
#        spin_ns(100)
#    end
#end

#@time spin_ns(100000 / 3)
#@time spin_ns(10000 / 3)
#@time spin_ns(1000 / 3)

test(GPIOPin(18, output=true))

#using InteractiveUtils
#code_native(test, [GPIOPin], debuginfo=:none)

"""
2 * (spin( 80) + d) =  520ns
2 * (spin(100) + d) =  632ns
2 * (spin(150) + d) =  936ns
2 * (spin(200) + d) = 1240ns

ns = 3 x 2 * spin + 34

2 * spin(100): 632ns
2 * spin(100): 632-34 = 598ns
spin(100): 598 / 2 = 300 ns
spin = 3ns
"""

# 2.1MHz
#test2()


text:
using CUDAapi

using Test

@testset "CUDAapi" begin

@testset "library types" begin
    @test CUDAapi.PATCH_LEVEL == CUDAapi.libraryPropertyType(2)
    @test CUDAapi.C_32U == CUDAapi.cudaDataType(Complex{UInt32})
end

@testset "discovery" begin
    CUDAapi.find_binary([Sys.iswindows() ? "CHKDSK" : "true"])
    CUDAapi.find_library([Sys.iswindows() ? "NTDLL" : "c"])

    dirs = find_toolkit()
    @test !isempty(dirs)

    if haskey(ENV, "CI")
        # CI deals with plenty of CUDA versions, which makes discovery tricky.
        # dump a relevant tree of files to help debugging
        function traverse(dir, level=0)
            for entry in readdir(dir)
                print("  "^level)
                path = joinpath(dir, entry)
                if isdir(path)
                    println("└ $entry")
                    traverse(path, level+1)
                else
                    println("├ $entry")
                end
            end
        end
        for dir in dirs
            println("File tree of toolkit directory $dir:")
            traverse(dir)
        end
    end

    if VERSION < v"1.1.0-DEV.472"
        isnothing(::Any) = false
        isnothing(::Nothing) = true
    else
        isnothing = Base.isnothing
    end

    @testset "CUDA tools and libraries" begin
        @test !isnothing(find_cuda_binary("ptxas", dirs))
        ptxas = find_cuda_binary("ptxas", dirs)
        ver = parse_toolkit_version(ptxas)
        @test !isnothing(find_cuda_library("cudart", dirs, [ver]))
        @test !isnothing(find_cuda_library("cupti", dirs, [ver]))   # in extras
        @test !isnothing(find_cuda_library("nvtx", dirs, [v"1"]))   # somewhere else
        @test !isnothing(find_libdevice(dirs))
        @test !isnothing(find_libcudadevrt(dirs))
    end
end

@testset "availability" begin
    @test isa(has_cuda(), Bool)
    @test isa(has_cuda_gpu(), Bool)
    @test isa(usable_cuda_gpus(), Int)
end

@testset "call" begin
    # ccall throws if the lib doesn't exist, even if not called
    foo(x) = (x && ccall((:whatever, "nonexisting"), Cvoid, ()); 42)
    if VERSION < v"1.4.0-DEV.653"
        @test_throws ErrorException foo(false)
    else
        foo(false)
    end

    # @runtime_ccall prevents that
    bar(x) = (x && @runtime_ccall((:whatever, "nonexisting"), Cvoid, ()); 42)
    @test bar(false) == 42
    # but should still error nicely if actually calling the library
    @test_throws ErrorException bar(true)

    # ccall also doesn't support non-constant arguments
    lib = Ref("libjulia")
    baz() = ccall((:jl_getpid, lib[]), Cint, ())
    @test_throws TypeError baz()

    # @runtime_ccall supports that
    qux() = @runtime_ccall((:jl_getpid, lib[]), Cint, ())
    @test qux() == getpid()

    # decoding ccall/@runtime_ccall
    @test decode_ccall_function(:(ccall((:fun, :lib)))) == "fun"
    @test decode_ccall_function(:(@runtime_ccall((:fun, :lib)))) == "fun"
end

@testset "@enum_without_prefix" begin
    mod = @eval module $(gensym())
        using ..CUDAapi
        @enum MY_ENUM MY_ENUM_VALUE
        @enum_without_prefix MY_ENUM MY_
    end

    @test mod.ENUM_VALUE == mod.MY_ENUM_VALUE
end

@testset "@checked" begin
    mod = @eval module $(gensym())
        using ..CUDAapi

        const checks = Ref(0)
        macro check(ex)
            esc(quote
                $checks[] += 1
                $ex
            end)
        end

        @checked function foo()
            ccall(:jl_getpid, Cint, ())
        end
    end

    @test mod.checks[] == 0
    @test mod.foo() == getpid()
    @test mod.checks[] == 1
    @test mod.unsafe_foo() == getpid()
    @test mod.checks[] == 1
end

end


text:
mutable struct ADM{U, O} <: AbstractSubspaceOptimiser
    λ::U
    R::O
end

function ADM(λ::T = 0.1) where T <: Real
    f = NormL1(λ)
    return ADM(λ, f)
end

function fit!(q::AbstractArray{T, 1}, Y::AbstractArray, opt::ADM; maxiter::Int64= 10) where T <: Real
    normalize!(q)
    x = Y*q
    for k in 1:maxiter
        prox!(x, opt.R, Y*q)
        mul!(q, Y', x/norm(Y'*x, 2))
    end
end

function fit!(q::AbstractArray{T, 2}, Y::AbstractArray, opt::ADM; maxiter::Int64= 10) where T <: Real
    @inbounds for i in 1:size(q, 2)
        fit!(view(q, :, i), Y, opt, maxiter = maxiter)
    end
end

# Test

#using LinearAlgebra
#using SparseArrays
#using ProximalOperators

#opt = ADM(1e-2)
#A = Matrix(sprandn(20, 40, 0.1))
#U, S, V = svd(A, full = true)
#
#
#N = nullspace(A, rtol = 0.99)
#X, Y = qr(randn(40,40))
#M = X*N
#rank(M)
#Q = M[:,3]
#norm(Q, 0)
#fit!(Q, N', opt, maxiter = 10000)
#norm(Q, 0)
#fit!(M, N', opt, maxiter = 10000)


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

JLLWrappers.@generate_main_file_header("Antlr4CppRuntime")
JLLWrappers.@generate_main_file("Antlr4CppRuntime", UUID("642bff37-c5db-510f-b67c-24692d406d7f"))
end  # module Antlr4CppRuntime_jll


text:
__precompile__()

module DiffEqBiological

using Compat, DataStructures, MacroTools, Parameters, Reexport, SparseArrays, SymEngine
using DiffEqBase, DiffEqJump
@reexport using DiffEqBase, DiffEqJump

import Base: (==)

const ExprValues = Union{Expr,Symbol,Float64,Int}

include("expression_utils.jl")
include("reaction_network.jl")
include("maketype.jl")
include("network_properties.jl")
include("massaction_jump_utils.jl")
include("problem.jl")

# reaction network macro
export @reaction_network, @reaction_func, @min_reaction_network, @empty_reaction_network

# functions to query network properties
export species, params, speciesmap, paramsmap, numspecies, numreactions, numparams
export oderhsfun, jacfun, paramjacfun, odefun, noisefun, sdefun, jumps, regularjumps
export odeexprs, jacobianexprs, noiseexprs, jumpexprs, rateexpr, oderatelawexpr, ssaratelawexpr
export substratestoich, productstoich, netstoich, ismassaction, dependants, dependents, substrates, products, substratesymstoich, productsymstoich
export rxtospecies_depgraph, speciestorx_depgraph, rxtorx_depgraph

# functions to extend empty_ and min_ reaction_networks
export addspecies!, addparam!, add_scale_noise_param!, addreaction!

# functions to add mathematical equations to the network
export addodes!, addsdes!, addjumps!

# problems that can be solved from the network
export ODEProblem, SDEProblem, DiscreteProblem, JumpProblem, SteadyStateProblem

end # module


text:
using BenchmarkTools
using Rotations
using StaticArrays
using Statistics
using Random
using Test

Rotations.params(q::Quat) = @SVector [q.w, q.x, q.y, q.z]
Rotations.params(spq::SPQuat) = @SVector [spq.x, spq.y, spq.z]

const comp = BenchmarkGroup()
comp["composition"] = BenchmarkGroup()
comp["rotation"] = BenchmarkGroup()
Random.seed!(1)

# ~~~~~~~~~~~~~~~~~~~~ Create two of each type of Quaternion ~~~~~~~~~~~~~~~~~~~~ #
q1 = rand(Quat)
q2 = rand(Quat)

u1 = UnitQuaternion(q1);
u2 = UnitQuaternion(q2);

# Test Equality
@test q1 == u1
@test u1 == q1

# Test composition
@test q1*q2 ≈ u1*u2
comp["composition"]["Quat"]           = @benchmarkable for i = 1:100; $q1 = $q1*$q2; end
comp["composition"]["UnitQuaternion"] = @benchmarkable for i = 1:100; $u1 = $u1*$u2; end

# Test rotation
r = @SVector rand(3)
@test q1*r ≈ u1*r
comp["rotation"]["Quat"]           = @benchmarkable for i = 1:100; $r = $q1*$r; end
comp["rotation"]["UnitQuaternion"] = @benchmarkable for i = 1:100; $r = $u1*$r; end

# ~~~~~~~~~~~~~~~~~~~~ MRP vs SPQuat ~~~~~~~~~~~~~~~~~~~~ #
p1 = rand(MRP)
p2 = rand(MRP)
s1 = SPQuat(params(p1)...)
s2 = SPQuat(params(p2)...)
@test p1 ≈ s1
@test p2 ≈ s2
@test params(p1) ≈ params(s1)
@test params(p2) ≈ params(s2)

# Test composition
@test p2*p1 ≈ s2*s1
comp["composition"]["SPQuat"] = @benchmarkable for i = 1:100; $s1 = $s1*$s2; end
comp["composition"]["MRP"]    = @benchmarkable for i = 1:100; $p1 = $p1*$p2; end

# Test rotation
r = @SVector rand(3)
@test p1*r ≈ s1*r
comp["rotation"]["SPQuat"]  = @benchmarkable for i = 1:100; $r = $s1*$r; end
comp["rotation"]["MRP"]     = @benchmarkable for i = 1:100; $r = $p1*$r; end

tune!(comp, verbose=true)
res = run(comp, verbose=true)

for (r1,r2) in ((Quat,UnitQuaternion),(SPQuat,MRP))
    println(r1, " vs ", r2)
    for key in keys(res)
        b1 = res[key][string(r1)]
        b2 = res[key][string(r2)]
        b = judge(median(b2), median(b1))
        println("  ", key, ": ", b)
    end
end


text:
using OrderedCollections

# Initialize model object
m = AnSchorfheide(testing = true)
m <= Setting(:forecast_horizons, 12)

estroot = joinpath(dirname(@__FILE__), "..", "reference")
overrides = forecast_input_file_overrides(m)
overrides[:mode] = joinpath(estroot, "optimize.h5")

# Zero out measurement error (or else filtering shocks won't work exactly)
for para in [:e_y, :e_π, :e_R]
    m[para].value = 0
end

# Simulate dummy scenarios
def = Scenario(:defscen, "Test Default Scenario", [:obs_gdp, :obs_cpi], Symbol[], "REF")
alt = Scenario(:altscen, "Test Alternative Scenario", [:obs_gdp, :obs_cpi], [:g_sh, :rm_sh], "REF")

forecast_scenario(m, def, verbose = :none)
forecast_scenario(m, alt, verbose = :none)

expect = JLD2.jldopen(get_scenario_input_file(m, alt), "r") do file
    read(file, "arr")
end

JLD2.jldopen(get_scenario_output_files(m, alt, [:forecastobs])[:forecastobs], "r") do file
    observable_indices = read(file, "observable_indices")
    global inds = map(var -> observable_indices[var], def.target_names)
    global actual = read(file, "arr")[:, inds, :]
end

@testset "Test dummy scenarios" begin
    @test expect ≈ actual
end

# Simulate always switching
probs_enter = vcat([1], zeros(11))
probs_exit  = zeros(12)
allalt = SwitchingScenario(:allalt, alt, def, probs_enter, probs_exit)
simulate_switching(m, allalt, verbose = :none)

original_draws = JLD2.jldopen(get_scenario_output_files(m, alt, [:forecastobs])[:forecastobs]) do file
    read(file, "arr")
end
actual = JLD2.jldopen(get_scenario_output_files(m, allalt, [:forecastobs])[:forecastobs]) do file
    read(file, "arr")
end
@testset "Test switching scenarios" begin
    @test original_draws ≈ actual
end

# Simulate never switching
probs_enter = zeros(12)
probs_exit  = zeros(12)
alldef = SwitchingScenario(:alldef, alt, def, probs_enter, probs_exit)
simulate_switching(m, alldef, verbose = :none)

default_draws = JLD2.jldopen(get_scenario_output_files(m, def, [:forecastobs])[:forecastobs]) do file
    read(file, "arr")
end

actual = JLD2.jldopen(get_scenario_output_files(m, alldef, [:forecastobs])[:forecastobs]) do file
    read(file, "arr")
end

@testset "Test never switching scenarios" begin
    for i = 1:10 # Check that the ith draw of actual matches some draw j from default_draws
        matches_one = false
        for j = 1:10
            (actual[i, :, :] ≈ default_draws[j, :, :]) && (matches_one = true)
        end
        @test matches_one
    end
end

# Simulate switching in at t = 5 and out at t = 9
probs_enter = vcat(zeros(4), [1], zeros(7))
probs_exit  = vcat(zeros(8), [1], zeros(3))
somealt = SwitchingScenario(:somealt, alt, def, probs_enter, probs_exit)
simulate_switching(m, somealt, verbose = :none)

@testset "Test switching scenarios, where the switch happens at particular points in time" begin
    global actual = JLD2.jldopen(get_scenario_output_files(m, somealt, [:forecastobs])[:forecastobs]) do file
        read(file, "arr")
    end
    for i = 1:10
        @test actual[i, :, 5:8] ≈ original_draws[i, :, 1:4]
        matches_one = false
        for j = 1:10
            (actual[i, :, [1:4;9:12]] ≈ default_draws[j, :, [1:4;9:12]]) && (matches_one = true)
        end
        @test matches_one
    end
end

# Transform single scenarios
scenario_means_bands(m, alt, verbose = :none)
scenario_means_bands(m, allalt, verbose = :none)
global actual = dropdims(mean(original_draws, dims = 1), dims=1)
@testset "Test single scenarios" begin
    mb1 = read_scenario_mb(m, alt, :forecastutobs)
    mb2 = read_scenario_mb(m, allalt, :forecastutobs)
    for mb in [mb1, mb2]
        for (var, ind) in mb.metadata[:indices]
            @test mb.means[!, var] ≈ actual[ind, :]
        end
    end

    mb1 = read_scenario_mb(m, alt, :forecastobs)
    mb2 = read_scenario_mb(m, allalt, :forecastobs)
    for mb in [mb1, mb2]
        @test mb.means[!, :obs_gdp] ≈ quartertoannual(actual[1, :])
        @test mb.means[!, :obs_cpi] ≈ quartertoannual(actual[2, :])
        @test mb.means[!, :obs_nominalrate] ≈ actual[3, :]
    end

    mb1 = read_scenario_mb(m, alt, :forecast4qobs)
    mb2 = read_scenario_mb(m, allalt, :forecast4qobs)
    for mb in [mb1, mb2]
        @test mb.means[!, :obs_gdp] ≈ loggrowthtopct_4q_approx(actual[1, :], zeros(3))
        @test mb.means[!, :obs_cpi] ≈ loggrowthtopct_4q_approx(actual[2, :], zeros(3))
        @test mb.means[!, :obs_nominalrate] ≈ actual[3, :]
    end
end
# Aggregate scenarios equally
aggall = ScenarioAggregate(:aggall, "Test Scenario Aggregate", AbstractScenario[def, alt],
                           [0.5, 0.5], 20, false, "REF")
scenario_means_bands(m, aggall, verbose = :none)
global actual = dropdims(mean(cat(default_draws, original_draws, dims = 1), dims = 1), dims=1)

@testset "Test scenario aggregation" begin
    mb = read_scenario_mb(m, aggall, :forecastutobs)
    for (var, ind) in mb.metadata[:indices]
        @test mb.means[!, var] ≈ actual[ind, :]
    end
    # Aggregate scenarios to only have alternative scenaro
    aggalt = ScenarioAggregate(:aggalt, "Test Scenario Aggregate", AbstractScenario[def, alt],
                               [0.0, 1.0], 20, false, "REF")
    scenario_means_bands(m, aggalt, verbose = :none)
    global actual = dropdims(mean(repeat(original_draws, outer = [2, 1, 1]), dims = 1), dims=1)

    mb = read_scenario_mb(m, aggalt, :forecastutobs)
    for (var, ind) in mb.metadata[:indices]
        @test mb.means[!, var] ≈ actual[ind, :]
    end
end

# Aggregate scenarios, drawing with replacement
aggrep1 = ScenarioAggregate(:aggrep1, "Test Scenario Aggregate", AbstractScenario[def, alt],
                            [0.5, 0.5], 20, true, "REF")
aggrep2 = ScenarioAggregate(:aggrep2, "Test Scenario Aggregate", AbstractScenario[def, alt],
                            [0.0, 1.0], 20, true, "REF")

# Simulate dummy scenario with shock scaling
scale_scenario = Scenario(:altscen, "Test Shock Scaling Scenario", [:obs_gdp, :obs_cpi], [:g_sh, :rm_sh], "REF",
                 shock_scaling = 2.0)
forecast_scenario(m, scale_scenario, verbose = :none)

expect = JLD2.jldopen(get_scenario_input_file(m, scale_scenario), "r") do file
    2.0 * read(file, "arr")
end

JLD2.jldopen(get_scenario_output_files(m, alt, [:forecastobs])[:forecastobs], "r") do file
    observable_indices = read(file, "observable_indices")
    global inds = map(var -> observable_indices[var], def.target_names)
    global actual = read(file, "arr")[:, inds, :]
end
@testset "Test scenario aggregation, drawing with replacement" begin
    @test expect ≈ actual
end

# Delete outputted files
#=
raw_fn_root  = "../reference/output_data/an_schorfheide/ss0/scenarios/raw/"
work_fn_root = "../reference/output_data/an_schorfheide/ss0/scenarios/work/"
rm(joinpath(dirname(@__FILE__), raw_fn_root, "forecastobs_scen=altscen_svin=REF_test.jld2"))
rm(joinpath(dirname(@__FILE__), raw_fn_root, "forecastpseudo_scen=altscen_svin=REF_test.jld2"))
rm(joinpath(dirname(@__FILE__), work_fn_root, "mbforecast4qobs_scen=altscen_svin=REF_test.jld2"))
rm(joinpath(dirname(@__FILE__), work_fn_root, "mbforecast4qpseudo_scen=altscen_svin=REF_test.jld2"))
rm(joinpath(dirname(@__FILE__), work_fn_root, "mbforecastobs_scen=altscen_svin=REF_test.jld2"))
rm(joinpath(dirname(@__FILE__), work_fn_root, "mbforecastpseudo_scen=altscen_svin=REF_test.jld2"))
rm(joinpath(dirname(@__FILE__), work_fn_root, "mbforecastutobs_scen=altscen_svin=REF_test.jld2"))
rm(joinpath(dirname(@__FILE__), work_fn_root, "mbforecastutpseudo_scen=altscen_svin=REF_test.jld2"))
=#
nothing


text:
"""
**Both latitude and longitude are present**

This filter will only retain occurrences that have *both* a latitude and a
longitude field.
"""
function have_both_coordinates(o::GBIFRecord)
  !(ismissing(o.latitude)&ismissing(o.latitude))
end

"""
**Neither latitude nor longitude are 0.0**

This filter will only retain occurrences that are *not* at the (0,0) coordinate.
"""
function have_neither_zero_coordinates(o::GBIFRecord)
  !((o.latitude==0.0)&(o.longitude==0.0))
end

"""
**At most one of latitude and longitude is 0.0**

This filter will only retain occurrences that have *at most* one coordinate
being exactly 0.0.
"""
function have_no_zero_coordinates(o::GBIFRecord)
  !((o.latitude==0.0)|(o.longitude==0.0))
end

"""
**The date is not missing**

This filter will remove observations with a missing date.
"""
function have_a_date(o::GBIFRecord)
  !ismissing(o.date)
end

"""
**No known issues at all**

This filter will retain no observation. At least not in theory, but it is very
rare to have GBIF records with absolutely no issues. Its use is discouraged.
"""
function have_no_issues(o::GBIFRecord)
  length(o.issues) == 0
end

"""
**No coordinates issues except rounding and WGS84 assumption**

This filter will retain observations that have no issues, with the exception of
assuming WGS84 projection (which a large number of GBIF records do), and rounded
coordinates (same thing). It is a reasonable filter for most use cases.
"""
function have_ok_coordinates(o::GBIFRecord)
  ok = true
  for i in o.issues
    if !(i in [:COORDINATE_ROUNDED, :GEODETIC_DATUM_ASSUMED_WGS84])
      ok = false
    end
  end
  return ok
end

"""
**Filters a series of records**

This function will take the filter function `f` and use it to *mask* the records
that do not satisfy it. Note that if a record is *already* masked due to the
application of a previous filter, its status will *not* be modified. The
application of filters is therefore cumulative.

Note that while the usual `filter!` function *removes* objects, this one will
only *mask* them. The objects can be unmasked with `allrecords!`. This design
choice was made because `GBIFRecords` track the position in the API pages, and
it seemed safer to keep all information. Note that when calling `length` or
iterating over a `GBIFRecords` object, only the *unmasked* records will be
shown.

It is possible to apply multiple filters at once, using the following syntax:

    filter!.([f1, f2, f3], records)

"""
function Base.filter!(f, o::GBIFRecords)
  keep = map(f, o.occurrences);
  o.show = o.show .& keep;
end

"""
**Make records broadcastable**
"""
Base.broadcastable(x::GBIFRecords) = Ref(x)

"""
**Cleans a search output**

This function is deprecated and will be removed in a later release, use
`filter!` instead.
"""
function qualitycontrol!(o::GBIFRecords; filters::Array{T,1}=[have_no_issues], verbose::Bool=true) where {T<:Function}
  @warn "The qualitycontrol! function is deprecated and will be removed in a future release -- use filter! instead."
  keep = ones(Bool, length(o.occurrences))
  if verbose
    @info "Starting quality control with $(length(o.occurrences)) records"
  end
  for f in filters
    keep_f = map(f, o.occurrences)
    keep = keep .* keep_f
    if verbose
      @info "$(count(keep)) records left after $(f)"
    end
  end
  o.show = keep
end

"""
**Show all occurrences**

This function reverses the action of `filter!`. It will unmask all
records that have been removed under the current filters.
"""
function showall!(o::GBIFRecords)
  @warn "The showall! function is deprecated and will be removed in a future release -- use allrecords! instead."
  allrecords!(o)
end

"""
**Show all occurrences**

This function reverses the action of `filter!`. It will unmask all
records that have been removed under the current filters.
"""
function allrecords!(o::GBIFRecords)
  o.show = ones(Bool, length(o.occurrences))
end


text:
using IncludeGuards

println("TestModuleWorking.jl, module = ", @__MODULE__)

@includeonce "ModB.jl"
@includeonce "ModA.jl"

module TestModuleWorking

using ..A
using ..B

p = A.Parameters(1)

b = B.BStruct(params = p)

B.g(b)

end # module

max_stars_repo_name: tgflynn/IncludeGuards

max_issues_repo_name: tgflynn/IncludeGuards

max_forks_repo_name: tgflynn/IncludeGuards


text:
"""
    EOS

Mark the end of the string to parse (helps with corner cases where a token
ends a document without being followed by a space).
"""
const EOS = '\0'

"""
    SPACE_CHARS

List of characters that correspond to a `\\s` regex + EOS.

Ref: https://github.com/JuliaLang/julia/blob/master/base/strings/unicode.jl.
"""
const SPACE_CHAR = [' ', '\r', '\n', '\t', '\f', '\v', EOS]

"""
    NUM_CHAR

Convenience list of characters corresponding to digits.
"""
const NUM_CHAR = ['0':'9'...]

"""
    ALPHA_LATIN

Convenience list of characters corresponding to letters a-zA-Z.
"""
const ALPHA_LATIN = ['a':'z'..., 'A':'Z'...]

"""
    ALPHANUM_LATIN

Convenience list of characters corresponding to a-zA-Z0-9.
"""
const ALPHANUM_LATIN = vcat(ALPHA_LATIN, NUM_CHAR)

"""
    ALPHA_ALL

All 10_000 first characters.
"""
const ALPHA_ALL = [Char(i) for i in 1:10_000 if isletter(Char(i))]

"""
    ALPHANUM_ALL

ALPHA_ALL and digits.
"""
const ALPHANUM_ALL = vcat(ALPHA_ALL, NUM_CHAR)

"""
    Chomp

Structure to encapsulate rules around a token such as whether it's fine at the
end of a string, what are allowed following characters and, in the greedy
case, what characters are allowed.
"""
struct Chomp
    # fixed style
    refstring::String
    ok_at_eos::Bool
    is_followed::Bool
    next_chars::Vector{Char}
    # greedy style
    head_chars::Vector{Vector{Char}}
    tail_chars::Vector{Char}
end
function Chomp(;
            refstring::String="",
            ok_at_eos::Bool=true,
            is_followed::Bool=true,
            next_chars::Vector{Char}=Char[],
            head_chars::Vector{Vector{Char}}=Vector{Vector{Char}}(),
            tail_chars::Vector{Char}=Char[])
    Chomp(
        refstring, ok_at_eos, is_followed,
        next_chars, head_chars, tail_chars)
end

"""
    TokenFinder

Structure to find a token keeping track of how many characters should be seen,
some rules with respect to positioning or following chars (see Chomp) and
possibly a validator that checks whether a candidate respects a rule.
"""
struct TokenFinder
    steps::Int
    chomp::Chomp
    check::Regex
end
TokenFinder(s::Int, c::Chomp) = TokenFinder(s, c, r"")

"""
    fixed_lookahead(tokenfinder, candidate, at_eos)

Applies a fixed lookahead step corresponding to a token finder.
This is used as a helper function in `find_tokens`.
"""
function fixed_lookahead(tf::TokenFinder, candidate::SS, at_eos::Bool)
    c = tf.chomp
    matches = false
    if at_eos
        matches  = (candidate == c.refstring)
        matches &= c.ok_at_eos
    elseif isempty(c.next_chars)
        matches  = candidate == c.refstring
    else
        matches  = chop(candidate) == c.refstring
        matches &= !xor(c.is_followed, candidate[end] ∈ c.next_chars)
    end
    offset = Int(!isempty(c.next_chars) & !at_eos)
    return matches, offset
end

"""
    greedy_lookahead(tokenfinder, nchars, probe_char)

Applies a greedy lookahead step corresponding to a token finder.
This is used as a helper function in `find_tokens`.
"""
function greedy_lookahead(tf::TokenFinder, nchars::Int, probe_char::Char)
    c = tf.chomp
    (nchars > length(c.head_chars)) && return (probe_char in c.tail_chars)
    return (probe_char in c.head_chars[nchars])
end

"""
    check(tokenfinder, ss)

Check whether a substring verifies the regex of a token finder.
"""
function check(tf::TokenFinder, ss::SS)
    return match(tf.check, ss) !== nothing
end

"""
    forward_match(refstring, next_chars, is_followed)

Return a TokenFinder corresponding to a forward lookup checking if a sequence
of characters matches a `refstring` and is followed (or not followed if
`is_followed==false`) by a char out of a list of chars (`next_chars`).
"""
@inline function forward_match(
            refstring::String,
            next_chars::Vector{Char}=Char[],
            is_followed::Bool=true
            )::TokenFinder
    # steps keeps track of the number of character to consume
    steps = lastindex(refstring)
    steps = ifelse(isempty(next_chars), prevind(refstring, steps), steps)

    # check whether it would be allowed to be at the end of the string
    ok_at_eos::Bool = isempty(next_chars) || !is_followed || EOS ∈ next_chars

    # form the validator and return the tokenfinder
    chomp = Chomp(; refstring, ok_at_eos, is_followed, next_chars)
    return TokenFinder(steps, chomp)
end

"""
    greedy_match(head_chars, tail_chars, check)

Lazily accept the next char and stop as soon as it fails to verify `λ(c)`.
"""
@inline function greedy_match(;
            head_chars::Vector{Vector{Char}}=Vector{Vector{Char}}(),
            tail_chars::Vector{Char}=Char[],
            check::Regex=r"")
    TokenFinder(-1, Chomp(; head_chars, tail_chars), check)
end


text:
module HELICS

__precompile__(true)

import Libdl

# Module initialization function
function __init__()
    atexit(helicsCloseLibrary)
end

module Lib

    using CEnum
    using HELICS_jll

    const HELICS_EXPORT = nothing
    const HELICS_NO_EXPORT = nothing

    include("ctypes.jl")

    include("common.jl")

    include("lib.jl")

end

include("wrapper.jl")

include("utils.jl")

include("enums.jl")
include("api.jl")

const HELICS_LIBRARY_VERSION = VersionNumber(split(helicsGetVersion())[1])
const HELICS_PACKAGE_VERSION = VersionNumber(match(r"version\s*=\s*\"(.+)\"", read(joinpath(@__DIR__, "..", "Project.toml"), String)).captures[1])

end # module


text:
abstract type AbstractLogger end

"Return true if `logger` is instantiated."
isdefined(::AbstractLogger) = false

"Reinitialize logger."
empty!(::AbstractLogger) = nothing

struct NoLogger <: AbstractLogger end

"""
Store evolution of objective, gradient and time across iterations.

"""
mutable struct Tracer{T <: Real} <: AbstractLogger
    obj::Vector{T}
    grad::Vector{T}
    time::Vector{T}
end
Tracer{T}() where T = Tracer(T[], T[], T[])

function empty!(tracer::Tracer{T}) where T
    tracer.obj = T[]
    tracer.grad = T[]
    tracer.time = T[]
end

isdefined(t::Tracer) = true


text:
using BenchmarkTools
using SequentialRaytrace
import SequentialRaytrace: gen_result, update_result!
using TimerOutputs
using StaticArrays

# 12.881 us AbstractVector
# 12.796 us Vector
# 1.74 SVector OpticalComponent
function testlens()
    Lens("", Object(Air, 200.0),
        SVector(
        #  [
             OpticalComponent(:a, Nothing, nothing,
            # [
             SVector(
                 Sphere(1/50, nothing, Silica, 15.0, :first_surface),
                 EvenAsphere(-1/50, 0.0, 0.0, 0.0, 0.0, nothing, Air, 65.0),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 ),
                 Plano(nothing, Silica, 15.0 ),
                 Plano(nothing, Air, 15.0 )
            )
            # ]
            )
        #  ]
        )
         )
end

function testray()
    x = 0.0
    y = 20.0
    z = 0.0
    X = 0.1
    Y = -0.1

    Ray(x, y, z, X, Y)
end

function test()
    x = 0.0
    y = 20.0
    z = 0.0
    X = 0.1
    Y = -0.1

    Z = sqrt(1.0 - X^2 - Y^2)
    ray = make_ray(x, y, z, X, Y, Z)
    lens = testlens()
    result = gen_result(Vector{Ray}, lens)
    r1 = trace!(result, lens, ray, 1.0)
    change_lens(lens)
    result = gen_result(Vector{Ray}, lens)
    r2 = trace!(result, lens, ray, 1.0)
    (r1, r2)
    # (unwrap(r1)[end], unwrap(r2)[end])
end

const to = TimerOutput()

function warmup()
    lens = testlens()
    ray = testray()
    result = gen_result(Vector{Ray}, lens)
    trace!(result, lens, ray, 1.0)
end

function timeit()
    lens = testlens()
    ray = testray()
    result = gen_result(Vector{Ray}, lens)
    @timeit to "trace!" trace!(result, lens, ray, 1.0)
end

function bench1()
    lens = testlens()
    ray = testray()
    result = gen_result(Vector{typeof(ray)}, lens)
    @benchmark trace!($result, $lens, $ray, 1.0)
end

function run2(lens, ray)
    result = gen_result(Vector{typeof(ray)}, lens)
    trace!(result, lens, ray, 1.0)
end

function bench2()
    lens = testlens()
    ray = testray()
    @benchmark run2($lens, $ray)
end

function timeit_xx()
    @timeit to "testlens" lens = testlens()
    @timeit to "testray" ray = testray()
    @timeit to "genresult" result = gen_result(Vector{typeof(ray)}, lens)
    @timeit to "trace!" trace!(result, lens, ray, 1.0)
end

function run2b()
    lens = testlens()
    ray = testray()
    result = gen_result(Vector{typeof(ray)}, lens)
    trace!(result, lens, ray, 1.0)
end

function bench2b()
    @benchmark run2b()
end

function gen_result(::Type{Array{T, 2}}, lens) where T
    Array{T, 2}(undef, 6, 2 + SequentialRaytrace.n_surfaces(lens))
    # Array{T, 2}(undef, 2 + SequentialRaytrace.n_surfaces(lens), 6)
end

function update_result!(result :: Array{T, 2}, index, _symbol, ray) where T
    result[1, index] = ray.x
    result[2, index] = ray.y
    result[3, index] = ray.z
    result[4, index] = ray.cx
    result[5, index] = ray.cy
    result[6, index] = ray.cz
    # result[index, 1] = ray.x
    # result[index, 2] = ray.y
    # result[index, 3] = ray.z
    # result[index, 4] = ray.cx
    # result[index, 5] = ray.cy
    # result[index, 6] = ray.cz
    return result
end

function bench3()
    lens = testlens()
    ray = testray()
    result = gen_result(Array{typeof(ray.x), 2}, lens)
    @benchmark trace!($result, $lens, $ray, 1.0)
end

function run4(lens, ray)
    result = gen_result(Array{typeof(ray.x), 2}, lens)
    trace!(result, lens, ray, 1.0)
end

function bench4()
    lens = testlens()
    ray = testray()
    @benchmark run4($lens, $ray)
end

struct OnlyLast{T}
    d :: Array{T, 1}
end

function gen_result(::Type{OnlyLast{T}}, lens) where T
    OnlyLast(Array{T, 1}(undef, 6))
end

function update_result!(result :: OnlyLast{T}, _index, _symbol, ray) where T
    result.d[1] = ray.x
    result.d[2] = ray.y
    result.d[3] = ray.z
    result.d[4] = ray.cx
    result.d[5] = ray.cy
    result.d[6] = ray.cz
    return result
end

function benchb()
    lens = testlens()
    ray = testray()
    result = gen_result(OnlyLast{typeof(ray.x)}, lens)
    @benchmark for i in range(1, stop=100000) trace!($result, $lens, $ray, 1.0) end
end

function timeit4b()
    @timeit to "testlens" lens = testlens()
    @timeit to "testray" ray = testray()
    @timeit to "genresult" result = gen_result(OnlyLast{typeof(ray.x)}, lens)
    @timeit to "trace!" trace!(result, lens, ray, 1.0)
end

function change_lens(lens)
    lens.components[1].surfaces[1] = sphere(1/40, nothing, Silica, 15.0)
end


text:
module benchmark

__precompile__(true)
reload("j4pr")
using DataArrays
import j4pr

function benchmark_training(a::T where T<:j4pr.CellData, w::S where S<:j4pr.CellFun)

	a[:,1:2] |> w

	@printf("* %s:", w.tinfo)
	@time out = a|>w
end

function benchmark_execution(a::T where T<:j4pr.CellData, w::S where S<:j4pr.AbstractCell)

	out1 = [];
	out2 = [];
	out3 = [];

	adc = deepcopy(a)
	ad = DataArray(deepcopy(a))
	aa = deepcopy(j4pr.getx(a))

	@printf("* %s ...\n", w.tinfo)

	# Test datacell interface
	@printf("\tDataCell |> :")
	#try
		deepcopy(adc[:,1:10]) |> w;
		@time out1 = adc |> w
	#catch
	#	@printf("FAIL\n")
	#end

	# Test DataArray interface
	#@printf("\tDataArray |> :")
	#try
	#	deepcopy(ad[:,1:10]) |>w
	#	@time out2 = ad |> w
	#catch
	#	@printf("FAIL\n")
	#end

	# Test Array interface
	@printf("\tArray |> :")
	try
		deepcopy(aa[:,1:10]) |>w
		@time out3 = aa |> w
	catch
		@printf("FAIL\n")
	end

	# Test returned values
	try
		all(out1.x.==out3) ? @printf("\tValues check OK.") : begin
			@printf("\tValues check failed.")
			# Specify which version(s) disagree ...
		end
	catch
		@printf("\tCould not execute values check.")
	end

	@printf("\n\n")
	return out1, out2, out3
end

function benchmark_execution_pipe(a::T where T<:j4pr.CellData, w::S where S<:j4pr.PipeGeneric)
	@printf("\n\n=============== Testing %s =========================\n",w.tinfo)
	@printf("Compiling pipe elements ...");
	a[:,1] |> w;

	@printf("Passing %s data through pipe ...", typeof(a));
	@time a |> w

end

function run_benchmarks_basic()

	@printf("\n\nRunning some basic tests ...");

	# Get some dataset
	#A= dgenr("datasets", "iris")

	#A = datacell( rand(10000,100), round(rand(10000)*20) )
	#A=datacell( [[1,2,1,1,1] [1,1,0,0,0.0] ["a","v","x","v","a"]] )

	M=100;
	N=100;
	C = 3;
	A = j4pr.datacell( rand(M, N), round.((C-1)*rand(N)) );

	# Execution only cells
	W = [		j4pr.cslice([0,1]),
      			j4pr.filterdomain!(Dict(i=>[0 10] for i in size(A,1))),
			j4pr.sample(0.5),
	     ]
	# Trainable cells
	Wt = [		j4pr.scaler!("mean"),
			j4pr.filterg("mean"),
			j4pr.ohenc("binary"),
			j4pr.ohenc("integer"),
	     ]

	# Benchmark training for trainable cells
	@printf("\n\n=============== Training [Basic] =========================\n")
	out = Any[]
	for (i,w) in enumerate(Wt)
		push!(out, benchmark_training(A,w))
	end

	# Benchmark execution for executable only cells
	@printf("\n\n=============== Testing [Basic] =========================\n")
	for w in W
		#@show w.m
		benchmark_execution(A,w)
	end

	# Benchmark execution for trainable cells
	for w in out
		benchmark_execution(A,w)
	end

end

function run_benchmarks_pipe()

	@printf("\n\nRunning some more advanced tests ...");

	# Get some dataset
	#A= dgenr("datasets", "iris")

	#A = datacell( rand(100,10000), round.(rand(10000)*20) )
	#A=datacell( [[1,2,1,1,1] [1,1,0,0,0.0] ["a","v","x","v","a"]] )
	M=100;
	N=10000;
	C = 3;
	A = 	j4pr.datacell( round.(50*rand(M, N)), round.((C-1)*rand(N)) ); # training dataset
	At = 	j4pr.datacell( round.(60*rand(M, N)), round.((C-1)*rand(N)) ); # testing dataset

	#Build pipe

	Wclass =	j4pr.cslice([0,1])
	Wdom = 		j4pr.filterdomain!(Dict(i=>[0 1] for i in 1:50))
	Wprocmv = 	j4pr.filterg("majority")
	Wremove =	j4pr.FunctionCell(getindex, (setdiff(1:100,collect(1:5)),:),"Data remover")
	Wsample =	j4pr.sample(0.5)
	Wuscale = 	j4pr.scaler!(Dict(i=>"mean" for i in 1:50))
	Wbu =		j4pr.ohenc("binary")
	Wiu	=	j4pr.ohenc("integer")

	# Train pipe
	@printf("\nPreliminary training ...\n")
	Wb = A |> Wbu
	PIPE = deepcopy(A) |> (deepcopy(A) |> (Wb+Wremove+Wdom+Wprocmv+Wuscale+Wsample))

	benchmark_execution(At,PIPE)
	return A, At, PIPE
end

run_benchmarks_basic();

a,at,pipe = run_benchmarks_pipe()

end


text:
@doc raw"""
    Ising{T<:Int} <: AbstractGraphPointProcess{T}

The [Ising model](https://en.wikipedia.org/wiki/Ising_model) characterizes a point process defined on the vertices of a graph ``(V, E)`` with joint density proportional to

```math
    \mathbb{P}\!\left[ \mathcal{X}=X \right]
    \propto
    \prod_{i \in V}
        \exp(h_i x_i)
    \prod_{\{i, j\} \in E}
        \exp(J x_i x_j)
```

where ``(h_i)_{i\in V}`` are called magnetization paremeters and ``J`` the interaction coefficient (``J \gtreqless 0`` characterizes ferro/antiferro magnetic interactions).

# Example

A realization from a ``5\times 5`` grid graph, with ``h = 0.2, J = 0.1``.

![assets/ising_grid_prs.png](assets/ising_grid_prs.png)
"""
struct Ising{T<:Int} <: AbstractGraphPointProcess{T}
    graph::LG.SimpleGraph{T}
    "Interaction"
    J::Float64
    "Magnetization"
    h::Union{Float64,Vector{Float64}}
end

function Base.show(io::IO, pp::Ising{T}) where {T}
    print(io, "Ising{$T}\n- graph = $(pp.graph)\n- J = $(pp.J) (interaction)\n- h = $(pp.h) (magnetization)")
end

"""
    Ising(
        dims::Vector{T},
        J::Real,
        h::Union{Real,Vector{Real}}=0;
        periodic::Bool=true
    ) where {T<:Int}

Construct [`PRS.Ising`](@ref) on a grid graph with dimension `dims`, interaction parameter `J` and magnetization `h`.

Periodic boundary conditions on the grid graph are set according to `periodic`.

```jldoctest; output = true
using PartialRejectionSampling

dims = [5, 5]
J, h = 0.01, 0
PRS.Ising(dims, J, h; periodic=true)

# output

Ising{Int64}
- graph = {25, 50} undirected simple Int64 graph
- J = 0.01 (interaction)
- h = 0.0 (magnetization)
```
"""
function Ising(
    dims::Vector{T},
    J::Real,
    h::Union{Real,Vector{Real}}=0;
    periodic::Bool=true
) where {T<:Int}
    g = LG.grid(dims; periodic=periodic)
    return Ising(g, J, float.(h))
end

"""
    Ising(
        graph::LG.SimpleGraph{T},
        J::Real,
        h::Union{Real,Vector{Real}}=0
    ) where {T<:Int}

Construct [`PRS.Ising`](@ref) on `graph` with interaction parameter `J` and magnetization `h`.

```jldoctest; output = true
using PartialRejectionSampling
using LightGraphs; const LG = LightGraphs

graph = LG.grid([5, 5])
J, h = 0.01, 0
PRS.Ising(graph, J, h)

# output

Ising{Int64}
- graph = {25, 40} undirected simple Int64 graph
- J = 0.01 (interaction)
- h = 0.0 (magnetization)
```
"""
function Ising(
    graph::LG.SimpleGraph{T},
    J::Real,
    h::Union{Real,Vector{Real}}=0
) where {T<:Int}
    h isa Vector && @assert length(h) == LG.nv(graph)
    return Ising{T}(graph, J, float.(h))
end

## Sampling

@doc raw"""
    generate_sample_prs(
        [rng::Random.AbstractRNG,]
        pp::Ising
    )

Generate an exact sample form `pp` using Partial Rejection Sampling (PRS), see Section 4.2 of [GuJeLi19](@cite).

Default sampler is [`PRS.generate_sample_grid_prs`](@ref).

**See also**

- [FeViYi19](@cite),
- [FeGuYi19](@cite) [`PRS.generate_sample_gibbs_perfect`](@ref).
"""
generate_sample_prs(rng::Random.AbstractRNG, pp::Ising) = generate_sample_grid_prs(rng, pp)
# generate_sample_prs(pp::Ising) = generate_sample_prs(Random.default_rng(), pp)

@doc raw"""
    generate_sample(
        [rng::Random.AbstractRNG,]
        pp::Ising{T},
        idx::T
    )::T where {T<:Int}

Generate an exact sample from the marginal distribution of state ``i=`` `idx` of `pp`.

More specifically,

```math
    x_i
        \sim
        \operatorname{Bernoulli}_{-1, 1}
            (\sigma(h_i)),
```

where ``\sigma`` denotes the [sigmoid function](https://en.wikipedia.org/wiki/Sigmoid_function).
"""
function generate_sample(
    rng::Random.AbstractRNG,
    pp::Ising{T},
    idx::T
)::T where {T<:Int}
    hᵢ = pp.h isa Real ? pp.h : pp.h[idx]
    return rand(rng) < sigmoid(hᵢ) ? one(T) : -one(T)
end

# function generate_sample(
#     pp::Ising{T},
#     idx::T
# )::T where {T<:Int}
#     return generate_sample(Random.default_rng(), pp, idx)
# end

@doc raw"""
    generate_sample!(
        rng::Random.AbstractRNG,
        state::AbstractVector{T},
        indices,
        ising::Ising{T}
    ) where {T<:Int}

Generate an exact sample from the marginal distribution of each states of the [`PRS.Ising`](@ref) model `ising` at the prescribed `indices`.
"""
function generate_sample!(
    rng::Random.AbstractRNG,
    state::AbstractVector{T},
    indices,
    ising::Ising{T}
) where {T<:Int}
    for i in indices
        state[i] = generate_sample(rng, ising, i)
    end
end

@doc raw"""
    generate_sample_conditional!(
        rng::Random.AbstractRNG,
        state::AbstractVector{T},
        i::T,
        ising::Ising{T}
    ) where {T<:Int}

Generate an exact sample from the conditional distribution of the state ``x_i`` given its neighboring states in `ising.graph`.
More specifically,
```math
    x_i \mid x_{N(i)}
        \sim
        \operatorname{Bernoulli}_{-1, 1}
            (\sigma(h_i + J \sum_{j \in N(i)} x_j)),

where ``\sigma`` denotes the [sigmoid function](https://en.wikipedia.org/wiki/Sigmoid_function).
```
"""
function generate_sample_conditional!(
    rng::Random.AbstractRNG,
    state::AbstractVector{T},
    i::T,
    ising::Ising{T}
) where {T<:Int}
    hᵢ = ising.h isa Number ? ising.h : ising.h[i]
    proba = sigmoid(hᵢ + ising.J * sum(state[j] for j in LG.neighbors(ising.graph, i)))
    state[i] = rand(rng) < proba ? 1 : -1
end

include("ising_sampling_gibbs_perfect.jl")
include("ising_sampling_grid_prs.jl")


text:
using DelimitedFiles, Random, ScoreDrivenModels, Statistics

# Load historical Consumer Price Index data
y = vec(readdlm("../test/data/cpichg.csv"))

# Set RNG seed to guarantee consistent results
Random.seed!(123)

# Specify GAS model: a student's t model with location scale transformation
# (see /src/distributions/non_native_dists.jl in the repository)
gas = Model(1, 1, TDistLocationScale, 0.0, time_varying_params=[1, 2])

# Fit specified model to historical data
f = fit!(gas, y)

# Next, we show examples using nondefault keyword arguments
# Note that we need to re-define `gas` since its parameters were populated by `fit!`
gas = Model(1, 1, TDistLocationScale, 0.0, time_varying_params=[1, 2])

# In this example, we set verbose to 2, prompting the package to output optimization results
f = fit!(gas, y; verbose=2)

# Re-define `gas` once again for another example
gas = Model(1, 1, TDistLocationScale, 0.0, time_varying_params=[1, 2])

# In this example, we set the optimization method to LBFGS with 5 initial points
f = fit!(gas, y; opt_method=LBFGS(gas, 5))

# Print estimation statistics
fit_stats(f)

# Perform forecast via simulations for 12 time periods ahead
forec = forecast(y, gas, 12)

# We can access the parameter forecasts in `forec.parameter_forecasts`
forec.parameter_forecast

# Similarly, we can access the simulated observation scenarios
forec.observation_scenarios


text:
"""
    replace_newlines(s::AbstractString, replacement = " ")

Replace newlines with something other string. This is useful when you want to
"wrap" multiple lines into a single line. Hence the default replacement is just
a blank space.
"""
function replace_newlines(s::AbstractString, replacement = " ")
    return replace(s, r"\n" => replacement)
end

"""
    @discord(ex)

Evaluate any Discord.jl function that is expected to return a future.
The response is captured and logged at Debug level before the
underlying value is returned.
"""
macro discord(ex)
    return quote
        future = $(esc(ex))
        response = fetch(future)
        @debug "Response received" response
        response.val
    end
end

"""
    audit(source, user_id, name, discriminator, message; audit_dir, audit_file)

Create an audit record in the system for a given event.

# Arguments
- `source`: unique string representing where the audit log comes from
- `user_id`: Discord user id
- `name`: Discord user name
- `discriminator`: Discord user discriminator
- `audit_dir`: audit log directory, defaults to `data/audit`
- `audit_file`: audit log file, defaults to `\$audit_dir/\$source.log`
"""
function audit(
    source::AbstractString,
    user_id::UInt64,
    name::AbstractString,
    discriminator::AbstractString,
    message::AbstractString;
    audit_dir = joinpath("data", "audit"),
    audit_file = joinpath(audit_dir, "$source.log")
)
    open(ensurepath!(audit_file); write = true, append = true) do io
        println(io, "$(now(tz"UTC"))\t$name#$discriminator\t$user_id\t$message")
    end
end

"""
    load_docs(filename)::Dict{String, Dict{String, Vector{String}}}

Return the documentation info contained in the given JSON file, which should
be in the form of a Dict, where the keys are the list of recorded packages,
and where each value is another Dict, which contains all the names in the
package as keys, with the values being a vector, where the first element
indicates whether the name is "exported" or "nonexported" from the package,
and the second element contains the corresponding docstring.
"""
function load_docs(filename)::Dict{String, Dict{String, Vector{String}}}
    all_docs = JSON.parsefile(filename)
    all_docs = convert(Dict{String, Dict{String, Vector{String}}}, all_docs)
    return all_docs
end

"""
    load_names(filename::String)::Dict{String,Dict{String,String}}

Return a Dict, read from the given JSON file, that contains the list of names
and the associated info from the packages each name appears in. The info
is another Dict, where the keys are the packages where the name appears
in, and the values indicate whether the name is "exported" from the packaged
or "nonexported".
"""
function load_names(filename::String)::Dict{String,Dict{String,String}}
    all_names = JSON.parsefile(filename)
    all_names = convert(Dict{String,Dict{String,String}}, all_names)
    return all_names
end

function update_names_count(
        source::AbstractString,
        name::AbstractString,
        channel_id::UInt64,
        channel_name::AbstractString,
        count_docs_dir = joinpath("data", "docs"),
        count_names_file = joinpath(count_docs_dir, "$source.json")
    )

    isfile(count_names_file) || write(ensurepath!(count_names_file), "{}")
    namescount = JSON.parsefile(count_names_file)

    info = [
        Dict(
            "when" => """$(now(tz"UTC"))""",
            "channel" => "$channel_name",
            "channel_id" => "$channel_id"
        )
    ]

    if name in keys(namescount)
        namescount[name]["count"] += 1
        append!(
            namescount[name]["info"], info
        )
    else
        push!(
            namescount,
                name => Dict(
                    "count" => 1,
                    "info" => info
                )
        )
    end
    write(count_names_file, JSON.json(namescount, 4))
    return namescount[name]["count"]
end

function stats_namescount(
        source::AbstractString;
        place::AbstractString = "top",
        number::Int = 20,
        name::AbstractString = "",
        count_docs_dir = joinpath("data", "docs"),
        count_names_file = joinpath(count_docs_dir, "$source.json")
    )

    isfile(count_names_file) || return "Sorry, no stats found."
    namescount = JSON.parsefile(count_names_file)
    name = strip(name)
    if length(name) > 0
        if name in keys(namescount)
            return "The docstring for `$name` has been looked upon $(namescount[name]["count"]) times!"
        else
            return "The docstring for `$name` has not been looked upon yet."
        end
    end
    rev = place == "top" ? true : false
    stats = sort(collect(namescount), by=x->x[2]["count"], rev=rev)[1:min(end, number, 100)]
    header = ifelse(place == "top", "**Top ", "**Bottom ") * "$(min(length(stats), number, 100)) stats**\n"
    header *= "≡"^(length(header)-8) * "\n__Count__\t__Name__\n"
    contents = prod(map(s -> "$(s[2]["count"])\t\t\t`$(s[1])`\n", stats))
    return header * contents
end
# Convenient Discord functions

"""
    delete_emoji(c::Client, channel_id::UInt64, message_id::UInt64, user_id::UInt64, emoji::Emoji)

Delete the specified user's emoji from a Discord message.
"""
function delete_emoji(c::Client, channel_id::UInt64, message_id::UInt64, user_id::UInt64, emoji::Emoji)
    message = Message(; id = message_id, channel_id = channel_id)
    user = User(; id = user_id)
    return @discord delete(c, emoji, message, user)
end

"""
    update_message(c::Client, channel_id::UInt64, message_id::UInt64, content::AbstractString)

Update a Discord message with new content.
"""
function update_message(c::Client, channel_id::UInt64, message_id::UInt64, content::AbstractString)
    message = Message(; id = message_id, channel_id = channel_id)
    return @discord update(c, message; content)
end

"""
    is_real_user(u::User)

Returns true if the user is not a bot.
"""
is_real_user(u::User) = ismissing(u.bot) || u.bot === false

"""
    is_bot(u::User)

Returns true if the user is a bot.
"""
is_bot(u::User) = !is_real_user(u)

# Generic flipping functions

"""
    previous

Flip state to the previous item. If it passes the beginning of the list,
reset the state to the last item.
"""
function previous end

"""
    next

Flip state to the next item. If it passes the end of the list,
reset the state to the first item.
"""
function next end

"""
    flipper(f::Union{typeof(previous), typeof(next)}, T::Type)

Return a flip function, either `previous` or `next`, that expects
argument of type `T`.
"""
flipper(f::Union{typeof(previous), typeof(next)}, T::Type) = TypedFunction{T}(f)

"""
    extract_command(command::AbstractString, s::AbstractString)

Extract command by removing the command prefix. For example:

```
julia> HoJBot.extract_command("j", ",j doc sin")
"doc sin"
```
"""
function extract_command(command::AbstractString, s::AbstractString)
    prefix = COMMAND_PREFIX * command
    return strip(replace(s, Regex("^" * prefix * " *") => ""))
end

"""
    ensurepath!(fileorpath::AbstractString)

Ensure that the path exists in a way that writing to path does not error.
Returns the argument afterwards for composability.
"""
function ensurepath!(fileorpath::AbstractString)
    mkpath(dirname(fileorpath))
    return fileorpath
end

# TODO move to constants.jl later
const Optional{T} = Union{T, Nothing}


text:
#=
Copyright (c) 2016, Intel Corporation
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.
=#

function pattern_match_call_kmeans(f::GlobalRef, cluster_out::RHSVar, arr::RHSVar,
                                   num_clusters::Union{RHSVar,Int}, iterNum::Union{RHSVar,Int}, start::LHSVar, count::LHSVar,
                                   col_size::Union{RHSVar,Int,Expr}, tot_row_size::Union{RHSVar,Int,Expr}, linfo)
    s = ""
    if f.name==:Kmeans_dist
        c_arr = ParallelAccelerator.CGen.from_expr(arr, linfo)
        c_num_clusters = ParallelAccelerator.CGen.from_expr(num_clusters, linfo)
        c_iters = ParallelAccelerator.CGen.from_expr(iterNum, linfo)
        c_col_size = ParallelAccelerator.CGen.from_expr(col_size, linfo)
        c_tot_row_size = ParallelAccelerator.CGen.from_expr(tot_row_size, linfo)
        c_cluster_out = ParallelAccelerator.CGen.from_expr(cluster_out, linfo)
        c_start = ParallelAccelerator.CGen.from_expr(start, linfo)
        c_count = ParallelAccelerator.CGen.from_expr(count, linfo)
        c_threads = ParallelAccelerator.CGen.USE_OMP==1 ? "omp_get_max_threads()" : "1"

        s *= """
        services::Environment::getInstance()->setNumberOfThreads($c_threads);
        byte   *nodeCentroids;
        size_t CentroidsArchLength;
        services::SharedPtr<NumericTable> centroids;
        InputDataArchive centroidsDataArch;
        int nIterations = $c_iters;
        int mpi_root = 0;
        int rankId = __hpat_node_id;

        HomogenNumericTable<double>* dataTable = new HomogenNumericTable<double>((double*)$c_arr.getData(), $c_col_size, $c_count);
        services::SharedPtr<NumericTable> dataTablePointer(dataTable);
        kmeans::init::Distributed<step1Local,double,kmeans::init::randomDense>
                       localInit($c_num_clusters, $c_tot_row_size, $c_start);
        localInit.input.set(kmeans::init::data, dataTablePointer);

        /* Compute k-means */
        localInit.compute();

        /* Serialize partial results required by step 2 */
        services::SharedPtr<byte> serializedData;
        InputDataArchive dataArch;
        localInit.getPartialResult()->serialize( dataArch );
        size_t perNodeArchLength = dataArch.getSizeOfArchive();

        /* Serialized data is of equal size on each node if each node called compute() equal number of times */
        if (rankId == mpi_root)
        {
            serializedData = services::SharedPtr<byte>( new byte[ perNodeArchLength * __hpat_num_pes ] );
        }

        byte *nodeResults = new byte[ perNodeArchLength ];
        dataArch.copyArchiveToArray( nodeResults, perNodeArchLength );

        /* Transfer partial results to step 2 on the root node */
        MPI_Gather( nodeResults, perNodeArchLength, MPI_CHAR, serializedData.get(), perNodeArchLength, MPI_CHAR, mpi_root,
                MPI_COMM_WORLD);

        delete[] nodeResults;

        if(rankId == mpi_root)
        {
            /* Create an algorithm to compute k-means on the master node */
            kmeans::init::Distributed<step2Master, double, kmeans::init::randomDense> masterInit($c_num_clusters);

            for( size_t i = 0; i < __hpat_num_pes ; i++ )
            {
                /* Deserialize partial results from step 1 */
                OutputDataArchive dataArch( serializedData.get() + perNodeArchLength * i, perNodeArchLength );

                services::SharedPtr<kmeans::init::PartialResult> dataForStep2FromStep1 = services::SharedPtr<kmeans::init::PartialResult>(
                                                                               new kmeans::init::PartialResult() );
                dataForStep2FromStep1->deserialize(dataArch);

                /* Set local partial results as input for the master-node algorithm */
                masterInit.input.add(kmeans::init::partialResults, dataForStep2FromStep1 );
        }

        /* Merge and finalizeCompute k-means on the master node */
        masterInit.compute();
        masterInit.finalizeCompute();

        centroids = masterInit.getResult()->get(kmeans::init::centroids);
        }

        for(int iter=0; iter<nIterations; iter++)
        {

            if(rankId == mpi_root)
            {
                /*Retrieve the algorithm results and serialize them */
                centroids->serialize( centroidsDataArch );
                CentroidsArchLength = centroidsDataArch.getSizeOfArchive();
            }

             /* Get partial results from the root node */
             MPI_Bcast( &CentroidsArchLength, sizeof(size_t), MPI_CHAR, mpi_root, MPI_COMM_WORLD );

              nodeCentroids = new byte[ CentroidsArchLength ];

            if(rankId == mpi_root)
            {
                centroidsDataArch.copyArchiveToArray( nodeCentroids, CentroidsArchLength );
            }

            MPI_Bcast( nodeCentroids, CentroidsArchLength, MPI_CHAR, mpi_root, MPI_COMM_WORLD );

            /* Deserialize centroids data */
            OutputDataArchive centroidsDataArch( nodeCentroids, CentroidsArchLength );

            centroids = services::SharedPtr<NumericTable>( new HomogenNumericTable<double>() );

            centroids->deserialize(centroidsDataArch);

            /* Create an algorithm to compute k-means on local nodes */
            kmeans::Distributed<step1Local> localAlgorithm($c_num_clusters);

            /* Set the input data set to the algorithm */
            localAlgorithm.input.set(kmeans::data,           dataTablePointer);
            localAlgorithm.input.set(kmeans::inputCentroids, centroids);

            /* Compute k-means */
            localAlgorithm.compute();

            /* Serialize partial results required by step 2 */
            services::SharedPtr<byte> serializedData;
            InputDataArchive dataArch;
            localAlgorithm.getPartialResult()->serialize( dataArch );
            size_t perNodeArchLength = dataArch.getSizeOfArchive();

            /* Serialized data is of equal size on each node if each node called compute() equal number of times */
            if (rankId == mpi_root)
            {
                serializedData = services::SharedPtr<byte>( new byte[ perNodeArchLength * __hpat_num_pes ] );
            }
            byte *nodeResults = new byte[ perNodeArchLength ];
            dataArch.copyArchiveToArray( nodeResults, perNodeArchLength );

            /* Transfer partial results to step 2 on the root node */
            MPI_Gather( nodeResults, perNodeArchLength, MPI_CHAR, serializedData.get(), perNodeArchLength, MPI_CHAR, mpi_root,
                        MPI_COMM_WORLD);

            delete[] nodeResults;

            if(rankId == mpi_root)
            {
               /* Create an algorithm to compute k-means on the master node */
               kmeans::Distributed<step2Master> masterAlgorithm($c_num_clusters);

               for( size_t i = 0; i < __hpat_num_pes ; i++ )
                {
                    /* Deserialize partial results from step 1 */
                    OutputDataArchive dataArch( serializedData.get() + perNodeArchLength * i, perNodeArchLength );

                    services::SharedPtr<kmeans::PartialResult> dataForStep2FromStep1 = services::SharedPtr<kmeans::PartialResult>(
                                                                               new kmeans::PartialResult() );
                    dataForStep2FromStep1->deserialize(dataArch);

                    /* Set local partial results as input for the master-node algorithm */
                    masterAlgorithm.input.add(kmeans::partialResults, dataForStep2FromStep1 );
                }

                /* Merge and finalizeCompute k-means on the master node */
                masterAlgorithm.compute();
                masterAlgorithm.finalizeCompute();

                /* Retrieve the algorithm results */
                centroids = masterAlgorithm.getResult()->get(kmeans::centroids);
            }
            delete[] nodeCentroids;
        }

        BlockDescriptor<double> block;

        //std::cout<<centroids->getNumberOfRows()<<std::endl;
        //std::cout<<centroids->getNumberOfColumns()<<std::endl;

        centroids->getBlockOfRows(0, $c_num_clusters, readOnly, block);
        double* out_arr = block.getBlockPtr();
        //std::cout<<"output ";
        //for(int i=0; i<$c_col_size*$c_num_clusters; i++)
        //{
        //    std::cout<<" "<<out_arr[i];
        //}
        //std::cout<<std::endl;
        int64_t res_dims[] = {$c_col_size,$c_num_clusters};
        double* out_data = new double[$c_col_size*$c_num_clusters];
        memcpy(out_data, block.getBlockPtr(), $c_col_size*$c_num_clusters*sizeof(double));
        j2c_array<double> kmeans_out(out_data,2,res_dims);
        $c_cluster_out = kmeans_out;
    """

    end
    return s
end

function pattern_match_call_kmeans(f::ANY, cluster_out::ANY, arr::ANY, num_clusters::ANY, start::ANY, count::ANY, cols::ANY, rows::ANY, linfo)
    return ""
end

function pattern_match_call_linear_regression(f::GlobalRef, coeff_out::RHSVar, points::RHSVar,
                                   responses::RHSVar, start_points::LHSVar, count_points::LHSVar,
                                   col_size_points::Union{RHSVar,Int,Expr}, tot_row_size_points::Union{RHSVar,Int,Expr},
                                   start_responses::LHSVar, count_responses::LHSVar,
                                   col_size_responses::Union{RHSVar,Int,Expr}, tot_row_size_responses::Union{RHSVar,Int,Expr}, linfo)
    s = ""
    if f.name==:LinearRegression_dist
        c_points = ParallelAccelerator.CGen.from_expr(points, linfo)
        c_responses = ParallelAccelerator.CGen.from_expr(responses, linfo)
        c_col_size_points = ParallelAccelerator.CGen.from_expr(col_size_points, linfo)
        c_tot_row_size_points = ParallelAccelerator.CGen.from_expr(tot_row_size_points, linfo)
        c_col_size_responses = ParallelAccelerator.CGen.from_expr(col_size_responses, linfo)
        c_tot_row_size_responses = ParallelAccelerator.CGen.from_expr(tot_row_size_responses, linfo)
        c_coeff_out = ParallelAccelerator.CGen.from_expr(coeff_out, linfo)
        c_start_points = ParallelAccelerator.CGen.from_expr(start_points, linfo)
        c_count_points = ParallelAccelerator.CGen.from_expr(count_points, linfo)
        c_start_responses = ParallelAccelerator.CGen.from_expr(start_responses, linfo)
        c_count_responses = ParallelAccelerator.CGen.from_expr(count_responses, linfo)
        c_threads = ParallelAccelerator.CGen.USE_OMP==1 ? "omp_get_max_threads()" : "1"

        s = """
            assert($c_tot_row_size_points==$c_tot_row_size_responses);
            int mpi_root = 0;
            int rankId = __hpat_node_id;
            services::Environment::getInstance()->setNumberOfThreads($c_threads);

            HomogenNumericTable<double>* dataTable = new HomogenNumericTable<double>((double*)$c_points.getData(), $c_col_size_points, $c_count_points);
            HomogenNumericTable<double>* responseTable = new HomogenNumericTable<double>((double*)$c_responses.getData(), $c_col_size_responses, $c_count_responses);
            services::SharedPtr<NumericTable> trainData(dataTable);
            services::SharedPtr<NumericTable> trainDependentVariables(responseTable);
            services::SharedPtr<linear_regression::training::Result> trainingResult;

            /* Create an algorithm object to train the multiple linear regression model based on the local-node data */
            linear_regression::training::Distributed<step1Local, double, linear_regression::training::qrDense> localAlgorithm;

            /* Pass a training data set and dependent values to the algorithm */
            localAlgorithm.input.set(linear_regression::training::data, trainData);
            localAlgorithm.input.set(linear_regression::training::dependentVariables, trainDependentVariables);

            /* Train the multiple linear regression model on local nodes */
            localAlgorithm.compute();

            /* Serialize partial results required by step 2 */
            services::SharedPtr<byte> serializedData;
            InputDataArchive dataArch;
            localAlgorithm.getPartialResult()->serialize( dataArch );
            size_t perNodeArchLength = dataArch.getSizeOfArchive();

            /* Serialized data is of equal size on each node if each node called compute() equal number of times */
            if (rankId == mpi_root)
            {
                serializedData = services::SharedPtr<byte>( new byte[ perNodeArchLength * __hpat_num_pes] );
            }

            byte *nodeResults = new byte[ perNodeArchLength ];
            dataArch.copyArchiveToArray( nodeResults, perNodeArchLength );

            /* Transfer partial results to step 2 on the root node */
            MPI_Gather( nodeResults, perNodeArchLength, MPI_CHAR, serializedData.get(), perNodeArchLength, MPI_CHAR, mpi_root,
                        MPI_COMM_WORLD);

            delete[] nodeResults;
            services::SharedPtr<NumericTable> trainingCoeffsTable;
            if(rankId == mpi_root)
            {
                /* Create an algorithm object to build the final multiple linear regression model on the master node */
                linear_regression::training::Distributed<step2Master, double, linear_regression::training::qrDense> masterAlgorithm;

                for( size_t i = 0; i < __hpat_num_pes; i++ )
                {
                    /* Deserialize partial results from step 1 */
                    OutputDataArchive dataArch( serializedData.get() + perNodeArchLength * i, perNodeArchLength );

                    services::SharedPtr<linear_regression::training::PartialResult> dataForStep2FromStep1 = services::SharedPtr<linear_regression::training::PartialResult>
                                                                               ( new linear_regression::training::PartialResult() );
                    dataForStep2FromStep1->deserialize(dataArch);

                    /* Set the local multiple linear regression model as input for the master-node algorithm */
                    masterAlgorithm.input.add(linear_regression::training::partialModels, dataForStep2FromStep1);
                }

                /* Merge and finalizeCompute the multiple linear regression model on the master node */
                masterAlgorithm.compute();
                masterAlgorithm.finalizeCompute();

                /* Retrieve the algorithm results */
                trainingResult = masterAlgorithm.getResult();
                // printNumericTable(trainingResult->get(linear_regression::training::model)->getBeta(), "Linear Regression coefficients:");
                trainingCoeffsTable = trainingResult->get(linear_regression::training::model)->getBeta();

                BlockDescriptor<double> block;

            //std::cout<<trainingCoeffsTable->getNumberOfRows()<<std::endl;
            //std::cout<<trainingCoeffsTable->getNumberOfColumns()<<std::endl;

                trainingCoeffsTable->getBlockOfRows(0, $c_col_size_responses, readOnly, block);
                double* out_arr = block.getBlockPtr();

            // assuming intercept is required
            int64_t coeff_size = $c_col_size_points+1;
            //std::cout<<"output ";
            //for(int i=0; i<coeff_size*$c_col_size_responses; i++)
            //{
            //    std::cout<<" "<<out_arr[i];
            //}
            //std::cout<<std::endl;

                int64_t res_dims[] = {coeff_size,$c_col_size_responses};
                double* out_data = new double[coeff_size*$c_col_size_responses];
                memcpy(out_data, block.getBlockPtr(), coeff_size*$c_col_size_responses*sizeof(double));
                j2c_array<double> linear_regression_out(out_data,2,res_dims);
                $c_coeff_out = linear_regression_out;
            }
            """
    end
    return s
end

function pattern_match_call_linear_regression(f::ANY, coeff_out::ANY, arr::ANY, num_clusters::ANY,
          start::ANY, count::ANY, cols::ANY, rows::ANY, start2::ANY, count2::ANY, cols2::ANY, rows2::ANY, linfo)
    return ""
end

function pattern_match_call_naive_bayes(f::GlobalRef, coeff_out::RHSVar, points::RHSVar,
                                   labels::RHSVar, num_classes::Union{RHSVar,Int,Expr}, start_points::LHSVar, count_points::LHSVar,
                                   col_size_points::Union{RHSVar,Int,Expr}, tot_row_size_points::Union{RHSVar,Int,Expr},
                                   start_labels::LHSVar, count_labels::LHSVar,
                                   col_size_labels::Union{RHSVar,Int,Expr}, tot_row_size_labels::Union{RHSVar,Int,Expr}, linfo)
    s = ""
    if f.name==:NaiveBayes_dist
        c_points = ParallelAccelerator.CGen.from_expr(points, linfo)
        c_labels = ParallelAccelerator.CGen.from_expr(labels, linfo)
        c_col_size_points = ParallelAccelerator.CGen.from_expr(col_size_points, linfo)
        c_tot_row_size_points = ParallelAccelerator.CGen.from_expr(tot_row_size_points, linfo)
        c_col_size_labels = ParallelAccelerator.CGen.from_expr(col_size_labels, linfo)
        c_tot_row_size_labels = ParallelAccelerator.CGen.from_expr(tot_row_size_labels, linfo)
        c_coeff_out = ParallelAccelerator.CGen.from_expr(coeff_out, linfo)
        c_num_classes = ParallelAccelerator.CGen.from_expr(num_classes, linfo)
        c_start_points = ParallelAccelerator.CGen.from_expr(start_points, linfo)
        c_count_points = ParallelAccelerator.CGen.from_expr(count_points, linfo)
        c_start_labels = ParallelAccelerator.CGen.from_expr(start_labels, linfo)
        c_count_labels = ParallelAccelerator.CGen.from_expr(count_labels, linfo)
        c_threads = ParallelAccelerator.CGen.USE_OMP==1 ? "omp_get_max_threads()" : "1"

        s = """
            assert($c_tot_row_size_points==$c_tot_row_size_labels);
            int mpi_root = 0;
            int rankId = __hpat_node_id;
            services::Environment::getInstance()->setNumberOfThreads($c_threads);

            HomogenNumericTable<double>* dataTable = new HomogenNumericTable<double>((double*)$c_points.getData(), $c_col_size_points, $c_count_points);
            HomogenNumericTable<double>* responseTable = new HomogenNumericTable<double>((double*)$c_labels.getData(), $c_col_size_labels, $c_count_labels);
            services::SharedPtr<NumericTable> trainData(dataTable);
            services::SharedPtr<NumericTable> trainGroundTruth(responseTable);
            services::SharedPtr<multinomial_naive_bayes::training::Result> trainingResult;

            /* Create an algorithm object to train the Na__ve Bayes model based on the local-node data */
            multinomial_naive_bayes::training::Distributed<step1Local> localAlgorithm($c_num_classes);

            /* Pass a training data set and dependent values to the algorithm */
            localAlgorithm.input.set(classifier::training::data,   trainData);
            localAlgorithm.input.set(classifier::training::labels, trainGroundTruth);

            /* Train the Na__ve Bayes model on local nodes */
            localAlgorithm.compute();

            /* Serialize partial results required by step 2 */
            services::SharedPtr<byte> serializedData;
            InputDataArchive dataArch;
            localAlgorithm.getPartialResult()->serialize(dataArch);
            size_t perNodeArchLength = dataArch.getSizeOfArchive();

            /* Serialized data is of equal size on each node if each node called compute() equal number of times */
            if (rankId == mpi_root)
            {
                serializedData = services::SharedPtr<byte>(new byte[perNodeArchLength * __hpat_num_pes]);
            }

            byte *nodeResults = new byte[perNodeArchLength];
            dataArch.copyArchiveToArray( nodeResults, perNodeArchLength );

            /* Transfer partial results to step 2 on the root node */
            MPI_Gather( nodeResults, perNodeArchLength, MPI_CHAR, serializedData.get(), perNodeArchLength, MPI_CHAR, mpi_root,
                        MPI_COMM_WORLD);

            delete[] nodeResults;

            if(rankId == mpi_root)
            {
                /* Create an algorithm object to build the final Na__ve Bayes model on the master node */
                multinomial_naive_bayes::training::Distributed<step2Master> masterAlgorithm($c_num_classes);

                for(size_t i = 0; i < __hpat_num_pes ; i++)
                {
                    /* Deserialize partial results from step 1 */
                    OutputDataArchive dataArch(serializedData.get() + perNodeArchLength * i, perNodeArchLength);

                    services::SharedPtr<classifier::training::PartialResult> dataForStep2FromStep1 = services::SharedPtr<classifier::training::PartialResult>
                                                                               (new classifier::training::PartialResult());
                    dataForStep2FromStep1->deserialize(dataArch);

                    /* Set the local Na__ve Bayes model as input for the master-node algorithm */
                    masterAlgorithm.input.add(classifier::training::partialModels, dataForStep2FromStep1);
                }

                /* Merge and finalizeCompute the Na__ve Bayes model on the master node */
                masterAlgorithm.compute();
                masterAlgorithm.finalizeCompute();

                /* Retrieve the algorithm results */
                trainingResult = masterAlgorithm.getResult();

                services::SharedPtr<NumericTable> trainingLogpTable = trainingResult->get(classifier::training::model)->getLogP();
                services::SharedPtr<NumericTable> trainingLogThetaTable = trainingResult->get(classifier::training::model)->getLogTheta();

                BlockDescriptor<double> block1, block2;

                //std::cout<<trainingLogThetaTable->getNumberOfRows()<<std::endl;
                //std::cout<<trainingLogThetaTable->getNumberOfColumns()<<std::endl;
                //std::cout<<trainingLogpTable->getNumberOfRows()<<std::endl;
                //std::cout<<trainingLogpTable->getNumberOfColumns()<<std::endl;

                trainingLogpTable->getBlockOfRows(0, $c_num_classes, readOnly, block1);
                trainingLogThetaTable->getBlockOfRows(0, $c_num_classes, readOnly, block2);
                double* out_arr1 = block1.getBlockPtr();
                double* out_arr2 = block2.getBlockPtr();

            //std::cout<<"output ";
            //for(int i=0; i<coeff_size*$c_col_size_labels; i++)
            //{
            //    std::cout<<" "<<out_arr[i];
            //}
            //std::cout<<std::endl;

                int64_t res_dims[] = {$c_num_classes,$c_num_classes+1};
                double* out_data = new double[$c_num_classes*($c_num_classes+1)];
                memcpy(out_data, block1.getBlockPtr(), $c_num_classes*sizeof(double));
                memcpy(out_data+$c_num_classes, block2.getBlockPtr(), $c_num_classes*$c_num_classes*sizeof(double));
                j2c_array<double> naive_bayes_out(out_data,2,res_dims);
                $c_coeff_out = naive_bayes_out;
            }
            """
    end
    return s
end

function pattern_match_call_naive_bayes(f::ANY, coeff_out::ANY, arr::ANY, arr2::ANY, numClass::Any,
          start::ANY, count::ANY, cols::ANY, rows::ANY, start2::ANY, count2::ANY, cols2::ANY, rows2::ANY, linfo)
    return ""
end


text:
function jump_mimo(solver, seed, n; verbose = false, test = false)

    # n = 3
    m = 10n
    s, H, y, L = mimo_data(seed, m, n)

    nvars = ProxSDP.sympackedlen(n + 1)

    model = Model(with_optimizer(solver))
    @variable(model, X[1:n+1, 1:n+1], PSD)
    for j in 1:(n+1), i in j:(n+1)
        @constraint(model, X[i, j] <=  1.0)
        @constraint(model, X[i, j] >= -1.0)
    end
    @objective(model, Min, sum(L[i, j] * X[i, j] for j in 1:n+1, i in 1:n+1))
    @constraint(model, ctr[i in 1:n+1], X[i, i] == 1.0)

    teste = @time optimize!(model)

    XX = value.(X)

    if test
        for j in 1:n+1, i in 1:n+1
            @test 1.01 > abs(XX[i,j]) > 0.99
        end
    end

    verbose && mimo_eval(s,H,y,L,XX)

    objval = objective_value(model)
    stime = MOI.get(model, MOI.SolveTime())

    rank = -1
    try
        @show rank = model.moi_backend.optimizer.model.optimizer.sol.final_rank
    catch
    end
    status = 0
    if JuMP.termination_status(model) == MOI.OPTIMAL
        status = 1
    end

    # SDP constraints
    max_spd_violation = minimum(eigen(XX).values)

    # test violations of linear constraints
    max_lin_viol = 0.0
    for j in 1:(n+1), i in j:(n+1)
        val = abs(XX[i,j]) - 1.0
        if val > 0.0
            if val > max_lin_viol
                max_lin_viol = val
            end
        end
    end
    for j in 1:(n+1)
        val = abs(XX[j,j] - 1.0)
        if val > 0.0
            if val > max_lin_viol
                max_lin_viol = val
            end
        end
    end

    return (objval, stime, rank, status, max_lin_viol, max_spd_violation)
end


text:
const res_data_key = "data";

function load_data(tree::ResDirTree)
    FileIO.load(datafile(tree))[res_data_key];
end

function save_data(tree::ResDirTree, data)
    FileIO.save(datafile(tree), res_data_key, data);
end

max_stars_repo_name: josePereiro/MyResultsManager

max_issues_repo_name: josePereiro/MyResultsManager

max_forks_repo_name: josePereiro/MyResultsManager


text:
# The sensor_status_t has a field named `type`, which is a reserved
# keyword in Julia v0.7 and below. To get around that, we can just
# construct the expression directly instead of relying on the parser.
eval(Expr(:struct, true,
    Expr(:(<:), :sensor_status_t, :LCMType),
    Expr(:block,
        Expr(:(::), :utime, :Int64),
        Expr(:(::), :sensor_name, :String),
        Expr(:(::), :rate, :Float64),
        Expr(:(::), :type, :Int16),
)))

@lcmtypesetup(sensor_status_t)