text:
module LE11NAFEMS_examples
using FinEtools

function LE11NAFEMS_H20()
    ## Solid cylinder/taper/sphere—-temperature loading; quadratic brick mesh
    #

    ## Description
    #
    # The solid cylinder/taper/sphere axially-symmetric part represented in
    # Figure 1 is exposed to linearly varying temperature in the plane of the
    # cross-section. The temperature in the coordinates $r$  (the coordinate)
    # and $z$ (the axial ccoordinate)  is given as $T=r+z$. The goal is to find
    # the mechanical stress at the point A induced by the thermal expansion.
    #

    ##
    # The part is constrained against axial expansion along the faces of HIH'I'
    # and ABA'B'. The Young's modulus is 210 GPa, the Poisson's ratio is .3,
    # and the coefficient of thermal expansion is 2.3e-4/degree Celsius.

    ##
    # This is a test recommended by the National Agency for Finite Element
    # Methods and Standards (U.K.): Test LE11 from NAFEMS Publication TNSB,
    # Rev. 3, “The Standard NAFEMS Benchmarks,” October 1990.
    #
    # Target solution: Compressive  axial stress $\sigma_z$  = –105 MPa along
    # the circle passing through point A.

    ##
    # The toolkit has a helpful physical-units facility.  The function phun()
    # allows use of basic  units and basic
    # multipliers (for instance, mega).

    ##
    # Set the material properties.
    Ea = 210000*phun("MEGA*PA");# Young's modulus
    nua = 0.3;# Poisson ratio
    alphaa = 2.3e-4;# coefficient of thermal expansion

    ##
    # This is the target stress value.
    sigmaA = -105*phun("MEGA*PA");

    ##
    # The mesh  will be created in a very coarse representation from the
    # key points in the drawing. The first coordinate is radial, the second coordinate is axial.
    rz=[1.     0.;#A
    1.4    0.;#B
    0.995184726672197   0.098017140329561;
    1.393258617341076 0.137223996461385;
    0.980785  0.195090;#
    1.37309939 0.27312645;
    0.956940335732209   0.290284677254462
    1.339716470025092 0.406398548156247
    0.9238795  0.38268;#C
    1.2124  0.7;#D
    0.7071  0.7071;#E
    1.1062  1.045;#F
    0.7071  (0.7071+1.79)/2;#(E+H)/2
    1.      1.39;#G
    0.7071  1.79;#H
    1.      1.79;#I
    ]*phun("M")
    tolerance =1.e-6*phun("M")

    ##
    # Note that the material object needs to be created with the proper
    # model-dimension reduction in mind.  In this case that is the fully three-dimensional solid.
    MR = DeforModelRed3D

    # This is the quadrilateral mesh of the cross-section.   It will be modified and
    # refined as  we go.
    fens = FENodeSet(rz);
    fes = FESetQ4([1 2 4 3; 3 4 6 5; 5 6 8 7; 7 8 10 9; 9 10 12 11; 11 12 14 13; 13 14 16 15]);

    ##
    # If needed, the initial mesh  can be refined by bisection.  Just set
    # `nref` greater than zero.  Note that  the nodes located along the
    # edges are moved onto the  spherical surface when they _should be_ on
    # the spherical surface.  This is important in order to ensure
    # convergence to the proper value of the stress.  Just refining  the
    # initial mesh without repositioning of the nodes onto the spherical surface would mean that the
    # refinement would preserve a concave corner where in reality there is
    # none.  The stress would be artificially raised and convergence would
    # not be guaranteed.

    nref = 0;
    for ref = 1:nref
        fens,fes = Q4refine(fens,fes);
        list = selectnode(fens, distance=1.0+0.1/2^nref, from=[0. 0.], inflate=tolerance);
        fens.xyz[list,:] = FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:],1.0);
    end

    ##
    # The mesh is extruded by sweeping around the axis of symmetry.
    # Only a single layer of elements is generated of internal angle
    # |angslice|.
    nLayers = 7;
    angslice  = 5*pi/16;

    ##
    # First the mesh is extruded to a block whose third dimension
    # represents the angular coordinate.
    fens,fes = H8extrudeQ4(fens, fes, nLayers,
    (rz,k)->[rz[1],rz[2],0.0]-(k)/nLayers*[0.,0.,angslice]);

    ##
    # The mesh is now converted to the serendipity 20-node elements.
    # We will reposition the nodes later.
    fens,fes = H8toH20(fens,fes);

    ##
    # The boundary of the block is extracted and the faces of the mesh on
    # the bounding cross-sections are identified. Recall that this is just
    # about the topology (connectivity), the geometry does not matter at
    # this point.
    bfes = meshboundary(fes);
    f1l = selectelem(fens, bfes, box=[-Inf,Inf,-Inf,Inf,0.0,0.0], inflate=tolerance);
    f2l = selectelem(fens, bfes, box=[-Inf,Inf,-Inf,Inf,-angslice,-angslice],
    inflate=tolerance);

    ##
    # The block is now converted  to the axially symmetric geometry by using the
    # third (angular) coordinate  to sweep out  an axially symmetric domain. The
    # ccoordinates of the nodes at this point are |rza|,  radial distance,
    # Z-coordinate, angle.
    sweep(rza) = [-rza[1]*sin(rza[3]+angslice/2.0), rza[1]*cos(rza[3]+angslice/2.0), rza[2]]
    for j=1:size(fens.xyz,1)
        fens.xyz[j,:] = sweep(fens.xyz[j,:])
    end

    ##
    # The nodes within the radial distance of 1.0 of the origin (i. e.
    # those on the spherical surface)  are repositioned one more time to be
    # located on the spherical surface for sure. (Recall  that we have
    # inserted additional nodes at the midpoints of the edges when the mesh
    # was converted to quadratic elements.)
    list = selectnode(fens,distance=1.0+0.1/2^nref,
    from=[0. 0. 0.], inflate=tolerance);
    fens.xyz[list,:]= FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:], 1.0);

    ##
    # We are ready to create the  finite element model machine and to use
    # it to construct  the global system for the displacements.
    ##
    # The material is created from the property object.  Note that the
    # |alpha| attribute is the thermal expansion coefficient.

    # Create isotropic elastic material
    material = MatDeforElastIso(MR, 1.0, Ea, nua, alphaa)

    ##
    # The finite element  model machine puts together the material, the
    # finite elements,  and the integration rule. The Gauss quadrature with
    # 3x3x3 points  gives good accuracy in this case. Compare it with 2x2x2
    # quadrature to appreciate the difference.

    femm = FEMMDeforLinear(MR, IntegData(fes, GaussRule(3, 3)), material)

    ##
    # The geometry nodal field is created from the node set.   The
    # displacement field is created by cloning the geometry and then
    # zeroing out the nodal parameters.
    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field
    nnodes(geom)
    ##
    # The EBCs are applied  next.  Only the axial (Z) degrees of freedom at
    # the bottom and top are fixed to zero.
    l1 = selectnode(fens, box=[-Inf Inf -Inf Inf 0.0 0.0], inflate=tolerance)
    setebc!(u, l1, true, 3, zeros(size(l1)))
    l1 = selectnode(fens, box=[-Inf Inf -Inf Inf 1.79  1.79], inflate=tolerance)
    setebc!(u, l1, true, 3, zeros(size(l1)))
    applyebc!(u)
    numberdofs!(u)

    ##
    # The restraints of the nodes on the bounding cross-sections in the direction
    # of the normal to the plane of the cross-section  in the
    # circumferential direction are introduced using a penalty formulation.
    # For that purpose we introduce  a finite element model machine for the
    # surface  finite elements on the cross-sections.
    springcoefficient =1.0 / ((abs(sigmaA)/1.0e12)/Ea)
    fl = vcat(f1l, f2l)
    xsfemm = FEMMDeforWinkler(IntegData(subset(bfes,fl), GaussRule(2, 3)))

    ##
    # We create the temperature field using the formula $T=r+z$.
    dT = NodalField(reshape(sqrt.(fens.xyz[:,1].^2+fens.xyz[:,2].^2)+fens.xyz[:,3],size(fens.xyz,1),1));

    ##
    # And we are ready to assemble the system matrix. Both the elastic stiffness of
    # the hexahedral elements ...
    K = stiffness(femm, geom, u)
    # ...  and the elastic stiffness    of the springs on the contact surfaces of the cross-sections.
    H = surfacenormalspringstiffness(xsfemm,  geom, u, springcoefficient)

    ##
    # The mechanical loads are computed from the thermal strains.
    F = thermalstrainloads(femm, geom, u, dT)

    ##
    # And  the solution for the free degrees of freedom is obtained.
    U=  (K+H)\F
    scattersysvec!(u, U[:])

    ##
    # The stress  is recovered from the stress calculated at the
    # integration points.

    fld= fieldfromintegpoints(femm, geom, u, dT, :Cauchy, 3)

    ##
    # Now that we have the nodal field  for the axial stress, we can plot
    # the axial stress painted on the deformed geometry.

    File =  "LE11NAFEMS_H20_sigmaz.vtk"
    vtkexportmesh(File, fens, fes;
    scalars=[("sigmaz", fld.values)], vectors=[("u", u.values)])
    @async run(`"paraview.exe" $File`)
    # File =  "LE11NAFEMS_H20_dT.vtk"
    # vtkexportmesh(File, fens, fes; scalars=dT.values,scalars_name ="dT", vectors=u.values,vectors_name="u")

    ##
    # The  computed stress at the node that is located at the point A  is
    # going to be now extracted from the nodal field for the stress.
    # Nodes at level Z=0.0
    l1 =selectnode(fens,box=FFlt[-Inf  Inf -Inf  Inf 0.0 0.0],inflate=tolerance);
    l2 =selectnode(fens,distance=1.0+0.1/2^nref,from=FFlt[0.0 0.0 0.0],inflate=tolerance);
    nA=intersect(l1,l2);
    sA = mean(fld.values[nA])/phun("MEGA*Pa")
    sAn = mean(fld.values[nA])/sigmaA
    println("Stress at point A: $(sA) i. e.  $( sAn*100  )% of reference value")

    ## Discussion
    #
    ##
    # The 3-D solution corresponds well to the 2-D axially symmetric model.
    # We also see good correspondence to other published solutions for
    # comparable finite element models.  For instance, Abaqus 6.11
    # Benchmark manual lists very similar numbers.

end # LE11NAFEMS_H20

function allrun()
    println("#####################################################")
    println("# LE11NAFEMS_H20 ")
    LE11NAFEMS_H20()
    return true
end # function allrun

end # module LE11NAFEMS_examples


text:
using DashBootstrapComponents, DashHtmlComponents

button = html_div([
    dbc_button(
        "Download",
        href = "/static/data_file.txt",
        download = "my_data.txt",
        external_link = true,
        color = "primary",
    ),
])

