text:
#unum-properties
#functions that assess properties of unum values.

#generally speaking, the "unum library form" of a test will have the form
#is_XXXXX, for some of these, they are equivalent to a julia form from floating
#point tests that are overloaded - e.g. isnan == is_nan.  A separate file
#unum-teoe-func.jl provides aliasing implementations for all of the function
#in appendix a of "The End of Error"

@universal decode_exp(x::Unum) = decode_exp(x.esize, x.exponent)

doc"""
  `Unums.@signof(x)` extracts the 16-bit unsigned "sign" component of the flag.
"""
macro signof(x)
  :($x.flags & UNUM_SIGN_MASK)
end

doc"""`Unums.@ubitof(x)` extracts the 16-bit unsigned "ubit" component of the flag."""
macro ubitof(x)
  :($x.flags & UNUM_UBIT_MASK)
end

doc"""`Uunums.@flipsign(x)` flips the sign of the the unum x"""
macro flipsign(x)
  :($x.flags $= UNUM_SIGN_MASK)
end

@universal uflag!(x::Unum, f::UInt16) = (x.flags = f; x)

#some really dumb ones, but we'll put these in for legibility.
@universal is_ulp(x::Unum)      = (@ubitof(x) != 0)
@universal is_exact(x::Unum)    = (@ubitof(x) == 0)
@universal is_negative(x::Unum) = (@signof(x) != 0)
@universal is_positive(x::Unum) = (@signof(x) == 0)
@universal is_neg_def(x::Unum)  = (!is_zero(x)) && is_negative(x)
@universal is_pos_def(x::Unum)  = (!is_zero(x)) && is_positive(x)
export is_ulp, is_exact, is_negative, is_positive

#a couple of testing conditions
@universal function __is_nan_or_inf(x::Unum)
  (x.fsize == max_fsize(FSS)) || return false
  (x.esize == max_esize(ESS)) || return false
  (x.exponent == max_biased_exponent(ESS)) || return false

  if FSS < 7
    x.fraction == mask_top(max_fsize(FSS))
  else
    is_all_ones(x.fraction)
  end
end

@universal is_nan(x::Unum) = is_ulp(x) && __is_nan_or_inf(x)
@universal Base.isnan(x::Unum) = is_ulp(x) && __is_nan_or_inf(x)
export is_nan

#isinf matches the julia definiton and triggers on either positive or negative
#infinity.  is_pos_inf and is_neg_inf both are Unum-specific functions that detect
#the expected values.
@universal is_inf(x::Unum) = is_exact(x) && __is_nan_or_inf(x)
@universal is_pos_inf(x::Unum) = is_positive(x) && is_exact(x) && __is_nan_or_inf(x)
@universal is_neg_inf(x::Unum) = is_negative(x) && is_exact(x) && __is_nan_or_inf(x)
export is_inf, is_pos_inf, is_neg_inf
#aliasing to base definition
@universal Base.isinf(x::Unum) = is_exact(x) && __is_nan_or_inf(x)

@universal function is_finite(x::Unum)
  #record the maximum esize and fsize values.  Any value less than this and
  #it's finite.
  x.esize < max_esize(ESS) && return true
  x.fsize < max_fsize(FSS) && return true
  x.exponent < max_biased_exponent(ESS) && return true

  if (FSS < 7)
    x.fraction < mask_top(max_fsize(FSS))
  else
    is_not_ones(x.fraction)
  end
end
@universal Base.isfinite(x::Unum) = is_finite(x)

#NB:  The difference between "is_exp_zero" and "issubnormal" - is_exp_zero admits
#zero as a solution; issubnormal is in compliance with the standard julia
#issubnormal function and does not admit zero as a true result.
@universal is_subnormal(x::Unum) = (x.exponent == z64) && (is_ulp(x) || is_not_zero(x.fraction))
@universal is_exp_zero(x::Unum) = x.exponent == z64

doc"""
  Unums.is_strange_subnormal(x) returns if the unum value x is a "strange subnormal".
  A "strange subnormal" is a subnormal where the esize is not maxed out.  These
  subnormal values can be "corrected" to normal rerpresentations.
"""
@universal function is_strange_subnormal(x::Unum)
  (x.esize < max_esize(ESS)) && is_subnormal(x)
end

@universal Base.issubnormal(x::Unum) = (x.exponent == z64) && is_not_zero(x.fraction) #alias the unum-form to the julia-compliant form.
@universal is_frac_zero(x::Unum) = is_all_zero(x.fraction)
@universal is_zero(x::Unum) = (x.exponent == z64) && is_exact(x) && is_frac_zero(x)

doc"""
  `is_unit(x::Unum{ESS,FSS})` tests if the value in a unum is +/- 1.  Because of
  asymmetric exponents, this is slightly more laborious than you might expect:
  all unums ESS != 0, there are two canonical representations for the value
  1, one being exponent 0; fraction 0 and the other being subnormal exponent,
  fraction 0b1000...0000.  When ESS == 0 there is only access to the subnormal
  form.
"""
@universal function is_unit(x::Unum)
  #asymmetric exponents make this slightly more laborious than might be expected
  #one is not an ulp, it is exact.
  is_ulp(x) && return false
  #case one:  An exponent of zero and nothing in the fraction.
  #note that when x.esize == 0, then decode_exp 0 is subnormal.
  (x.esize != 0) && (decode_exp(x) == 0) && is_all_zero(x.fraction) && return true
  is_subnormal(x) && (x.esize == 0) && (is_top(x.fraction)) && return true
  return false
end

@universal is_one(x::Unum) = is_positive(x) && is_unit(x)
@universal is_neg_one(x::Unum) = is_negative(x) && is_unit(x)
#checks if the value is sss ("smaller than small subnormal")
@universal function is_sss(x::Unum)
  is_ulp(x) && (x.esize == max_esize(ESS)) && (x.exponent == z64) && (x.fsize == max_fsize(FSS)) && is_all_zero(x.fraction)
end
@universal is_pos_sss(x::Unum) = is_positive(x) && is_sss(x)
@universal is_neg_sss(x::Unum) = is_negative(x) && is_sss(x)

################################################################################
##  MMR CHECKS
doc"""
  `is_mmr_frac(::Unum)` overloads the is_mmr_frac(::ArrayNum) procedure to allow
  for transparent checking of whether or not a Unum has the look of MMR, whether
  or not the Unum has a big or small structure
"""
function is_mmr_frac{ESS,FSS}(x::UnumSmall{ESS,FSS})
  is_mmr_frac(x.fraction, Val{FSS})
end
function is_mmr_frac{ESS,FSS}(x::UnumLarge{ESS,FSS})
  is_mmr_frac(x.fraction)
end

doc"""
  `is_mmr(::Unum)` sign-agnostically checks to see if the passed unum is the positive
  'more than maxreal' interval, or the open bound ±(maxreal, ∞)
"""
@universal function is_mmr(x::Unum)
  is_ulp(x) || return false
  x.esize == max_esize(ESS) || return false
  x.fsize == max_fsize(FSS) || return false
  x.exponent == max_biased_exponent(ESS)|| return false
  is_mmr_frac(x)
end

doc""" `is_pos_mmr(::Unum)` checks to see if the passed unum is the positive 'more than maxreal' interval, (maxreal, ∞)"""
@universal is_pos_mmr(x::Unum) = is_positive(x) && is_mmr(x)

doc""" `is_neg_mmr(::Unum)` checks to see if the passed unum is the negative 'more than maxreal' interval, (-∞, -maxreal)"""
@universal is_neg_mmr(x::Unum) = is_negative(x) && is_mmr(x)

export is_subnormal, is_exp_zero
export is_frac_zero, is_zero, is_sss, is_pos_sss, is_neg_sss
export is_mmr, is_pos_mmr, is_neg_mmr

doc"""
  `Unums.is_magnitude_less_than_one(::Unum)` checks to see if the magnitude (absolute
  value) is less than one.  This is important for multiplication, in particular
  checking for smaller than smallest subnormal and more than maxreal
  calculations.  As is evident by this awkward syntax, you're really not supposed
  to use this outside of the Unums namespace.
"""
@universal function is_magnitude_less_than_one(x::Unum)
  #checking subnormal status.
  if (is_exp_zero(x))
    #check to see if the top bit is zero
    return (x.esize > 0) | is_top_frac_bit_zero(x)
  else
    #check to see if the top bit of the exponent is zero.
    return (decode_exp(x) < 0)
  end
end

################################################################################
# detecting if something is a zero ulp
@universal function is_zero_ulp(x::Unum)
  is_ulp(x) || return false
  (x.exponent == z64) || return false
  is_all_zero(x.fraction)
end

################################################################################
# detecting if something is an inf ulp
@universal function is_inf_ulp(x::Unum)
  is_ulp(x) || return false
  (decode_exp(x) == max_exponent(ESS)) || return false

  if x.fsize == max_fsize(FSS)
    is_mmr_frac(x)
  else
    is_all_ones(x.fraction, x.fsize)
  end
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using LinearAlgebra

function get_module_3_excersise_4_author_name()
    return "YOUR NAME HERE"
end

function problem_4(a)

    A = [
    a 1 0 0 0 1
    1 a 1 0 0 0
    0 1 a 1 0 0
    0 0 1 a 1 0
    0 0 0 1 a 1
    1 0 0 0 1 a
    ]

    eigen_vals = []
    eigen_vectors = []
    return eigen_vals, eigen_vectors
end

ext: jl


text:
"""
Copyright Marcus Greiff 2020

Tests for the matrix_solve_posdef() function
"""

# Include modules and recompile if running test standalone
if !isdefined(Main, :runAllTests)
    using Test, AerialVehicleControl, LinearAlgebra
    recompile();
end

# Testing fixture
function fixture(m,n)
    A = rand(m,m)
    A = A + A'
    A = A + (abs(minimum(eigvals(A))) + 0.01) * Matrix{Float64}(I, m, m)
    Ac = deepcopy(A)
    Am = matrix_double_t(A)
    B = rand(m,n)
    Bc = deepcopy(B)
    Bm = matrix_double_t(B)
    return A, Ac, Am, B, Bc, Bm
end

# Compact call for the matrix_double_addition function
function test_call_matrix_double_solve_posdef(Am, Bm)
    status = ccall((:matrix_double_solve_posdef, AerialVehicleControl.CONT_LIB_PATH),
        Cint,
        (Ref{matrix_double_t},Ref{matrix_double_t}),
        Am,
        Bm)
    return status
end

@verbose(1, "Testing matrix PSD solver...")
@testset "Matrix solve PD linear systems" begin
    @testset "Good inputs" begin
        for m = 1:10
            for n = 1:10
                # Fixture
                A, Ac, Am, B, Bc, Bm = fixture(m, n)
                # Solve system in Julia
                X = A\B;
                # Call
                status = test_call_matrix_double_solve_posdef(Am, Bm)
                # Test
                @test isapprox(status, 1)
                @test !isapprox(A, Ac)
                @test !isapprox(B, Bc)
                @test isapprox(X,B)
                @test isapprox(Ac*B,Bc)
            end
        end
    end
end


text:
#studenti.jl
using Printf

struct Student
    ime::String
    prezime::String
    brojIndeksa::Int
end

studenti = []
push!(studenti, Student("Marko", "Markovic", 13))
push!(studenti, Student("Nikola", "Tomic", 8))
push!(studenti, Student("Srdjan", "Knezevic", 5))
push!(studenti, Student("Vanja", "Kovacevic", 17))
push!(studenti, Student("Nenad", "Nikolic", 16))
push!(studenti, Student("Ivana", "Ivanovic", 2))
push!(studenti, Student("Branislav", "Tomic", 9))
push!(studenti, Student("Stefan", "Zivkovic", 3))

function ispisStudenata(studenti)
  for i in 1:length(studenti)
     @printf("%2i. %-10s %-10s %3i\n", i, studenti[i].ime,
              studenti[i].prezime, studenti[i].brojIndeksa)
  end
end


text:
juliacore_typeof(line::Ref{UInt}, ns::Namespace, arg::AetherlangObject)::AetherlangObject = AetherlangObject(string(typeof(arg.dataref)))
juliacore_eval(line::Ref{UInt}, ns::Namespace, exp::AetherlangObject)::AetherlangObject = AetherlangObject(eval(Meta.parse(exp.dataref)))

juliacore_namespace_modify = Namespace(
    "typeof" => AetherlangObject(juliacore_typeof),
    "eval" => AetherlangObject(juliacore_eval)
)

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
##### DegeneracyHunter.jl
# Created by Alexander Dowling

module DegeneracyHunter

export IrreducibleDegenerateSet, DegenSettings, DegenData, printBound, printVariablesInEquation,
	printIDS, degeneracyHunter, printInfeasibleEquations, printInactiveEquations,
	checkVarBounds, printRows, ProblemStats, assembleProblemStats, stringProblemStats, printProblemStats,
	printVariableDiagnostics, round_to_bounds!, checkEquationScaling

importall JuMP
import MathProgBase

include("types.jl")
include("model_diagnostics.jl")
include("classic_degeneracy_hunter.jl")
include("problem_stats.jl")
include("auxillary.jl")

end # module


text:
using ..Utils: LDS, CLCDS

using LazySets: CachedMinkowskiSumArray,
                isdisjoint

import LazySets.Approximations: overapproximate

"""
    reach(problem, options)

Interface to reachability algorithms for an LTI system.

### Input

- `problem`   -- initial value problem
- `options`   -- additional options

### Output

A sequence of [`SparseReachSet`](@ref)s.

### Notes

A dictionary with available algorithms is available via
`Reachability.available_algorithms`.

The numeric type of the system's coefficients and the set of initial states
is inferred from the first parameter of the system (resp. lazy set), ie.
`NUM = first(typeof(problem.s).parameters)`.
"""
function reach(problem::Union{IVP{<:CLDS{NUM}, <:LazySet{NUM}},
                              IVP{<:CLCDS{NUM}, <:LazySet{NUM}}},
               options::TwoLayerOptions
              ) where {NUM <: Real}
    # optional matrix conversion
    problem = matrix_conversion(problem, options)

    # list containing the arguments passed to any reachability function
    args = []

    # coefficients matrix
    A = problem.s.A
    push!(args, A)

    # determine analysis mode (sparse/dense) for lazy_expm mode
    if A isa SparseMatrixExp
        push!(args, Val(options[:assume_sparse]))
    end

    n = statedim(problem)
    blocks = options[:blocks]
    partition = convert_partition(options[:partition])
    T = options[:T]
    N = (T == Inf) ? nothing : ceil(Int, T / options[:δ])

    # Cartesian decomposition of the initial set
    if length(partition) == 1 && length(partition[1]) == n &&
            options[:block_options_init] == LinearMap
        info("- Skipping decomposition of X0")
        Xhat0 = LazySet{NUM}[problem.x0]
    else
        info("- Decomposing X0")
        @timing begin
            Xhat0 = array(decompose(problem.x0, options[:partition],
                                    options[:block_options_init]))
        end
    end

    # compute dimensions
    dimensions = compute_dimensions(partition, blocks)

    # determine output function: linear map with the given matrix
    output_function = options[:output_function] != nothing ?
        (x -> options[:output_function] * x) :
        nothing

    # preallocate output vector
    if output_function == nothing
        res_type = SparseReachSet{CartesianProductArray{NUM, LazySet{NUM}}}
    else
        # by default, this algorithm uses box overapproximation
        res_type = SparseReachSet{Hyperrectangle{NUM, Vector{NUM}, Vector{NUM}}}
    end
    res = (N == nothing) ? Vector{res_type}() : Vector{res_type}(undef, N)

    # shortcut if only the initial set is required
    if N == 1
        res[1] = SparseReachSet(CartesianProductArray(Xhat0[blocks]), zero(NUM), options[:δ], dimensions)
        return Flowpipe(res)
    end
    push!(args, Xhat0)

    # inputs
    if !options[:assume_homogeneous] && inputdim(problem) > 0
        U = inputset(problem)
    else
        U = nothing
    end
    push!(args, U)

    # overapproximation function for states
    block_options_iter = options[:block_options_iter]
    if block_options_iter isa AbstractVector ||
            block_options_iter isa Dict{Int, Any}
        # individual overapproximation options per block
        overapproximate_fun = (i, X) -> overapproximate(X, block_options_iter[i])
    else
        # uniform overapproximation options for each block
        overapproximate_fun = (i, X) -> overapproximate(X, block_options_iter)
    end
    push!(args, overapproximate_fun)

    # overapproximate function for inputs
    lazy_inputs_interval = options[:lazy_inputs_interval]
    if lazy_inputs_interval == lazy_inputs_interval_always
        overapproximate_inputs_fun = (k, i, x) -> overapproximate_fun(i, x)
    else
        # first set in a series
        function _f(k, i, x::LinearMap{NUM}) where {NUM}
            @assert k == 1 "a LinearMap is only expected in the first iteration"
            return CachedMinkowskiSumArray(LazySet{NUM}[x])
        end
        # further sets of the series
        function _f(k, i, x::MinkowskiSum{NUM, <:CachedMinkowskiSumArray}) where NUM
            if has_constant_directions(block_options_iter, i)
                # forget sets if we do not use epsilon-close approximation
                forget_sets!(x.X)
            end
            push!(array(x.X), x.Y)
            if lazy_inputs_interval(k)
                # overapproximate lazy set
                y = overapproximate_fun(i, x.X)
                return CachedMinkowskiSumArray(LazySet{NUM}[y])
            end
            return x.X
        end
        function _f(k, i, x)
            # other set types
            if lazy_inputs_interval(k)
                # overapproximate lazy set
                return overapproximate_fun(i, x.X)
            end
            return x
        end
        overapproximate_inputs_fun = _f
    end
    push!(args, overapproximate_inputs_fun)

    # ambient dimension
    push!(args, n)

    # number of computed sets
    push!(args, N)

    # output function
    push!(args, output_function)

    # add mode-specific block(s) argument
    push!(args, blocks)
    push!(args, partition)
    push!(args, dimensions)

    # time step
    push!(args, options[:δ])

    # termination function
    invariant = stateset(problem.s)
    termination = get_termination_function(N, invariant)
    push!(args, termination)

    info("- Computing successors")

    # progress meter
    progress_meter = (N != nothing) ?
        Progress(N, 1, "Computing successors ") :
        nothing
    push!(args, progress_meter)

    # choose algorithm backend
    algorithm = options[:algorithm]
    if algorithm == "explicit"
        algorithm_backend = "explicit_blocks"
    elseif algorithm == "wrap"
        algorithm_backend = "wrap"
    else
        error("Unsupported algorithm: ", algorithm)
    end
    push!(args, res)

    # call the adequate function with the given arguments list
    @timing begin
        index, skip = available_algorithms[algorithm_backend]["func"](args...)
    end

    # shrink result array
    if skip || index < N
        if N != nothing
            info("termination after only $index of $N steps")
        end
        deleteat!(res, (skip ? index : index + 1):length(res))
    end

    # return the result
    return Flowpipe(res)
end

function reach(problem::Union{IVP{<:CLCS{NUM}, <:LazySet{NUM}},
                              IVP{<:CLCCS{NUM}, <:LazySet{NUM}}},
               options::TwoLayerOptions
              ) where {NUM <: Real}
    info("Time discretization...")
    Δ = @timing discretize(problem, options[:δ],
        algorithm=options[:discretization], exp_method=options[:exp_method],
        sih_method=options[:sih_method])
    return reach(Δ, options)
end

function get_termination_function(N::Nothing, invariant::Universe)
    termination_function = (k, set, t0) -> termination_unconditioned()
    warn("no termination condition specified; the reachability analysis will " *
         "not terminate")
    return termination_function
end

function get_termination_function(N::Int, invariant::Universe)
    return (k, set, t0) -> termination_N(N, k, t0)
end

function get_termination_function(N::Nothing, invariant::LazySet)
    return (k, set, t0) -> termination_inv(invariant, set, t0)
end

function get_termination_function(N::Int, invariant::LazySet)
    return (k, set, t0) -> termination_inv_N(N, invariant, k, set, t0)
end

function termination_unconditioned()
    return (false, false)
end

function termination_N(N, k, t0)
    return (k >= N, false)
end

function termination_inv(inv, set, t0)
    if isdisjoint(set, inv)
        return (true, true)
    else
        return (false, false)
    end
end

function termination_inv_N(N, inv, k, set, t0)
    if k >= N
        return (true, false)
    elseif isdisjoint(set, inv)
        return (true, true)
    else
        return (false, false)
    end
end

function overapproximate(X::LazySet, pair::Pair)
    return overapproximate(X, pair[1], pair[2])
end

function overapproximate(X::LazySet, ::Nothing)
    return X
end

function overapproximate(msa::MinkowskiSumArray, ::Nothing)
    return MinkowskiSumArray(copy(array(msa)))
end

function has_constant_directions(block_options::AbstractVector, i::Int)
    return has_constant_directions(block_options[i], i)
end

function has_constant_directions(block_options::Dict{<:UnionAll, <:Real},
                                 i::Int)
    return has_constant_directions(block_options[i], i)
end

function has_constant_directions(block_options::Pair{<:UnionAll, <:Real},
                                 i::Int)
    return has_constant_directions(block_options[2], i)
end

function has_constant_directions(block_options::Real, i::Int)
    return ε == Inf
end

function has_constant_directions(block_options, i::Int)
    return true
end


text:
module FESolversTests

using Test
using Gridap.Arrays
using Gridap.Algebra
using Gridap.TensorValues
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Integration
using Gridap.Fields
using Gridap.FESpaces
using Gridap.CellData

domain =(0,1,0,1,0,1)
partition = (3,3,3)
model = CartesianDiscreteModel(domain,partition)
trian = get_triangulation(model)
degree = 4
quad = CellQuadrature(trian,degree)

order = 2
grid_topology = get_grid_topology(model)
polytopes = get_polytopes(grid_topology)
reffes = [LagrangianRefFE(Float64,p,order) for p in polytopes]

dirichlet_tags = [1,10]
V = GradConformingFESpace(reffes,model,dirichlet_tags)

U = TrialFESpace(V)

f(x) = x[2]

v = get_cell_basis(V)
u = get_cell_basis(U)

cellmat = integrate(∇(v)⊙∇(u),trian,quad)
cellvec = integrate(v⊙f,trian,quad)
cellids = collect(1:num_cells(trian))

assem = SparseMatrixAssembler(U,V)
matdata = ([cellmat],[cellids],[cellids])
vecdata = ([cellvec],[cellids])
A =  assemble_matrix(assem,matdata)
b =  assemble_vector(assem,vecdata)
x = A \ b
x0 = zeros(length(x))

op = AffineFEOperator(U,V,A,b)
solver = LinearFESolver()
test_fe_solver(solver,op,x0,x)
uh = solve(solver,op)
@test get_free_values(uh) ≈ x
uh = solve(op)
@test get_free_values(uh) ≈ x

solver = NonlinearFESolver()
test_fe_solver(solver,op,x0,x)
uh = solve(solver,op)
@test get_free_values(uh) ≈ x

end # module


text:
# This file includes testset for ODESystem