ext: jl

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
@testset "Constructors for core methods" begin
    @testset "From Float vector" begin
        v = [1.0, 0]
        @testset "$T" for T in OPT_TYPES
            @test T(v) isa T
        end
    end

    @testset "From Int vector" begin
        v = [1, 0]
        @testset "$T" for T in OPT_TYPES
            @test T(v) isa T
        end
    end

    @testset "From Rational vector" begin
        v = [1//1, 0]
        @testset "$T" for T in OPT_TYPES
            @test T(v) isa T
        end
    end
end


text:
using Gadfly, Distributions

set_default_plot_size(7inch, 3inch)

x = sort!(rand(Poisson(20), 1000))
p1 = plot(x=x, Scale.x_discrete, Geom.histogram)
p2 = plot(x=x, Scale.x_discrete, Geom.histogram(density=true))
hstack(p1, p2)


text:
# using Base:splat
using JSON
using SQLite

splat_db = SQLite.DB("../splatalogue_v3.db")
SQLite.tables(splat_db)

freq_start = 688.2213591803346
freq_end = 689.3824263880999

strSQL = "SELECT * FROM lines WHERE frequency>=$freq_start AND frequency<=$freq_end;"
println(strSQL)

has_molecules = false
resp = IOBuffer()
write(resp, "{\"molecules\" : [")

for row in SQLite.DBInterface.execute(splat_db, strSQL)
    global has_molecules = true
    json = JSON.json(row)

    write(resp, json)
    write(resp, ",")
end

json = String(take!(resp))

if !has_molecules
    json = "{\"molecules\" : []}"
else
    # remove the last character (comma) from json
    json = chop(json, tail=1) * "]}"
end

println(json)


text:
Ints = Union{Array{Int8, 1}, Array{Int16, 1}, Array{Int32, 1}, Array{Int64, 1}, Array{Int128, 1}, Array{UInt8, 1}, Array{UInt16, 1}, Array{UInt32, 1}, Array{UInt64, 1}, Array{UInt128, 1}}

function count_sort!(data::Ints, output::Ints, c::Array{Int32, 1}, bits::Int32, offset::Int32)
	combs = 2^bits

	for i in 1:combs
		c[i] = 0
	end

	for i in data
		c[(i >> offset) & (combs - 1) + 1] += 1
	end

	for i in 2:combs
		c[i] += c[i - 1]
	end

	index = 0

	for i in length(data) : -1 : 1
		index = (data[i] >> offset) & (combs - 1) + 1
		output[c[index]] = data[i]
		c[index] -= 1
	end
end

function count_sort_sign!(data::Ints, output::Ints, c::Array{Int32, 1})
	c[1] = 0
	c[2] = 0

	bit = sizeof(data[1]) * 8

	for i in data
		c[~(i >> (bit - 1)) & 1 + 1] += 1
	end

	c[2] += c[1]

	index = 0

	for i in length(data) : -1 : 1
		index = ~(data[i] >> (bit - 1)) & 1 + 1
		output[c[index]] = data[i]
		c[index] -= 1
	end
end

function radix_sort!(data::Ints, log_bits::Int=3)
	if length(data) == 0
		return
	end
	T = typeof(data[1])
	mirror = zeros(T, length(data))
	bits::Int32 = 2^log_bits
	count = zeros(Int32, 2^bits)
	for i in 1:fld(sizeof(T) * 8, bits)
		count_sort!(data, mirror, count, bits, Int32(((i - 1) * 2) * bits))
		count_sort!(mirror, data, count, bits, Int32(((i - 1) * 2 + 1) * bits))
	end
	if !(T <: Unsigned)
		count_sort_sign!(data, mirror, count)
		for i in 1:length(data)
			data[i] = mirror[i]
		end
	end
end

function radix_sort!(data::Array{Float16, 1}, log_bits::Int)
	intdata::Array{Int16, 1} = unsafe_wrap(Array, convert(Ptr{Int16}, pointer(data)), length(data))

	for i in 1:length(intdata)
		if intdata[i] < 0
			intdata[i] = ~intdata[i] + Int16(-32768)
		end
	end

	radix_sort!(intdata, log_bits)

	for i in 1:length(intdata)
		if intdata[i] < 0
			intdata[i] = ~intdata[i] + Int16(-32768)
		end
	end
end

function radix_sort!(data::Array{Float32, 1}, log_bits::Int)
	intdata::Array{Int32, 1} = unsafe_wrap(Array, convert(Ptr{Int32}, pointer(data)), length(data))

	for i in 1:length(intdata)
		if intdata[i] < 0
			intdata[i] = ~intdata[i] + Int32(-2147483648)
		end
	end

	radix_sort!(intdata, log_bits)

	for i in 1:length(intdata)
		if intdata[i] < 0
			intdata[i] = ~intdata[i] + Int32(-2147483648)
		end
	end
end

function radix_sort!(data::Array{Float64, 1}, log_bits::Int)
	intdata::Array{Int64, 1} = unsafe_wrap(Array, convert(Ptr{Int64}, pointer(data)), length(data))

	for i in 1:length(intdata)
		if intdata[i] < 0
			intdata[i] = ~intdata[i] + Int64(-9223372036854775808)
		end
	end

	radix_sort!(intdata, log_bits)

	for i in 1:length(intdata)
		if intdata[i] < 0
			intdata[i] = ~intdata[i] + Int64(-9223372036854775808)
		end
	end
end


text:
using Revise, CSV, DataFrames
using Tables

using MLJ
@load DecisionTreeClassifier
tree_model = DecisionTreeClassifier()

function RandomForestClassifierTS(X, y; n_trees::Int=200, min_interval::Int=3)
    transform_xt = InvFeatureGen(X, n_trees=n_trees, min_interval=min_interval)
    forest = Array{Machine{DecisionTreeClassifier},1}()
    for i in range(1, stop=n_trees)
        tree = machine(tree_model, transform_xt[i], y)
        mdl = deepcopy(tree)
        fit!(mdl)
        push!(forest, mdl)
    end
    forest
end

function InvFeatureGen(X; n_trees::Int=200, min_interval::Int=3)
    n_samps, series_length = size(X)
    transform_xt = Array{Array{Float64,2},1}()
    n_intervals = floor(Int, sqrt(series_length))
    intervals = zeros(Int, n_trees, n_intervals, 2)
    for i in range(1, stop = n_trees)
       transformed_x = Array{Float64,2}(undef, 3*n_intervals, n_samps)
       for j in range(1, stop = n_intervals)
           intervals[i,j,1] = rand(1:(series_length - min_interval))
           len = rand(1:(series_length - intervals[i,j,1]))
           if len < min_interval
               len = min_interval
           end
           intervals[i,j,2] = intervals[i,j,1] + len
           x = Array(1:len+1)
           Y = X[:, intervals[i,j,1]:intervals[i,j,2]]
           means = mean(Y, dims=2)
           stds =  std(Y, dims=2)
           slope = (mean(transpose(x).*Y, dims=2) -
                    mean(x)*mean(Y, dims=2)) / (mean(x.*x) - mean(x)^2)
           transformed_x[3*j-2,:] =  means
           transformed_x[3*j-1,:] =  stds
           transformed_x[3*j,:]   =  slope
       end
           push!(transform_xt, transpose(transformed_x))
    end
    return transform_xt
end

function proba_predict(forest, features)
    n_trees = length(forest)
    n_samps, _ = size(features[1])
    class_ = zeros(Float64, n_samps, 2)
    for i in range(1, stop=n_trees)
        for j in range(1, stop=n_samps)
            class_[j,1] = class_[j,1] + predict(forest[i], features[i])[j].prob_given_class[1]
        end
    end
    class_[:,2] = n_samps .- class_[:,1]
    class_
end

file  = "/Users/aa25desh/TimeSerise/GunPoint/GunPoint_TRAIN.ts"
df = DataFrame(CSV.File(file))
mat = Tables.matrix(df)
X = float.(mat[:, 1:149])
y = [Meta.parse(mat[i, 150]).args[3] for i in range(1, stop=49)]
y = CategoricalArray(y)

forest = RandomForestClassifierTS(X, y)

file1  = "/Users/aa25desh/TimeSerise/GunPoint/GunPoint_TEST.ts"
df1 = DataFrame(CSV.File(file1))
mat1 = Tables.matrix(df1)
X1 = float.(mat1[:, 1:149])
y1 = [Meta.parse(mat1[i, 150]).args[3] for i in range(1, stop=149)]

features = InvFeatureGen(X1)

print(proba_predict(forest, features))


text:
module Rematch

import MacroTools
import MacroTools: @capture
import Compat: occursin, fieldcount

macro splice(iterator, body)
  @assert iterator.head == :call
  @assert iterator.args[1] == :in
  Expr(:..., :(($(esc(body)) for $(esc(iterator.args[2])) in $(esc(iterator.args[3])))))
end

struct MatchFailure
    value
end

function assert_num_fields(::Type{T}, matched::Integer) where {T}
    actual = fieldcount(T)
    @assert actual == matched "Tried to match $matched fields of $actual field struct $T"
end

"""
Statically get the fieldcount of a type. Useful to avoid runtime calls to
fieldcount.
"""
@generated function evaluated_fieldcount(t::Type{T}) where T
    # TODO(nathan.daly) Does this need to be @generated, or is this function even necessary?
    fieldcount(T)
end

"""
    handle_destruct_fields(value, pattern, subpatterns, len, get, bound, asserts; allow_splat=true)

Destruct a tuple, vector, or struct `value` with the given `pattern`.
Returns a boolean expression that evaluates to true if the pattern matches.
Variables bound in the pattern are added to the `bound` set.
Assertions are added to the `asserts` vector.

The value is indexed by `get` from 1 to `len` and element i is matched against `subpattern` i.
If `allow_splat` is true, one `...` is allowed among the subpatterns.
"""
function handle_destruct_fields(
    value::Symbol,
    pattern,
    subpatterns,
    len,
    get::Symbol,
    bound::Set{Symbol},
    asserts::Vector{Expr};
    allow_splat=true
)
    # NOTE we assume `len` is cheap
    fields = []
    seen_splat = false

    for (i,subpattern) in enumerate(subpatterns)
        if (subpattern isa Expr) && (subpattern.head == :(...))
            @assert allow_splat && !seen_splat "Too many ... in pattern $pattern"
            @assert length(subpattern.args) == 1
            seen_splat = true
            push!(fields, (:($i:($len-$(length(subpatterns)-i))), subpattern.args[1]))
        elseif seen_splat
            push!(fields, (:($len-$(length(subpatterns)-i)), subpattern))
        elseif (subpattern isa Expr) && (subpattern.head == :kw)
            field_symbol = Meta.quot(subpattern.args[1])
            push!(fields, (field_symbol, subpattern.args[2]))
        else
            push!(fields, (i, subpattern))
        end
    end

    return Expr(:&&,
        if seen_splat
            :($len >= $(length(subpatterns)-1))
        else
            :($len == $(length(subpatterns)))
        end,
        @splice (i, (field, subpattern)) in enumerate(fields) quote
            $(Symbol("$(value)_$i")) = $get($value, $field)
            $(handle_destruct(Symbol("$(value)_$i"), subpattern, bound, asserts))
        end)
end

"""
handle_destruct(value, pattern, bound, asserts)

Destruct `value` with the given `pattern`.

The pattern is compiled to a boolean expression which evaluates to `true` if the
pattern matches. Variables bound in the pattern are added to the `bound` set.
Assertions are added to the `asserts` vector.
"""
function handle_destruct(value::Symbol, pattern, bound::Set{Symbol}, asserts::Vector{Expr})
    if pattern == :(_)
        # wildcard
        return true
    elseif !(pattern isa Expr || pattern isa Symbol) ||
           pattern == :nothing ||
           @capture(pattern, _quote_macrocall) ||
           @capture(pattern, Symbol(_))
        # constant
        return quote
            $value == $pattern
        end
    elseif (pattern isa Expr && pattern.head == :$)
        # interpolated value
        return quote
            $value == $(esc(pattern.args[1]))
        end
    elseif @capture(pattern, subpattern_Symbol)
        # variable

        # if the pattern doesn't match, we don't want to set the variable
        # so for now just set a temp variable
        our_sym = Symbol("variable_$pattern")

        if pattern in bound
            # already bound, check that this value matches
            return quote
                $our_sym == $value
            end
        else
            # bind
            push!(bound, pattern)
            return quote
                $our_sym = $value;
                true
            end
        end
    elseif @capture(pattern, subpattern1_ || subpattern2_) ||
          (@capture(pattern, f_(subpattern1_, subpattern2_)) && f == :|)
        # disjunction

        # need to only bind variables which exist in both branches
        bound1 = copy(bound)
        bound2 = copy(bound)

        body1 = handle_destruct(value, subpattern1, bound1, asserts)
        body2 = handle_destruct(value, subpattern2, bound2, asserts)
        union!(bound, intersect(bound1, bound2))

        return quote
            $body1 || $body2
        end
    elseif @capture(pattern, subpattern1_ && subpattern2_) ||
          (@capture(pattern, f_(subpattern1_, subpattern2_)) && f == :&)
        # conjunction

        body1 = handle_destruct(value, subpattern1, bound, asserts)
        body2 = handle_destruct(value, subpattern2, bound, asserts)

        return quote
            $body1 && $body2
        end
    elseif @capture(pattern, _where)
        # guard
        @assert length(pattern.args) == 2

        subpattern = pattern.args[1]
        guard = pattern.args[2]

        return quote
            $(handle_destruct(value, subpattern, bound, asserts)) &&
            let $(bound...)
                # bind variables locally so they can be used in the guard
                $(@splice variable in bound quote
                    $(esc(variable)) = $(Symbol("variable_$variable"))
                end)
                $(esc(guard))
            end
        end
    elseif @capture(pattern, ~p_) && @capture(p, f_(subpatterns__))
        # Extractor call.
        # Extractor calls must begin with `~`.

        result = gensym("unapply")
        len = length(subpatterns)

        function replace_call(f)
            if f isa Symbol
                return Symbol("unapply_$f")
            elseif @capture(f, g_(args__))
                h = replace_call(g)
                if h == nothing
                    return nothing
                end
                return Expr(:call, replace_call(g), args...)
            else
                return nothing
            end
        end

        f = replace_call(f)

        if f === nothing
            error("Extractor pattern could not be parsed: should be a call to a named function.")
        end

        # The function should take one argument and return a `Union{Tuple{T,...}, Nothing}`.
        # If the argument does not match, nothing should be returned.
        # If the argument does match, a tuple should be returned, which is then matched against
        # the subpatterns.
        if len == 0
            # If there are no subpatterns, the result value is just checked for
            # nothingness.
            return quote
                $(esc(f))($value) !== nothing
            end
        else
            # If there are subpatterns, the result value is matched against a tuple pattern.
            return quote
                begin
                    $result = $(esc(f))($value)
                    $result !== nothing &&
                    ($result isa Tuple) &&
                    $(handle_destruct_fields(
                        result,
                        pattern,
                        subpatterns,
                        :(length($result)),
                        :getindex,
                        bound,
                        asserts;
                        allow_splat=true
                    ))
                end
            end
        end
    elseif @capture(pattern, T_(subpatterns__))
        # Struct.

        # Since extractors and struct patterns have similar syntax, we require structs to
        # be either symbols or qualified names. We then eval the struct name to check if
        # it is a DataType. Extractors should be Functions instead, and can have arbitrary
        # syntax.

        len = length(subpatterns)

        named_fields = [pat.args[1] for pat in subpatterns
                                    if (pat isa Expr) && pat.head == :kw]
        named_count = length(named_fields)

        @assert named_count == length(unique(named_fields))
            "Pattern $pattern has duplicate named arguments ($(named_fields))."
        @assert named_count == 0 || named_count == len
            "Pattern $pattern mixes named and positional arguments."

        if named_count == 0
            # Pattern uses positional arguments to refer to fields e.g. Foo(0, true)
            expected_fieldcount = gensym("$(T)_expected_fieldcount")
            actual_fieldcount = gensym("$(T)_actual_fieldcount")
            push!(asserts, quote
                t = $(esc(T))
                if typeof(t) <: Function
                    throw(LoadError("Attempted to match on a function",
                                    @__LINE__, AssertionError("Incorrect match usage")))
                end
                if !(isstructtype(typeof(t)) || issabstracttype(typeof(t)))
                    throw(LoadError("Attempted to match on a type that is not a struct",
                                    @__LINE__, AssertionError("Incorrect match usage")))
                end
                # This assertion is necessary:
                # If $expected_fieldcount < $actual_fieldcount, to catch missing fields.
                # If $expected_fieldcount > $actual_fieldcount, to avoid a BoundsError.
                $expected_fieldcount = evaluated_fieldcount(t)
                $actual_fieldcount = $(esc(len))
                if $expected_fieldcount != $actual_fieldcount
                    error("Pattern field count is $($actual_fieldcount), ",
                            "expected $($expected_fieldcount)")
                end

            end)
        else
            # Pattern uses named arguments to refer to fields e.g. Foo(x=0, y=true)
            # Could assert that the expected field names are a subset of the actual
            # field names. Can omit the assertion because if the field doesn't exist
            # getfield() will fail with "type $T has no field $field".
        end

        return quote
            # I would prefer typeof($value) == $(esc(T)) but this doesn't convey type
            # information in Julia 0.6
            $value isa $(esc(T)) &&
            $(handle_destruct_fields(
                value,
                pattern,
                subpatterns,
                length(subpatterns),
                :getfield,
                bound,
                asserts;
                allow_splat=false
            ))
        end
    elseif @capture(pattern, (subpatterns__,))
        # tuple
        return quote
            ($value isa Tuple) &&
            $(handle_destruct_fields(
                value,
                pattern,
                subpatterns,
                :(length($value)),
                :getindex,
                bound,
                asserts;
                allow_splat=true
            ))
        end
    elseif @capture(pattern, [subpatterns__])
        # array
        return quote
            ($value isa AbstractArray) &&
            $(handle_destruct_fields(
                value,
                pattern,
                subpatterns,
                :(length($value)),
                :getindex,
                bound,
                asserts;
                allow_splat=true
            ))
        end
    elseif @capture(pattern, subpattern_::T_)
        # typeassert
        return quote
            ($value isa $(esc(T))) &&
            $(handle_destruct(value, subpattern, bound, asserts))
        end
    else
        error("Unrecognized pattern syntax: $pattern")
    end
end

function handle_match_eq(expr)
    if @capture(expr, pattern_ = value_)
        asserts = Expr[]
        bound = Set{Symbol}()
        body = handle_destruct(:value, pattern, bound, asserts)
        return quote
            $(asserts...)
            value = $(esc(value))
            $body || throw(MatchFailure(value))
            $(@splice variable in bound quote
                $(esc(variable)) = $(Symbol("variable_$variable"))
            end)
            value
        end
    else
        error("Unrecognized match syntax: $expr")
    end
end

function handle_match_case(value, case, tail, asserts)
    if @capture(case, pattern_ => result_)
        bound = Set{Symbol}()
        body = handle_destruct(:value, pattern, bound, asserts)

        return quote
            if $body
                let $(bound...)
                    # export bindings
                    $(@splice variable in bound quote
                      $(esc(variable)) = $(Symbol("variable_$variable"))
                      end)
                    $(esc(result))
                end
            else
                $tail
            end
        end
    else
        error("Unrecognized case syntax: $case")
    end
end

function handle_match_cases(value, match)
    if @capture(match, begin cases__ end)
        tail = :(throw(MatchFailure(value)))
        asserts = Expr[]

        for case in reverse(cases)
            tail = handle_match_case(:value, case, tail, asserts)
        end

        return quote
            $(asserts...)
            value = $(esc(value))
            $tail
        end
    else
        error("Unrecognized match syntax: $value $match")
    end
end

"""
    @match pattern = value

If `value` matches `pattern`, bind variables and return `value`.
Otherwise, throw `MatchFailure`.
"""
macro match(expr)
    handle_match_eq(expr)
end

"""
    @match value begin
        pattern1 => result1
        pattern2 => result2
        ...
    end

Return `result` for the first matching `pattern`.
If there are no matches, throw `MatchFailure`.
"""
macro match(value, cases)
    handle_match_cases(value, cases)
end

"""
Patterns:

  * `_` matches anything
  * `foo` matches anything, binds value to `foo`
  * `~Foo(x,y,z)` calls the extractor function `unapply_Foo(value)` which returns a tuple
    matching `(x,y,z)`
  * `Foo(x,y,z)` matches structs of type `Foo` with fields matching `x,y,z`
  * `Foo(y=1)` matches structs of type `Foo` whose `y` field equals `1`
  * `[x,y,z]` matches `AbstractArray`s with 3 entries matching `x,y,z`
  * `(x,y,z)` matches `Tuple`s with 3 entries matching `x,y,z`
  * `[x,y...,z]` matches `AbstractArray`s with at least 2 entries, where `x` matches
    the first entry, `z` matches the last entry and `y` matches the remaining entries.
  * `(x,y...,z)` matches `Tuple`s with at least 2 entries, where `x` matches the first
    entry, `z` matches the last entry and `y` matches the remaining entries.
  * `_::T` matches any subtype (`isa`) of T
  * `x || y` matches values which match either `x` or `y` (only variables which exist in
    both branches will be bound)
  * `x && y` matches values which match both `x` and `y`
  * `x where condition` matches only if `condition` is true (`condition` may use any
    variables that occur earlier in the pattern eg `(x, y, z where x + y > z)`)
  * Anything else is treated as a constant and tested for equality
  * Expressions can be interpolated in as constants via standard interpolation
    syntax `\$(x)`

Patterns can be nested arbitrarily.

Repeated variables only match if they are equal (`==`).
For example `(x,x)` matches `(1,1)` but not `(1,2)`.

Extractors:

An extractor function must take one argument--the value to be matched against--and should
return either a tuple of values or `nothing`. Returning `nothing` indicates the extractor
does not match. For example to destruct an array into its head and tail:

    function unapply_Cons(xs)
        if isempty(xs)
            return nothing
        else
            return ([xs[1], xs[2:end]])
        end
    end

    @match [1,2,3] begin
        ~Cons(x, xs) => @assert x == 1 && xs == [2,3]
    end

Or here's one to extract the polar coordinates of a point:

    function unapply_Polar(p)
        @match p begin
            (x, y) =>
                begin
                    r = sqrt(x^2+y^2)
                    theta = atan(y, x)
                    return (r, theta)
                end
            _ => return nothing
        end
    end

    @match (1,1) begin
        ~Polar(r, theta) => @assert r == sqrt(2) && theta == pi/4
    end
"""
:(@match)

export @match, @pattern

end


text:
struct Scatterer
    m::ComplexF64
    x::Float64
end

function Scatterer(m::ComplexF64, radius::Float64, wavelength::Float64, n_medium::Float64 = 1.00)
    m /= n_medium
    wavelength /= n_medium  # effective wavelength in the environment

    k = 2π / wavelength  # free space wavenumber in the environment
    x = k .* radius
    return Scatterer(m, x)
end

scattering(s::Scatterer, m) = scattering(s.m, s.x, m)
absorption(s::Scatterer, m) = absorption(s.m, s.x, m)
extinction(s::Scatterer, m) = extinction(s.m, s.x, m)
backscattering(s::Scatterer, m) = backscattering(s.m, s.x, m)

@doc raw"""
mx(m, radius, wavelength, n_medium=1.00)

Conversion of parameters for calculation of light scattering by a small sphere based on Mie theory

# Arguments
- `m`: ``n_{\mathrm{material}}``; refractive index of the material of the sphere
- `radius`: radius of the sphere
- `wavelength`: target vacuum wavelength (range)
- `n_medium`: n_environement; refractive index of the environment
  around the sphere (default: air)

# Return values
- `m`: relative refractive index (``n_\mathrm{material} / n_\mathrm{environment}``)
- `x`: size parameter (wavenumber * radius)
"""
function mx(m, radius, wavelength, n_medium = 1.00)
    m /= n_medium
    wavelength /= n_medium  # effective wavelength in the environment

    k = 2π / wavelength  # free space wavenumber in the environment
    x = k .* radius
    return m, x
end

max_stars_repo_name: fekad/jlmie

max_issues_repo_name: fekad/jlmie

max_forks_repo_name: fekad/jlmie


text:
using MagSim
using Plots;
gr();
using DataFrames
using CSV
using LaTeXStrings

# %%

βs = range(1e-3, stop = 1, length = 201)
Ls = [4, 8, 16, 32, 64];
prefix = "2020-5-14-simulating-magnetism-part-4";

q = 2
@time sims = simulate_Potts_parallel(q, metropolis!, Ls, βs);
df = sim_dict_to_df(sims);
CSV.write(prefix * "_q=$(q).csv", df);

# %%

function plot_data(df, Os, y_limits, y_ticks)
    plots = Dict()
    labels = latexstring.(Os)
    Ls = unique(df.L)

    for (i, O) in enumerate(Os)
        O ∈ [:e] ? legend_bool = true : legend_bool = false
        p = plot(
            dpi = 200,
            xaxis = (latexstring("\\beta")),
            yaxis = (labels[i], y_limits[i], y_ticks[i]),
            legend = legend_bool,
        )

        colors = [:firebrick1, :orange, :forestgreen, :cyan3, :dodgerblue]
        for (i, L) in enumerate(Ls)
            df_ = df[(df.L.==L), :]
            x = df_[:, :β]
            y = df_[:, O]
            plot!(
                x,
                y,
                marker = (:star4, colors[i], 4, stroke(0)),
                line = (colors[i], 2),
                label = latexstring("L = $(string(L))"),
            )
        end

        plots[O] = p
    end

    return plots
end

Os = [:e, :c]
y_limits = [(-2.1, -0.9), (-0.1, 2.6)]
y_ticks = [-2:0.2:-1, 0:0.5:2.5]

plots = plot_data(df, Os, y_limits, y_ticks);
plot_e_c = plot(plots[:e], plots[:c], layout = (2, 1), size = (600, 500));
plot_e_c
savefig(plot_e_c, "../images/" * prefix * "/potts_2_metropolis_e_c.png")


text:
#Almost exact copy of the sparse_int_set in DataStructures.jl
const INT_PER_PAGE = div(ccall(:jl_getpagesize, Clong, ()), sizeof(Int))
# we use this to mark pages not in use, it must never be written to.
const NULL_INT_PAGE = Vector{Int}()

mutable struct Indices
    packed ::Vector{Int}
    reverse::Vector{Vector{Int}}
    counters::Vector{Int}  # counts the number of real elements in each page of reverse.
end

Indices() = Indices(Int[], Vector{Int}[], Int[])

Indices(indices) = union!(Indices(), indices)

Base.eltype(::Type{Indices}) = Int

Base.empty(::Indices) = Indices()

function Base.empty!(s::Indices)
    empty!(s.packed)
    empty!(s.reverse)
    empty!(s.counters)
    return s
end

Base.isempty(s::Indices) = isempty(s.packed)

Base.copy(s::Indices) = copy!(Indices(), s)

function Base.copy!(to::Indices, from::Indices)
    to.packed = copy(from.packed)
    #we want to keep the null pages === NULL_INT_PAGE
    resize!(to.reverse, length(from.reverse))
    for i in eachindex(from.reverse)
        page = from.reverse[i]
        if page === NULL_INT_PAGE
            to.reverse[i] = NULL_INT_PAGE
        else
            to.reverse[i] = copy(from.reverse[i])
        end
    end
    to.counters = copy(from.counters)
    return to
end

Base.lastindex(s::Indices) = s.packed[end]

function pageid_offset(s::Indices, i)
    pageid = div(i - 1, INT_PER_PAGE) + 1
    return pageid, (i - 1) & (INT_PER_PAGE - 1) + 1
end

@inline function Base.in(i, s::Indices)
    pageid, offset = pageid_offset(s, i)
    if pageid > length(s.reverse)
        return false
    else
        page = @inbounds s.reverse[pageid]
        return page !== NULL_INT_PAGE &&  @inbounds page[offset] != 0
    end
end

Base.length(s::Indices) = length(s.packed)

@inline function Base.getindex(s::Indices, i::Integer)
    pageid, offset = pageid_offset(s, i)

    @boundscheck if pageid > length(s.reverse)
        throw(BoundsError(s, i))
    end
    page = @inbounds s.reverse[pageid]

    @boundscheck if page === NULL_INT_PAGE
        throw(BoundsError(s, i))
    end
    id = @inbounds page[offset]
    @boundscheck if id === 0
        throw(BoundsError(s, i))
    end
    return id
end

@inline function Base.push!(s::Indices, i::Integer)
    i <= 0 && throw(DomainError("Only positive Ints allowed."))

    pageid, offset = pageid_offset(s, i)
    pages = s.reverse
    plen = length(pages)

    if pageid > plen
        # Create new null pages up to pageid and fresh (zero-filled) one at pageid
        sizehint!(pages, pageid)
        sizehint!(s.counters, pageid)
        for i in 1:pageid - plen - 1
            push!(pages, NULL_INT_PAGE)
            push!(s.counters, 0)
        end
        push!(pages, zeros(Int, INT_PER_PAGE))
        push!(s.counters, 0)
    elseif pages[pageid] === NULL_INT_PAGE
        #assign a page to previous null page
        pages[pageid] = zeros(Int, INT_PER_PAGE)
    end
    page = pages[pageid]
    if page[offset] == 0
        @inbounds page[offset] = length(s) + 1
        @inbounds s.counters[pageid] += 1
        push!(s.packed, i)
        return s
    end
    return s
end

@inline function Base.push!(s::Indices, is::Integer...)
    for i in is
        push!(s, i)
    end
    return s
end

@inline Base.@propagate_inbounds function Base.pop!(s::Indices)
    if isempty(s)
        throw(ArgumentError("Cannot pop an empty set."))
    end
    id = pop!(s.packed)
    pageid, offset = pageid_offset(s, id)
    @inbounds s.reverse[pageid][offset] = 0
    @inbounds s.counters[pageid] -= 1
    cleanup!(s, pageid)
    return id
end

@inline Base.@propagate_inbounds function Base.pop!(s::Indices, id::Integer)
    id < 0 && throw(ArgumentError("Int to pop needs to be positive."))

    @boundscheck if !in(id, s)
        throw(BoundsError(s, id))
    end
    @inbounds begin
        packed_endid = s.packed[end]
        from_page, from_offset = pageid_offset(s, id)
        to_page, to_offset = pageid_offset(s, packed_endid)

        packed_id = s.reverse[from_page][from_offset]
        s.packed[packed_id] = packed_endid
        s.reverse[to_page][to_offset] = s.reverse[from_page][from_offset]
        s.reverse[from_page][from_offset] = 0
        s.counters[from_page] -= 1
        pop!(s.packed)
    end
    cleanup!(s, from_page)
    return id
end

@inline function cleanup!(s::Indices, pageid::Int)
    if s.counters[pageid] == 0
        s.reverse[pageid] = NULL_INT_PAGE
    end
end

@inline function Base.pop!(s::Indices, id::Integer, default)
    id < 0 && throw(ArgumentError("Int to pop needs to be positive."))
    return in(id, s) ? (@inbounds pop!(s, id)) : default
end
Base.popfirst!(s::Indices) = pop!(s, first(s))

@inline Base.iterate(set::Indices, args...) = iterate(set.packed, args...)

Base.last(s::Indices) = isempty(s) ? throw(ArgumentError("Empty set has no last element.")) : last(s.packed)

Base.union(s::Indices, ns) = union!(copy(s), ns)
function Base.union!(s::Indices, ns)
    for n in ns
        push!(s, n)
    end
    return s
end

Base.intersect(s1::Indices) = copy(s1)
Base.intersect(s1::Indices, ss...) = intersect(s1, intersect(ss...))
function Base.intersect(s1::Indices, ns)
    s = Indices()
    for n in ns
        n in s1 && push!(s, n)
    end
    return s
end

Base.intersect!(s1::Indices, ss...) = intersect!(s1, intersect(ss...))

#Is there a more performant way to do this?
Base.intersect!(s1::Indices, ns) = copy!(s1, intersect(s1, ns))

Base.setdiff(s::Indices, ns) = setdiff!(copy(s), ns)
function Base.setdiff!(s::Indices, ns)
    for n in ns
        pop!(s, n, nothing)
    end
    return s
end

function Base.:(==)(s1::Indices, s2::Indices)
    length(s1) != length(s2) && return false
    return all(in(s1), s2)
end

issubset(a::Indices, b::Indices) = isequal(a, intersect(a, b))

Base.:(<)(a::Indices, b::Indices) = ( a<=b ) && !isequal(a, b)
Base.:(<=)(a::Indices, b::Indices) = issubset(a, b)

function findfirst_packed_id(i, s::Indices)
    pageid, offset = pageid_offset(s, i)
    if pageid > length(s.counters) || s.counters[pageid] == 0
        return 0
    end
    @inbounds id = s.reverse[pageid][offset]
    return id
end

Base.collect(s::Indices) = copy(s.packed)

function Base.permute!(s::Indices, p)
    permute!(s.packed)
    @inbounds for (i, eid) in s.packed
        p[i] == i && continue #nothing was changed
        pageid, offset = pageid_offset(s, eid)
        s.reverse[page][offset] = i
    end
end

struct IndicesIterator{I, T<:Function}
    shortest::I
    test::T
    len::Int
end

Base.length(it::IndicesIterator) = it.len

@inline indices(i::Indices) = i

@inline function Base.iterate(it::IndicesIterator, state=1)
    it_length = length(it)
    for i=state:it_length
        id = indices(it.shortest).packed[i]
        if it.test(id)
            return id, i+1
        end
    end
end

macro indices_in(indices_expr)
    expr, t_sets, t_orsets = expand_indices_bool(indices_expr)
    return esc(quote
        sets = $(Expr(:tuple, t_sets...))
        orsets = $(Expr(:tuple, t_orsets...))

        if isempty(sets)
            minlen, minid = findmin(map(length, orsets))
            t_shortest = orsets[minid]
        else
            minlen, minid = findmin(map(length, sets))
            t_shortest = sets[minid]
        end
        if $(!isempty(t_orsets))
            shortest = deepcopy(t_shortest)
            for s in orsets
                union!(shortest, s)
            end
        else
            shortest = t_shortest
        end
        Overseer.IndicesIterator(shortest, x -> $expr, length(shortest))
    end)
end

function expand_indices_bool(expr)
    if expr isa Symbol || expr.head == :ref || (expr.head == :call && expr.args[1] != :!)
        return Expr(:call, :in, :x, expr), [expr], Symbol[]
    end
    # if !in(expr.head, (:||, :&&)) && !(expr.head == :call && expr.args[1] in == :!)
    #     error("Can only expand expressions with ||, && and !")
    # end
    sets = Union{Symbol, Expr}[]
    orsets = Union{Symbol, Expr}[]
    if expr.args[1] == :!
        nothing
    elseif isa(expr.args[1], Symbol)
        if expr.head != :||
            push!(sets, expr.args[1])
        else
            push!(orsets, expr.args[1])
        end
        expr.args[1] = Expr(:call, :in, :x, expr.args[1])
    else
        expr_, sets_, orsets_ = expand_indices_bool(expr.args[1])
        append!(sets,  sets_)
        append!(orsets,  orsets_)
        expr.args[1] = expr_
    end
    if isa(expr.args[2], Symbol)
        if expr.args[1] != :!
            if expr.head != :||
                push!(sets, expr.args[2])
            else
                push!(orsets, expr.args[2])
            end
        end
        expr.args[2] = Expr(:call, :in, :x, expr.args[2])
    else
        expr_, sets_, orsets_ = expand_indices_bool(expr.args[2])
        append!(sets, sets_)
        append!(orsets, orsets_)
        expr.args[2] = expr_
    end
    return expr, sets, orsets
end

#This is assuming that it's all inbounds
function set_packed_id!(ids::Indices, reverse_id::Int, new_packed_id::Int)
    @inbounds begin
        from_pageid, from_offset = pageid_offset(ids, reverse_id)
        from_packed_id = ids.reverse[from_pageid][from_offset]
        ids.reverse[from_pageid][from_offset] = new_packed_id

        to_pageid, to_offset = pageid_offset(ids, ids.packed[new_packed_id])
        ids.reverse[to_pageid][to_offset] = from_packed_id

        ids.packed[from_packed_id], ids.packed[new_packed_id] =
            ids.packed[new_packed_id], ids.packed[from_packed_id]
    end
end


text:
# Autogenerated wrapper script for LLVMCompilerRT_jll for aarch64-linux-gnu
JLLWrappers.@generate_wrapper_header("LLVMCompilerRT")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
##################################################
## Deflation
## when a Q[k] matrix become a "D" matrix, we deflate. This is checked by the sine term being basically 0.
function check_deflation(state::FactorizationType{T,St,P, Tw}, tol = eps(T)) where {T,St,P,Tw}
    for k in state.ctrs.stop_index:-1:state.ctrs.start_index
        if abs(vals(state.Q[k])[2]) <= tol
            deflate(state, k)
            return
        end
    end
end

# deflate a term
# turn on `show_status` to view sequence
function deflate(state::FactorizationType{T,Val{:DoubleShift},P, Tw}, k) where {T,P,Tw}

    # make a D matrix for Q. C could be +/- 1
    c,s = vals(state.Q[k])
    c = c/norm(c)
    state.Q[k] = RealRotator(c, zero(T), state.Q[k].i)

    # shift zero counter
    state.ctrs.zero_index = k      # points to a matrix Q[k] either RealRotator(-1, 0) or RealRotator(1, 0)
    state.ctrs.start_index = k + 1

    # reset counter
    state.ctrs.it_count = 1
end

# deflate a term
# deflation for ComplexReal is different, as
# we replace Qi with I and move diagonal part into D
function deflate(state::FactorizationType{T,Val{:SingleShift},P, Val{:NotTwisted}}, k) where {T, P}

    # when we deflate here we want to leave Q[k] = I and

    # move Dk matrix over to merge with D
    # Qi       Qi            Qi           Qi
    #   Qj   ->   I Dj  -->    I Dj   -->   I Dj
    #     Qk        QK           Qk Dj        Qk Dj    fuse Dj's with D
    #       QL        QL           QL           Ql Dj

    # then the Dk's are collected into D

    alpha, s = vals(state.Q[k])
    state.Q[k] = Rotator(one(Complex{T}), zero(T), idx(state.Q[k]))

    cascade(state.Q, state.D, alpha, k, state.ctrs.stop_index)

    # shift zero counter
    state.ctrs.zero_index = k      # points to a matrix Q[k] either RealRotator(-1, 0) or RealRotator(1, 0)
    state.ctrs.start_index = k + 1

    # reset counter
    state.ctrs.it_count = 1
end

##################################################


text:
const parconv = 4.57mol*W^-1*s^-1

struct ManualTemperature{T}
    airtemperature::T
    soiltemperature::T
end

apply_environment!(plant, organs, u, t) =
    apply_environment!(plant, environment(plant), organs, u, t)
apply_environment!(plant, ::Nothing, organs, u, t) = nothing
apply_environment!(plant, env::ManualTemperature, organs, u, t) = begin
    envpos = ustrip(calc_envtime(plant, t))
    shoot, root = organs
    update_temp!(shoot, env.airtemperature)
    update_temp!(root, env.soiltemperature)
end

apply_environment!(plant, env::AbstractMicroclimate, organs, u, t) = begin
    envpos = ustrip(calc_envtime(plant, t))
    shoot, root = organs
    # Get the climate at a particular height and time
    # Use half the height for the median temperature/windspeed/humidity etc.
    # relating temperature to the distribution of biomass may be more accurate.
    shootenv = MicroclimInstant(env, height(root)/2, envpos)
    rootenv = MicroclimInstant(env, depth(root)/2, envpos)
    # Update any environment dependent variables
    apply_environment!(assimilation_pars(shoot), shoot, u, shootenv, rootenv)
    apply_environment!(assimilation_pars(root), root, u, rootenv, shootenv)
    nothing
end

apply_environment!(plant, env::MicroclimControl, organs, u, t) = begin
    shoot, root = organs
    apply_environment!(assimilation_pars(shoot), shoot, u, env, env)
    apply_environment!(assimilation_pars(root), root, u, env, env)
    nothing
end

apply_environment!(a::FvCBPhotosynthesis, o, u, shootenv, rootenv) = begin
    v = a.vars

    v.tair = temp = airtemperature(shootenv)
    v.windspeed = windspeed(shootenv)
    v.rh = relhumidity(shootenv)
    v.rnet = radiation(shootenv)
    v.par = radiation(shootenv) * parconv
    # v.soilmoist = mean_soilwatercontent(rootenv)

    update_temp!(o, temp)

    if is_germinated(o, u)
        run_enbal!(assimilation_pars(o).photoparams, v)
    else
        v.tleaf = v.tair
    end
end

apply_environment!(a::AbstractAssim, o, u, shootenv, rootenv) =
    update_temp!(o, airtemperature(shootenv))
apply_environment!(a::AbstractNAssim, o, u, rootenv, shootenv) =
    update_temp!(o, soiltemperature(rootenv))

apply_environment!(a::AbstractKooijmanPhoto, o, u, shootenv, rootenv) = begin
    update_temp!(o, airtemperature(shootenv))
    a.vars.J_L_F = radiation(shootenv) * parconv
    # Use the mean soil water potential between zero and the full root depth
    a.vars.soilwaterpotential = mean_soilwaterpotential(rootenv.microclimate, depth(o), rootenv.t)
end

update_temp!(o, temp) = begin
    set_temp!(o, temp)
    set_tempcorrection!(o, tempcorr(tempcorr_pars(o), temp))
end


text:
function nonLinChem(dy,y,p,t)
 dy[1] = -y[1]
 dy[2] = y[1]-(y[2])^2
 dy[3] = (y[2])^2
end
y0 = [1.0;0.0;0.0]
tspan = (0.0,20)
nlc_analytic(u0,p,t) = [exp(-t);
    (2sqrt(exp(-t))besselk(1,2sqrt(exp(-t)))-2besselk(1,2)/besseli(1,2)*sqrt(exp(-t))besseli(1,2sqrt(exp(-t))))/(2besselk(0,2sqrt(exp(-t)))+(2besselk(1,2)/besseli(1,2))besseli(0,2sqrt(exp(-t))));
    -exp(-t)+1+(-2sqrt(exp(-t))*besselk(1,2sqrt(exp(-t)))+sqrt(exp(-t))*besseli(1,2sqrt(exp(-t)))*2besselk(1,2)/besseli(1,2))/(2besselk(0,2sqrt(exp(-t)))+2besselk(1,2)/besseli(1,2)*besseli(0,2sqrt(exp(-t))))]
nonLinChem_f = ODEFunction(nonLinChem,analytic = nlc_analytic)
prob_ode_nonlinchem = ODEProblem(nonLinChem,y0,tspan)


text:
include("../src/Rand.jl")
include("testing.jl")

print("checking rand interface ")

bytes = Rand.cryptographic_rand(10)
@assert_equal(10, length(bytes))

c = copy(bytes)
@assert_equal(bytes, c)
Rand.cryptographic_rand(c)
@assert_not_equal(bytes, c)
println()

println()


text:
function prune_system_devices!(sys::PSY.System, prune_dict::Dict{Type{<:PSY.Component}, Array{AbstractString}})
    for (type, device_names) in prune_dict
        for name in device_names
            device = PSY.get_component(type, sys, name)
            PSY.remove_component!(sys, device)
        end
    end
    return"""
    THis function removes the pruned devices from PSY System.
    """
    function prune_system_devices!(sys::PSY.System, prune_dict::Dict{Type{<:PSY.Component}, Array{AbstractString}})
        for (type, device_names) in prune_dict
            for name in device_names
                device = PSY.get_component(type, sys, name)
                PSY.remove_component!(sys, device)
            end
        end
        return
    end

    """
    Use this function to specify which devices to remove from the RTS at the start of the simulation.
    """
    function specify_pruned_units()
        pruned_unit = Dict{Type{<:PSY.Component}, Array{AbstractString}}()
        pruned_unit[PSY.ThermalStandard] = ["115_STEAM_1", "115_STEAM_2", "315_STEAM_1", "315_STEAM_2", "315_STEAM_3", "315_STEAM_4", "315_STEAM_5",
                                         "101_CT_1", "101_CT_2", "102_CT_1", "102_CT_2",
                                         "201_CT_1", "201_CT_2", "202_CT_1", "202_CT_2",
                                         "301_CT_1", "301_CT_2", "302_CT_1", "302_CT_2",
                                         "207_CT_1", "307_CT_1", "101_STEAM_4"]
        pruned_unit[PSY.RenewableFix] = ["308_RTPV_1", "313_RTPV_1", "313_RTPV_2", "313_RTPV_3", "313_RTPV_4", "313_RTPV_5", "313_RTPV_6", "313_RTPV_7",
                                        "313_RTPV_8", "313_RTPV_9", "313_RTPV_10", "313_RTPV_11", "313_RTPV_12", "320_RTPV_1", "320_RTPV_2", "320_RTPV_3",
                                        "313_RTPV_13", "320_RTPV_4", "320_RTPV_5", "118_RTPV_1", "118_RTPV_2", "118_RTPV_3", "118_RTPV_4", "118_RTPV_5",
                                        "118_RTPV_6", "320_RTPV_6", "118_RTPV_7", "118_RTPV_8", "118_RTPV_9", "118_RTPV_10", "213_RTPV_1"]
        pruned_unit[PSY.RenewableDispatch] = ["309_WIND_1", "212_CSP_1"]
        pruned_unit[PSY.GenericBattery] = ["313_HEAD_STORAGE"]

        return pruned_unit
    end

    """
    This function creates a PSY System for Unit Commitment.
    """
    function create_rts_sysUC(rts_dir::String,
                            base_power::Float64)
        rawsys1 = PSY.PowerSystemTableData(rts_dir, base_power, joinpath(rts_dir,"../FormattedData/SIIP/user_descriptors.yaml"),
        generator_mapping_file = joinpath(rts_dir,"../FormattedData/SIIP/generator_mapping.yaml"),
        timeseries_metadata_file = joinpath(rts_dir,"../FormattedData/SIIP/timeseries_pointers.json"))
        sys_UC = PSY.System(rawsys1; forecast_resolution = Dates.Hour(1))

        pruned_unit = specify_pruned_units()
        device_list = Dict{Type{<:PSY.Component}, Array{AbstractString}}()
        prune_system_devices!(sys_UC, pruned_unit)

        return sys_UC
    end

    """
    This function creates a PSY System for Economic Dispatch.
    """
    function create_rts_sysED(rts_dir::String,
                            base_power::Float64)
        rawsys1 = PSY.PowerSystemTableData(rts_dir, base_power, joinpath(rts_dir,"../FormattedData/SIIP/user_descriptors.yaml"),
        generator_mapping_file = joinpath(rts_dir,"../FormattedData/SIIP/generator_mapping.yaml"),
        timeseries_metadata_file = joinpath(rts_dir,"../FormattedData/SIIP/timeseries_pointers.json"))
        sys_ED = PSY.System(rawsys1; forecast_resolution = Dates.Hour(1))

        pruned_unit = specify_pruned_units()

        device_list = Dict{Type{<:PSY.Component}, Array{AbstractString}}()

        prune_system_devices!(sys_ED, pruned_unit)

        return sys_ED
    end

end

function specify_pruned_units()
    pruned_unit = Dict{Type{<:PSY.Component}, Array{AbstractString}}()
    pruned_unit[PSY.ThermalStandard] = ["115_STEAM_1", "115_STEAM_2", "315_STEAM_1", "315_STEAM_2", "315_STEAM_3", "315_STEAM_4", "315_STEAM_5",
                                     "101_CT_1", "101_CT_2", "102_CT_1", "102_CT_2",
                                     "201_CT_1", "201_CT_2", "202_CT_1", "202_CT_2",
                                     "301_CT_1", "301_CT_2", "302_CT_1", "302_CT_2",
                                     "207_CT_1", "307_CT_1", "101_STEAM_4"]
    pruned_unit[PSY.RenewableFix] = ["308_RTPV_1", "313_RTPV_1", "313_RTPV_2", "313_RTPV_3", "313_RTPV_4", "313_RTPV_5", "313_RTPV_6", "313_RTPV_7",
                                    "313_RTPV_8", "313_RTPV_9", "313_RTPV_10", "313_RTPV_11", "313_RTPV_12", "320_RTPV_1", "320_RTPV_2", "320_RTPV_3",
                                    "313_RTPV_13", "320_RTPV_4", "320_RTPV_5", "118_RTPV_1", "118_RTPV_2", "118_RTPV_3", "118_RTPV_4", "118_RTPV_5",
                                    "118_RTPV_6", "320_RTPV_6", "118_RTPV_7", "118_RTPV_8", "118_RTPV_9", "118_RTPV_10", "213_RTPV_1"]
    pruned_unit[PSY.RenewableDispatch] = ["309_WIND_1", "212_CSP_1"]
    pruned_unit[PSY.GenericBattery] = ["313_HEAD_STORAGE"]

    return pruned_unit
end

function create_rts_sysUC(rts_dir::String,
                        base_power::Float64)
    rawsys1 = PSY.PowerSystemTableData(rts_dir, base_power, joinpath(rts_dir,"../FormattedData/SIIP/user_descriptors.yaml"),
    generator_mapping_file = joinpath(rts_dir,"../FormattedData/SIIP/generator_mapping.yaml"),
    timeseries_metadata_file = joinpath(rts_dir,"../FormattedData/SIIP/timeseries_pointers.json"))
    sys_UC = PSY.System(rawsys1; forecast_resolution = Dates.Hour(1))

    pruned_unit = specify_pruned_units()
    device_list = Dict{Type{<:PSY.Component}, Array{AbstractString}}()
    prune_system_devices!(sys_UC, pruned_unit)

    return sys_UC
end

function create_rts_sysED(rts_dir::String,
                        base_power::Float64)
    rawsys1 = PSY.PowerSystemTableData(rts_dir, base_power, joinpath(rts_dir,"../FormattedData/SIIP/user_descriptors.yaml"),
    generator_mapping_file = joinpath(rts_dir,"../FormattedData/SIIP/generator_mapping.yaml"),
    timeseries_metadata_file = joinpath(rts_dir,"../FormattedData/SIIP/timeseries_pointers.json"))
    sys_ED = PSY.System(rawsys1; forecast_resolution = Dates.Hour(1))

    pruned_unit = specify_pruned_units()

    device_list = Dict{Type{<:PSY.Component}, Array{AbstractString}}()

    prune_system_devices!(sys_ED, pruned_unit)

    return sys_ED
end


text:
export BlockMap
export remoteIDList, lid, gid
export minAllGID, maxAllGID, minMyGID, maxMyGID, minLID, maxLID
export numGlobalElements, myGlobalElements
export uniqueGIDs, globalIndicesType, sameBlockMapDataAs, sameAs
export linearMap, myGlobalElementIDs
export myGID, myLID, distributedGlobal, numMyElements

# methods and docs based straight off Epetra_BlockMap to match Comm

# ignoring indexBase methods and sticking with 1-based indexing
# ignoring elementSize methods since, type information is carried anyways
# ignoring point-related code, since elementSize is ignored

# TODO figure out expert users and developers only functions

"""
A type for partitioning block element vectors and matrices
"""
struct BlockMap{GID <: Integer, PID <:Integer, LID <: Integer}
    data::BlockMapData{GID, PID, LID}

    function BlockMap{GID, PID, LID}(data::BlockMapData) where {GID <: Integer, PID <:Integer, LID <: Integer}
        new(data)
    end
end

"""
    BlockMap(numGlobalElements, comm)

Constructor for petra-defined uniform linear distribution of elements
"""
function BlockMap(numGlobalElements::Integer, comm::Comm{GID, PID, LID}) where GID <: Integer where PID <: Integer where LID <: Integer
    BlockMap(GID(numGlobalElements), comm)
end

function BlockMap(numGlobalElements::GID, comm::Comm{GID, PID, LID}) where GID <: Integer where PID <: Integer where LID <: Integer
    if numGlobalElements < 0
        throw(InvalidArgumentError("NumGlobalElements = $(numGlobalElements).  Should be >= 0"))
    end

    data = BlockMapData(numGlobalElements, comm)
    map = BlockMap{GID, PID, LID}(data)

    numProcVal = numProc(comm)
    data.linearMap = true

    myPIDVal = myPid(comm) - 1

    data.numMyElements = floor(typeof(data.numGlobalElements),
        data.numGlobalElements/numProcVal)
    remainder = data.numGlobalElements % numProcVal
    startIndex = myPIDVal * (data.numMyElements+1)

    if myPIDVal < remainder
        data.numMyElements += 1
    else
        startIndex -= (myPIDVal - remainder)
    end

    data.minAllGID = 1
    data.maxAllGID = data.minAllGID + data.numGlobalElements - 1
    data.minMyGID = startIndex + 1
    data.maxMyGID = data.minMyGID + data.numMyElements - 1
    data.distributedGlobal = isDistributedGlobal(map, data.numGlobalElements,
        data.numMyElements)

    EndOfConstructorOps(map)
    map
end

"""
    BlockMap(numGlobalElements, numMyElements, comm)

Constructor for user-defined linear distribution of elements
"""
function BlockMap(numGlobalElements::Integer, numMyElements::Integer, comm::Comm{GID, PID, LID}) where GID <: Integer where PID <: Integer where LID <: Integer
    BlockMap(GID(numGlobalElements), LID(numMyElements), comm)
end

function BlockMap(numGlobalElements::GID, numMyElements::LID, comm::Comm{GID, PID, LID}) where GID <: Integer where PID <: Integer where LID <: Integer
    if numGlobalElements < -1
        throw(InvalidArgumentError("NumGlobalElements = $(numGlobalElements).  Should be >= -1"))
    end
    if numMyElements < 0
        throw(InvalidArgumentError("NumMyElements = $(numMyElements). Should be >= 0"))
    end

    data = BlockMapData(numGlobalElements, comm)
    map = BlockMap{GID, PID, LID}(data)

    data.numMyElements = numMyElements
    data.linearMap = true

    data.distributedGlobal = isDistributedGlobal(map, numGlobalElements, numMyElements)

    #Local Map and uniprocessor case: Each processor gets a complete copy of all elements
    if !data.distributedGlobal || numProc(comm) == 1
        data.numGlobalElements = data.numMyElements

        data.minAllGID = 1
        data.maxAllGID = data.minAllGID + data.numGlobalElements - 1
        data.minMyGID = 1
        data.maxMyGID = data.minMyGID + data.numMyElements - 1
    else
        tmp_numMyElements = data.numMyElements
        data.numGlobalElements = sumAll(data.comm, tmp_numMyElements)

        data.minAllGID = 1
        data.maxAllGID = data.minAllGID + data.numGlobalElements - 1

        tmp_numMyElements = data.numMyElements
        data.maxMyGID = scanSum(data.comm, tmp_numMyElements)

        startIndex = data.maxMyGID - data.numMyElements
        data.minMyGID = startIndex + 1
        data.maxMyGID = data.minMyGID + data.numMyElements - 1
    end
    checkValidNGE(map, numGlobalElements)

    EndOfConstructorOps(map)
    map
end

"""
    BlockMap(numGlobalElements, myGlobalElements, comm)

Constructor for user-defined arbitrary distribution of elements
"""
function BlockMap(numGlobalElements::Integer, myGlobalElements::AbstractArray{<:Integer}, comm::Comm{GID, PID,LID}
        ) where GID <: Integer where PID <: Integer where LID <: Integer
    BlockMap(numGlobalElements, Array{GID}(myGlobalElements), comm)
end

function BlockMap(numGlobalElements::Integer, myGlobalElements::UnitRange{GID}, comm::Comm{GID, PID, LID}
        ) where {GID <: Integer, PID <: Integer, LID <: Integer}

    numMyElements = LID(length(myGlobalElements))

    data = BlockMapData(GID(0), comm)
    map = BlockMap{GID, PID, LID}(data)

    data.numMyElements = numMyElements

    data.myGlobalElements = collect(myGlobalElements)
    data.minMyGID = first(myGlobalElements)
    data.maxMyGID = last(myGlobalElements)

    #TODO this doesn't check if there is overlap between processors
    # call the reduce to allow mixing this version with the abstract version
    data.linearMap = Bool(minAll(data.comm, 1))

    if numProc(comm) == 1
        data.numGlobalElements = data.numMyElements
        data.minAllGID = data.minMyGID
        data.maxAllGID = data.maxMyGID
    else
        tmp_send = [
            -((data.numMyElements > 0) ?
                data.minMyGID : Inf)
            , data.maxMyGID]

        tmp_recv = maxAll(data.comm, tmp_send)

        @assert typeof(tmp_recv[1]) <: Integer "Result type is $(typeof(tmp_recv[1])), should be subtype of Integer"

        data.minAllGID = -tmp_recv[1]
        data.maxAllGID =  tmp_recv[2]

        if numGlobalElements != -1
            data.numGlobalElements = numGlobalElements
        else
            if data.linearMap
                data.numGlobalElements = sumAll(data.comm, data.numMyElements)
            else
                #if 1+ GIDs shared between processors, need to total that correctly
                allIDs = gatherAll(data.comm, myGlobalElements)

                indexModifier = 1 - data.minAllGID
                maxGID = data.maxAllGID

                count = 0
                arr = falses(maxGID + indexModifier)
                for id in allIDs
                    if !arr[GID(id + indexModifier)]
                        arr[GID(id + indexModifier)] = true
                        count += 1
                    end
                end
                data.numGlobalElements = count
            end
        end
    end

    data.distributedGlobal = isDistributedGlobal(map, data.numGlobalElements, numMyElements)

    EndOfConstructorOps(map)
    map
end

function BlockMap(numGlobalElements::Integer, myGlobalElements::AbstractArray{GID}, comm::Comm{GID, PID,LID}
        ) where GID <: Integer where PID <: Integer where LID <: Integer
    numMyElements = LID(length(myGlobalElements))

    data = BlockMapData(GID(0), comm)
    map = BlockMap{GID, PID, LID}(data)

    data.numMyElements = numMyElements

    linear = 1
    if numMyElements > 0
        data.myGlobalElements = Array{GID, 1}(undef, numMyElements)

        data.myGlobalElements[1] = myGlobalElements[1]
        data.minMyGID = myGlobalElements[1]
        data.maxMyGID = myGlobalElements[1]

        for i = 2:numMyElements
            data.myGlobalElements[i] = myGlobalElements[i]
            data.minMyGID = min(data.minMyGID, myGlobalElements[i])
            data.maxMyGID = max(data.maxMyGID, myGlobalElements[i])

            if myGlobalElements[i] != myGlobalElements[i-1] + 1
                linear = 0
            end
        end
    else
        data.minMyGID = 1
        data.maxMyGID = 0
    end

    #TODO this doesn't check if there is overlap between processors
    data.linearMap = Bool(minAll(data.comm, linear))

    if numProc(comm) == 1
        data.numGlobalElements = data.numMyElements
        data.minAllGID = data.minMyGID
        data.maxAllGID = data.maxMyGID
    else
        tmp_send = [
            -((data.numMyElements > 0) ?
                data.minMyGID : Inf)
            , data.maxMyGID]

        tmp_recv = maxAll(data.comm, tmp_send)

        @assert typeof(tmp_recv[1]) <: Integer "Result type is $(typeof(tmp_recv[1])), should be subtype of Integer"

        data.minAllGID = -tmp_recv[1]
        data.maxAllGID =  tmp_recv[2]

        if numGlobalElements != -1
            data.numGlobalElements = numGlobalElements
        else
            if data.linearMap
                data.numGlobalElements = sumAll(data.comm, data.numMyElements)
            else
                #if 1+ GIDs shared between processors, need to total that correctly
                allIDs = gatherAll(data.comm, myGlobalElements)

                indexModifier = 1 - data.minAllGID
                maxGID = data.maxAllGID

                count = 0
                arr = falses(maxGID + indexModifier)
                for id in allIDs
                    if !arr[GID(id + indexModifier)]
                        arr[GID(id + indexModifier)] = true
                        count += 1
                    end
                end
                data.numGlobalElements = count
            end
        end
    end

    data.distributedGlobal = isDistributedGlobal(map, data.numGlobalElements, numMyElements)

    EndOfConstructorOps(map)
    map
end

"""
    BlockMap(numGlobalElements, numMyElements, myGlobalElements, isDistributedGlobal, minAllGID, maxAllGID, comm)

Constructor for user-defined arbitrary distribution of elements with all information on globals provided by the user
"""
function BlockMap(numGlobalElements::Integer, numMyElements::Integer,
        myGlobalElements::AbstractArray{GID}, userIsDistributedGlobal::Bool,
        userMinAllGID::Integer, userMaxAllGID::Integer, comm::Comm{GID, PID, LID}) where GID <: Integer where PID <: Integer where LID <: Integer
    BlockMap(GID(numGlobalElements), LID(numMyElements), Array{GID}(myGlobalElements), userIsDistributedGlobal,
        GID(userMinAllGID), GID(userMaxAllGID), comm)
end

function BlockMap(numGlobalElements::GID, numMyElements::LID,
        myGlobalElements::AbstractArray{GID}, userIsDistributedGlobal::Bool,
        userMinAllGID::GID, userMaxAllGID::GID, comm::Comm{GID, PID, LID}) where GID <: Integer where PID <: Integer where LID <: Integer
    if numGlobalElements < -1
        throw(InvalidArgumentError("NumGlobalElements = $(numGlobalElements).  Should be >= -1"))
    end
    if numMyElements < 0
        throw(InvalidArgumentError("NumMyElements = $(numMyElements). Should be >= 0"))
    end
    if userMinAllGID < 1
        throw(InvalidArgumentError("Minimum global element index = $(data.minAllGID).  Should be >= 1"))
    end

    data = BlockMapData(numGlobalElements, comm)
    map = BlockMap{GID, PID, LID}(data)

    data.numMyElements = numMyElements

    linear = 1
    if numMyElements > 0
        data.myGlobalElements = Vector{GID}(undef, numMyElements)

        data.myGlobalElements[1] = myGlobalElements[1]
        data.minMyGID = myGlobalElements[1]
        data.maxMyGID = myGlobalElements[1]

        for i = 2:numMyElements
            data.myGlobalElements[i] = myGlobalElements[i]
            data.minMyGID = min(data.minMyGID, myGlobalElements[i])
            data.maxMyGID = max(data.maxMyGID, myGlobalElements[i])

            if myGlobalElements[i] != myGlobalElements[i-1] + 1
                linear = 0
            end
        end

    else
        data.minMyGID = 1
        data.maxMyGID = 0
    end

    data.linearMap = Bool(minAll(comm, linear))

    data.distributedGlobal = userIsDistributedGlobal

     if !data.distributedGlobal || numProc(comm) == 1
        data.numGlobalElements = data.numMyElements
        checkValidNGE(map, numGlobalElements)
        data.minAllGID = data.minMyGID
        data.maxAllGID = data.maxMyGID
    else
        if numGlobalElements == -1
            data.numGlobalElements = sumAll(data.comm, data.numMyElements)
        else
            data.numGlobalElements = numGlobalElements
        end
        checkValidNGE(data.numGlobalELements)

        data.minAllGID = userMinAllGID
        data.maxAllGID = userMaxAllGID
    end
    EndOfConstructorOps(map)
    map
end

##### internal construction methods #####
function isDistributedGlobal(map::BlockMap{GID, PID, LID}, numGlobalElements::GID,
        numMyElements::LID) where GID <: Integer where PID <: Integer where LID <: Integer
    data = map.data
    if numProc(data.comm) > 1
        localReplicated = numGlobalElements == numMyElements
        !Bool(minAll(data.comm, localReplicated))
    else
        false
    end
end

function EndOfConstructorOps(map::BlockMap)
    map.data.minLID = 1
    map.data.maxLID = max(map.data.numMyElements, 1)

    GlobalToLocalSetup(map);
end

function GlobalToLocalSetup(map::BlockMap)
    data = map.data
    numMyElements = data.numMyElements
    myGlobalElements = data.myGlobalElements

    if data.linearMap || numMyElements == 0
        return map
    end
    if length(data.numGlobalElements) == 0
        return map
    end

    val = myGlobalElements[1]
    i = 1
    while i < numMyElements
        if val+1 != myGlobalElements[i+1]
            break
        end
        val += 1
    end

    data.lastContiguousGIDLoc = i
    if data.lastContiguousGIDLoc <= 1
        data.lastContiguousGID = myGlobalElements[1]
    else
        data.lastContiguousGID = myGlobalElements[data.lastContiguousGIDLoc]
    end

    if i < numMyElements
        data.lidHash = empty!(data.lidHash)

        sizehint!(data.lidHash, numMyElements - i + 2)

        for i = i:numMyElements
            data.lidHash[myGlobalElements[i]] = i
        end
    end
    map
end

function checkValidNGE(map::BlockMap{GID, PID, LID}, numGlobalElements::GID) where GID <: Integer where PID <: Integer where LID <: Integer
    if (numGlobalElements != -1) && (numGlobalElements != map.data.numGlobalElements)
        throw(InvalidArgumentError("Invalid NumGlobalElements.  "
              * "NumGlobalElements = $(numGlobalElements)"
              * ".  Should equal $(map.data.numGlobalElements)"
              * ", or be set to -1 to compute automatically"))
    end
end

##### external methods #####

getComm(map::BlockMap) = map.data.comm

"""
    myGID(map::BlockMap, gidVal::Integer)

Return true if the GID passed in belongs to the calling processor in this
map, otherwise returns false.
"""
myGID(map::BlockMap, gidVal) = lid(map, gidVal) != 0

"""
    myLID(map::BlockMap, lidVal::Integer)

Return true if the LID passed in belongs to the calling processor in this
map, otherwise returns false.
"""
@inline myLID(map::BlockMap, lidVal) = gid(map, lidVal) != 0

"""
    distributedGlobal(map::BlockMap)

Return true if map is defined across more than one processor
"""
distributedGlobal(map::BlockMap) = map.data.distributedGlobal

"""
    numMyElements(map::BlockMap{GID, PID, LID})::LID

Return the number of elements across the calling processor
"""
numMyElements(map::BlockMap) = map.data.numMyElements

"""
    minMyGID(map::BlockMap{GID, PID, LID})::GID

Return the minimum global ID owned by this processor
"""
minMyGID(map::BlockMap) = map.data.minMyGID

"""
    maxMyGID(map::BlockMap{GID, PID, LID})::GID

Return the maximum global ID owned by this processor
"""
maxMyGID(map::BlockMap) = map.data.maxMyGID

"""
    getLocalMap(::BlockMap{GID, PID, LID})::BlockMap{GID, PID, LID}

Creates a copy of the given map that doesn't support any inter-process actions
"""
function getLocalMap(map::BlockMap{GID, PID, LID})::BlockMap{GID, PID, LID} where {GID, PID, LID}
    oldData = map.data
    data = BlockMapData(oldData.numGlobalElements, LocalComm(oldData.comm))

    data.directory = nothing
    data.lid = copy(oldData.lid)
    #maps shouldn't be modified anyways, may as well share array
    #data.myGlobalElements =  copy(oldData.myGlobalElements)
    data.myGlobalElements = oldData.myGlobalElements
    data.numMyElements = oldData.numMyElements
    data.minAllGID = oldData.minAllGID
    data.maxAllGID = oldData.maxAllGID
    data.minMyGID = oldData.minMyGID
    data.maxMyGID = oldData.maxMyGID
    data.minLID = oldData.minLID
    data.maxLID = oldData.maxLID
    data.linearMap = oldData.linearMap
    data.distributedGlobal = oldData.distributedGlobal
    data.oneToOneIsDetermined = oldData.oneToOneIsDetermined
    data.oneToOne = oldData.oneToOne
    data.lastContiguousGID = oldData.lastContiguousGID
    data.lastContiguousGIDLoc = oldData.lastContiguousGIDLoc
    data.lidHash = copy(oldData.lidHash)

    BlockMap{GID, PID, LID}(data)
end

##local/global ID accessor methods##

"""
    remoteIDList(map::BlockMap{GID, PID, LID}, gidList::AbstractArray{<: Integer}::Tuple{AbstractArray{PID}, AbstractArray{LID}}

Return the processor ID and local index value for a given list of global indices.
The returned value is a tuple containing
1. an Array of processors owning the global ID's in question
2. an Array of local IDs of the global on the owning processor
"""
function remoteIDList(map::BlockMap{GID, PID, LID}, gidList::AbstractArray{<:Integer}
        )::Tuple{AbstractArray{PID}, AbstractArray{LID}} where GID <: Integer where PID <: Integer where LID <: Integer
    remoteIDList(map, Vector{GID}(gidList))
end

function remoteIDList(map::BlockMap{GID, PID, LID}, gidList::AbstractArray{GID}
        )::Tuple{AbstractArray{PID}, AbstractArray{LID}} where GID <: Integer where PID <: Integer where LID <: Integer
    data = map.data
    if data.directory == nothing
        data.directory = createDirectory(data.comm, map)
    end

    getDirectoryEntries(data.directory, map, gidList)
end

"""
    lid(map::BlockMap{GID, PID, LID}, gid::Integer)::LID

Return local ID of global ID, or 0 if not found on this processor
"""
@inline function lid(map::BlockMap{GID, PID, LID}, gid::Integer) where GID <: Integer where PID <: Integer where LID <: Integer
    data = map.data
    if (gid < data.minMyGID) || (gid > data.maxMyGID)
         LID(0)
    elseif data.linearMap
        LID(gid - data.minMyGID + 1)
    elseif gid >= data.myGlobalElements[1] && gid <= data.lastContiguousGID
        LID(gid - data.myGlobalElements[1] + 1)
    elseif haskey(data.lidHash, GID(gid))
        data.lidHash[gid]
    else
        LID(0)
    end
end

"""
    gid(map::BlockMap{GID, PID, LID}, lid::Integer)::GID

Return global ID of local ID, or 0 if not found on this processor
"""
@inline function gid(map::BlockMap{GID, PID, LID}, lid::Integer) where GID <: Integer where PID <: Integer where LID <: Integer
    data = map.data
    if (data.numMyElements == LID(0)) || (lid < data.minLID) || (lid > data.maxLID)
        GID(0)
    elseif data.linearMap
        GID(lid + data.minMyGID - 1)
    else
        GID(data.myGlobalElements[lid])
    end
end

"""
    minAllGID(map::BlockMap{GID, PID, LID})::GID

Return the minimum global ID across the entire map
"""
minAllGID(map::BlockMap) = map.data.minAllGID

"""
    maxAllGID(map::BlockMap{GID, PID, LID})::GID

Return the maximum global ID across the entire map
"""
maxAllGID(map::BlockMap) = map.data.maxAllGID

"""
    minLID(map::BlockMap{GID, PID, LID})::LID

Return the mimimum local index value on the calling processor
"""
minLID(map::BlockMap) = map.data.minLID

"""
    maxLID(map::BlockMap{GID, PID, LID})::LID

Return the maximum local index value on the calling processor
"""
maxLID(map::BlockMap) = map.data.maxLID

##size/dimension accessor functions##

"""
    numGlobalElements(map::BlockMap{GID, PID, LID})::GID

Return the number of elements across all processors
"""
numGlobalElements(map::BlockMap) = map.data.numGlobalElements

"""
    myGlobalElements(map::BlockMap{GID, PID, LID})::AbstractArray{GID}

Return a list of global elements on this processor
"""
function myGlobalElements(map::BlockMap{GID})::AbstractArray{GID} where GID <: Integer
    data = map.data

    if length(data.myGlobalElements) == 0
        myGlobalElements = Vector{GID}(undef, data.numMyElements)
        @inbounds for i = GID(1):GID(data.numMyElements)
            myGlobalElements[i] = data.minMyGID + i - 1
        end
        data.myGlobalElements = myGlobalElements
    else
        data.myGlobalElements
    end
end

##Miscellaneous boolean tests##

"""
    uniqueGIDs(map::BlockMap)::Bool

Return true if each map GID exists on at most 1 processor
"""
uniqueGIDs(map::BlockMap) = isOneToOne(map)

"""
    globalIndicesType(map::BlockMap{GID, PID, LID})::Type{GID}

Return the type used for global indices in the map
"""
globalIndicesType(map::BlockMap{GID}) where GID <: Integer = GID

"""
    sameBlockMapDataAs(this::BlockMap, other::BlockMap)::Bool

Return true if the maps have the same data
"""
sameBlockMapDataAs(this::BlockMap, other::BlockMap) = this.data == other.data

"""
    sameAs(this::BlockMap, other::BlockMap)::Bool

Return true if this and other are identical maps
"""
sameAs(this::BlockMap, other::BlockMap) = false
# behavior by specification

function sameAs(this::BlockMap{GID, PID, LID}, other::BlockMap{GID, PID, LID})::Bool where GID <: Integer where PID <: Integer where LID <: Integer
    tData = this.data
    oData = other.data
    if tData == oData
        return true
    end

    if ((tData.minAllGID != oData.minAllGID)
        || (tData.maxAllGID != oData.maxAllGID)
        || (tData.numGlobalElements != oData.numGlobalElements))
        return false
    end

    mySameMap = 1

    if tData.numMyElements != oData.numMyElements
        mySameMap = 0
    end

    if tData.linearMap && oData.linearMap
        # For linear maps, just need to check whether lower bound is the same
        if tData.minMyGID != oData.minMyGID
            mySameMap = 0
        end
    else
        for i = 1:tData.numMyElements
            if gid(this, i) != gid(other, i)
                mySameMap = 0
                break
            end
        end
    end

    Bool(minAll(tData.comm, mySameMap))
end

"""
    linearMap(map::BlockMap)::Bool

Return true if the global ID space is contiguously divided (but
not necessarily uniformly) across all processors
"""
linearMap(map::BlockMap) = map.data.linearMap

##Array accessor functions##

"""
    myGlobalElementsIDs map::BlockMap{GID, PID, LID})::AbstractArray{GID}

Return list of global IDs assigned to the calling processor
"""
function myGlobalElementIDs(map::BlockMap{GID})::AbstractArray{GID} where GID <: Integer
    data = map.data
    if length(data.myGlobalElements) == 0
        base = 0:data.numMyElements-1
        rng = data.minMyGID .+ base
        collect(rng)
    else
        copy(data.myGlobalElements)
    end
end

function isOneToOne(map::BlockMap)::Bool
    data = map.data
    if !(data.oneToOneIsDetermined)
        data.oneToOne = determineIsOneToOne(map)
        data.oneToOneIsDetermined = true
    end
    data.oneToOne
end

function determineIsOneToOne(map::BlockMap)::Bool
    data = map.data
    if numProc(data.comm) < 2
        true
    else
        if data.directory == nothing
            data.directory = createDirectory(data.comm, map)
        end
       gidsAllUniquelyOwned(data.directory)
    end
end


text:
# # TODO: is this nunecessary? Maybe i'll need once I start doing unregularized...
# """Takes the indices of igrams to ignore, then shrinks down the full size B matrix and
# timediffs array to remove zeros
# """
# function cut_igram_dates(A, B, timediffs, ignore_int_indices)
#     row_masks = find_row_masks(size(A, 1), ignore_int_indices)
#     column_masks = find_column_masks(A[row_masks, :])
#
#     timediffs_cut = collapse_timediffs(timediffs, column_masks)
#     A_cut = A[row_masks, column_masks]
#     return timediffs_cut, B_cut
# end
#
# function find_row_masks(num_rows, rows_to_delete)
#     row_masks = trues(num_rows)
#     row_masks[rows_to_delete] .= false
#     return row_masks
# end
#
#
# """Take the row-deleted A, find columns that lost all igrams to delete"""
# function find_column_masks(A_short)
#     # Now find any columns that are all zero to delete, make it 1D
#     column_sums = vec(sum(abs2, A_short, dims=1))
#
#     # Trues are good columns (nonzero), false is a bad column
#     column_masks = .!iszero.(column_sums)
# end
#
#
# """Using the missing columns, recalculate the time diff vector"""
# function collapse_timediffs(timediffs, column_masks)
#     out_tds = Array{Int}(undef, 0)
#     cur_td = 0
#     for (idx, td) in enumerate(timediffs)
#         cur_td += td
#         if column_masks[idx]
#             append!(out_tds, cur_td)
#             cur_td = 0
#         end
#     end
#     return out_tds
# end
#


text:
using ModelParams, Test

mdl = ModelParams;

function param_test()
    @testset "Param" begin
        pValue = [1.1 2.2; 3.3 4.4]

        # Simple constructor
        p1 = Param(:p1, "param1", "\$p_{1}\$", pValue);
        show(p1);
        @test p1.value == pValue
        @test size(p1.lb) == size(pValue)
        calibrate!(p1)
        @test p1.isCalibrated == true
        fix!(p1)
        @test p1.isCalibrated == false
        validate(p1)

        # Full constructor
        pValue2 = 1.23
        lb2 = -2.8
        ub2 = 7.3
        p2 = Param(:p2, "param2", "\$p_{2}\$", pValue2, pValue2 + 0.5,
            lb2, ub2, true)
        show(p2)
        @test p2.value == pValue2
        @test size(p2.ub) == size(pValue2)
        newValue = 9.27
        set_value!(p2, newValue)
        @test p2.value ≈ newValue
        str1 = ModelParams.short_string(p2)
        @test startswith(str1,  "p2: 9.27")

        update!(p2, value = 98.2)
        @test p2.value == 98.2
        @test p2.lb == lb2
        update!(p2, lb = -2.0, ub = 8.0)
        validate(p2)
        @test p2.lb == -2.0
        @test p2.value == 98.2
    end
end

function vec1_test()
    @testset "Vector of length 1" begin
        pValue = fill(1.2, 1);
        p1 = Param(:p1, "param1", "\$p_{1}\$", pValue, pValue .+ 0.1,
            fill(0.1, 1), fill(5.0, 1), true);
        show(p1);
        @test isequal(pValue, mdl.value(p1))
    end
end

function array_test(p)
    @testset "Array param $p" begin
        @test validate(p);
        @test size(value(p)) == size(default_value(p));
        @test size(mdl.lb(p)) == size(mdl.ub(p)) == size(value(p));
        closeToLb = mdl.close_to_lb(p; rtol = 0.1);
        closeToUb = mdl.close_to_ub(p; rtol = 0.1);
        @test mdl.close_to_bounds(p; rtol = 0.1) == (closeToLb || closeToUb);
        isCal = is_calibrated(p);
        mdl.set_calibration_status!(p, !isCal);
        @test is_calibrated(p) == !isCal;

        lbnd = mdl.lb(p) .+ 0.1;
        ubnd = mdl.ub(p) .+ 0.2;
        set_bounds!(p; lb = lbnd, ub = ubnd);
        @test mdl.lb(p) ≈ lbnd;
        @test mdl.ub(p) ≈ ubnd;
    end
end

function make_test_array_param()
    v = [1.0 2.0; 3.0 4.0; 5.0 6.0];
    p = Param(:x, "x", "x", v, v .+ 0.1, v .- 2.0, v .+ 2.0, true);
    return p
end

@testset "Parameters" begin
    param_test();
    vec1_test();
    for p in (
        make_test_array_param(),
        mdl.make_test_cal_array(:x, 2),
        mdl.make_test_bvector(:x; increasing = :increasing),
        mdl.make_test_bvector(:x; increasing = :decreasing),
        )
        array_test(p);
    end
end

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


text:
# License for this file: MIT (expat)
# Copyright 2020, DLR Institute of System Dynamics and Control

"""
    module ModiaMathTest_ode_pendulum2

ODE-model of a mass point attached via a rod to a revolute joint with 2 states.

This test is used to evalute how to use DifferentialEquations.jl in ModiaMath.
The difference to ModiaMathTest_ode_pendulum.jl is that the model just uses a Float64 type.
"""
module ModiaMathTest_ode_pendulum2

import DifferentialEquations
using PyPlot        # For plotting

using StaticArrays  # Test of StaticArrays with parameterized data

include("ModiaMathTest.jl")
using .ModiaMathTest

#--------------------------------------- Generic (model independent) definitions

mutable struct Variable
    typ
    start
    fixed::Union{Bool,Nothing}
    info::AbstractString

    Variable(;typ=Real, start=nothing, fixed=nothing, info="") = new(typ, start, fixed, info)
end
const VariableDict = Dict{Symbol, Variable}

mutable struct ModelResult
    t::Vector{Float64}
    model::AbstractVector
    variables::VariableDict

    ModelResult(m) = new(Float64[], typeof(m.model)[], m.variables)
end

function addToResult!(result::ModelResult, t, m)::Nothing
    push!(result.t, deepcopy(t))
    push!(result.model, deepcopy(m.model))
    return nothing;
end

#---------------------------------------- Specific to a particular model (here: Pendulum)

mutable struct Model
    # Parameters
    L::Float64
    m::Float64
    d::Float64
    g::Float64

    # Variables during integration
    phi::Float64
    w::Float64
    der_w::Float64

    # Variables at communication points
    r::SVector{2,Float64}   # absolute position of end point

    function Model(;L=1.0, m=1.0, d=0.1, g=9.81)
        @assert(L > 0)
        @assert(m > 0)
        @assert(d >= 0)
        new(L, m, d, g, 0, 0, 0, SVector{2,Float64}(0,0))
    end
end

mutable struct Pendulum
    # Required:
    simulationState::SimulationState{Float64,Float64}

    # Defined by model
    model::Model               # Data structure of model variable values
    variables::VariableDict    # Data structure of model variable definition
end

function ModelVariables()::VariableDict
    var = VariableDict()

    # Parameters
    var[:L] = Variable(info="Angle")
    var[:m] = Variable(info="Mass")
    var[:d] = Variable(info="Damping factor")
    var[:g] = Variable(info="Gravity acceleration")

    # Variables during integration
    var[:phi]   = Variable(start=π/2, info="Angle")
    var[:w]     = Variable(start=0.0, info="Angular Velocity (= der(phi))")
    var[:der_w] = Variable(info="Angular Acceleration (= der(w))")

    # Variables at communication points
    var[:r]     = Variable(typ=SVector{2,Real}, info="Absolute position of pendulum tip")

    return var
end

function InitialStates(model::Model, var::VariableDict)::Vector{Float64}
    x_start = [Float64(var[:phi].start),
               Float64(var[:w].start)]
    return x_start
end

function setStates!(obj::Pendulum, x)::Nothing
    m = obj.model
    m.phi = x[1]
    m.w   = x[2]
    return nothing
end

function getDerivatives!(obj::Pendulum, der_x)::Nothing
    m = obj.model
    der_x[1] = m.w
    der_x[2] = m.der_w
    return nothing
end

function evaluate!(obj::Pendulum)::Nothing
    m   = obj.model
    sim = obj.simulationState
    m.der_w = (-m.m * m.g * m.L * sin(m.phi) - m.d * m.w) / (m.m * m.L^2)

    if sim.storeResult
        m.r = SVector{2,Float64}(m.L*sin(m.phi), -m.L*cos(m.phi))
    end
    return nothing
end

"Pendulum: Constructor for Pendulum SimulationModel"
function Pendulum(;kwargs...)
    model     = Model(;kwargs...)
    variables = ModelVariables()
    x_start   = InitialStates(model, variables)

    simulationState = SimulationState{Float64,Float64}("Pendulum", x_start, setStates!, getDerivatives!, evaluate!, ModelResult, addToResult!)
    simulationModel = Pendulum(simulationState, model, variables)
end

#------------------------------------------------------------------------------------------------
#---------------------------------------- Simulation of the pendulum model
#------------------------------------------------------------------------------------------------

function standardPlot(solution, fig1, fig2)
    result = solution.prob.p.simulationState.result
    Float64 = string( typeof( solution.prob.p.simulationState.x_start[1] ) )
    figure(fig1)
    clf()
    t   = result.t
    phi = getfield.(result.model, :phi)
    plot(t, phi, label="\$\\varphi\$")
    grid(true)
    legend()
    title("SimplePendulum with " * Float64)

    figure(fig2)
    clf()
    r  = getfield.(result.model, :r)
    r1 = getindex.(r,1)
    r2 = getindex.(r,2)
    plot(t, r1, label="r[1]")
    plot(t, r2, label="r[2]")
    grid(true)
    legend()
    title("SimplePendulum with " * Float64)
end

# ODE integrator with Float64
pendulum = Pendulum(m=1.1)
      solution = solve!(pendulum, DifferentialEquations.Tsit5(), stopTime=7.0, tolerance=1e-4)
@time solution = solve!(pendulum, DifferentialEquations.Tsit5(), stopTime=7.0, tolerance=1e-4)
standardPlot(solution, 1, 2)

# DAE integrator with FLoat64
      solution = solve!(pendulum, DifferentialEquations.IDA(), stopTime=7.0, tolerance=1e-4)
@time solution = solve!(pendulum, DifferentialEquations.IDA(), stopTime=7.0, tolerance=1e-4)
standardPlot(solution, 3, 4)

end


text:
module DropoutM
    using LoopVectorization, VectorizedRNG

    mutable struct Dropout
        save_layer::Any
        load_layer::Any
        activate::Any
        initialize::Any
        update::Any

        input_shape::Tuple
        output_shape::Tuple

        rate::Float64
        table::Array{Float32}

        output::Array{Float32}
        δ::Array{Float32}

        function Dropout(;input_shape::Tuple, rate::Float64)
            new(save_Dropout, load_Dropout, activate_Dropout, init_Dropout, update_Dropout, input_shape, input_shape, rate)
        end
    end

    function init_Dropout(layer::Dropout, mini_batch::Int64)
        layer.table = zeros(Float32, layer.output_shape..., mini_batch)
        layer.output = zeros(Float32, layer.output_shape..., mini_batch)
        layer.δ = zeros(Float32, layer.input_shape..., mini_batch)
    end

    function activate_Dropout(layer::Dropout, input::Array{Float32})
        rand!(local_rng(), layer.table)
        @avxt for i in eachindex(input)
            layer.output[i] = ifelse(layer.table[i]<=layer.rate, 0.0f0, input[i])
        end
    end

    function update_Dropout(layer::Dropout, optimizer::String, input::Array{Float32}, δₗ₊₁::Array{Float32}, α::Float64, parameters::Tuple, direction::Int64=1)
        @avxt for i in eachindex(δₗ₊₁)
            layer.δ[i] = ifelse(layer.table[i]<=layer.rate, 0.0f0, δₗ₊₁[i])
        end
    end

    function save_Dropout(layer::Dropout, file::Any, id::String)
        write(file, id, "Dropout")
        write(file, id*"input_shape", collect(layer.input_shape))
        write(file, id*"rate", layer.rate)
    end

    function load_Dropout(layer::Dropout, file::Any, id::String)
    end

    function get_args(file::Any, id::String)
        input_shape = tuple(read(file, id*"input_shape")...)
        rate = read(file, id*"rate")
        return (input_shape=input_shape, rate=rate)
    end
end


text:
using SeisMain
using Test

#From USGS, National Petroleum Reserve - Alaska Data Archive
#https://energy.usgs.gov/GeochemistryGeophysics/SeismicDataProcessingInterpretation/NPRASeismicDataArchive.aspx#3862174-data-

download("http://certmapper.cr.usgs.gov/data/NPRA/seismic/1976/125_76/PROCESSED/125_76_PT2_PR.SGY","125_76.SGY")

SegyToSeis("125_76.SGY","125_76")
d,h,ext = SeisRead("125_76")
imx = SeisMain.ExtractHeader(h,"imx")
println("size(d)=",size(d))
println("imx=",imx)
SeisWrite("125_76_copy",d,h,ext)


text:
using LinearAlgebra

include("./BEMUtils.jl")

hwlist = haarCube(2)
chargeList = points(1.0 .+ rand(3, 100))

# BEMatrix = assembleMatrix(Φ, hwlist, hwlist)

@time constructMatrixElement(Φ, hwlist[1], hwlist[3])


text:
@testset "parse single session" begin
    ses_path = joinpath(
        @__DIR__,
        "data",
        "bids_root",
        "sub-subtest",
        "ses-1"
    )

    my_ses = Session(ses_path)

    @test my_ses.path == ses_path
    @test my_ses.identifier == "1"
    @test my_ses.scans_detail == DataFrame()
    @test total_files(my_ses) == 1
    @test sprint(show, my_ses) == """
        Session:
            identifier = 1
            total files = 1"""
end


text:
# Shape inference

#### preceding or succeeding length

function prec_length(s::Dims, d::Int)
    d == 1 ? 1 :
    d == 2 ? s[1] :
    d == 3 ? s[1] * s[2] : prod(s[1:d-1])
end

function succ_length{N}(s::NTuple{N,Int}, d::Int)
    d == N ? 1 :
    d == N-1 ? s[N] : prod(s[d+1:N])
end

### Shape of mapping

# on size tuples

mapshape(s::Dims) = s

mapshape{N}(s1::NTuple{N,Int}, s2::NTuple{N,Int}) = (s1 == s2 || error("Argument dimensions are not map-compatible."); s1)

function mapshape{N1,N2}(s1::NTuple{N1,Int}, s2::NTuple{N2,Int})
    if N1 < N2
        for i = 1 : N1
            s1[i] == s2[i] || error("Argument dimensions are not map-compatible.")
        end
        for i = N1+1 : N2
            s2[i] == 1 || error("Argument dimensions are not map-compatible.")
        end
        return s2
    else
        for i = 1 : N2
            s1[i] == s2[i] || error("Argument dimensions are not map-compatible.")
        end
        for i = N2+1 : N1
            s1[i] == 1 || error("Argument dimensions are not map-compatible.")
        end
        return s1
    end
end

mapshape(s1::(), s2::()) = ()
mapshape(s1::Dims, s2::()) = s1
mapshape(s1::(), s2::Dims) = s2

mapshape(s1::(Int,Int), s2::(Int,)) = ((s1[1] == s2[1] && s1[2] == 1) || error("Argument dimensions are not map-compatible."); s1)
mapshape(s1::(Int,), s2::(Int,Int)) = ((s1[1] == s2[1] && s2[2] == 1) || error("Argument dimensions are not map-compatible."); s2)

mapshape(s1::Dims, s2::Dims, s3::Dims, rs::Dims...) = mapshape(mapshape(s1, s2), mapshape(s3, rs...))

# on arrays / numbers

mapshape(a1::Number) = ()
mapshape(a1::AbstractArray) = size(a1)

mapshape(a1::Number, a2::Number) = ()
mapshape(a1::AbstractArray, a2::Number) = size(a1)
mapshape(a1::Number, a2::AbstractArray) = size(a2)
mapshape(a1::AbstractArray, a2::AbstractArray) = mapshape(size(a1), size(a2))

mapshape(a1::Number, a2::Number, a3::Number, ra::Number...) = ()
mapshape(a1::Union(AbstractArray,Number),
         a2::Union(AbstractArray,Number),
         a3::Union(AbstractArray,Number),
         ra::Union(AbstractArray,Number)...) = mapshape(mapshape(a1, a2), mapshape(a3, ra...))

### Length of mapping

maplength(a1::Number) = 1
maplength(a1::AbstractArray) = length(a1)

maplength(a1::Number, a2::Number) = 1
maplength(a1::AbstractArray, a2::Number) = length(a1)
maplength(a1::Number, a2::AbstractArray) = length(a2)
maplength(a1::AbstractArray, a2::AbstractArray) = prod(mapshape(a1, a2))

maplength(a1::Number, a2::Number, a3::Number, ra::Number...) = 1
maplength(a1::Union(AbstractArray,Number),
          a2::Union(AbstractArray,Number),
          a3::Union(AbstractArray,Number),
          ra::Union(AbstractArray,Number)...) = prod(mapshape(a1, a2, a3, ra...))