@testset "ODESystemTestSet" begin
    @info "Running ODESystemTestSet ..."

    # ODESystem construction
    sfunc1(dx, x, u, t) = (dx .= -x)
    ofunc1(x, u, t) = x
    ds = ODESystem(righthandside = sfunc1, readout=ofunc1, state=[1.], solverkwargs=(dt=0.1,), input=nothing, output=Outport())
    ds = ODESystem(righthandside = sfunc1, readout=ofunc1, state=[1.], solverkwargs=(dt=0.1, dense=true), input=nothing, output=Outport())
    ds = ODESystem(righthandside = sfunc1, readout=ofunc1, state=[1.], input=nothing, output=Outport())
    @test typeof(ds.trigger) == Inpin{Float64}
    @test typeof(ds.handshake) == Outpin{Bool}
    @test ds.input === nothing
    @test typeof(ds.output) == Outport{Outpin{Float64}}
    @test length(ds.output) == 1
    @test ds.state == [1.]
    @test ds.t == 0.
    @test ds.integrator.sol.prob.p === nothing

    function sfunc2(dx, x, u, t)
        dx[1] = x[1] + u[1](t)
        dx[2] = x[2] - u[2](t)
        dx[3] = x[3] + sin(u[1](t))
    end
    ofunc2(x, u, t) = x
    ds = ODESystem(righthandside=sfunc2, readout=ofunc2, state=ones(3), input=Inport(2), output=Outport(3))
    @test isa(ds.input, Inport)
    @test isa(ds.output, Outport)
    @test length(ds.input) == 2
    @test length(ds.output) == 3

    ds = ODESystem(righthandside=sfunc2, readout=nothing, state=ones(3), input=Inport(2), output=nothing)
    @test isa(ds.input, Inport)
    @test length(ds.input) == 2
    @test ds.readout === nothing
    @test ds.output === nothing

    # Driving ODESystem
    sfunc3(dx, x, u, t) = (dx .= -x)
    ofunc3(x, u, t) = x
    ds = ODESystem(righthandside=sfunc3, readout=ofunc3, state=[1.], input=nothing, output=Outport())
    iport = Inport()
    trg = Outpin()
    hnd = Inpin{Bool}()
    connect!(ds.output, iport)
    connect!(trg, ds.trigger)
    connect!(ds.handshake, hnd)
    tsk = launch(ds)
    tsk2 = @async while true
        all(take!(iport) .=== NaN) && break
    end
    for t in 1. : 10.
        put!(trg, t)
        take!(hnd)
        @test ds.t == t
        @test [read(pin.link.buffer) for pin in iport] == ds.state
    end
    put!(trg, NaN)
    sleep(0.1)
    @test istaskdone(tsk)
    put!(ds.output, NaN * ones(length(ds.output)))
    sleep(0.1)
    @test istaskdone(tsk2)

    # LinaerSystem tests
    ds = ContinuousLinearSystem(input=nothing, output=Outport(1))
    @test ds.A == fill(-1, 1, 1)
    @test ds.B == fill(1., 1, 1)
    @test ds.C == fill(1, 1, 1)
    @test ds.D == fill(0, 1, 1)
    @test_throws Exception ds.γ == 1.
    ds = ContinuousLinearSystem(input=nothing, output=Outport(2), A=ones(2,2), C=ones(2,2))
    @test ds.input === nothing
    @test isa(ds.output, Outport)
    @test length(ds.output) == 2
    ds = ContinuousLinearSystem(input=Inport(2), output=Outport(3), A=ones(4,4), B=ones(4,2), C=ones(3,4), D=ones(3, 2), state=zeros(4))
    @test ds.t == 0.
    @test ds.state == zeros(4)

    ds = ContinuousLinearSystem(input=Inport(2), output=Outport(3), A=ones(4,4), B=ones(4,2), C=ones(3,4), D=ones(3, 2))
    @test typeof(ds.integrator.sol.prob.p) <: Interpolant
    @test size(ds.integrator.sol.prob.p.timebuf) == (3,)
    @test size(ds.integrator.sol.prob.p.databuf) == (2, 3)
    @test isa(ds.output, Outport)
    @test length(ds.input) == 2
    @test length(ds.output) == 3
    @test length(ds.state) == 4
    oport = Outport(2)
    iport = Inport(3)
    trg = Outpin()
    hnd = Inpin{Bool}()
    connect!(ds.output, iport)
    connect!(ds.handshake, hnd)
    connect!(trg, ds.trigger)
    connect!(oport, ds.input)
    tsk = launch(ds)
    tsk2 = @async while true
        all(take!(iport) .=== NaN) && break
    end
    for t in 1. : 10.
        put!(trg, t)
        u = [sin(t), cos(t)]
        put!(oport, u)
        take!(hnd)
        @test ds.t == t
        @test [read(pin.link.buffer) for pin in iport] == ds.C * ds.state + ds.D * u
    end
    put!(trg, NaN)
    sleep(0.1)
    @test istaskdone(tsk)
    put!(ds.output, NaN * ones(length(ds.output)))
    sleep(0.1)
    @test istaskdone(tsk2)

    # Other 3-Dimensional AbstractODESystem tests
    for (DSystem, defaults) in zip(
        [
            LorenzSystem,
            ChenSystem,
            ChuaSystem,
            RosslerSystem
        ], [
            (σ = 10., β = 8/3, ρ = 28., γ = 1.),
            (a = 35., b = 3., c = 28., γ = 1.),
            (diode = Causal.PiecewiseLinearDiode(), α = 15.6, β = 28., γ = 1.),
            (a = 0.38, b = 0.3, c = 4.82, γ = 1.)
        ]
        )
        ds = DSystem(input=nothing, output=Outport(3); defaults...)  # System with key-value pairs with no input and bus output.
        ds = DSystem(input=nothing, output=Outport(3))  # System with no input
        @test ds.input === nothing
        @test isa(ds.output, Outport)
        @test length(ds.output) == 3
        @test length(ds.state) == 3
        iport = Inport(3)
        trg = Outpin()
        hnd = Inpin{Bool}()
        connect!(ds.output, iport)
        connect!(trg, ds.trigger)
        connect!(ds.handshake, hnd)
        tsk = launch(ds)
        tsk2 = @async while true
            all(take!(iport) .=== NaN) && break
        end
        for t in 1. : 10.
            put!(trg, t)
            take!(hnd)
            @test ds.t == t
            @test [read(pin.link.buffer) for pin in iport] == ds.state
        end
        put!(trg, NaN)
        sleep(0.1)
        @test istaskdone(tsk)
        put!(ds.output, NaN * ones(length(ds.output)))
        sleep(0.1)
        @test istaskdone(tsk2)
    end

    for (DSystem, defaults) in zip(
        [
            ForcedLorenzSystem,
            ForcedChenSystem,
            ForcedChuaSystem,
            ForcedRosslerSystem
        ], [
            (σ = 10., β = 8/3, ρ = 28., γ = 1.),
            (a = 35., b = 3., c = 28., γ = 1.),
            (diode = Causal.PiecewiseLinearDiode(), α = 15.6, β = 28., γ = 1.),
            (a = 0.38, b = 0.3, c = 4.82, γ = 1.)
        ]
        )
        ds = DSystem(input=Inport(3), output=Outport(3), state=rand(3), cplmat=[1. 0 0; 0 1 0; 0 0 0])
        @test typeof(ds.integrator.sol.prob.p) <: Interpolant
        @test size(ds.integrator.sol.prob.p.timebuf) == (3,)
        @test size(ds.integrator.sol.prob.p.databuf) == (3,3)
        @test isa(ds.input, Inport)
        @test isa(ds.output, Outport)
        @test length(ds.input) == 3
        @test length(ds.output) == 3
        iport = Inport(3)
        oport = Outport(3)
        trg = Outpin()
        hnd = Inpin{Bool}()
        connect!(ds.output, iport)
        connect!(oport, ds.input)
        connect!(trg, ds.trigger)
        connect!(ds.handshake, hnd)
        tsk = launch(ds)
        tsk2 = @async while true
            all(take!(iport) .=== NaN) && break
        end
        for t in 1. : 10.
            put!(trg, t)
            u = [sin(t), cos(t), log(t)]
            put!(oport, u)
            take!(hnd)
            @test ds.t == t
        end
        put!(trg, NaN)
        sleep(0.1)
        @test istaskdone(tsk)
        put!(ds.output, NaN * ones(length(ds.output)))
        sleep(0.1)
        @test istaskdone(tsk2)
    end

    # VanderpolSystem tests.
    ds = VanderpolSystem(input=nothing, output=Outport(2))  # System with no input
    @test ds.input === nothing
    @test isa(ds.output, Outport)
    @test length(ds.output) == 2
    @test length(ds.state) == 2
    @test ds.μ == 5.
    @test ds.γ == 1.
    iport = Inport(2)
    trg = Outpin()
    hnd = Inpin{Bool}()
    connect!(ds.output, iport)
    connect!(trg, ds.trigger)
    connect!(ds.handshake, hnd)
    tsk = launch(ds)
    tsk2 = @async while true
        all(take!(iport) .=== NaN) && break
    end
    for t in 1. : 10.
        put!(trg, t)
        take!(hnd)
        @test ds.t == t
        @test [read(pin.link.buffer) for pin in iport] == ds.state
    end
    put!(trg, NaN)
    sleep(0.1)
    @test istaskdone(tsk)
    put!(ds.output, NaN * ones(length(ds.output)))
    sleep(0.1)
    @test istaskdone(tsk2)

    ds = ForcedVanderpolSystem(input=Inport(2), output=Outport(2), cplmat=[1 0; 0 0])  # Add input values to 1. state
    @test isa(ds.input, Inport)
    @test isa(ds.output, Outport)
    @test length(ds.input) == 2
    @test length(ds.output) == 2
    iport = Inport(2)
    oport = Outport(2)
    trg = Outpin()
    hnd = Inpin{Bool}()
    connect!(oport, ds.input)
    connect!(ds.output, iport)
    connect!(trg, ds.trigger)
    connect!(ds.handshake, hnd)
    tsk = launch(ds)
    tsk2 = @async while true
        all(take!(iport) .=== NaN) && break
    end
    for t in 1. : 10.
        put!(trg, t)
        u = [sin(t), cos(t)]
        put!(oport, u)
        take!(hnd)
        @test ds.t == t
        @test [read(pin.link.buffer) for pin in iport] == ds.state
    end
    put!(trg, NaN)
    sleep(0.1)
    @test istaskdone(tsk)
    put!(ds.output, NaN * ones(length(ds.output)))
    sleep(0.1)
    @test istaskdone(tsk2)

    # Test defining new ODE systems
    # The type must be mutabe
    @test_throws Exception @eval @def_ode_system struct ODESystem{RH, RO, ST, IP, OP}
        righthandside::RH
        readout::RO
        state::ST
        input::IP
        output::OP
    end

    # The type must be a subtype of AbstractODESystem
    @test_throws Exception @eval @def_ode_system mutable struct ODESystem{RH, RO, ST, IP, OP}
        righthandside::RH
        readout::RO
        state::ST
        input::IP
        output::OP
    end

    # The type must be subtype of AbstractODESystem
    @test_throws Exception @eval @def_ode_system mutable struct ODESystem{RH, RO, ST, IP, OP} <: MyDummyyAbstractODESystem
        righthandside::RH
        readout::RO
        state::ST
        input::IP
        output::OP
    end

    @info "Done ODESystemTestSet."
end  # testset


text:
libpath = "/home/hessammehr/harminv-1.4.1/.libs"

const C = Complex128
const F = Float64

struct HarminvNative
    c::Ptr{C}
    n::Int
    K::Int
    J::Int
    nfreqs::Int
    fmin::F
    fmax::F
    z::Ptr{C}
    U0::Ptr{C}
    U1::Ptr{C}
    G0::Ptr{C}
    G0_M::Ptr{C}
    D0::Ptr{C}
    B::Ptr{C}
    u::Ptr{C}
    amps::Ptr{C}
    errs::Ptr{F}
end

if !(libpath in Base.LD_LOAD_PATH)
    push!(Base.LD_LOAD_PATH, libpath)
end

ccall((:harminv_data_create, "libharminv"), )


text:
struct Multiplication{T<:Sequence{<:SequenceSpace}}
    sequence :: T
end

(ℳ::Multiplication)(a::Sequence) = *(ℳ, a)
Base.:*(ℳ::Multiplication, a::Sequence) = *(ℳ.sequence, a)

Base.:+(ℳ::Multiplication) = Multiplication(+(ℳ.sequence))
Base.:-(ℳ::Multiplication) = Multiplication(-(ℳ.sequence))
Base.:^(ℳ::Multiplication, n::Int) = Multiplication(ℳ.sequence ^ n)

for f ∈ (:+, :-, :*)
    @eval begin
        Base.$f(ℳ₁::Multiplication, ℳ₂::Multiplication) = Multiplication($f(ℳ₁.sequence, ℳ₂.sequence))
        Base.$f(a::Number, ℳ::Multiplication) = Multiplication($f(a, ℳ.sequence))
        Base.$f(ℳ::Multiplication, a::Number) = Multiplication($f(ℳ.sequence, a))
    end
end

Base.:/(ℳ::Multiplication, a::Number) = Multiplication(/(ℳ.sequence, a))
Base.:\(a::Number, ℳ::Multiplication) = Multiplication(\(a, ℳ.sequence))

LinearAlgebra.opnorm(ℳ::Multiplication, d::Norm) = norm(ℳ.sequence, d)

function project(ℳ::Multiplication, domain::SequenceSpace, codomain::SequenceSpace, ::Type{T}=eltype(ℳ.sequence)) where {T}
    _iscompatible(domain, codomain) & _iscompatible(space(ℳ.sequence), domain) || return throw(DimensionMismatch)
    C = LinearOperator(domain, codomain, zeros(T, dimension(codomain), dimension(domain)))
    _project!(C, ℳ)
    return C
end

function project!(C::LinearOperator{<:SequenceSpace,<:SequenceSpace}, ℳ::Multiplication)
    domain_C = domain(C)
    _iscompatible(domain_C, codomain(C)) & _iscompatible(space(ℳ.sequence), domain_C) || return throw(DimensionMismatch)
    coefficients(C) .= zero(eltype(C))
    _project!(C, ℳ)
    return C
end

#

function _project!(C::LinearOperator{<:TensorSpace,<:TensorSpace}, ℳ::Multiplication)
    space_ℳ = space(ℳ.sequence)
    @inbounds for β ∈ _mult_domain_indices(domain(C)), α ∈ indices(codomain(C))
        if _isvalid(space_ℳ, α, β)
            C[α,_extract_valid_index(space_ℳ, β, zero.(α))] += ℳ.sequence[_extract_valid_index(space_ℳ, α, β)]
        end
    end
    return C
end

_mult_domain_indices(s::TensorSpace) = Base.Iterators.ProductIterator(map(_mult_domain_indices, spaces(s)))

_isvalid(s::TensorSpace{<:NTuple{N,BaseSpace}}, α::NTuple{N,Int}, β::NTuple{N,Int}) where {N} =
    @inbounds _isvalid(s[1], α[1], β[1]) & _isvalid(Base.tail(s), Base.tail(α), Base.tail(β))
_isvalid(s::TensorSpace{<:Tuple{BaseSpace}}, α::Tuple{Int}, β::Tuple{Int}) =
    @inbounds _isvalid(s[1], α[1], β[1])

_extract_valid_index(s::TensorSpace{<:NTuple{N,BaseSpace}}, α::NTuple{N,Int}, β::NTuple{N,Int}) where {N} =
    @inbounds (_extract_valid_index(s[1], α[1], β[1]), _extract_valid_index(Base.tail(s), Base.tail(α), Base.tail(β))...)
_extract_valid_index(s::TensorSpace{<:Tuple{BaseSpace}}, α::Tuple{Int}, β::Tuple{Int}) =
    @inbounds (_extract_valid_index(s[1], α[1], β[1]),)

# Taylor

function _project!(C::LinearOperator{Taylor,Taylor}, ℳ::Multiplication)
    order_codomain = order(codomain(C))
    ord = order(ℳ.sequence)
    @inbounds for j ∈ indices(domain(C)), i ∈ j:min(order_codomain, ord+j)
        C[i,j] = ℳ.sequence[i-j]
    end
    return C
end

_mult_domain_indices(s::Taylor) = indices(s)
_isvalid(s::Taylor, i::Int, j::Int) = 0 ≤ i-j ≤ order(s)
_extract_valid_index(::Taylor, i::Int, j::Int) = i-j

# Fourier

function _project!(C::LinearOperator{<:Fourier,<:Fourier}, ℳ::Multiplication)
    order_codomain = order(codomain(C))
    ord = order(ℳ.sequence)
    @inbounds for j ∈ indices(domain(C)), i ∈ max(-order_codomain, -ord+j):min(order_codomain, ord+j)
        C[i,j] = ℳ.sequence[i-j]
    end
    return C
end

_mult_domain_indices(s::Fourier) = indices(s)
_isvalid(s::Fourier, i::Int, j::Int) = abs(i-j) ≤ order(s)
_extract_valid_index(::Fourier, i::Int, j::Int) = i-j

# Chebyshev

function _project!(C::LinearOperator{Chebyshev,Chebyshev}, ℳ::Multiplication)
    order_codomain = order(codomain(C))
    ord = order(ℳ.sequence)
    @inbounds for j ∈ indices(domain(C)), i ∈ indices(codomain(C))
        if abs(i-j) ≤ ord
            if j == 0
                C[i,j] = ℳ.sequence[i]
            else
                C[i,j] = ℳ.sequence[abs(i-j)]
                idx2 = i+j
                if idx2 ≤ ord
                    C[i,j] += ℳ.sequence[idx2]
                end
            end
        end
    end
    return C
end

_mult_domain_indices(s::Chebyshev) = -order(s):order(s)
_isvalid(s::Chebyshev, i::Int, j::Int) = abs(i-j) ≤ order(s)
_extract_valid_index(::Chebyshev, i::Int, j::Int) = abs(i-j)


text:
using Rimu
using BenchmarkTools

const SUITE = @benchmarkset "Rimu" begin
    @case "(10, 20) Mom space with projected energy and initiator" begin
        add = BoseFS(ntuple(i -> ifelse(i == 10, 10, 0), 20))
        ham = HubbardMom1D(add, u=6.0)
        dv = InitiatorDVec(add => 1.0; style=IsDynamicSemistochastic())
        post_step = ProjectedEnergy(ham, dv)
        s_strat = DoubleLogUpdate(targetwalkers=20_000)

        lomc!(ham, dv; s_strat, post_step, dτ=1e-4, laststep=4000)
    end seconds=150
    @case "(4+1, 11) 2C Mom space with G2Correlators" begin
        add = BoseFS2C(ntuple(i -> ifelse(i == 5, 4, 0), 11), ntuple(==(5), 11))
        ham = BoseHubbardMom1D2C(add, v=0.1)
        dv = DVec(add => 1.0f0; style=IsDynamicSemistochastic{Float32}())
        s_strat = DoubleLogUpdate(targetwalkers=10_000)
        replica = AllOverlaps(2, ntuple(i -> G2Correlator(i - 1), 7))

        lomc!(ham, dv; s_strat, replica, laststep=2000)
    end seconds=150
    @case "(50, 50) Real space" begin
        add = near_uniform(BoseFS{50,50})
        ham = HubbardReal1D(add, u=6.0)
        dv = DVec(add => 1.0; style=IsDynamicSemistochastic())
        s_strat = DoubleLogUpdate(targetwalkers=50_000)

        lomc!(ham, dv; s_strat, dτ=1e-4, laststep=2000)
    end seconds=150
    @case "(5+3, 12) Real space fermions in a 4×3 lattice with initiator" begin
        add = CompositeFS(
            near_uniform(FermiFS{5,12}),
            near_uniform(FermiFS{3,12}),
        )
        ham = HubbardRealSpace(add, geometry=PeriodicBoundaries(4, 3))
        dv = InitiatorDVec(add => 1.0, style=IsDynamicSemistochastic())
        s_strat = DoubleLogUpdate(targetwalkers=10_000)

        lomc!(ham, dv; s_strat, laststep=15_000, dτ=1e-2)
    end seconds=150
end


text:
module Day4
export solvePart1, solvePart2

function getSortedInput()
    inputs = Array{String}(undef, 0)
    open("day4input.txt") do f
        for line in eachline(f)
            push!(inputs, string(line))
        end
    end
    return sort(inputs)
end

function extractMinute(line)
    return parse(Int64, line[16:17])
end

function extractGuardId(lines)
    split(lines[1])[4][2:end]
end

function partitionByGuard(sortedInput)
    guards = []
    guard = []
    for line in sortedInput
        if occursin("Guard", line) && length(guard) > 0
            push!(guards, guard)
            guard = []
        end
        push!(guard, line)
    end

    dict = Dict()
    for guard in guards
        id = extractGuardId(guard)
        dict[id] = append!(get(dict, id, []), guard)
    end
    return collect(values(dict))
end

function lineContainsSleepOrWakeTime(line)
    occursin(r"asleep|wakes", line)
end

function sumTimeAsleepForGuard(lines)
    id = extractGuardId(lines)
    sleepAndWakeMinutes = map(extractMinute, filter(lineContainsSleepOrWakeTime, lines))
    minutesAsleep = Iterators.partition(sleepAndWakeMinutes, 2) |>
        iter -> map(xy -> xy[2] - xy[1], iter) |>
        iter -> reduce(+, iter)
    return (id, minutesAsleep)
end

function findGuardAsleepMost(inputs)
    guards = partitionByGuard(inputs)
    guardIdToTimeAsleep = map(sumTimeAsleepForGuard, guards)
    max = 0
    maxId = "N/A"
    for (id, timeAsleep) in guardIdToTimeAsleep
        if timeAsleep > max
            max = timeAsleep
            maxId = id
        end
    end
    return (maxId, max)
end

function getGuardLinesForId(id, inputs)
    guards = partitionByGuard(inputs)
    for guard in guards
        if occursin("Guard #$id", guard[1])
            return guard
        end
    end
end

function findMinuteGuardIsAsleepMost(lines)
    sleepAndWakeMinutes = map(extractMinute, filter(lineContainsSleepOrWakeTime, lines))
    daysAsleepByMinute = Dict()
    for (asleep, awake) in Iterators.partition(sleepAndWakeMinutes, 2)
        for m in asleep:awake-1
            daysAsleepByMinute[m] = get(daysAsleepByMinute, m, 0) + 1
        end
    end
    max = 0
    maxMinute = -1
    for (minute, timesAsleep) in daysAsleepByMinute
        if timesAsleep > max
            max = timesAsleep
            maxMinute = minute
        end
    end
    return (max, maxMinute)
end

function solvePart1()
    inputs = getSortedInput()
    (id, minutes) = findGuardAsleepMost(inputs)
    guard = getGuardLinesForId(id, inputs)
    (max, minute) = findMinuteGuardIsAsleepMost(guard)
    println("Guard #$id is asleep the most with $minutes minutes asleep.")
    println("Guard #$id is asleep most frequently in minute $minute: $max times")
    return parse(Int, id) * minute
end

function solvePart2()
    inputs = getSortedInput()
    guards = partitionByGuard(inputs)
    max = 0
    maxMinute = -1
    maxId = "n/a"
    for guard in guards
        (timesAsleep, minute) = findMinuteGuardIsAsleepMost(guard)
        if timesAsleep >= max
            max = timesAsleep
            maxMinute = minute
            maxId = extractGuardId(guard)
        end
    end
    println("Guard #$maxId is asleep most in one minute, with $max times in minute $maxMinute")
    return parse(Int, maxId) * maxMinute
end

end


text:
module NumericalIntegration

using LinearAlgebra
using Logging
using Base.Iterators
using Interpolations

export integrate, cumul_integrate
export Trapezoidal, TrapezoidalEven, TrapezoidalFast, TrapezoidalEvenFast
export SimpsonEven, SimpsonEvenFast
export RombergEven
export IntegrationMethod

abstract type IntegrationMethod end

struct Trapezoidal         <: IntegrationMethod end
struct TrapezoidalEven     <: IntegrationMethod end
struct TrapezoidalFast     <: IntegrationMethod end
struct TrapezoidalEvenFast <: IntegrationMethod end
struct SimpsonEven         <: IntegrationMethod end # https://en.wikipedia.org/wiki/Simpson%27s_rule#Alternative_extended_Simpson.27s_rule
struct SimpsonEvenFast     <: IntegrationMethod end
struct RombergEven{T<:AbstractFloat} <: IntegrationMethod
    acc::T
end # https://en.wikipedia.org/wiki/Romberg%27s_method
RombergEven() = RombergEven(1e-12)

const HALF = 1//2

#documentation

"""
    integrate(x,y...)

Compute numerical integral of y(x) from x=x[1] to x=x[end]. Return a scalar of the same type as the input. If no method is supplied, use Trapezdoial.
"""
function integrate end

"""
    cumul_integrate(x,y...)

Compute cumulative numerical integral of y(x) from x=x[1] to x=x[end]. Return a vector with elements of the same type as the input. If no method is supplied, use Trapezdoial.
"""
function cumul_integrate end

#implementation

"""
    integrate(x::AbstractVector, y::AbstractVector, ::Trapezoidal)

Use Trapezoidal rule. This is the default when no method is supplied.
"""
function integrate(x::AbstractVector, y::AbstractVector, ::Trapezoidal)
    n = length(x)
    length(x) == length(y) || error("x and y vectors must be of the same length!")
    n ≥ 2 || error("vectors must contain at least two elements")

    return integrate(x, y, TrapezoidalFast())
end
function integrate(x::AbstractRange, y::AbstractVector, ::Trapezoidal)
    return integrate(x, y, TrapezoidalEven())
end

"""
    integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEven)

Use Trapezoidal rule, assuming evenly spaced vector x.
"""
function integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEven)
    n = length(x)
    length(y) == n || error("x and y vectors must be of the same length!")
    n ≥ 2 || error("vectors must contain at least two elements")

    return integrate(x, y, TrapezoidalEvenFast())
end

"""
    integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalFast)

Use Trapezoidal rule. Unsafe method: no bound checking.
"""
function integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalFast)
    @inbounds retval = (x[2] - x[1]) * (y[1] + y[2])
    @inbounds @fastmath @simd for i in 2:(length(y) - 1)
        retval += (x[i+1] - x[i]) * (y[i] + y[i+1])
    end
    return HALF * retval
end
function integrate(x::AbstractRange, y::AbstractVector, ::TrapezoidalFast)
    return integrate(x, y, TrapezoidalEvenFast())
end

"""
    integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEvenFast)

Use Trapezoidal rule, assuming evenly spaced vector x. Unsafe method: no bound checking.
"""
function integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEvenFast)
    if length(x) < 3
        @inbounds return HALF * (x[2] - x[1]) * (y[1] + y[2])
    else
        @inbounds retval = y[2]
        @inbounds @fastmath @simd for i in 3:(length(y) - 1)
            retval += y[i]
        end
        @inbounds return (x[2] - x[1]) * (retval + HALF * (y[1] + y[end]))
    end
end

"""
    integrate(x::AbstractVector, y::AbstractVector, ::SimpsonEven)

Use Simpson's rule, assuming evenly spaced vector x.
"""
function integrate(x::AbstractVector, y::AbstractVector, ::SimpsonEven)
    length(x) == length(y) || error("x and y vectors must be of the same length!")
    length(x) ≥ 4 || error("vectors must contain at least 4 elements")

    return integrate(x, y, SimpsonEvenFast())
end

"""
    integrate(x::AbstractVector, y::AbstractVector, ::SimpsonEven)

Use Simpson's rule, assuming evenly spaced vector x.  Unsafe method: no bound checking.
"""
function integrate(x::AbstractVector, y::AbstractVector, ::SimpsonEvenFast)
    @inbounds retval = (17 * (y[1] + y[end]) + 59 * (y[2] + y[end-1]) +
                        43 * (y[3] + y[end-2]) + 49 * (y[4] + y[end-3])) / 48
    @fastmath @inbounds for i in 5:(length(y) - 4)
        retval += y[i]
    end
    @inbounds return (x[2] - x[1]) * retval
end

"""
    integrate(x::AbstractVector, y::AbstractMatrix, method; dims=2)

When y is an array, compute integral along dimension specified by dims (default 2: columns).
"""
function integrate(x::AbstractVector, y::AbstractMatrix, M::IntegrationMethod; dims=2)
    out = [integrate(x,selectdim(y,dims,j),M) for j=1:size(y,dims)]
    return out
end

function integrate(x::AbstractVector, y::AbstractVector, m::RombergEven)
    @assert length(x) == length(y) "x and y vectors must be of the same length!"
    @assert ((length(x) - 1) & (length(x) - 2)) == 0 "Need length of vector to be 2^n + 1"
    maxsteps::Integer = Int(log2(length(x)-1))
    @inbounds h = x[end] - x[1]
    @inbounds v = (y[1] + y[end]) * h * HALF
    rombaux = Matrix{typeof(v)}(undef, maxsteps, 2)
    rombaux[1,1] = v
    prevcol = 1
    currcol = 2
    # Precomputed values for norm check
    acc_compare = m.acc * oneunit(norm(v,Inf))
    minsteps = maxsteps ÷ 3
    @inbounds for i in 1 : (maxsteps-1)
        h *= HALF
        npoints = 1 << (i-1)
        jumpsize = div(length(x)-1, 2*npoints)
        c = zero(eltype(y))
        for j in 1 : npoints
            c += y[1 + (2*j-1)*jumpsize]
        end
        rombaux[1, currcol] = h*c + HALF*rombaux[1, prevcol]
        for j in 2 : (i+1)
            n_k = 4^(j-1)
            rombaux[j, currcol] = (n_k*rombaux[j-1, currcol] - rombaux[j-1, prevcol])/(n_k - 1)
        end

        # Clunky way to address Unitful compatibility, while also allowing for vectors
        normval = norm(rombaux[i, prevcol] - rombaux[i+1, currcol], Inf)
        if i > minsteps && normval < acc_compare
            return rombaux[i+1, currcol]
        end

        prevcol, currcol = currcol, prevcol
    end
    finalerr = norm(rombaux[maxsteps-1, currcol] - rombaux[maxsteps, prevcol], Inf)
    @warn "RombergEven :: final step reached, but accuracy not: $finalerr > $(m.acc)"
    @inbounds return rombaux[maxsteps, prevcol]
end

@inline function _midpoints(x::AbstractVector{T}) where T
    length(x) == 1 && return x
    retval = Vector{T}(undef, length(x)-1)
    @simd for i in eachindex(retval)
        retval[i] = HALF * (x[i] + x[i+1])
    end
    return retval
end
@inline function _midpoints(x::AbstractRange)
    length(x) == 1 && return x
    Δx = HALF*step(x)
    return range(first(x)+Δx, stop=last(x)-Δx, length=length(x)-1)
end

function integrate(X::NTuple{N,AbstractVector}, Y::AbstractArray{T,N}, ::Trapezoidal) where {T,N}
    @assert length.(X) == size(Y)
    return integrate(X, Y, TrapezoidalFast())
end

function integrate(X::NTuple{N,AbstractVector}, Y::AbstractArray{T,N}, ::TrapezoidalFast) where {T,N}
    midpnts = map(_midpoints, X)
    Δ(x::AbstractVector) = length(x) > 1 ? diff(x) : 1
    Δxs = map(Δ, X)
    interp = LinearInterpolation(X,Y)
    f((Δx,x)) = prod(Δx)*interp(x...)
    return sum(f, zip(product(Δxs...), product(midpnts...)))
end

function integrate(X::NTuple{N,AbstractVector}, Y::AbstractArray{T,N}, ::TrapezoidalEvenFast) where {T,N}
    midpnts = map(_midpoints, X)
    Δ(x::AbstractVector) = length(x) > 1 ? x[2] - x[1] : 1
    Δx = prod(Δ, X)
    interp = LinearInterpolation(X,Y)
    f(x) = interp(x...)
    return Δx*sum(f, product(midpnts...))
end

function integrate(X::NTuple{N,AbstractRange}, Y::AbstractArray{T,N}, ::TrapezoidalEvenFast) where {T,N}
    midpnts = map(_midpoints, X)
    Δ(x::AbstractVector) = length(x) > 1 ? step(x) : 1
    Δx = prod(Δ, X)
    interp = LinearInterpolation(X,Y)
    f(x) = interp(x...)
    return Δx*sum(f, product(midpnts...))
end

#function integrate(X::Tuple{AbstractVector}, Y::AbstractVector{T}, M::IntegrationMethod) :: T where {T}
#    return integrate(X[1], Y, M)
#end
#"""
#    integrate(X::NTuple{N,AbstractVector}, Y::AbstractArray{T,N}, method, cache=nothing)
#
#Given an n-dimensional grid of values, compute the total integral along each dim
#"""
#function integrate(X::NTuple{N,AbstractVector}, Y::AbstractArray{T,N}, M::IntegrationMethod) :: T where {T,N}
#    n = last(size(Y))
#    cache = Vector{T}(undef, n)
#    x = X[1:N-1]
#    @inbounds for i in 1:n
#        cache[i] = integrate(x, selectdim(Y,N,i), M)
#    end
#    return integrate(X[end], cache, M)
#end

# cumulative integrals

"""
    cumul_integrate(x::AbstractVector, y::AbstractVector, ::Trapezoidal)

Use Trapezoidal rule. This is the default when no method is supplied.
"""
function cumul_integrate(x::AbstractVector, y::AbstractVector, ::Trapezoidal)
    @assert length(x) == length(y) "x and y vectors must be of the same length!"
    length(x) ≥ 2 || error("vectors must contain at least two elements")

    return cumul_integrate(x, y, TrapezoidalFast())
end
function cumul_integrate(x::AbstractRange, y::AbstractVector, ::Trapezoidal)
    return cumul_integrate(x, y, TrapezoidalEven())
end

"""
    cumul_integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEven)

Use Trapezoidal rule, assuming evenly spaced vector x.
"""
function cumul_integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEven)
    @assert length(x) == length(y) "x and y vectors must be of the same length!"
    length(x) ≥ 2 || error("vectors must contain at least two elements")

    return cumul_integrate(x, y, TrapezoidalEvenFast())
end

"""
    cumul_integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalFast)

Use Trapezoidal rule. Unsafe method: no bound checking.
"""
function cumul_integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalFast)
    # compute initial value
    @inbounds init = (x[2] - x[1]) * (y[1] + y[2])
    n = length(x)
    retarr = Vector{typeof(init)}(undef, n)
    retarr[1] = init

    # for all other values
    @inbounds @fastmath for i in 2:n # not sure if @simd can do anything here
        retarr[i] = retarr[i-1] + (x[i] - x[i-1]) * (y[i] + y[i-1])
    end

    return HALF * retarr
end

"""
    cumul_integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEvenFast)

Use Trapezoidal rule, assuming evenly spaced vector x. Unsafe method: no bound checking.
"""
function cumul_integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEvenFast)
    @inbounds init = y[1] + y[2]
    n = length(x)
    retarr = Vector{typeof(init)}(undef, n)
    retarr[1] = init

    # for all other values
    @fastmath for i in 2 : length(y)
        @inbounds retarr[i] = retarr[i-1] + (y[i] + y[i-1])
    end
    @inbounds return (x[2] - x[1]) * HALF * retarr
end

"""
    cumul_integrate(x::AbstractVector, y::AbstractMatrix, method; dims=2)

When y is an array, compute integral along each dimension specified by dims (default 2: columns)
"""
function cumul_integrate(x::AbstractVector, y::AbstractMatrix, M::IntegrationMethod; dims=2)
    return hcat([cumul_integrate(x,selectdim(y,dims,j),M) for j=1:size(y,dims)]...)
end

#default behaviour
integrate(x, y; kwargs...) = integrate(x, y, Trapezoidal(); kwargs...)
cumul_integrate(x, y; kwargs...) = cumul_integrate(x, y, Trapezoidal(); kwargs...)

end # module


text:
include("task_func.jl")

function movements!(r, n)
    while !isborder(r, Ost)
        up = 0
        for _ = 1:n
            shag = goToSideWithShag!(r, Ost, n-1, save = true, checkBorder = true, mark = true)[1]
            goToSideWithShag!(r, West, shag)
            putmarker!(r)

            up += goToSideWithShag!(r, Nord, 1, save = true, checkBorder = true)[1]
        end
        goToSideWithShag!(r, Sud, up)
        goToSideWithShag!(r, Ost, 2*n, checkBorder = true)
    end

    goToSide!(r, West)
    goToSideWithShag!(r, Nord, n, checkBorder = true)
end

function draw!(r, n)
    listOfShags = []
    push!(listOfShags, goToSide!(r, Sud, save = true))
    push!(listOfShags, goToSide!(r, West, save = true))

    while !isborder(r, Nord)
        movements!(r, n)

        if isborder(r, Nord) && ismarker(r) break end
        goToSideWithShag!(r, Ost, n, checkBorder = true)

        movements!(r, n)
    end

    goToWestSud!(r)
    goToStart!(r, listOfShags)
end

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
"""
TEST FILE FOR VARIOUS TEMPORAL META GRAPH THINGS
"""

using QuNet
using LightGraphs, SimpleWeightedGraphs, GraphPlot, MetaGraphs, Plots, Colors, Statistics
using LaTeXStrings
using JLD
using Parameters

datafile = "data/timedepth_test"

num_trials = 100::Int64
timedepth = 5
num_pairs = 8::Int64
grid_size = 4::Int64
async_cost = 0.001

G = GridNetwork(grid_size, grid_size)
# Create a Temporal Graph from G with timedepth i
T = QuNet.TemporalGraph(G, timedepth, memory_costs = zero_costvector())
# Get data
p, p_e, pat, pat_e = net_performance(T, num_trials, num_pairs, max_paths=4,
async_cost=async_cost, src_layer=-1, dst_layer=3)

println("p")
println("$p")
println("pat")
println("$pat")

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
# plotting (ζ -> 2020's relative deviation from benchmark)
afunc_PctRelDeviaFrom1stElement( x::AbstractArray ) = (x ./ x[1] .- 1.0) .* 100
figure(figsize = (8,4))
      plot( afunc_PctRelDeviaFrom1stElement(res_ContriRatRise["zetalevels"]) ,
            afunc_PctRelDeviaFrom1stElement(res_ContriRatRise["gap/exp"]), "-" )
      plot( afunc_PctRelDeviaFrom1stElement(res_ContriRatRise["zetalevels"]) ,
            afunc_PctRelDeviaFrom1stElement(res_ContriRatRise["gap/in"]), "-.b" )
      # plot( afunc_PctRelDeviaFrom1stElement(res_ContriRatRise["zetalevels"]) ,
      #       afunc_PctRelDeviaFrom1stElement(res_ContriRatRise["gap/gdp"]), "--r" )
      # plot( afunc_PctRelDeviaFrom1stElement(res_ContriRatRise["zetalevels"]) ,
      #       afunc_PctRelDeviaFrom1stElement(res_ContriRatRise["gap/taxrev"]), ":m" )
      ylabel("Relative Deviation from Benchmark (%)");
      xlabel("ζ - Relative Deviation from Benchmark (%)");
      grid(true);
      legend([ "Gap / Pool Expenditure",
               # "Gap / Tax Revenues"
               # "Gap / GDP",
               "Gap / Pool Income", ], loc = "best");
      tight_layout()
# save figure
savefig( "./output/ContributionRateRise_Elasticity.pdf", format = "pdf" )

#


text:
## methods.jl : refactored methods from run.jl
#
# Implements refactored methods from run.jl.
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods (c) Pieterjan Robbe, 2019

#
# basics
#
ndims(::Estimator{<:AbstractIndexSet{d}}) where d = d

get_index_set(estimator::Estimator, sz) = get_index_set(estimator.index_set, sz)

get_tols(estimator::Estimator, tol::T) where T<:Real = estimator[:continuate] ? estimator[:continuation_mul_factor].^(estimator[:nb_of_tols]-1:-1:0)*tol : T[tol]

mse(estimator::Estimator) = varest(estimator) + bias(estimator)^2

rmse(estimator::Estimator) = sqrt(mse(estimator))

converged(estimator::Estimator, ϵ::Real, θ::Real) = ( bias(estimator)^2 ≤ (1-θ)*ϵ^2 || mse(estimator) ≤ ϵ^2 )

max_level_exceeded(estimator::Estimator) = sz(estimator) ≥ estimator[:max_index_set_param]

#
# inspector functions: mean, var, varest...
#
qoi_with_max_var(estimator::Estimator{<:AbstractIndexSet, <:MC}) = estimator[:qoi_with_max_var] == 0 ? argmax(map(n_qoi -> sum(var(samples_diff(estimator, n_qoi, index)) for index in keys(estimator)), 1:estimator[:nb_of_qoi])) : estimator[:qoi_with_max_var]

cost(estimator::Estimator, index::Index) = estimator[:cost_model] isa EmptyFunction ? time(estimator, index) : work(estimator, index)

for f in [:mean, :var]
    @eval begin

        $f(estimator::Estimator{<:AbstractIndexSet, <:MC}, index::Index) = $f(samples_diff(estimator, qoi_with_max_var(estimator), index))

        $(Symbol(f, 0))(estimator::Estimator{<:AbstractIndexSet, <:MC}, index::Index) = $f(samples(estimator, qoi_with_max_var(estimator), index))

        $f(estimator::Estimator{<:AbstractIndexSet, <:QMC}, index::Index) = mean($f(samples_diff(estimator, qoi_with_max_var(estimator), shift, index)) for shift in 1:estimator[:nb_of_shifts](index))

        $(Symbol(f, 0))(estimator::Estimator{<:AbstractIndexSet, <:QMC}, index::Index) = mean($f(samples(estimator, qoi_with_max_var(estimator), shift, index)) for shift in 1:estimator[:nb_of_shifts](index))

        $f(estimator::Estimator) = sum($f(estimator, index) for index in keys(estimator))

    end
end

varest(estimator::Estimator, index::Index) = var(estimator, index)/nb_of_samples(estimator, index)

varest(estimator::Estimator) =  sum(varest(estimator, index) for index in keys(estimator))

#
# rates
#
for (f, g, sgn) in zip([:α, :β, :γ], [:mean, :var, :cost], [-1, -1, 1])
    @eval begin

        $f(estimator::Estimator{<:SL}) = nothing

        $f(estimator::Estimator{<:AbstractIndexSet}) = $sgn.*getindex.(broadcast(i->$(Symbol("rates_", f))(estimator, i), 1:ndims(estimator)), 2)

        $(Symbol("rates_", f))(estimator::Estimator{<:AbstractML}) = $(Symbol("rates_", f))(estimator, 1)

        $(Symbol("rates_", f))(estimator::Estimator{<:AbstractIndexSet}, dir::Integer) = $(Symbol("rates_", f))(estimator, (maximum(getindex.(keys(estimator), dir)) + 1) * Index(ntuple(i -> i == dir, ndims(estimator))), dir)

        function $(Symbol("rates_", f))(estimator::Estimator{<:AbstractIndexSet}, idx::Index, dir::Integer)
            m = idx[dir] - 1
            if m < 2
                return (NaN, NaN)
            else
                x = m:-1:1
                y = map(xᵢ -> $g(estimator, idx - xᵢ * Index(ntuple(i -> i == dir, ndims(estimator)))), 1:m)
                idcs = .!(isnan.(y) .| iszero.(y))
                θ = interp1(view(x, idcs), log2.(abs.(view(y, idcs))))
                return tuple(θ...)
            end
        end

        function $(Symbol("_regress_", g))(estimator::Estimator{<:AbstractIndexSet, <:MC}, index::Index)
            p = broadcast(dir->$(Symbol("rates_", f))(estimator, index, dir), 1:ndims(estimator))
            estimates = broadcast(dir->2^(p[dir][1]+index[dir]*p[dir][2]), 1:ndims(estimator))
            estimate = mean(filter(!isnan, estimates))
            if isnan(estimate)
                p = interp($g, estimator)
                return 2^(p[1]+sum(p[2:end].*index.I))
            else
                return estimate
            end
        end
    end
end

#
# regression
#
function interp1(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})
    A = hcat(ones(eltype(x), length(x)), x)
    A\y
end

function interp(f::Function, estimator::Estimator)
	filter_fcn = estimator isa Estimator{<:AbstractIndexSet, QMC} ? i -> !isempty(samples(estimator)[1][i]) : i -> length(samples(estimator)[1][i]) > 1
    idx_set = filter(filter_fcn, CartesianIndices(size(samples(estimator)[1])))
    A = [i == 0 ? 1 : getindex(index - one(index), i) for index in idx_set, i in 0:ndims(estimator)]
    y = map(i -> log2(f(estimator, i - one(i))), idx_set)
    try
        return A \ y
    catch e
        return fill(NaN, length(y))
    end
end

regress_mean(estimator, index) = _regress_mean(estimator, index)
regress_var(estimator, index) = _regress_var(estimator, index)
regress_cost(estimator, index) = estimator[:cost_model] isa EmptyFunction ? _regress_cost(estimator, index) : estimator[:cost_model](index)

regress_nb_of_samples(estimator::Estimator{<:SL}, index_set, ϵ::Real, θ::Real, L::Integer) = Dict(Level(0) => estimator[:nb_of_warm_up_samples])

regress_nb_of_samples(estimator::Estimator{<:SL, <:QMC}, index_set, ϵ::Real, θ::Real, L::Integer) = Dict(Level(0) => estimator[:nb_of_warm_up_samples])

regress_nb_of_samples(estimator::Estimator{<:AbstractIndexSet, <:QMC}, index_set, ϵ::Real, θ::Real, L::Integer) = Dict(index => estimator[:nb_of_warm_up_samples] for index in index_set)

function regress_nb_of_samples(estimator::Estimator, index_set, ϵ::Real, θ::Real, L::Integer)
    if estimator[:do_regression] && L > 2
        return _regress_nb_of_samples(estimator, index_set, ϵ, θ)
    else
        return Dict(index => estimator[:nb_of_warm_up_samples] for index in index_set)
    end
end

function _regress_nb_of_samples(estimator::Estimator{<:AbstractIndexSet, <:MC}, index_set, ϵ::Real, θ::Real)
    vars = Dict(index => regress_var(estimator, index) for index in index_set)
    costs = Dict(index => regress_cost(estimator, index) for index in index_set)
    Σ_estimate = Σ(estimator)
    for index in keys(vars)
        Σ_estimate += sqrt(vars[index] * costs[index])
    end
    Dict(index => begin
             n_opt = optimal_nb_of_samples(ϵ, θ, vars[index], costs[index], Σ_estimate)
             max(2, min(n_opt, estimator[:nb_of_warm_up_samples]))
         end for index in index_set)
end

function compute_splitting(estimator::Estimator, ϵ::Real)
    bias_estimate = bias(estimator, max_sz(estimator))
    a = estimator[:min_splitting]
    b = estimator[:max_splitting]

    splitting = 1 - bias_estimate^2/ϵ^2
    isnan(splitting) ? a : min(b, max(a, splitting))
end

Σ(estimator::Estimator) = sum(sqrt.(map(index -> var(estimator, index) * cost(estimator, index), keys(estimator))))

optimal_nb_of_samples(estimator::Estimator, index::Index, ϵ::Real, θ::Real) = optimal_nb_of_samples(ϵ, θ, var(estimator, index), cost(estimator, index), Σ(estimator))

optimal_nb_of_samples(ϵ::Real, θ::Real, var_estimate::Real, cost_estimate::Real, Σ_estimate::Real) = ceil(Int, 1/(θ*ϵ^2) * sqrt(var_estimate/cost_estimate) * Σ_estimate)

#
# bias computation
#
boundary(estimator::Estimator, cntr::Integer) = setdiff(get_index_set(estimator, cntr), get_index_set(estimator, cntr-1))

new_index_set(estimator::Estimator, cntr::Integer) = boundary(estimator, cntr)

bias(estimator::Estimator{<:SL}) = 0.0

bias(estimator::Estimator) = bias(estimator, sz(estimator))

function bias(estimator::Estimator, sz::Integer)
    if !isempty(boundary(estimator, sz + 1) ∩ keys(estimator)) && !estimator[:robustify_bias_estimate]
        return abs(sum(broadcast(i -> mean(estimator, i), boundary(estimator, sz + 1))))
    else
        x = max(1, sz - 2):sz
        y = Float64[log2(abs(sum(broadcast(i -> mean(estimator, i), boundary(estimator, xᵢ))))) for xᵢ in x]
        idcs = isfinite.(y)
        p = interp1(x[idcs], y[idcs])
        return 2^(p[1]+(sz+1)*p[2])
    end
end

#
# adaptivity
#
profit(estimator::Estimator{<:AD}, index::Index) = abs(mean(estimator, index)) / sqrt(var(estimator, index) * cost(estimator, index))^estimator[:penalization]

max_level_exceeded(estimator::Estimator{<:AD}) = isempty(setdiff(get_index_set(estimator[:max_search_space], estimator[:max_index_set_param]), keys(estimator)))

function find_index_with_max_profit(estimator::Estimator{<:AD})
    indices = collect(active_set(estimator))
    profits = [profit(estimator, index) for index in indices]
    (max_profit, idx) = findmax(profits)
    if rand() < estimator[:acceptance_rate] || length(profits)==1
        max_index = indices[idx]
    else
        max_index = indices[rand(deleteat!(collect(1:length(profits)), idx))]
    end
    estimator[:verbose] && print_largest_profit(estimator, max_index, max_profit, indices, profits)
    return max_index
end

function new_index_set(estimator::Estimator{<:AD}, sz::Integer)
    d = ndims(estimator)
    if isempty(active_set(estimator))
        new_indices = Set{Index{d}}()
        max_index = Index(ntuple(i -> 0, d))
        add_to_active_set(estimator, max_index)
        push!(new_indices, max_index)
    else
        max_index = find_index_with_max_profit(estimator)
        add_to_old_set(estimator, max_index)
        remove_from_active_set(estimator, max_index)
        new_indices = Set{Index{d}}()
        push!(new_indices, max_index)
        for k in 1:d
            new_index = max_index + Index(ntuple(i -> i == k, d))
            if is_admissable(estimator, new_index)
                if new_index ∈ get_index_set(estimator[:max_search_space], estimator[:max_index_set_param])
                    add_to_active_set(estimator, new_index)
                    push!(new_indices, new_index)
                else
                    warn_max_index(estimator, new_index)
                    add_to_max_index_set(estimator, new_index)
                end
            end
        end
    end
    log_adaptive_index_set(estimator, max_index)
    return new_indices
end

function bias(estimator::Estimator{<:AD}, sz::Integer)
    b(indices) = abs(sum(broadcast(i -> mean(estimator, i), collect(indices))))
    if sz != max_sz(estimator)
        return b(active_set(estimator))
    else
        return min(b(active_set(estimator)), b(boundary(estimator)))
    end
end

#
# QMC related functions
#
qoi_with_max_var(estimator::Estimator{<:AbstractIndexSet, <:QMC}) = estimator[:qoi_with_max_var] == 0 ? argmax(map(n_qoi -> sum(mean(var(samples_diff(estimator, n_qoi, n_shift, index)) for n_shift in 1:estimator[:nb_of_shifts](index)) for index in keys(estimator)), 1:estimator[:nb_of_qoi])) : estimator[:qoi_with_max_var]

function next_number_of_samples(estimator, index)
    if estimator[:sample_mul_factor] == 2
        Dict(index => nextpow(2, nb_of_samples(estimator, index) + 1))
    elseif estimator[:sample_mul_factor] ≤ 1
        Dict(index => nb_of_samples(estimator, index) + 1)
    else
        Dict(index => ceil(Int, nb_of_samples(estimator, index) * estimator[:sample_mul_factor]))
    end
end

function find_index_with_max_var_over_cost(estimator::Estimator{<:AbstractIndexSet, <:QMC})
    indices = collect(keys(estimator))
    vars = [varest(estimator, index) / cost(estimator, index) for index in indices]
    (max_var, idx) = findmax(vars)
    indices[idx]
end

varest(estimator::Estimator{<:AbstractIndexSet, <:QMC}, index::Index) = var(mean(samples_diff(estimator, qoi_with_max_var(estimator), n_shift, index)) for n_shift in 1:estimator[:nb_of_shifts](index), corrected=true) / estimator[:nb_of_shifts](index)

#
# Unbiased estimation
#
qoi_with_max_var(estimator::Estimator{<:U, <:MC}) = estimator[:qoi_with_max_var] == 0 ? argmax(map(n_qoi -> var(accumulator(estimator, n_qoi)) / length(accumulator(estimator, n_qoi)), 1:estimator[:nb_of_qoi])) : estimator[:qoi_with_max_var]

qoi_with_max_var(estimator::Estimator{<:U, <:QMC}) = estimator[:qoi_with_max_var] == 0 ? argmax(map(n_qoi -> var(mean(accumulator(estimator, n_qoi, n_shift)) for n_shift in 1:size(accumulator(estimator), 2)) / estimator[:nb_of_qoi], 1:estimator[:nb_of_qoi])) : estimator[:qoi_with_max_var]

function next_number_of_samples(estimator::Estimator{<:U})
    n_total = sum(values(nb_of_samples(estimator)))
    n0 = max(estimator[:nb_of_warm_up_samples], n_total)
    n = estimator[:sample_mul_factor] ≤ 1 ? 1 : ceil(Int, n0 * (estimator[:sample_mul_factor] - 1))
    sample(pmf(estimator), n)
end

function sample(pmf::Dict{Index{d}, Float64}, n::Integer) where d
    u = rand(n)
    x = Vector{keytype(pmf)}(undef, n)
    psum = 0.
    for (index, p) in pmf
        for i in 1:length(x)
            if psum ≤ u[i] ≤ psum + p
                x[i] = index
            end
        end
        psum = psum + p
    end
    Dict(index => sum(x .== Ref(index)) for index in keys(pmf))
end

Geometric(p, k) = (1 - p)^k*p

function update_pmf(estimator::Estimator{<:U})
    f(estimator, index) = sqrt(var(estimator, index)/cost(estimator, index))
    p = interp(f, estimator)
    if !any(isnan.(p)) && all(p[2:end] .< 0)
        for index in keys(estimator)
            if isnan(f(estimator, index)) || isinf(f(estimator, index))
                set_pmf_key(estimator, index, 2^(p[1]+sum(p[2:end].*index.I)))
            else
                set_pmf_key(estimator, index, f(estimator, index))
            end
        end
        normalize!(pmf(estimator))
    end
end

function normalize!(pmf::Dict{<:Index, Float64})
    tot_sum = sum(values(pmf))
    for (key, val) in pmf
        pmf[key] /= tot_sum
    end
end

Prob(estimator::Estimator{<:U}, index::Index) = sum(val for (key, val) in pmf(estimator) if key ≥ index)

Prob(estimator::Estimator{<:U}) = Dict(index => Prob(estimator, index) for index in keys(estimator))

varest(estimator::Estimator{<:U, <:MC}) = var(accumulator(estimator, qoi_with_max_var(estimator))) / length(accumulator(estimator, qoi_with_max_var(estimator)))

varest(estimator::Estimator{<:U, <:QMC}) = var(mean(accumulator(estimator, qoi_with_max_var(estimator), n_shift)) for n_shift in 1:size(accumulator(estimator), 2)) / estimator[:nb_of_qoi]

bias(estimator::Estimator{<:U}) = 0.0

converged(estimator::Estimator{<:U}, ϵ::Real) = mse(estimator) ≤ ϵ^2

mean(estimator::Estimator{<:U}) = mean(accumulator(estimator, qoi_with_max_var(estimator)))

mean0(estimator::Estimator{<:U}) = mean0(accumulator(estimator, qoi_with_max_var(estimator)))

var(estimator::Estimator{<:U}) = var(accumulator(estimator, qoi_with_max_var(estimator)))

var0(estimator::Estimator{<:U}) = var0(accumulator(estimator, qoi_with_max_var(estimator)))


text:
using Revise
using AdFem
using PyPlot
using LinearAlgebra

mmesh = Mesh(10, 10, 0.1)
nu = ones(get_ngauss(mmesh))
u = rand(mmesh.ndof)
f = compute_fem_laplace_term1(u, nu, mmesh)
f1 = compute_fem_laplace_term1(constant(u), nu, mmesh)

K = compute_fem_laplace_matrix1(nu, mmesh)
f0 = K * u
@info norm(f0 - f)

sess = Session(); init(sess)
f1 = run(sess, f1)
@info norm(f0-f1)


text:
"""
    se_hessian(semfit::SemFit; hessian = :finitediff)

Return hessian based standard errors.

# Arguments
- `hessian`: how to compute the hessian. Options are
    - `:analytic`: (only if an analytic hessian for the model can be computed)
    - `:finitediff`: for finite difference approximation
"""
function se_hessian(sem_fit::SemFit; hessian = :finitediff)

    c = H_scaling(sem_fit.model)

    if hessian == :analytic
        par = solution(sem_fit)
        H = zeros(eltype(par), length(par), length(par))
        hessian!(H, sem_fit.model, sem_fit.solution)
    elseif hessian == :finitediff
        H = FiniteDiff.finite_difference_hessian(
                x -> objective!(sem_fit.model, x),
                sem_fit.solution
                )
    elseif hessian == :optimizer
        throw(ArgumentError("standard errors from the optimizer hessian are not implemented yet"))
    elseif hessian == :expected
        throw(ArgumentError("standard errors based on the expected hessian are not implemented yet"))
    else
        throw(ArgumentError("I dont know how to compute `$how` standard-errors"))
    end

    invH = c*inv(H)
    se = sqrt.(diag(invH))

    return se
end

# Addition functions -------------------------------------------------------------
H_scaling(model::AbstractSemSingle) =
    H_scaling(
        model,
        model.observed,
        model.imply,
        model.optimizer,
        model.loss.functions...)

H_scaling(model, obs, imp, optimizer, lossfun::Union{SemML, SemWLS}) =
    2/(n_obs(model)-1)

H_scaling(model, obs, imp, optimizer, lossfun::SemFIML) =
    2/n_obs(model)

H_scaling(model::SemEnsemble) =
    2/n_obs(model)


text:
using ANYmalCRobot
using Test

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


text:
module BinaryTraitsTestX

export CanFly, CannotFly, flytrait, FlyTrait

using BinaryTraits

@trait Fly
@implement CanFly by liftoff()

struct Swan end
@assign Swan with Fly

function diagnose()
    println("__init__: prefix_map = ", BinaryTraits.prefix_map)
    println("__init__: traits_map = ", BinaryTraits.traits_map)
end

function testme()
    diagnose()
    check(Swan)
end

__init__() = diagnose()

end # module


text:
isdefined(Base, :__precompile__) && __precompile__()

module MathProgComplex

using DataStructures, JuMP

# see https://invenia.github.io/Memento.jl/latest/faq/pkg-usage.html#Using-Memento-in-Julia-packages?-1
using Memento
using Compat: @__MODULE__
const LOGGER = getlogger(@__MODULE__)
__init__() = Memento.register(LOGGER)

# message hierarchy: debug < info < warn < error
# do not print efficiency-related warnings
setlevel!(LOGGER, "info")

# print efficiency-related warnings
# setlevel!(LOGGER, "debug")

import Base: ==, !=, <<, >>, isless, isconst, isreal, isnull, isequal
import Base: +, -, *, /, ^, conj, conj!, abs2, norm, real, imag
import Base: show, print, convert, copy, hash, merge
import Base: start, next, done, length, setindex!, getindex, haskey, keys, values, deepcopy

export AbstractPolynomial

abstract type  AbstractPolynomial end

include("utils_internal.jl")

include("variable.jl")
include("exponent.jl")
include("polynomial.jl")

include("constraint.jl")
include("problem.jl")
include("point.jl")

include("cplx2real.jl")
include("evaluate.jl")
include("iterators.jl")
include("utils_Poly.jl")
include("problem_accessors.jl")

## algebra
include(joinpath("algebra", "add_algebra.jl"))
include(joinpath("algebra", "mult_algebra.jl"))
include(joinpath("algebra", "order.jl"))
include(joinpath("algebra", "unaries.jl"))

## export dat
include(joinpath("export_dat", "utils_ampl.jl"))
# include(joinpath("export_dat", "utils_dat_compare.jl"))
include(joinpath("export_dat", "utils_dat_export.jl"))
include(joinpath("export_dat", "utils_dat_import.jl"))

## export JuMP
include(joinpath("export_JuMP", "utils_jump.jl"))

# SDPhierarchy function
include(joinpath("SDPhierarchy", "SDPhierarchy.jl"))

end


text:
using ProblemParser
using Utils

using Lazy
using Underscores

const file = find_input(@__FILE__)
const input = @>> file slurp parse(Lines())
@assert @_ map(length(_), input) |> Set |> length == 1

function part1(input)
    matrix = @_ input |>
        map(map(==('1'), collect(_)), __) |>
        reduce(hcat, __) |>
        transpose
    len = size(matrix, 1)

    gamma =
        @_ [sum(matrix[:,r]) > len ÷ 2 for r in axes(matrix, 2)] |>
        map(curry(+, '0'), __) |>
        String |>
        parse(Int, __; base=2)
    epsilon = gamma ⊻ (2^size(matrix, 2) - 1)
    gamma * epsilon
end
assertequal(part1(input), 693486)

function sortdown(digit1, input)
    input = sort(input)
    pos = 1
    while length(input) > 1
        first1 = @_ findfirst(_[pos] == '1', input)
        num0 = first1 - 1
        if digit1 ⊻ (num0 <= length(input) ÷ 2)
            input = input[begin:num0]
        else
            input = input[first1:end]
        end
        pos += 1
    end
    parse(Int, input[begin]; base=2)
end

function part2(input)
    oxygen = sortdown(true, input)
    co2 = sortdown(false, input)
    oxygen * co2
end
assertequal(part2(input), 3379326)

max_stars_repo_name: CmdQ/AoC


text:
include("IntervalSpace.jl")
include("PolynomialSpace.jl")

## Union

# union_rule dictates how to create a space that both spaces can be converted to
# in this case, it means
function union_rule(s1::PiecewiseSpace{S1},s2::PiecewiseSpace{S2}) where {S1<:Tuple{Vararg{PolynomialSpace}},
                    S2<:Tuple{Vararg{PolynomialSpace}}}
    PiecewiseSpace(map(Space,merge(domain(s1),domain(s2)).domains))
end

function union_rule(s1::PiecewiseSpace{S1},s2::PolynomialSpace) where S1<:Tuple{Vararg{PolynomialSpace}}
    PiecewiseSpace(map(Space,merge(domain(s1),domain(s2)).domains))
end

include("Chebyshev/Chebyshev.jl")
include("Ultraspherical/Ultraspherical.jl")
include("Jacobi/Jacobi.jl")
include("Hermite/Hermite.jl")
include("Laguerre/Laguerre.jl")
include("CurveSpace.jl")

## Heaviside

conversion_rule(sp::HeavisideSpace,sp2::PiecewiseSpace{NTuple{k,PS}}) where {k,PS<:PolynomialSpace} = sp

Conversion(a::HeavisideSpace,b::PiecewiseSpace{NTuple{kk,CC},DD,RR}) where {kk,CC<:PolynomialSpace,DD<:Domain{<:Number},RR<:Real} =
    ConcreteConversion(a,b)
bandwidths(::ConcreteConversion{HS,PiecewiseSpace{NTuple{kk,CC},DD,RR}}) where {HS<:HeavisideSpace,CC<:PolynomialSpace,DD<:Domain{<:Number},RR<:Real,kk} =
    0,0

getindex(C::ConcreteConversion{HS,PiecewiseSpace{NTuple{kk,CC},DD,RR}},k::Integer,j::Integer) where {HS<:HeavisideSpace,CC<:PolynomialSpace,DD<:Domain{<:Number},RR<:Real,kk} =
    k ≤ dimension(domainspace(C)) && j==k ? one(eltype(C)) : zero(eltype(C))


text:
using Images, TestImages, StaticArrays, OffsetArrays

#Lucas-Kanade Optical Flow tests

    # Testing constants
    test_image = "mandrill"
    number_test_pts = 500
    difference = 0.3
    max_error_points_percentage = 10
    max_allowed_error = 0.1
    max_lost_points_percentage = 40

function test_lk(number_test_pts::Int, flow::Array{SVector{2, Float64}, 1}, x_flow::Float64, y_flow::Float64, status::BitArray{1}, err::Array{Float64, 1}, difference::Float64)
    #Testing parameters
    error_pts = 0
    max_err = 0
    total_err = 0
    lost_points = 0

    for i = 1:number_test_pts
        if status[i]
            if abs(flow[i][2] - x_flow) > difference || abs(flow[i][1] - y_flow) > difference
                error_pts += 1
            end
            if err[i] > max_err
                max_err = err[i]
            end
            total_err += err[i]
        else
            lost_points += 1
        end
    end
    return error_pts, max_err, total_err, lost_points
end

@testset "lucas-kanade" begin

    println("Lucas-Kanade TestSet!")

    #Basic translations (Horizontal)
    img1 = Gray{Float64}.(testimage(test_image))
    img2 = similar(img1)
    for i = 1:size(img1)[1]
        for j = 4:size(img1)[2]
            img2[i,j] = img1[i,j-3]
        end
    end

    corners = imcorner(img1, method=shi_tomasi)
    y, x = findn(corners)
    a = map((yi, xi) -> SVector{2}(yi, xi), y, x)

    srand(9876)
    pts = rand(a, (number_test_pts,))

    flow, status, err = optical_flow(img1, img2, LK(pts, [SVector{2}(0.0,0.0)], 11, 4, false, 20))

    error_pts, max_err, total_err, lost_points = test_lk(number_test_pts, flow, 3.0, 0.0, status, err, difference)

    println("Horizontal Motion")
    println("Error Points Percentage = ", (error_pts/length(pts))*100)
    @test ((error_pts/length(pts))*100) < max_error_points_percentage
    println("Maximum Error = ", max_err)
    @test max_err < max_allowed_error
    println("Lost Points Percentage = ", (lost_points/length(pts))*100)
    @test ((lost_points/length(pts))*100) < max_lost_points_percentage

    # Basic translations (Vertical)
    img2 = OffsetArray(img1, 3, 0)

    flow, status, err = optical_flow(img1, img2, LK(pts, [SVector{2}(0.0,0.0)], 11, 4, false, 20))

    error_pts, max_err, total_err, lost_points = test_lk(number_test_pts, flow, 0.0, 3.0, status, err, difference)

    println("Vertical Motion")
    println("Error Points Percentage = ", (error_pts/length(pts))*100)
    @test ((error_pts/length(pts))*100) < max_error_points_percentage
    println("Maximum Error = ", max_err)
    @test max_err < max_allowed_error
    println("Lost Points Percentage = ", (lost_points/length(pts))*100)
    @test ((lost_points/length(pts))*100) < max_lost_points_percentage

    # Basic translations (Both)
    img2 = OffsetArray(img1, 3, 1)

    flow, status, err = optical_flow(img1, img2, LK(pts, [SVector{2}(0.0,0.0)], 11, 4, false, 20))

    error_pts, max_err, total_err, lost_points = test_lk(number_test_pts, flow, 1.0, 3.0, status, err, difference)

    println("Combined Motion")
    println("Error Points Percentage = ", (error_pts/length(pts))*100)
    @test ((error_pts/length(pts))*100) < max_error_points_percentage
    println("Maximum Error = ", max_err)
    @test max_err < max_allowed_error
    println("Lost Points Percentage = ", (lost_points/length(pts))*100)
    @test ((lost_points/length(pts))*100) < max_lost_points_percentage
end

#Farneback Optical Flow tests

    # Testing constants
    test_image = "mandrill"
    difference = 0.3
    max_error_points_percentage = 7.5

function test_fb(dims::SVector{2, Int}, flow::Array{SVector{2, Float64}, 2}, x_flow::Float64, y_flow::Float64, difference::Float64)
    #Testing parameters
    error_pts = 0
    max_err = 0

    for i = 1:dims[2]
        for j = 1:dims[1]
            if abs(flow[j,i][1] - x_flow) > difference || abs(flow[j,i][2] - y_flow) > difference
                    error_pts += 1
            end
            err = sqrt((abs(flow[j,i][1] - x_flow))^2 + (abs(flow[j,i][2] - y_flow))^2)
            if err > max_err
                max_err = err
            end
        end
    end
    return error_pts, max_err
end

@testset "farneback" begin

    println("Farneback TestSet!")

    #Basic translations (Horizontal)
    img1 = Gray{Float64}.(testimage(test_image))
    img2 = similar(img1)
    for i = 1:size(img1)[1]
        for j = 4:size(img1)[2]
            img2[i,j] = img1[i,j-3]
        end
    end

    flow = optical_flow(img1, img2, Farneback(rand(SVector{2,Float64},2,2), 7, 39, 6.0, 11, 1.5, false, true))

    error_pts, max_err = test_fb(SVector{2}(size(img1)), flow, 3.0, 0.0, difference)

    println("Horizontal Motion")
    println("Error Points Percentage = ", (error_pts/(size(img1)[1]*size(img1)[2]))*100)
    @test ((error_pts/(size(img1)[1]*size(img1)[2]))*100) < max_error_points_percentage
    println("Maximum Error = ", max_err)

    # Basic translations (Vertical)
    img2 = similar(img1)
    for i = 4:size(img1)[1]
        for j = 1:size(img1)[2]
            img2[i,j] = img1[i-3,j]
        end
    end

    flow = optical_flow(img1, img2, Farneback(rand(SVector{2,Float64},2,2), 7, 39, 6.0, 11, 1.5, false, true))

    error_pts, max_err = test_fb(SVector{2}(size(img1)), flow, 0.0, 3.0, difference)

    println("Horizontal Motion")
    println("Error Points Percentage = ", (error_pts/(size(img1)[1]*size(img1)[2]))*100)
    @test ((error_pts/(size(img1)[1]*size(img1)[2]))*100) < max_error_points_percentage
    println("Maximum Error = ", max_err)

    # Basic translations (Both)
    img2 = similar(img1)
    for i = 4:size(img1)[1]
        for j = 2:size(img1)[2]
            img2[i,j] = img1[i-3,j-1]
        end
    end

    flow = optical_flow(img1, img2, Farneback(rand(SVector{2,Float64},2,2), 7, 39, 6.0, 11, 1.5, false, true))

    error_pts, max_err = test_fb(SVector{2}(size(img1)), flow, 1.0, 3.0, difference)

    println("Horizontal Motion")
    println("Error Points Percentage = ", (error_pts/(size(img1)[1]*size(img1)[2]))*100)
    @test ((error_pts/(size(img1)[1]*size(img1)[2]))*100) < max_error_points_percentage
    println("Maximum Error = ", max_err)
end


text:
using Printf

function main()

    infile = "xy_hayabusa"
    znum = 5
    dz = 2.0e-2

    xnum, ynum, x, y = read_wing(infile)

    outdir = "grid"
    make_dir(outdir)
    result = "result"
    make_dir(result)
    result = "post_result"
    make_dir(result)
    nodes, xnum_max, ynum_max, znum_max = mk_gird(xnum, ynum, znum, x, y, dz, outdir)
    vecA(nodes, xnum_max, ynum_max, znum_max, outdir)
end

function read_wing(infile)

    fff=[]
    open(infile, "r") do f  # 全行格納
        fff = read(f,String)
    end
    fff = split(fff,"\n",keepempty=false)   # 改行分割(\n)

    for i in 1:length(fff)
        fff[i] = replace(fff[i]," \r" => "")  # 改行削除(\r)
    end

    temp   = split(fff[1]," ")
    filter!(e->e≠"",temp)

    xnum = Int(parse(Float64,temp[1]))
    ynum = Int(parse(Float64,temp[2]))

    x = zeros(xnum, ynum)
    y = zeros(xnum, ynum)
    skip = 1
    fnum = length(fff)-skip
    for i in 1:fnum
        temp = split(fff[i+skip]," ")
        filter!(e->e≠"",temp)

        tx = parse(Int64,temp[1])
        ty = parse(Int64,temp[2])
        x[tx,ty] = parse(Float64,temp[3])
        y[tx,ty] = parse(Float64,temp[4])
    end

    return xnum, ynum, x, y
end

function mk_gird(xnum, ynum, znum, x, y, dz, outdir)
    """
    nodes[i,j,k,l]
    i   : x,r方向の番号
    j   : y,theta方向の番号
    k   : z方向の番号
    l=1 : 点のx座標
    l=2 : 点のy座標

    nodes[1,:]      : x方向境界
    nodes[:,1]      : y方向境界
    nodes[xnum+2,:] : x方向境界
    nodes[:,ynum+2] : y方向境界
    """
    icell = 2
    xnum_max = xnum + 2*icell
    ynum_max = ynum + 2*icell
    znum_max = znum + 2*icell
    nodes = zeros(xnum_max, ynum_max, znum_max, 3)

    for j in 1+icell:ynum_max-icell
        for i in 1+icell:xnum_max-icell
            nodes[i,j,1,1] = x[i-icell,j-icell]
            nodes[i,j,1,2] = y[i-icell,j-icell]
        end
    end

    #=
    仮想セルの作成：現在あるセル境界線を延長して作成
    そのまま延長した際にセルが交差するのを防ぐため、仮想セルの大きさを100^-1倍(n)にする.
    ベクトルで考えれば下記のような計算になる．たぶん
    また、[1,1]等の角にはダミー値が入っている.
    =#
    n = 1.0
    for j in 1:ynum_max
        for l in 1:2
            nodes[2,j,1,l]        =  (1.0+n)*nodes[3,j,1,l] - n*nodes[4,j,1,l]
            nodes[1,j,1,l]        =  (1.0+n)*nodes[2,j,1,l] - n*nodes[3,j,1,l]
            nodes[xnum_max-1,j,1,l] =  (1.0+n)*nodes[xnum_max-2,j,1,l] - n*nodes[xnum_max-3,j,1,l]
            nodes[xnum_max,j,1,l] =  (1.0+n)*nodes[xnum_max-1,j,1,l] - n*nodes[xnum_max-2,j,1,l]
        end
    end
    for i in 1:xnum_max
        for l in 1:2
            nodes[i,2,1,l]        =  (1.0+n)*nodes[i,3,1,l] - n*nodes[i,4,1,l]
            nodes[i,1,1,l]        =  (1.0+n)*nodes[i,2,1,l] - n*nodes[i,3,1,l]
            nodes[i,ynum_max-1,1,l] =  (1.0+n)*nodes[i,ynum_max-2,1,l] - n*nodes[i,ynum_max-3,1,l]
            nodes[i,ynum_max,1,l] =  (1.0+n)*nodes[i,ynum_max-1,1,l] - n*nodes[i,ynum_max-2,1,l]
        end
    end

    for k in 1:znum_max
        for j in 1:ynum_max
            for i in 1:xnum_max
                nodes[i,j,k,1] = nodes[i,j,1,1]
                nodes[i,j,k,2] = nodes[i,j,1,2]
                nodes[i,j,k,3] = dz*(k-1)
            end
        end
    end

    fff = outdir*"/nodes"
    open(fff,"w") do f
        write(f,"nodes: xnum, ynum, znum, x, y, z\n")
        for i in 1:xnum_max
            for j in 1:ynum_max
                for k in 1:znum_max
                    x = @sprintf("%8.8e", nodes[i,j,k,1])
                    y = @sprintf("%8.8e", nodes[i,j,k,2])
                    z = @sprintf("%8.8e", nodes[i,j,k,3])
                    write(f,string(i)*" "*string(j)*" "*string(k)*" "*x*" "*y*" "*z*"\n")
                end
            end
        end
    end
    println("write "*fff)

    # nodes_forvtk
    nodes_num = zeros(Int, xnum_max, ynum_max, znum_max)
    fff=outdir*"/nodes_forvtk"
    open(fff,"w") do f
        write(f,"nodes: xnum, ynum, znum, x, y, z\n")
        ite=1
        for i in 1+icell:xnum_max-icell
            for j in 1+icell:ynum_max-icell
                for k in 1+icell:znum_max-icell
                x = @sprintf("%8.8e", nodes[i,j,k,1])
                y = @sprintf("%8.8e", nodes[i,j,k,2])
                z = @sprintf("%8.8e", nodes[i,j,k,3])
                write(f,string(k)*" "*x*" "*y*" "*z*"\n")
                nodes_num[i,j,k] = ite
                ite = ite+1
                end
            end
        end
    end
    println("write "*fff)

    fff=outdir*"/nodesnum"
    open(fff,"w") do f
        write(f,"nodesnum: xnum_max, ynum_max, znum_max\n")
        write(f,string(xnum_max) * " "*string(ynum_max) * " "*string(znum_max)*"\n")
    end

    ### element ###
    fff=outdir*"/element_forvtk"
    open(fff,"w") do f
        write(f,"elements:cell_xnum, lup,rup,ldown,rdown, lup,rup,ldown,rdown \n")
        k=1
        for i in 1+icell:xnum_max-icell-1
            for j in 1+icell:ynum_max-icell-1
                for k in 1+icell:znum_max-icell-1
                    d1 = @sprintf("%1.0f", nodes_num[  i,  j,  k])
                    d2 = @sprintf("%1.0f", nodes_num[  i,j+1,  k])
                    d3 = @sprintf("%1.0f", nodes_num[i+1,j+1,  k])
                    d4 = @sprintf("%1.0f", nodes_num[i+1,  j,  k])
                    d5 = @sprintf("%1.0f", nodes_num[  i,  j,k+1])
                    d6 = @sprintf("%1.0f", nodes_num[  i,j+1,k+1])
                    d7 = @sprintf("%1.0f", nodes_num[i+1,j+1,k+1])
                    d8 = @sprintf("%1.0f", nodes_num[i+1,  j,k+1])
                    write(f, string(k)*" "*d1*" "*d2*" "*d3*" "*d4)
                    write(f, " "*d5*" "*d6*" "*d7*" "*d8*"\n")
                    k = k+1
                end
            end
        end
    end
    println("write "*fff)

    return  nodes, xnum_max, ynum_max, znum_max
end

function vecA(nodes,xnum_max,ynum_max,znum_max,outdir)

    vecAx = zeros(xnum_max, ynum_max-1, znum_max-1, 3)
    for k in 1:znum_max-1
        for j in 1:ynum_max-1
            for i in 1:xnum_max
                # 3dim Ax=(, , 0)
                a1 = nodes[  i,j+1,k+1,1] - nodes[  i,  j,  k,1]
                a2 = nodes[  i,j+1,k+1,2] - nodes[  i,  j,  k,2]
                a3 = nodes[  i,j+1,k+1,3] - nodes[  i,  j,  k,3]
                b1 = nodes[  i,  j,k+1,1] - nodes[  i,j+1,  k,1]
                b2 = nodes[  i,  j,k+1,2] - nodes[  i,j+1,  k,2]
                b3 = nodes[  i,  j,k+1,3] - nodes[  i,j+1,  k,3]
                vecAx[i,j,k,1] = 0.5*(a2*b3 - a3*b2)
                vecAx[i,j,k,2] = 0.5*(a3*b1 - a1*b3)
                vecAx[i,j,k,3] = 0.5*(a1*b2 - a2*b1)
            end
        end
    end

    vecAy = zeros(xnum_max-1, ynum_max, znum_max-1, 3)
    for k in 1:znum_max-1
        for j in 1:ynum_max
            for i in 1:xnum_max-1
                a1 = nodes[i+1,  j,k+1,1] - nodes[  i,  j,  k,1]
                a2 = nodes[i+1,  j,k+1,2] - nodes[  i,  j,  k,2]
                a3 = nodes[i+1,  j,k+1,3] - nodes[  i,  j,  k,3]
                b1 = nodes[i+1,  j,  k,1] - nodes[  i,  j,k+1,1]
                b2 = nodes[i+1,  j,  k,2] - nodes[  i,  j,k+1,2]
                b3 = nodes[i+1,  j,  k,3] - nodes[  i,  j,k+1,3]
                vecAy[i,j,k,1] = 0.5*(a2*b3 - a3*b2)
                vecAy[i,j,k,2] = 0.5*(a3*b1 - a1*b3)
                vecAy[i,j,k,3] = 0.5*(a1*b2 - a2*b1)
            end
        end
    end

    vecAz = zeros(xnum_max-1, ynum_max-1, znum_max, 3)
    for k in 1:znum_max
        for j in 1:ynum_max-1
            for i in 1:xnum_max-1
                a1 = nodes[i+1,j+1,  k,1] - nodes[  i,  j,  k,1]
                a2 = nodes[i+1,j+1,  k,2] - nodes[  i,  j,  k,2]
                a3 = nodes[i+1,j+1,  k,3] - nodes[  i,  j,  k,3]
                b1 = nodes[  i,j+1,  k,1] - nodes[i+1,  j,  k,1]
                b2 = nodes[  i,j+1,  k,2] - nodes[i+1,  j,  k,2]
                b3 = nodes[  i,j+1,  k,3] - nodes[i+1,  j,  k,3]
                vecAz[i,j,k,1] = 0.5*(a2*b3 - a3*b2)
                vecAz[i,j,k,2] = 0.5*(a3*b1 - a1*b3)
                vecAz[i,j,k,3] = 0.5*(a1*b2 - a2*b1)
            end
        end
    end

    fff = outdir*"/vecAx"
    open(fff,"w") do f
        write(f,"vecAx: xnum, ynum, znum, x vec, y vec, z vec\n")
        for i in 1:xnum_max
            for j in 1:ynum_max-1
                for k in 1:znum_max-1
                    x = @sprintf("%8.8f", vecAx[i,j,k,1])
                    y = @sprintf("%8.8f", vecAx[i,j,k,2])
                    z = @sprintf("%8.8f", vecAx[i,j,k,3])
                    write(f, string(i)*" "*string(j)*" "*string(k)*" "*x*" "*y*" "*z*"\n")
                end
            end
        end
    end
    println("write "*fff)

    fff = outdir*"/vecAy"
    open(fff,"w") do f
        write(f,"vecAy: xnum, ynum, znum, x vec, y vec, z vec\n")
        for i in 1:xnum_max-1
            for j in 1:ynum_max
                for k in 1:znum_max-1
                    x = @sprintf("%8.8f", vecAy[i,j,k,1])
                    y = @sprintf("%8.8f", vecAy[i,j,k,2])
                    z = @sprintf("%8.8f", vecAy[i,j,k,3])
                    write(f, string(i)*" "*string(j)*" "*string(k)*" "*x*" "*y*" "*z*"\n")
                end
            end
        end
    end
    println("write "*fff)

    fff = outdir*"/vecAz"
    open(fff,"w") do f
        write(f,"vecAz: xnum, ynum, znum, x vec, y vec, z vec\n")
        for i in 1:xnum_max-1
            for j in 1:ynum_max-1
                for k in 1:znum_max
                    x = @sprintf("%8.8f", vecAz[i,j,k,1])
                    y = @sprintf("%8.8f", vecAz[i,j,k,2])
                    z = @sprintf("%8.8f", vecAz[i,j,k,3])
                    write(f, string(i)*" "*string(j)*" "*string(k)*" "*x*" "*y*" "*z*"\n")
                end
            end
        end
    end
    println("write "*fff)
end

function make_dir(outdir)
    k = 0
    try rm(outdir,recursive=true)
    catch
        mkdir(outdir)
        k = 1
    end

    if k == 0
        mkdir(outdir)
    end
end

# ---------------------------------
main()


text:
"""
mapをparalellで実行してidを各プロセスで表示する
"""

using Base.Iterators
using Distributed

@everywhere function miseru(x, y)
    @show x y myid()
    sleep(1)
    return x^2 + y
end

function main()
    N::Int64 = 8
    # tuple -> list
    vmiseru = pmap(x -> miseru(x...), product(1:N, 1:N))
    @show vmiseru
    @show size(vmiseru)
    @show typeof(vmiseru)
end

main()

max_stars_repo_name: Soyukke/julia-paralell-examples

max_issues_repo_name: Soyukke/julia-paralell-examples

max_forks_repo_name: Soyukke/julia-paralell-examples