text:
@testset "Broadcast" begin
    for T in (Float32, Float64)
        v = rand(Vec{3, T})
        A = rand(Mat{2, 3, T})
        vs = [rand(Vec{3, T}) for _ in 1:4]
        As = [rand(Mat{2, 3, T}) for _ in 1:4]
        # with array (tensors behave like scalar)
        @test (@inferred(v .+ vs))::Vector{Vec{3, T}} ≈ map(y -> v + y, vs)
        @test (@inferred(v .- vs))::Vector{Vec{3, T}} ≈ map(y -> v - y, vs)
        @test (@inferred(v .+ vs .- vs .+ v))::Vector{Vec{3, T}} ≈ map(y -> v + y - y + v, vs)
        @test (@inferred(broadcast(+, v, vs, vs, v)))::Vector{Vec{3, T}} ≈ map(y -> v + y + y + v, vs)
        @test (vs .= v)::Vector{Vec{3, T}} ≈ map(y -> v, vs)
        @test (Array(v) .= v)::Vector{T} ≈ v
        @test (@inferred(A .+ As))::Vector{Mat{2, 3, T, 6}} ≈ map(y -> A + y, As)
        @test (@inferred(A .- As))::Vector{Mat{2, 3, T, 6}} ≈ map(y -> A - y, As)
        @test (@inferred(A .+ As .- As .+ A))::Vector{Mat{2, 3, T, 6}} ≈ map(y -> A + y - y + A, As)
        @test (@inferred(broadcast(+, A, As, As, A)))::Vector{Mat{2, 3, T, 6}} ≈ map(y -> A + y + y + A, As)
        @test (As .= A)::Vector{Mat{2, 3, T, 6}} ≈ map(y -> A, As)
        @test (Array(A) .= A)::Matrix{T} ≈ A
        # with scalar
        @test (@inferred(v .+ 1))::Vec{3, T} ≈ map(y -> y + 1, v)
        @test (@inferred(A .+ 1))::Mat{2, 3, T} ≈ map(y -> y + 1, A)
        @test (@inferred(v .+ 1 .+ 2 .+ v))::Vec{3, T} ≈ map(y -> y + 1 + 2 + y, v)
        @test (@inferred(A .+ 1 .+ 2 .+ A))::Mat{2, 3, T} ≈ map(y -> y + 1 + 2 + y, A)
        # with tuple
        @test (@inferred(v .+ (1,2,3)))::Vec{3, T} ≈ map(+, v, (1,2,3))
        @test (@inferred(v .- (1,2,3)))::Vec{3, T} ≈ map(-, v, (1,2,3))
        # others
        @test (@inferred(broadcast(sqrt, v)))::Vec{3, T} ≈ map(sqrt, v)
        @test (@inferred(broadcast(sqrt, A)))::Mat{2, 3, T} ≈ map(sqrt, A)
    end
end


text:
function solve(pep, μ, δ, β)
    Tstar, wstar = oracle(pep, μ);
    ⋆ = istar(pep, μ);

    # lower bound
    kl = (1-2δ)*log((1-δ)/δ);
    lbd = Tstar*kl;

    # more practical lower bound with the employed threshold β
    practical = binary_search(t -> t-Tstar*β(t), max(1, lbd), 1e10);

    Tstar, wstar, ⋆, lbd, practical;
end

function dump_stats(pep, μ, δs, βs, srs, datas)

    for i in 1:length(δs)
        δ = δs[i];
        β = βs[i];
        data = getindex.(datas, i);

        Tstar, wstar, ⋆, lbd, practical = solve(pep, μ, δ, β)

        rule = repeat("-", 60);

        println("");
        println(rule);
        println("$pep at δ = $δ");
        println(@sprintf("%-30s", "Arm"),
                join(map(k -> @sprintf("%6s", k), 1:length(μ))), " ",
                @sprintf("%7s", "total"), "  ",
                @sprintf("%7s", "err"), "  ",
                @sprintf("%7s", "time"));
        println(@sprintf("%-30s", "μ"),
                join(map(x -> @sprintf("%6.2f", x), μ)));
        println(@sprintf("%-30s", "w⋆"),
                join(map(x -> @sprintf("%6.2f", x), wstar)));
        println(rule);
        println(@sprintf("%-30s", "oracle"),
                join(map(w -> @sprintf("%6.0f", lbd*w), wstar)), " ",
                @sprintf("%7.0f", lbd));
        println(@sprintf("%-30s", "practical"),
                join(map(w -> @sprintf("%6.0f", practical*w), wstar)), " ",
                @sprintf("%7.0f", practical));
        println(rule);

        for r in eachindex(srs)
            Eτ = sum(x->sum(x[2]), data[r,:])/N;
            err = sum(x->x[1].!=⋆, data[r,:])/N;
            tim = sum(x->x[3],     data[r,:])/N;
            println(@sprintf("%-30s", long(srs[r])),
                    join(map(k -> @sprintf("%6.0f", sum(x->x[2][k], data[r,:])/N), eachindex(μ))), " ",
                    @sprintf("%7.0f", Eτ), "  ",
                    @sprintf("%7.5f", err), "  ",
                    @sprintf("%7.5f", tim/1e6)
                    );
            if err > δ
                @warn "too many errors for $(srs[r])";
            end
        end
        println(rule);
    end
end

function dump_stats(pep, μ, δs, βs, srs, datas, repeats)
    K = narms(pep, µ)

    for i in 1:length(δs)
        δ = δs[i];
        β = βs[i];
        data = getindex.(datas, i);

        #Tstar, wstar, ⋆, lbd, practical = solve(pep, μ, δ, β)
        ⋆ = istar(pep, µ)

        rule = repeat("-", 60);

        println("");
        println(rule);
        println("$pep at δ = $δ");
        println(@sprintf("%-30s", "Arm"),
                join(map(k -> @sprintf("%6s", k), 1:K)), " ",
                @sprintf("%7s", "total"), "  ",
                @sprintf("%7s", "err"), "  ",
                @sprintf("%7s", "time"));
        println(@sprintf("%-30s", "a\'μ"),
                join(map(x -> @sprintf("%6.2f", x'µ), pep.arms)));
        #println(@sprintf("%-30s", "w⋆"),
        #        join(map(x -> @sprintf("%6.2f", x), wstar)));
        println(rule);
        #println(@sprintf("%-30s", "oracle"),
        #        join(map(w -> @sprintf("%6.0f", lbd*w), wstar)), " ",
        #        @sprintf("%7.0f", lbd));
        #println(@sprintf("%-30s", "practical"),
        #        join(map(w -> @sprintf("%6.0f", practical*w), wstar)), " ",
        #        @sprintf("%7.0f", practical));
        println(rule);

        for r in eachindex(srs)
            Eτ = sum(x->sum(x[2]), data[r,:])/repeats;
            err = sum(x->x[1].!=⋆, data[r,:])/repeats;
            tim = sum(x->x[4],     data[r,:])/repeats;
            println(@sprintf("%-30s", long(srs[r])),
                    join(map(k -> @sprintf("%6.0f", sum(x->x[2][k], data[r,:])/repeats), 1:K)), " ",
                    @sprintf("%7.0f", Eτ), "  ",
                    @sprintf("%7.5f", err), "  ",
                    @sprintf("%7.5f", tim/1e6)
                    );
            if err > δ
                @warn "too many errors for $(srs[r])";
            end
        end
        println(rule);
    end
end

function τhist(pep, μ, δ, β, srs, data)
    Tstar, wstar, ⋆, lbd, practical = solve(pep, μ, δ, β)

    stephist(map(x -> sum(x[2]), data)', label=permutedims(collect(abbrev.(srs))));
    vline!([lbd], label="lower bd");
    vline!([practical], label="practical");
end

function boxes(pep, μ, δ, β, srs, data)
    #Tstar, wstar, ⋆, lbd, practical = solve(pep, μ, δ, β)

    xs = permutedims(collect(abbrev.(srs)));

    means = sum(sum.(getindex.(data,2)),dims=2)/repeats;

    #plot([lbd practical], seriestype=:hline, label=["lower bd" "practical"], legend=:top)
    plot(legend=:top)

    if true
        # violin(
        #     xs,
        #     map(x -> sum(x[2]), data)',
        #     label=""
        # )
        boxplot!(
            xs,
            map(x -> sum(x[2]), data)',
            label="",
            notch=true,
            outliers=true)

        plot!(xs,
              means', marker=(:star4,10,:black), label="");
        #top = max(4*practical, 2*maximum(means));
        # hard coded xlims are a dirty hack; I cannot figure out how
        # to draw the first two hlines without upsetting the x-axis
        #plot!(xlims=(1.2,15.4), ylims=(-.025*top, top));
    else
        bar(collect(abbrev.(srs)),
            means,
            label="");
    end

end

function randmax(vector, rank = 1)
   # returns an integer, not a CartesianIndex
    vector = vec(vector)
    Sorted = sort(vector, rev = true)
    m = Sorted[rank]
    Ind = findall(x -> x == m, vector)
    index = Ind[floor(Int, length(Ind) * rand())+1]
    return index
end

function randmin(vector, rank = 1)
   # returns an integer, not a CartesianIndex
    vector = vec(vector)
    Sorted = sort(vector, rev = false)
    m = Sorted[rank]
    Ind = findall(x -> x == m, vector)
    index = Ind[floor(Int, length(Ind) * rand())+1]
    return index
end

ext: jl


text:
# Find the slope of two point
find_slope(x1, y1, x2, y2) = (y2 - y1) / (x2 - x1)

# Find the closest whole numbers on both sides of a x
neighbors(x) = round(x), ceil(x)

# Find the predicted value of y given a value x
function find_y(x, known_x)
    # Check if the given value is exactly one in the known list
    if round(x) == x
        return known_x[x]
    end

    left_x, right_x = neighbors(x)

    left_y = known_x[left_x]
    right_y = known_x[right_x]

    slope = find_slope(left_x, left_y, right_x, right_y)

    return left_y + (slope * (x - left_x))
end

println(find_slope(2, 3, 5, 6))
println(neighbors(4.5))
println(find_y(2, [3, 4, 5, 9]))

max_stars_repo_name: JakeRoggenbuck/colorgradient-julia

max_issues_repo_name: JakeRoggenbuck/colorgradient-julia

max_forks_repo_name: JakeRoggenbuck/colorgradient-julia


text:
import Base: getindex, setindex!
export get_nv, getvertex, edgeparse, StingerEdge, StingerVertex

"""
    getindex(x::Stinger, field::StingerFields)

Obtain value of a field from the Stinger data structure.
For `batch_time` and `update_time`, use `get_batchtime`
and `get_updatetime` respectively.
"""
function getindex(x::Stinger, field::StingerFields)
    idx = Int(field)
    if field == batch_time || field == update_time
        error("For $field use get_$field()")
    end
    basepointer = convert(Ptr{Int64}, x.handle)
    unsafe_load(basepointer, idx)
end

function get_batchtime(x::Stinger)
    basepointer = convert(Ptr{Float64}, x.handle)
    unsafe_load(basepointer, Int(batch_time))
end

function get_updatetime(x::Stinger)
    basepointer = convert(Ptr{Float64}, x.handle)
    unsafe_load(basepointer, Int(update_time))
end

"""
    setindex!(x::Stinger, val, field::StingerFields)

Set value of a field from the Stinger data structure.
"""
function setindex!(x::Stinger, val, field::StingerFields)
    idx = Int(field)
    ftype = fieldtype(StingerGraph, idx)
    @assert isa(val, ftype)
    basepointer = convert(Ptr{ftype}, x.handle)
    unsafe_store!(basepointer,val,idx)
end

"""
    get_nv(x::Stinger)

Returns number of active vertices in the graph. This is based on the largest
vertex ID which has a non-zero indegree or outdegree.
"""
function get_nv(x::Stinger)
    nv = ccall(
        dlsym(stinger_core_lib, "stinger_max_active_vertex"),
        Int64,
        (Ptr{Void},),
        x
    )
    nv==0 && return nv
    nv+1
end

"""
    storageptr(s::Stinger)

Get a pointer to the storage array of the STINGER data structure
"""
storageptr(s::Stinger) = s.handle + sizeof(StingerGraph) + 5*sizeof(UInt64)

"""
The STINGER Vertex representation.
"""
immutable StingerVertex
    vtype::Int64
    weight::Int64
    indegree::Int64
    outdegree::Int64
    degree::Int64
    edges::Int64
end

const NUMEDGEBLOCKS = 14

"""
The STINGER Edge representation.
"""
immutable StingerEdge
    neighbor::Int64
    weight::Int64
    timefirst::Int64
    timerecent::Int64
end

immutable StingerEdgeBlock
    next::UInt64
    etype::Int64
    vertexid::Int64
    numedges::Int64
    high::Int64
    smallstamp::Int64
    largestamp::Int64
    cache_pad::Int64
end

"""
    edgeparse(edge::StingerEdge)

Parse the direction and neighbor given a `StingerEdge`.
The first 2 bits of the `neighbor` field of the edge denotes the direction.
1 - in, 2 - out, 3 - both
"""
function edgeparse(edge::StingerEdge)
    direction = edge.neighbor >> 61 #The first 2 bits denote the direction, 1 - in, 2 - out, 3 - both
    neighbor = ~(7 << 61) & edge.neighbor
    return direction, neighbor
end

"""
    getvertex(s::Stinger, v::Int64)

Load the specified vertex from the STINGER datastructure.
"""
function getvertex(s::Stinger, v::Int64)
    vertices = convert(Ptr{StingerVertex}, storageptr(s) + sizeof(Int64)) #Read the StingerVertex array
    vertex = unsafe_load(vertices, v+1)
    vertex
end

function getvertexedgesoffset(s::Stinger, v::Int64)
    vertex = stinger_vertex_get(s,v)
    vertex.edges
end


text:
module aoc_17
    using ImageFiltering
    using Test

    # Replaced by NeighborhoodFilter type below
    const nhood_filter = setindex!( centered( ones(3,3,3) ) ,0,0,0,0 )
    const nhood_filter_4d = setindex!( centered( ones(3,3,3,3) ) ,0,0,0,0,0 )

    function read_input(filename, nd = 3)
        lines = readlines(filename)
        universe = permutedims(hcat(Vector{Char}.(lines)...)) .== '#'
        reshape(universe, ( size(universe)..., ntuple(i->1, nd-2)... ) )
    end

    remain_active(neighbors) = (neighbors .== 2) .| (neighbors .== 3)
    become_active(neighbors) = neighbors .== 3

    function compact(universe)
        e = extrema(findall(universe))

        # 3D
        # v = view(universe, e[1][1]:e[2][1], e[1][2]:e[2][2], e[1][3]:e[2][3])
        # 4D
        # v = view(universe, e[1][1]:e[2][1], e[1][2]:e[2][2], e[1][3]:e[2][3], e[1][4]:e[2][4])

        # ND
        v = view(universe, [a:b for (a,b) in zip(Tuple.(e)...)]...)
        copy(v)
    end

    function part1(universe = read_input("input.txt",3) )
        simulate(universe, 6)[1]
    end

    function test1()
        @test simulate( read_input("demo.txt",3) , 6 )[1] == 112
    end

    function answer1()
        @test part1() == 322
    end

    # Part 2

    # === Begin NeighborhoodFilter type === #

    struct NeighborhoodFilter{N} <: AbstractArray{Float64,N}
        linear_center::UInt
        center::NTuple{N,UInt}
    end
    function NeighborhoodFilter{N}() where N
        linear_center = ceil(UInt,3^N/2)
        center = ntuple(i->UInt(0),N)
        NeighborhoodFilter{N}( linear_center, center )
    end

    import Base: size, getindex, axes
    size(f::NeighborhoodFilter{N}) where N = ntuple(i->3, N)
    getindex(f::NeighborhoodFilter{N}, i::Int) where N = f.linear_center == i ? 0 : 1
    getindex(f::NeighborhoodFilter{N}, I::Vararg{Int,N}) where N = I == f.center ? 0 : 1
    axes(f::NeighborhoodFilter{N}) where {N} = ntuple(i->-1:1,N)

    # === End NeighborhoodFilter type === #

    # n-Dimensional step
    function step(universe::BitArray{N}, make_compact = true) where {N}
        universe = padarray(universe, Fill(0, ntuple(i->1, N) ) )
        neighbors = round.(imfilter(universe, NeighborhoodFilter{N}(), Fill(0)))
        # If a cube is active and exactly 2 or 3 of its neighbors are also active, the cube remains active. Otherwise, the cube becomes inactive.
        active_state = copy(universe)
        inactive_state = .~active_state
        active = @view universe[ active_state ]
        inactive = @view universe[ inactive_state ]
        active .= remain_active( neighbors[ active_state ] )
        inactive .= become_active( neighbors[ inactive_state ] )
        if make_compact
            compact(universe)
        else
            universe
        end
    end

    function test2()
        @test simulate( read_input("demo.txt", 4), 6 )[1] == 848
    end

    function part2(universe = read_input("input.txt", 4); cycles = 6)
        simulate(universe, cycles)[1]
    end

    function answer2()
        @test part2() == 2000
    end

    function simulate(universe, cycles = 6)
        for i = 1:cycles
            @debug i
            universe = step(universe)
        end
        sum(universe), universe
    end

    function test()
        test1()
        answer1()
        test2()
        answer2()
    end

end


text:
using Sinograms: bdd_2d
using Test: @test, @testset, @inferred
# using JLD: load

@testset "bdd_2d" begin
    deg = 1
    geo = (DSD = 100000, DSO = 99700, pSize = 1, dSize = 0.5, nPix = 256, nDet = 1024,
          theta = deg2rad.(0:deg:360-deg), isoX = 0, isoY = 0)
    phantomImg = shepp_logan(geo.nPix, SheppLoganToft())
    sinogramB = projection(phantomImg',geo)
    imageB = backprojection(sinogramB,geo)

    # sinogramB_correct = load("sinomat.jld")["sinomat"]
    # imageB_correct = load("imagemat.jld")["imagemat"]

    sino = @inferred projection(phantomImg',geo)
    @test sino isa Matrix
    @test size(sinogramB) == (360, 1024)
    image = @inferred backprojection(sinogramB,geo)
    @test image isa Matrix
    @test size(imageB) == (256, 256)

    # @test sinogramB == sinogramB_correct
    # @test imageB == imageB_correct

    #=
    using MIRTjim: jim

    p1 = jim(1:geo.nDet, 1:length(geo.theta), sinogramB'; aspect_ratio = :auto, clim=(0,68))
    p2 = jim(1:geo.nPix, 1:geo.nPix, imageB'; aspect_ratio = :auto, clim=(30,95))
    jim(p1,p2)
    =#

end


text:
# ---------------------------------------------------------------------------- #
#
#   basisFuncLineLag.jl
#
#   Get Legendre basis functions for Line
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

# function basisFuncLineLeg( ::Type{Val{0}}, s::Array{Float64} )
#
#   dim = 1
#   nx  = length( s )
#
#   # Constant
#   nphi = 1
#   phi  = Array{Float64}(nx, nphi)
#   dphi = Array{Float64}(nx, nphi, dim)
#
#   # value
#   phi[:,1] = 1.0
#
#   # derivative in x
#   dphi[:,1,1] = 0.0
#
#   # transpose
#   phi   = phi'
#   dphi  = dphi[:,:,1]'
#
#   return phi, dphi
#
# end

function basisFuncLineLeg( P::P1, s::Array{Float64} )

  dim = 1
  nx  = length( s )

  # Linear
  nphi = 2
  phi  = Array{Float64}(nx, nphi)
  dphi = Array{Float64}(nx, nphi, dim)

  t = 1 - s

  phi[:,1] = 1
  phi[:,2] = sqrt(3).*(s - t)

  dphi[:,1,1] =  0
  dphi[:,2,1] =  2.*sqrt(3)

  # transpose
  phi   = phi'
  dphi  = dphi[:,:,1]'

  return phi, dphi

end

function basisFuncLineLeg( P::P2, s::Array{Float64} )

  dim = 1
  nx  = length( s )

  # Quadratric
  nphi = 3
  phi  = Array{Float64}(nx, nphi)
  dphi = Array{Float64}(nx, nphi, dim)

  t = 1 - s

  phi[:,1] = 1
  phi[:,2] = sqrt(3).*(s - t)
  phi[:,3] = sqrt(5).*(1 - 6.*s.*t)

  dphi[:,1,1] =  0
  dphi[:,2,1] =  2.*sqrt(3)
  dphi[:,3,1] =  6.*sqrt(5).*(s - t)

  # transpose
  phi   = phi'
  dphi  = dphi[:,:,1]'

  return phi, dphi

end

function basisFuncLineLeg( P::P3, s::Array{Float64} )

  dim = 1
  nx  = length( s )

  # Cubic
  nphi = 4
  phi  = Array{Float64}(nx, nphi)
  dphi = Array{Float64}(nx, nphi, dim)

  t = 1 - s

  phi[:,1] = 1
  phi[:,2] = sqrt(3).*(s - t)
  phi[:,3] = sqrt(5).*(1 - 6.*s.*t)
  phi[:,4] = sqrt(7).*(1 - 10.*s.*t).*(s - t)

  dphi[:,1,1] =  0
  dphi[:,2,1] =  2.*sqrt(3)
  dphi[:,3,1] =  6.*sqrt(5).*(s - t)
  dphi[:,4,1] = 12.*sqrt(7).*(1 - 5.*s.*t)

  # transpose
  phi   = phi'
  dphi  = dphi[:,:,1]'

  return phi, dphi

end

function basisFuncLineLeg( P::P4, s::Array{Float64} )

  dim = 1
  nx  = length( s )

  # Quartic
  nphi = 5
  phi  = Array{Float64}(nx, nphi)
  dphi = Array{Float64}(nx, nphi, dim)

  t = 1 - s

  phi[:,1] = 1
  phi[:,2] = sqrt(3).*(s - t)
  phi[:,3] = sqrt(5).*(1 - 6.*s.*t)
  phi[:,4] = sqrt(7).*(1 - 10.*s.*t).*(s - t)
  phi[:,5] = 3.*(1 - 10.*s.*t.*(2 - 7.*s.*t))

  dphi[:,1,1] =  0
  dphi[:,2,1] =  2.*sqrt(3)
  dphi[:,3,1] =  6.*sqrt(5).*(s - t)
  dphi[:,4,1] = 12.*sqrt(7).*(1 - 5.*s.*t)
  dphi[:,5,1] = 60.*(1 - 7.*s.*t).*(s - t)

  # transpose
  phi   = phi'
  dphi  = dphi[:,:,1]'

  return phi, dphi

end

# function basisFuncLineLeg( P::P5, s::Array{Float64} )
#
#   dim = 1
#   nx  = length( s )
#
#   # Quintic
#   nphi = 6
#   phi  = Array{Float64}(nx, nphi)
#   dphi = Array{Float64}(nx, nphi, dim)
#
#   t = 1 - s
#
#   phi[:,1] = 1;
#   phi[:,2] = sqrt(3).*(s - t);
#   phi[:,3] = sqrt(5).*(1 - 6.*s.*t);
#   phi[:,4] = sqrt(7).*(1 - 10.*s.*t).*(s - t);
#   phi[:,5] = 3.*(1 - 10.*s.*t.*(2 - 7.*s.*t));
#   phi[:,6] = sqrt(11).*(1 - 14.*s.*t.*(s - 2.*t).*(2.*s - t)).*(s - t);
#
#
#   dphi[:,1,1] =  0
#   dphi[:,2,1] =  2.*sqrt(3)
#   dphi[:,3,1] =  6.*sqrt(5).*(s - t)
#   dphi[:,4,1] = 12.*sqrt(7).*(1 - 5.*s.*t)
#   dphi[:,5,1] = 60.*(1 - 7.*s.*t).*(s - t)
#   dphi[:,6,1] = 30.*sqrt(11).*(1 - 14.*s.*t.*(1 - 3.*s.*t))
#
#   # transpose
#   phi   = phi'
#   dphi  = dphi[:,:,1]'
#
#   return phi, dphi
#
# end


text:
using Random, Statistics, ProgressMeter

#2x2-J1-J2-checkboard
function e_dif(model::Ising22, array, lattice, x, y)
    r = model.r
    a , b = mod(x,2), mod(y,2)

    top     = array[x, y + 1 - lattice * (y==lattice)]
    bottom  = array[x,y - 1 + lattice * (y==1)]
    left    = array[x - 1 + lattice * (x==1), y]
    right   = array[x + 1 - lattice * (x==lattice), y]

    if a==1
        right *= r
    else
        left *= r
    end
    if b==1
        top *= r
    else
        bottom *= r
    end

    return 2 * array[x, y] * (top+bottom+left+right)
end

#3x3-J1-J2
function e_dif(model::Ising33, array, lattice, x, y)
    r = model.r
    a , b = mod(x,3), mod(y,3)

    top     = array[x, y + 1 - lattice * (y==lattice)]
    bottom  = array[x,y - 1 + lattice * (y==1)]
    left    = array[x - 1 + lattice * (x==1), y]
    right   = array[x + 1 - lattice * (x==lattice), y]

    if a == 1
        if b == 2
            right *= r
        end
    elseif a == 2
        if b == 1
            top *= r
        elseif b == 2
            top *= r
            bottom *= r
            left *= r
            right *= r
        else
            bottom *= r
        end
    else
        if b == 2
            left *= r
        end
    end

    return 2 * array[x, y] * (top+bottom+left+right)
end

function onestep!(model::MT, spin_array, lattice, β) where {MT <: HamiltonianModel}
    for i = 1:lattice
        for j = 1:lattice
            e = e_dif(model, spin_array, lattice, i, j)
            if e <= 0
                spin_array[i,j] = -spin_array[i,j];
            elseif exp(-e*β) > rand()
                spin_array[i,j] = -spin_array[i,j];
            end
        end
    end
    return spin_array
end

function energy(model::MT, spin_array, lattice) where {MT <: HamiltonianModel}
    tol_en = 0
    for i = 1:lattice
        for j = 1:lattice
            e = e_dif(model, spin_array, lattice, i, j)
            tol_en += e
        end
    end
    return -tol_en/4
end

function MCMC(model::MT,lattice,β,Sweeps_heat,Sweeps) where {MT <: HamiltonianModel}
#     spin_array = ones(lattice,lattice)
    spin_array = (bitrand(lattice,lattice).-0.5)*2
    for j = 1:Sweeps_heat
        spin_array = onestep!(model,spin_array,lattice ,β);
    end
    mag = zeros(1,Sweeps);
    ene = zeros(1,Sweeps)
    Threads.@threads for j = 1:Sweeps
        spin_array = onestep!(model,spin_array,lattice ,β)
        mag[j] = abs(sum(spin_array)/lattice^2)
        ene[j] = energy(model,spin_array, lattice)/lattice^2
    end
    mag_ave = sum(mag)/Sweeps;
    ene_ave = sum(ene)/Sweeps;
    return mag_ave, ene_ave
end

function mutiMC(model::MT,β,lattice,bins,Sweeps_heat,Sweeps) where {MT <: HamiltonianModel}
    eMC_bins = zeros(bins,1)
    Threads.@threads for j=1:bins
        eMC_bins[j] =  MCMC(model,lattice,β,Sweeps_heat,Sweeps)
    end
    eMC = mean(eMC_bins)
    return eMC
end

function MutiMC(model::MT,β,lattice,Bins,bins,Sweeps_heat,Sweeps) where {MT <: HamiltonianModel}
    eMC_Bins = zeros(1,Bins)
    p = Progress(Bins)
    for j = 1:Bins
        eMC_Bins[j] = mutiMC(model,β,lattice,bins,Sweeps_heat,Sweeps)
        update!(p,j)
    end
    eMC_m = mean(eMC_Bins)
    err = 1.96*sqrt(var(eMC_Bins)/Bins)
    return eMC_m,err
end


text:
using DelimitedFiles, ProgressMeter, PyPlot

include("./util.jl")
include("./parsers.jl")
include("./Algorithms/TEAVAR.jl")
include("./simulation.jl")

function topologyCutoffError(topology,
                             num_demands,
                             iterations,
                             cutoffs; paths="SMORE",
                             weibull_scale=.0001,
                             demand_downscale=5000,
                             plot=true,
                             dirname="./data/raw/cutoff_error/")
    env = Gurobi.Env()
    y_vals = []
    scenarios = []
    probs = []
    links, capacity, link_probs, nodes = readTopology("$(topology)")
    demand, flows = readDemand("$(topology)/demand", length(nodes), 1, scale=1, downscale=demand_downscale)
    T, Tf, k = parsePaths("$(topology)/paths/$(paths)", links, flows)

    dir = nextRun(dirname)
    for j in 1:length(cutoffs)
        writedlm("$(dir)/cutoff_$(cutoffs[j])_error", [["cvar_o", "cvar", "error_cvar"]])
    end

    progress = ProgressMeter.Progress(num_demands*iterations*length(cutoffs), .1, "Computing cutoff error...", 50)
    cutoff_errors = zeros(length(cutoffs))
    for d in 1:num_demands
        demand, flows = readDemand("$(topology)/demand", length(nodes), d, scale=1, downscale=demand_downscale)
        for i in 1:iterations
            weibull_probs = weibullProbs(length(link_probs), shape=.8, scale=weibull_scale)
            optimal_scenarios, optimal_probs = subScenarios(weibull_probs, (sum(weibull_probs)/length(weibull_probs))^2, first=true, last=true)
            beta = optimal_probs[1]
            var_o, cvar_o, a_o, max_u_o = TEAVAR(env, links, capacity, flows, demand, beta, k, T, Tf, optimal_scenarios, optimal_probs)
            for j in 1:length(cutoffs)
                scenarios, scenario_probs = subScenarios(weibull_probs, cutoffs[j], first=true, last=true)
                var, cvar, a, max_u = TEAVAR(env, links, capacity, flows, demand, beta, k, T, Tf, scenarios, scenario_probs)
                cvar_o = sum((1 .- max_u_o) .* optimal_probs)
                cvar = sum((1 .- max_u) .* scenario_probs)
                err = abs((cvar_o - cvar)/cvar_o)
                cutoff_errors[j] += err
                open("$(dir)/cutoff_$(cutoffs[j])_error", "a") do io
                    writedlm(io, transpose([cvar_o, cvar, err]))
                end
                ProgressMeter.next!(progress, showvalues = [(:topology,topology), (:demand,"$(d)/$(num_demands)"), (:iterations,"$(i)/$(iterations)"), (:cutoff, cutoffs[j]), (:error, err)])
            end
        end
    end
    y_vals = cutoff_errors ./ (num_demands * iterations)
    writedlm("$(dir)/cutoffs", cutoffs)
    writedlm("$(dir)/y_vals", y_vals)
    writedlm("$(dir)/params", [["topology", "num_demands", "iterations", "cutoffs", "paths", "weibull_scale", "demand_downscale"], [topology, num_demands, iterations, cutoffs, paths, weibull_scale, demand_downscale]])

    if plot
        PyPlot.clf()
        nbars = length(cutoffs)
        barWidth = 1/(nbars)
        for i in 1:length(y_vals)
            PyPlot.bar(barWidth * i, y_vals[i], alpha=0.8, width=barWidth)
        end
        PyPlot.ylabel("Percent Error", fontweight="bold")
        PyPlot.legend(cutoffs, loc="upper right")
        PyPlot.savefig("$(dir)/plot.png")
        PyPlot.show()
    end
end


text:
module SpringCollab2020PauseBadelineBossesTrigger

using ..Ahorn, Maple

@mapdef Trigger "SpringCollab2020/PauseBadelineBossesTrigger" PauseBadelineBossesTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight)

const placements = Ahorn.PlacementDict(
    "Pause Badeline Bosses (Spring Collab 2020)" => Ahorn.EntityPlacement(
        PauseBadelineBossesTrigger,
        "rectangle",
    ),
)

end


text:
module TestScenes

using ..GfxBase
using ..Scenes
using ..Materials
using ..WWUMeshes
using ..Cameras

make_diffuse(color) = Material(0.0, 0.0, color)

black = Vec3(0,0,0)
red = Vec3(1,0,0)
green = Vec3(0,1,0)
blue = Vec3(0,0,1)
white = Vec3(1,1,1)
purple = Vec3(1,0,1)

""" Return an Array of all Triangles belonging to the given mesh, assigning
each one the given material. """
function create_triangles(mesh::OBJMesh, material)
    triangles = []

    for t in mesh.triangles
        positions = [mesh.positions[t.positions[1]], mesh.positions[t.positions[2]], mesh.positions[t.positions[3]]]

        normals = nothing
        if t.normals !== nothing && length(t.normals) == 3
            normals = [mesh.normals[t.normals[1]], mesh.normals[t.normals[2]], mesh.normals[t.normals[3]]]
        end

        uvs = nothing
        if t.uvs !== nothing && length(t.uvs) == 3
            uvs = [mesh.uvs[t.uvs[1]], mesh.uvs[t.uvs[2]], mesh.uvs[t.uvs[3]]]
        end

        push!(triangles, Triangle(positions, normals, uvs, material))
    end

    triangles
end

function mesh_helper(mesh, material, scale=Vec3(1.0,1.0,1.0), translation=Vec3(0,0,0), rotation=Vec3(0, 0, 0))
    Rx = [1 0                0
          0 cos(rotation[1]) -sin(rotation[1])
          0 sin(rotation[1]) cos(rotation[1])]

    Ry = [cos(rotation[2])  0 sin(rotation[2])
          0                 1 0
          -sin(rotation[2]) 0 cos(rotation[2])]

    Rz = [cos(rotation[3]) -sin(rotation[3]) 0
          sin(rotation[3]) cos(rotation[3])  0
          0                0                 1]

    R = Rx * Ry * Rz

    for i in 1:length(mesh.positions)
        mesh.positions[i] = (R * mesh.positions[i]) .* scale + translation
    end

    create_triangles(mesh, material)
end

function camera_1(img_height, img_width)
    CanonicalCamera(img_height, img_width)
end

function camera_2(img_height, img_width)
    eye = Vec3(20, 4, 10)
    view = Vec3(-1, 0, -5) - eye
    up = Vec3(0, 1, 0)
    focal = 8.0
    Cameras.PerspectiveCamera(eye, view, up, focal, img_height, img_width)
end

function camera_3(img_height, img_width)
    Cameras.PerspectiveCamera(Vec3(-1, 0.8, -1.2), Vec3(1, -1, -1), Vec3(0, 1, 0), 0.3, img_height, img_width)
end

cameras = [camera_1, camera_2, camera_3]

function get_camera(i, img_height, img_width)
    cameras[i](img_height, img_width)
end

function bunny()
    bunny = read_obj("data/bunny.obj")
    bunny = mesh_helper(bunny, 1, 1.0, Vec3(0.2, 0, -5))

    cube = mesh_helper(cube_mesh(), 2, 10.0, Vec3(-11.2, 0, 0))

    lights = [ Light(Float32(1.3) * white, Vec3(1,2,-5), 1),
               Light(Float32(0.8) * white, Vec3(0,0,10), 3),
               Light(Float32(0.8) * white, Vec3(0,10,10), 3),
               Light(Float32(0.8) * white, Vec3(10,10,10), 3),
               Light(Float32(0.8) * white, Vec3(0,10,0), 3) ]

    materials = [Material(0.0, 0.0, Vec3(0.6, 0.5, 0.5)), Material(0.6, 0.0, white)]
    triangles = vcat(cube, bunny)
    Scene(black, materials, triangles, lights, bvh(triangles))
end

function teapot()
    teapot = read_obj("data/teapot.obj")
    estimate_normals(teapot)
    teapot = mesh_helper(teapot, 1, 0.5, Vec3(0, -1, -5))

    cube = mesh_helper(cube_mesh(), 1, 10.0, Vec3(0, -11, 0))

    lights = [ Light(4 * white, Vec3(1,2,-5), 1) ]

    materials = [Material(0.0, 0.0, white)]
    triangles = vcat(cube, teapot)
    Scene(black, materials, triangles, lights, bvh(triangles))
end

function cornell()
    cube_mat = Material(0, 0.0, .7 * white)

    cube_back = mesh_helper(cube_mesh(), 1, 2.0, Vec3(0, 0, -8))
    cube_left = mesh_helper(cube_mesh(), 2, 2.0, Vec3(-4, 0, -6))
    cube_right = mesh_helper(cube_mesh(), 3, 2.0, Vec3(4, 0, -6))
    cube_top = mesh_helper(cube_mesh(), 1, 2.0, Vec3(0, 4, -6))
    cube_bottom = mesh_helper(cube_mesh(), 1, 2.0, Vec3(0, -4, -6))

    cube_1 = mesh_helper(cube_mesh(), 1, Vec3(0.5, 1.0, 0.5), Vec3(-0.75, -1, -5), Vec3(0, pi / 8, 0))
    cube_2 = mesh_helper(cube_mesh(), 1, 0.5, Vec3(0.75, -1.5, -5), Vec3(0, -pi / 8, 0))

    lights = [ Light(4 * white, Vec3(0,1.49,-4), 0.5) ]

    materials = [Material(0, 0.0, .7 * white), Material(0, 0.0, .7 * red), Material(0, 0.0, .7 * green)]
    triangles = vcat(cube_back, cube_left, cube_right, cube_top, cube_bottom, cube_1, cube_2)
    Scene(.7 * white, materials, triangles, lights, bvh(triangles))
end

function cubes()
    cube_bottom = mesh_helper(cube_mesh(), make_diffuse(white), 6.0, Vec3(0, -9, -12))
    cube_back = mesh_helper(cube_mesh(), make_diffuse(blue), 6.0, Vec3(0, 0, -16))
    cube_opaque = mesh_helper(cube_mesh(), make_diffuse(green), 0.5, Vec3(-2, -1, -7))
    cube_semi_transparent = mesh_helper(cube_mesh(), Material(0.0, 0.75, white), 0.5, Vec3(0, -1, -7))

    lights = [ Light(4 * white, Vec3(0, 10, 10), 3)]

    triangles = vcat(cube_bottom, cube_back, cube_opaque, cube_semi_transparent)
    Scene(black, triangles, lights, bvh(triangles))
end

function glass()
    glass_mat = Material(0, 0.5, white)
    glass = read_obj("data/glass.obj")
    estimate_normals(glass)
    glass = mesh_helper(glass, glass_mat, 0.01, Vec3(0, 0.1, -3))

    cube_mat = Material(0.6, 0.0, white)
    cube = mesh_helper(cube_mesh(), cube_mat, 10.0, Vec3(0, -11, 0))

    lights = [ Light(4 * white, Vec3(1,2,-2), 1) ]

    triangles = vcat(cube, glass)
    Scene(black, triangles, lights, bvh(triangles))
end

function get_scene(name)
    scenes = Dict("bunny" => bunny(), "cornell" => cornell(), "teapot" => teapot())
    scenes[name]
end

end

max_stars_repo_name: aedans/raytracer

max_issues_repo_name: aedans/raytracer

max_forks_repo_name: aedans/raytracer


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Mesh{Dim,T}

A mesh embedded in a `Dim`-dimensional space with coordinates of type `T`.
"""
abstract type Mesh{Dim,T} <: Domain{Dim,T} end

"""
    topology(mesh)

Return the topological structure of the `mesh`.
"""
function topology end

"""
    vertices(mesh)

Return the vertices of the `mesh`.
"""
function vertices end

"""
    nvertices(mesh)

Return the number of vertices of the `mesh`.
"""
nvertices(m::Mesh) = nvertices(topology(m))

"""
    faces(mesh, rank)

Return an iterator with `rank`-faces of the `mesh`.

## Examples

Consider a mesh of tetrahedrons embedded in a 3D space. We can loop over
all 3-faces (a.k.a. elements) or over all 2-faces to handle the interfaces
(i.e. triangles) between adjacent elements:

```julia
tetrahedrons = faces(mesh, 3)
triangles    = faces(mesh, 2)
segments     = faces(mesh, 1)
```
"""
faces(m::Mesh, rank) = (materialize(f, vertices(m)) for f in faces(topology(m), rank))

"""
    nfaces(mesh, rank)

Return the number of `rank`-faces of the `mesh`.
"""
nfaces(m::Mesh, rank) = nfaces(topology(m), rank)

"""
    elements(mesh)

Return the top-faces (a.k.a. elements) of the `mesh`.

## Examples

The elements of a volume embedded in 3D space can be tetrahedrons, hexahedrons,
or any 3-face. The elements of a surface embedded in 3D space can be triangles,
quadrangles or any 2-face.
"""
elements(m::Mesh) = (element(m, ind) for ind in 1:nelements(m))

"""
    element(mesh, ind)

Return the element of the `mesh` at index `ind`.
"""
element(m::Mesh, ind) = materialize(element(topology(m), ind), vertices(m))

"""
    nelements(mesh)

Return the number of elements of the `mesh`.
"""
nelements(m::Mesh) = nelements(topology(m))

"""
    facets(mesh)

Return the (top-1)-faces (a.k.a. facets) of the `mesh`.
"""
facets(m::Mesh) = (facet(m, ind) for ind in 1:nfacets(m))

"""
    facet(mesh, ind)

Return the facet of the `mesh` at index `ind`.
"""
facet(m::Mesh, ind) = materialize(facet(topology(m), ind), vertices(m))

"""
    nfacets(mesh)

Return the number of facets of the `mesh`.
"""
nfacets(m::Mesh) = nfacets(topology(m))

"""
    topoconvert(T, mesh)

Convert underlying topology of the `mesh` to topology of type `T`.

## Examples

Convert underlying topology to [`HalfEdgeTopology`](@ref) for
efficient topological relations.

```julia
newmesh = topoconvert(HalfEdgeTopology, mesh)
```
"""
topoconvert(TP::Type{<:Topology}, m::Mesh) =
  SimpleMesh(vertices(m), convert(TP, topology(m)))

"""
    merge(mesh₁, mesh₂)

Merge `mesh₁` with `mesh₂`, i.e. concatenate the vertices
and adjust the connectivities accordingly.
"""
function Base.merge(m₁::Mesh, m₂::Mesh)
  t₁ = topology(m₁)
  t₂ = topology(m₂)
  v₁ = vertices(m₁)
  v₂ = vertices(m₂)
  e₁ = collect(elements(t₁))
  e₂ = collect(elements(t₂))

  # concatenate vertices
  points = [v₁; v₂]

  # concatenate connectivities
  pT = promote_type(eltype(e₁), eltype(e₂))
  connec = pT[]
  append!(connec, e₁)
  offset = length(v₁)
  for e in e₂
    c  = indices(e)
    c′ = ntuple(i -> c[i] + offset, length(c))
    push!(connec, connect(c′))
  end

  SimpleMesh(points, connec)
end

==(m₁::Mesh, m₂::Mesh) =
  vertices(m₁) == vertices(m₂) &&
  topology(m₁) == topology(m₂)

# -----------
# IO METHODS
# -----------

function Base.show(io::IO, ::MIME"text/plain", m::Mesh{Dim,T}) where {Dim,T}
  t = topology(m)
  verts = vertices(m)
  elems = elements(t)
  nvert = nvertices(m)
  nelms = nelements(m)
  println(io, m)
  println(io, "  $nvert vertices")
  println(io, io_lines(verts, "    "))
  println(io, "  $nelms elements")
  print(  io, io_lines(elems, "    "))
end

# ----------------
# IMPLEMENTATIONS
# ----------------

include("mesh/cartesiangrid.jl")
include("mesh/simplemesh.jl")


text:
# Checks to make sure the fields in the datatype are of type Parameter{T} or Parameter{U} where
# T<:AbstractFloat and U<:Integer
function fieldparameters(obj::DataType)
    fields = fieldnames(obj)
    0 < length(obj.parameters) <= 2 || error("Type must have one or two parameters (T & U)")
    for param in obj.parameters
        if param.name == :T
            param.ub == AbstractFloat || error("Parameter T must be T<:AbstractFloat")
        elseif param.name == :U
            param.ub == Integer || error("Parameter U must be U<:Integer")
        else
            error("Parameter $(param.name) not recognized; must be T<:AbstactFloat or U<:Integer")
        end
    end
    for field in fields
        field_type = fieldtype(obj, field)
        isa(field_type, DataType) || error("Field $field must be ::Parameter")
        if field_type.name.name != :HyperParameter
            error("Fields must consist of Parameter types")
        end
        length(field_type.parameters) == 1 || error("Parameter type must have one type parameter")
        field_param = field_type.parameters[1]
        if field_param.name == :T
            if !(field_param.ub == AbstractFloat)
                error("Field $field must be ::Parameter{T<:AbstractFloat}")
            end
        elseif field_param.name == :U
            if !(field_param.ub == Integer)
                error("Field $field must be ::Parameter{U<:Integer}")
            end
        else
            error("$field type must be ::Parameter{T<:AbstractFloat} or ::Parameter{U<:Integer}")
        end
    end
    field_parameters = Symbol[fieldtype(obj, field).parameters[1].name for field in fields]
    return (fields, field_parameters)
end

function generate_conversions(obj::DataType)
    fields, field_parameters = fieldparameters(obj)
    obj_sym = obj.name.name
    if length(obj.parameters) == 2
        # convert(::Type{obj{T,U}}, ::obj) = ...
        convert_target = Expr(:(::), Expr(:curly, :Type, Expr(:curly, obj_sym, :T, :U)))
        converted_arguments =
            Expr[:(Variable(convert($(field_parameters[i]), obj.$(fields[i]).value),
                   obj.$(fields[i]).isfixed)) for i in eachindex(fields)]
        conversion1 = Expr(:(=), Expr(:call, Expr(:curly, :convert, :T, :U), convert_target,
                                             Expr(:(::), :obj, obj_sym)),
                                 Expr(:call, obj_sym, converted_arguments...))
        # convert(::Type{obj{T}}, ::obj{_,U}) = ...
        convert_target = Expr(:(::), Expr(:curly, :Type, Expr(:curly, obj_sym, :T)))
        convert_source = Expr(:(::), Expr(:curly, :Type, Expr(:curly, obj_sym, :_, :U)))
        conversion2 = Expr(:(=), Expr(:call, Expr(:curly, :convert, :T, :_, :U), convert_target,
                                             Expr(:(::), :obj, Expr(:curly, obj_sym, :_, :U))),
                                 Expr(:call, obj_sym, converted_arguments...))
        return Expr(:block, conversion1, conversion2)
    elseif length(obj.parameters) == 1
        # convert(::Type{obj{T}}, ::obj) = ...
        convert_target      = Expr(:(::), Expr(:curly, :Type, Expr(:curly, obj_sym, :T)))
        converted_arguments = Expr[:(Variable(convert(T, obj.$(fields[i]).value),
                                     obj.$(fields[i]).isfixed)) for i in eachindex(fields)]
        return Expr(:(=), Expr(:call, Expr(:curly, :convert, :T), convert_target,
                                      Expr(:(::), :obj, obj_sym)),
                          Expr(:call, obj_sym, converted_arguments...))
    else
        error("Data type should have 1 or 2 parameters")
    end
end

# Notes:
#   Assumes arguments are organized by type parameter
#   Assumes Variable{} is the argument type
function generate_outer_constructor(obj::DataType, defaults::Tuple{Vararg{Real}})
    fields, parameters = fieldparameters(obj)
    length(defaults) == length(fields) || error("Default count does not match field count")
    first_idx  = Dict(findfirst(parameters, :T) => :Float64, findfirst(parameters, :U) => :Int64)
    # Produces [:Float64, T, T, ..., Int64, U, U, ...]
    default_parameters = Symbol[get(first_idx, i, parameters[i]) for i in eachindex(parameters)]
    # Produces [arg1::Argument{T} = convert(Float64, default1),
    #           arg2::Argument{T} = convert(T,default2)...]
    constructor_params = [Expr(:(<:), p.name, p.ub.name.name) for p in obj.parameters]
    constructor_args   = [Expr(:kw, :($(fields[i])::Argument{$(parameters[i])}),
                                    :(convert($(default_parameters[i]), $(defaults[i]))))
                          for i in eachindex(fields)]
    type_params = [p.name for p in obj.parameters]
    type_args   = [:(Variable($(fields[i]))) for i in eachindex(fields)]

    obj_sym = obj.name.name
    Expr(:(=), Expr(:call, Expr(:curly, obj_sym, constructor_params...), constructor_args...),
               Expr(:call, Expr(:curly, obj_sym, type_params...), type_args...))
end

macro outer_constructor(obj, defaults)
    eval(generate_outer_constructor(eval(obj), eval(defaults)))
    eval(generate_conversions(eval(obj)))
end


text:
# =======================================
# Flowpipe composition with a time-shift
# =======================================

"""
    ShiftedFlowpipe{FT<:AbstractFlowpipe, NT<:Number} <: AbstractFlowpipe

Type that lazily represents a flowpipe that has been shifted in time.

### Fields

- `F`  -- original flowpipe
- `t0` -- time shift

### Notes

This type can wrap any concrete subtype of `AbstractFlowpipe`, and the extra
field `t0` is such that the time spans of each reach-set in `F` are shifted
by the amount `t0` (which should be a subtype of `Number`).

A convenience constructor alias `Shift` is given.
"""
struct ShiftedFlowpipe{FT<:AbstractFlowpipe, NT<:Number} <: AbstractFlowpipe
    F::FT
    t0::NT
end

function ShiftedFlowpipe(vec::AbstractVector{<:AbstractLazyReachSet}, t0::Number)
    return ShiftedFlowpipe(Flowpipe(vec), t0)
end

# getter functions
@inline array(fp::ShiftedFlowpipe) = array(fp.F)
@inline flowpipe(fp::ShiftedFlowpipe) = fp.F
@inline time_shift(fp::ShiftedFlowpipe) = fp.t0

# time domain interface
@inline tstart(fp::ShiftedFlowpipe) = tstart(first(fp)) + fp.t0
@inline tend(fp::ShiftedFlowpipe) = tend(last(fp)) + fp.t0
@inline tspan(fp::ShiftedFlowpipe) = TimeInterval(tstart(fp), tend(fp))

@inline tstart(fp::ShiftedFlowpipe, i::Int) = tstart(fp.F[i]) + fp.t0
@inline tend(fp::ShiftedFlowpipe, i::Int) = tend(fp.F[i]) + fp.t0
@inline tspan(fp::ShiftedFlowpipe, i::Int) = TimeInterval(tstart(fp, i), tend(fp, i))

# TODO use interface?
project(fp::ShiftedFlowpipe, vars::AbstractVector) = project(fp, Tuple(vars))
project(fp::ShiftedFlowpipe; vars) = project(fp, Tuple(vars))

function project(fp::ShiftedFlowpipe, vars::NTuple{D, T}) where {D, T<:Integer}
    Xk = array(fp)
    # TODO: use projection of the reachsets
    if 0 ∈ vars # projection includes "time"
        # we shift the vars indices by one as we take the Cartesian prod with the time spans
        aux = vars .+ 1
        t0 = time_shift(fp)
        return map(X -> project(convert(Interval, tspan(X) + t0) × set(X), aux), Xk)
    else
        return map(X -> project(set(X), vars), Xk)
    end
end

# this method is analogue to project(::AbstractLazyReachSet, vars; check_vars=true)
# TODO add check_vars ?
function project(fp::ShiftedFlowpipe, i::Int, vars::NTuple{D, M}) where {D, M<:Integer}
    t0 = time_shift(fp)
    R = fp[i]
    if 0 ∈ vars
        # if the projection involves "time", we shift the vars indices by one as
        # we will take the Cartesian product of the reach-set with the time interval
        aux = vars .+ 1

        Δt = convert(Interval, tspan(R) + t0)
        proj =  project(Δt × set(R), aux)
    else
        proj = project(set(R), vars)
    end

    return SparseReachSet(proj, tspan(R) + t0, vars)
end

# TODO: improve using mutable ShiftedFlowpipe struct (so that we can modify t0->t0+t1)
function shift(fp::ShiftedFlowpipe, t1::Number)
    return ShiftedFlowpipe(shift(fp.F, t1), fp.t0)
end


text:
const __MAIN__ = MainCommand()

function make(main::String)
    if main != name(__MAIN__)
        __MAIN__.name = main
    end
end

function make(name::Symbol, cmd::T) where {T <: Union{MainCommand, NodeCommand}}
    name = string(name)
    leaf = findcmd(name, cmd.cmds)
    if leaf === nothing
        leaf = NodeCommand(name)
        push!(cmd.cmds, leaf)
    end
    return leaf
end

function make(name::NTuple{N, Symbol}, cmd::T) where {N, T <: Union{MainCommand, NodeCommand}}
    child = make(first(name), cmd)
    make(Base.tail(name), child)
end

function make(name::Tuple{}, cmd::T) where {T <: Union{MainCommand, NodeCommand}}
    return cmd
end

isparameter(expr::Expr) = expr.head == :parameters

# TODO: support user-defined short/long options
function getoption(expr::Expr)
    key = expr.args[1].args[1]
    dtype = eval(expr.args[1].args[2])
    ShortOption(dtype, key)
end

function gettype(expr::Expr)
    eval(expr.args[2])
end

function parameter(expr::Expr)
    kwds = AbstractOption[]
    for each in expr.args
        push!(kwds, getoption(each))
    end
    kwds
end

function arguments(expr::Vector)
    args = DataType[]
    for each in expr
        push!(args, gettype(each))
    end
    args
end

function signature(expr::Expr)
    args = DataType[]; kwds = AbstractOption[]

    if isa(expr.args[2], Expr) && isparameter(expr.args[2])
        kwds = parameter(expr.args[2])
        args = arguments(expr.args[3:end])
    else
        args = arguments(expr.args[2:end])
    end

    Signature(args, kwds)
end

import Base.Docs: Binding, DocStr, MultiDoc, parsedoc, catdoc, meta, getdoc, defined

function commanddoc(var::Symbol)
    results, groups = DocStr[], MultiDoc[]
    m = current_module()
    binding = Binding(m, var)

    if defined(binding)
        result = getdoc(binding)
        result === nothing || return result
    end

    dict = meta(m)
    if haskey(dict, binding)
        multidoc = dict[binding]
        push!(groups, multidoc)

        for msig in multidoc.order
            push!(results, multidoc.docs[msig])
        end
    end

    if isempty(groups)
        return "No Documents"
    end

    md = catdoc(map(parsedoc, results)...)

    if isa(md, Markdown.MD)
        md.meta[:results] = results
        md.meta[:binding] = binding
        md.meta[:typesig] = Union{}
    end
    return md
end

"""
    @command "..."

add docstring (description) to main command/node command

    @command

add a function to the CLI tree, like:

    @command function add(x, y)
        x + y
    end

or just

    @command add(x, y) = x + y

```julia
module math
@command "math operations"

module add
@command "add numbers"

@command function plus(x::Int, y::Int)
    # some code for plus
end

end

end
```

and when this CLI is executed in command line with `-h`/`--help`,
will looks like

```shell
> math --help
[SYNOPSIS]
	math [-h,--help] [-v,--version] <command>
[DESCRIPTION]
	math operations

[COMMANDS]
	add
	  add numbers
```
"""
:(@command)

function module_filename(m::Module)
    basename(m.eval(parse("@__FILE__")))
end

macro command(doc::String)
    m = current_module()
    nodes = fullname(m)

    if isempty(nodes)
        name = module_filename(m)
        make(name)
        __MAIN__.doc = doc
    else
        make(string(first(nodes)))
        tail = make(Base.tail(nodes), __MAIN__)
        tail.doc = doc
    end

end

import Core: @__doc__

macro command(expr::Expr)
    m = current_module()
    nodes = fullname(m)

    if !isempty(nodes)
        make(string(first(nodes)))
        tail = make(Base.tail(nodes), __MAIN__)
    else
        tail = __MAIN__
    end

    # insert leaf
    name = expr.args[1].args[1]

    sig = signature(expr.args[1])

    namestr = string(name)
    quote
        # f = $(esc(expr))
        @__doc__ $(esc(expr))
        f = $(esc(name))
        docstring = commanddoc(Symbol($namestr))
        push!($(tail.cmds), LeafCommand(Symbol($namestr), f, $sig, docstring))
    end
end

export @command


text:
################################################################################
#
#          Misc/BigIntInterface.jl : BigInt Handling and Conversions
#
################################################################################j

################################################################################
#
#  Inline Bitwise Assignment Methods for BigInt
#
################################################################################

@inline function assign_zero() end::BigInt
@inline function assign_one() end
@inline function assign(Union{::Int64, ::Int128,
    ::BigInt, ::Float32, ::Float64}) end

################################################################################
#
#  Conversions
#
################################################################################

# Convert BigInt to Hex Array
function int_to_bytes(x::BigInt)
    n_bytes_with_zeros = x.size*sizeof(Sys.WORD_SIZE)
    uint8_ptr = convert(Ptr{UInt8}, x.d)
    n_bytes_without_zeros = 1

    if ENDIAN_BOM == 0x04030201
        for i in n_bytes_with_zeros:-1:1
            if unsafe_load(uint8_ptr, i) != 0x00
                n_bytes_without_zeros = i
                break
            end
        end

        result = Array{UInt8}(undef, n_bytes_without_zeros)

        for i in 1:n_bytes_without_zeros
            @inbounds result[n_bytes_without_zeros + 1 - i] = unsafe_load(uint8_ptr, i)
        end
    else
        for i in 1:n_bytes_with_zeros
            if unsafe_load(uint8_ptr, i) != 0x00
                n_bytes_without_zeros = i
                break
            end
        end

        result = Array{UInt8}(undef, n_bytes_without_zeros)

        for i in 1:n_bytes_without_zeros
            @inbounds result[i] = unsafe_load(uint8_ptr, i)
        end
    end
    return result
end

# Convert Hex Array to BigInt
function bytes_to_big(x::Array{UInt8, 1})
    hex = bytes2hex(x)
    return parse(BigInt, hex, base=16)
end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Data: 05/03/2014

# make a string from value by "%g" format
# it was implemented since Julia did not support "%g" format.
# this function can be obsolete if Julia introduces the format.

function outputGFormatString(v::Real)

    num_str = string(v)

    temp = split(num_str, 'e')

    digits = ""
    exponent = ""

    if length(temp) == 1
        digits = temp[1]
    elseif length(temp) == 2
        digits = temp[1]
        exponent = temp[2]
    else
        error("not a valid number: ", v)
    end

    b_dot = false
    sig_digit = 0

    i = 1
    for i = 1:length(digits)
        if digits[i] == '.'
            if sig_digit >= 6
                i = i - 1
                break
            end

            b_dot = true
        elseif isdigit(digits[i])
            if sig_digit != 0 || digits[i] != '0'
                sig_digit += 1

                if b_dot && sig_digit == 6
                    break
                end
            end
        elseif digits[i] == '-'
        else
            error("not a valid number: ", v)
        end
    end

    end_ind = i

    if b_dot
        for j = 1:i
            k = i - j + 1

            if digits[k] == '.'
                end_ind = k - 1
                break
            elseif digits[k] != '0'
                end_ind = k
                break
            end

        end
    end

    if exponent == ""
        return digits[1:end_ind]
    else
        return digits[1:end_ind] * "e" * exponent
    end
end

# validate output_g_format_string()
#
# validation step
# 1. generate a sample file without output_g_format_string()
# 2. generate a sample file with output_g_format_string()
# 3. compare values of two files

function validate_output_g_format_string(sample_filename_1, sample_filename_2)

    f1 = open(sample_filename_1, "r")
    f2 = open(sample_filename_2, "r")

    readline(f1)
    readline(f2)

    for line1 = eachline(f1)
        line2 = readline(f2)

        values1 = float(split(chomp(line1)))
        values2 = float(split(chomp(line2)))

        for i = 1:length(values1)
            if abs((values1[i] - values2[i]) / values2[i] * 100) > 0.01
                println("line 1: ", chomp(line1))
                println("line 2: ", chomp(line2))
                println("i: ", i)
                println("value 1: ", values1[i])
                println("value 2: ", values2[i])
                error("Invalid")
            end
        end
    end

    close(f1)
    close(f2)
end


text:
using DrWatson
@quickactivate :FractalDimension # uses DynamicalSystems, PyPlot
using DynamicalSystems, PyPlot, Statistics
include(srcdir("style.jl"))

# Plot dependence on time and i.c.
function plot_return_dim(ds, εs, ic, e; Tmax = 20_000_000, ici=1)
    fig, axs = subplots(2, 1; sharex = true)
    mcts = []
    for (i, T) ∈ enumerate(ic)
        println("ic $i ...")
        u0 = trajectory(ds, T; dt = T, Ttr = 100)[end]
        n, c = mean_return_times(ds, u0, εs, Tmax)
        j = findfirst(isequal(0), c)
        j = isnothing(j) ? length(c) : j-1
        push!(mcts, n[1:j])
        o, d = linreg(log.(εs[1:j]), log.(n[1:j]))
        axs[1].plot(log.(εs), log.(n);
                    marker = "o", ls="-", lw = 1.0, ms = 5, color = "C$(i)",
                    label = "i.c. $(i), \$d = $(rdspl(abs(d)))\$"
        )
    end
    for (i, T) ∈ enumerate(e)
        println("Time $i ...")
        u0 = trajectory(ds, ic[ici]; Ttr = 100)[end]
        n, c = mean_return_times(ds, u0, εs, 10^T)
        j = findfirst(isequal(0), c)
        j = isnothing(j) ? length(c) : j-1
        o, d = linreg(log.(εs[1:j]), log.(n[1:j]))
        axs[2].plot(log.(εs), log.(n);
                    marker = "o", ls="-", lw = 1.0, ms = 5, color = lighten_color("C0", 1.2i),
                    label = "\$T=10^$(T), d = $(rdspl(abs(d)))\$"
        )
    end

    minlength = minimum(length.(mcts))
    average_mct = [mean(mcts[j][i] for j in 1:length(ic)) for i in 1:minlength]
    x, y = log.(εs[1:minlength]), log.(average_mct)
    _, d = linreg(x, -y)
    # _, d = linear_region(x, -y) # this doesn't work well. The line is too jagged.
    axs[1].plot(x, y; ls = "dashed", label = "avg., \$d=$(rdspl(d))\$", lw = 3.0)

    axs[1].legend(); axs[2].legend()
    axs[1].set_ylabel("log(⟨t⟩)")
    axs[2].set_ylabel("log(⟨t⟩)")
    axs[2].set_xlabel("log(ε)")
    fig.tight_layout()
    return fig, axs
end

# %% Henon map
# first plot random initial conditions
ds = Systems.henon()
ic = (500, 1000, 1500, 2000, 2500)
e = (3, 5, 7)
εs = sort!(ℯ .^ (-9:0.5:-3); rev = true)
fig, axs = plot_return_dim(ds, εs, ic, e)
axs[1].set_title("Henon map")

tr = trajectory(ds, 50000; Ttr = 100)
figure()
scatter(columns(tr)..., s = 2)
for (i, T) ∈ enumerate(ic)
    u0 = trajectory(ds, T; Ttr = 100)[end]
    scatter(u0[1], u0[2], color = "C$i")
end

# wsave(plotsdir("returntime", "henon"), fig)

# %% Towel map
# first plot random initial conditions
ds = Systems.towel()
εs = sort!(ℯ .^ (-7:0.5:-3); rev = true)
ic = (500, 1000, 1500)
e = (4, 6, 8)

tr = trajectory(ds, 50000; Ttr = 100)
using3D(); figure();
scatter3D(columns(tr)...; s = 2)
for (i, T) ∈ enumerate(ic)
    u0 = trajectory(ds, T; Ttr = 100)[end]
    scatter3D([u0[1]], [u0[2]], [u0[3]]; color = "C$i", s = 1000)
end

fig, axs = plot_return_dim(ds, εs, ic, e)
axs[1].set_title("Towel map")
# wsave(plotsdir("returntime", "towel"), fig)

# %% Roessler system
# first plot random initial conditions
ds = Systems.roessler([1,1,1.0])
εs = sort!(ℯ .^ (-4:0.5:-1); rev = true)
ic = (563.0, 1300.0, 1000.0, 1500.0, 222.0)
e = (4, 5, 6)

tr = trajectory(ds, 1000; Ttr = 100)
using3D(); figure();
plot3D(columns(tr)...; lw = 0.5)
for (i, T) ∈ enumerate(ic)
    u0 = trajectory(ds, T; Ttr = 100)[end]
    scatter3D([u0[1]], [u0[2]], [u0[3]]; color = "C$i", s = 1000)
end

# fig, axs = plot_return_dim(ds, εs, ic, e; Tmax = 10^6, ici=2)
# axs[1].set_title("Roessler system")
# wsave(plotsdir("returntime", "roessler"), fig)

# %% Roessler system
# first plot random initial conditions
ds = Systems.henonheiles([0.0, -0.25, 0.42081, 0.0])
εs = sort!(ℯ .^ (-4:0.5:-1); rev = true)
ic = (563.0, 1300.0, 1002.0, 1500.5, 222.0)
e = (4, 5, 6)

tr = trajectory(ds, 1000; Ttr = 100)
using3D(); figure();
x, y, z = columns(tr)
plot3D(x, y, z; lw = 0.5)
for (i, T) ∈ enumerate(ic)
    u0 = trajectory(ds, T; Ttr = 100)[end]
    scatter3D([u0[1]], [u0[2]], [u0[3]]; color = "C$i", s = 1000)
end

# fig, axs = plot_return_dim(ds, εs, ic, e; Tmax = 10^6, ici=2)
# axs[1].set_title("Roessler system")
# wsave(plotsdir("returntime", "roessler"), fig)

# %% Lorenz system
# first plot random initial conditions
ds = Systems.lorenz([1,10.0,1.0])
εs = sort!(ℯ .^ (-5:0.5:-0); rev = true)
ic = (563, 1000, 1500.5)
e = (4, 5, 6)

tr = trajectory(ds, 500; Ttr = 100, dt = 0.02)
using3D(); figure();
plot3D(columns(tr)...; lw = 1.0)
for (i, T) ∈ enumerate(ic)
    u0 = trajectory(ds, T; Ttr = 100)[end]
    scatter3D([u0[1]], [u0[2]], [u0[3]]; color = "C$i", s = 1000)
end

fig, axs = plot_return_dim(ds, εs, ic, e; Tmax = 10^5, ici=2)
axs[1].set_title("Lorenz system")
# wsave(plotsdir("returntime", "lorenz"), fig)

# %% Lorenz-96 system
D = 4
ds = Systems.lorenz96(D, range(0; length = D, step = 0.1); F = 24.0)
εs = sort!(ℯ .^ (-5:0.5:-0); rev = true)
ic = (26.0, 563.0, 876.8, 1000.0)
e = (4, 5, 6)

# TODO: Must increase max return time.

tr = trajectory(ds, 50; Ttr = 100.0, dt = 0.01)
using3D(); figure();
x, y, z = columns(tr)
plot3D(x, y, z; lw = 0.5)
for (i, T) ∈ enumerate(ic)
    u0 = trajectory(ds, T; Ttr = 100)[end]
    scatter3D([u0[1]], [u0[2]], [u0[3]]; color = "C$i", s = 1000)
end

fig, axs = plot_return_dim(ds, εs, ic, e; Tmax = 10^5, ici=2)
axs[1].set_title("Lorenz96, D = $D")
wsave(plotsdir("returntime", "lorenz96_D=$D"), fig)

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
struct MyOperator{T}
    n::Int
    kind::Symbol
end

Base.axes(O::MyOperator) = (Base.OneTo(O.n),Base.OneTo(O.n))
Base.axes(O::MyOperator,i) = axes(O)[i]
Base.size(O::MyOperator) = (O.n,O.n)
Base.eltype(::MyOperator{T}) where T = T

struct MyApplyStyle <: ApplyStyle end

@materialize function *(Ac::Adjoint{<:Any,<:AbstractMatrix},
                        O::MyOperator,
                        B::AbstractMatrix)
    MyApplyStyle
    T -> begin
        A = parent(Ac)

        if O.kind == :diagonal
            Diagonal(Vector{T}(undef, O.n))
        else
            Tridiagonal(Vector{T}(undef, O.n-1),
                        Vector{T}(undef, O.n),
                        Vector{T}(undef, O.n-1))
        end
    end
    dest::Diagonal -> begin
        dest.diag .= 1
    end
    dest::Tridiagonal{T} -> begin
        dest.dl .= -2
        dest.d .= 1
        dest.du .= 3
    end
end

@testset "Materialize DSL" begin
    o = ones(10)
    M = ones(10,10)
    D = MyOperator{Float64}(10, :diagonal)
    T = MyOperator{ComplexF64}(10, :tridiagonal)

    @test LazyArrays.ApplyStyle(*, typeof(M'), typeof(D), typeof(M)) == MyApplyStyle()
    @test LazyArrays.ApplyStyle(*, typeof(M'), typeof(T), typeof(M)) == MyApplyStyle()

    d = apply(*, M', D, M)
    @test d isa Diagonal{Float64}
    @test all(d.diag .== 1)

    t = apply(*, M', T, M)
    @test t isa Tridiagonal
    @test all(t.dl .== -2)
    @test all(t.d .== 1)
    @test all(t.du .== 3)

    M̃ = ones(11,11)
    @test_throws DimensionMismatch apply(*, M̃', D, M̃)
end


text:
# -------- Schut et al (2020):
struct GreedyGenerator <: AbstractGradientBasedGenerator
    loss::Symbol # loss function
    δ::AbstractFloat # perturbation size
    n::Int # maximum number of times any feature can be changed
end

"""
    GreedyGenerator(
        ;
        loss::Symbol=:logitbinarycrossentropy,
        δ::Union{AbstractFloat,Nothing}=nothing,
        n::Union{Int,Nothing}=nothing
    )

An outer constructor method that instantiates a greedy generator.

# Examples

```julia-repl
generator = GreedyGenerator()
```
"""
function GreedyGenerator(
    ;
    loss::Symbol=:logitbinarycrossentropy,
    δ::Union{AbstractFloat,Nothing}=nothing,
    n::Union{Int,Nothing}=nothing
)
    if all(isnothing.([δ, n]))
        δ = 0.1
        n = 10
    elseif isnothing(δ) && !isnothing(n)
        δ = 1/n
    elseif !isnothing(δ) && isnothing(n)
        n = 1/δ
    end

    generator = GreedyGenerator(loss,δ,n)

    return generator
end

"""
    ∇(generator::GreedyGenerator, counterfactual_state::CounterfactualState)

he default method to compute the gradient of the counterfactual search objective for a greedy generator. Since no complexity penalty is needed, this gradients just correponds to the partial derivative with respect to the loss function.

"""
∇(generator::GreedyGenerator, counterfactual_state::CounterfactualState) = ∂ℓ(generator, counterfactual_state)

"""
    generate_perturbations(generator::GreedyGenerator, counterfactual_state::CounterfactualState)

The default method to generate perturbations for a greedy generator. Only the most salient feature is perturbed.
"""
function generate_perturbations(generator::GreedyGenerator, counterfactual_state::CounterfactualState)
    𝐠ₜ = ∇(generator, counterfactual_state.M, counterfactual_state) # gradient
    𝐠ₜ[counterfactual_state.params[:mutability] .== :none] .= 0
    Δx′ = reshape(zeros(length(counterfactual_state.x′)), size(𝐠ₜ))
    iₜ = argmax(abs.(𝐠ₜ)) # choose most salient feature
    Δx′[iₜ] -= generator.δ * sign(𝐠ₜ[iₜ]) # counterfactual update
    return Δx′
end

"""
    mutability_constraints(generator::GreedyGenerator, counterfactual_state::CounterfactualState)

The default method to return search state dependent mutability constraints for a greedy generator. Features that have been perturbed `n` times already can no longer be perturbed.
"""
function mutability_constraints(generator::GreedyGenerator, counterfactual_state::CounterfactualState)
    mutability = counterfactual_state.params[:mutability]
    mutability[counterfactual_state.search[:times_changed_features] .>= generator.n] .= :none # constrains features that have already been exhausted
    return mutability
end

"""
    conditions_satisified(generator::GreedyGenerator, counterfactual_state::CounterfactualState)

If all features have been perturbed `n` times already, then the search terminates.
"""
function conditions_satisified(generator::GreedyGenerator, counterfactual_state::CounterfactualState)
    status = all(counterfactual_state.search[:times_changed_features].>=generator.n)
    return status
end


text:
"""
    FeMat{T,S}

Term with an explicit, constant matrix representation

# Fields
* `x`: matrix
* `wtx`: weighted matrix
* `piv`: pivot `Vector{Int}` for moving linearly dependent columns to the right
* `rank`: computational rank of `x`
* `cnames`: vector of column names

"""
mutable struct FeMat{T,S<:AbstractMatrix}
    x::S
    wtx::S
    piv::Vector{Int}
    rank::Int
    cnames::Vector{String}
end

"""
    FeMat(X::AbstractMatrix, cnms)

Convenience constructor for [`FeMat`](@ref) that computes the rank and pivot with unit weights.

See the vignette "[Rank deficiency in mixed-effects models](@ref)" for more information on the
computation of the rank and pivot.
"""
function FeMat(X::AbstractMatrix, cnms)
    T = eltype(X)
    if size(X,2) > 0
        st = statsqr(X)
        pivot = st.p
        rank = findfirst(<=(0), diff(st.p))
        rank = isnothing(rank) ? length(pivot) : rank
        Xp = pivot == collect(1:size(X, 2)) ? X : X[:, pivot]
        # single-column rank deficiency is the result of a constant column vector
        # this generally happens when constructing a dummy response, so we don't
        # warn.
        if rank < length(pivot) && size(X,2) > 1
            @warn "Fixed-effects matrix is rank deficient"
        end
    else
        # although it doesn't take long for an empty matrix,
        # we can still skip the factorization step, which gets the rank
        # wrong anyway
        pivot = Int[]
        Xp = X
        rank = 0
    end
    FeMat{T,typeof(X)}(Xp, Xp, pivot, rank, cnms[pivot])
end

"""
    FeMat(X::SparseMatrixCSC, cnms)

Convenience constructor for a sparse [`FeMat`](@ref) assuming full rank, identity pivot and unit weights.

Note: automatic rank deficiency handling may be added to this method in the future, as discused in
the vignette "[Rank deficiency in mixed-effects models](@ref)" for general `FeMat`.
"""
function FeMat(X::SparseMatrixCSC, cnms)
    @debug "Full rank is assumed for sparse fixed-effect matrices."
    rank = size(X,2)
    FeMat{eltype(X),typeof(X)}(X, X, 1:rank, rank, cnms)
end

function reweight!(A::FeMat{T}, sqrtwts::Vector{T}) where {T}
    if !isempty(sqrtwts)
        if A.x === A.wtx
            A.wtx = similar(A.x)
        end
        mul!(A.wtx, Diagonal(sqrtwts), A.x)
    end
    A
end

Base.adjoint(A::FeMat) = Adjoint(A)

Base.eltype(A::FeMat{T}) where {T} = T

fullrankwtx(A::FeMat) = rank(A) == size(A, 2) ? A.wtx : A.wtx[:, 1:rank(A)]

Base.length(A::FeMat) = length(A.wtx)

LinearAlgebra.rank(A::FeMat) = A.rank

Base.size(A::FeMat) = size(A.wtx)

Base.size(A::Adjoint{T,<:FeMat{T}}) where {T} = reverse(size(A.parent))

Base.size(A::FeMat, i) = size(A.wtx, i)

Base.copyto!(A::FeMat{T}, src::AbstractVecOrMat{T}) where {T} = copyto!(A.x, src)

*(adjA::Adjoint{T,<:FeMat{T}}, B::FeMat{T}) where {T} =
    fullrankwtx(adjA.parent)' * fullrankwtx(B)

LinearAlgebra.mul!(R::StridedVecOrMat{T}, A::FeMat{T}, B::StridedVecOrMat{T}) where {T} =
    mul!(R, A.x, B)

LinearAlgebra.mul!(C, adjA::Adjoint{T,<:FeMat{T}}, B::FeMat{T}) where {T} =
    mul!(C, fullrankwtx(adjA.parent)', fullrankwtx(B))

"""
    isfullrank(A::FeMat)

Does `A` have full column rank?
"""
isfullrank(A::FeMat) = A.rank == length(A.piv)


text:
#############################################
# hurdle
#############################################

abstract type TwoPartModel{Z<:RegressionModel,P<:RegressionModel} <: RegressionModel end

"Hurdle returned object"
mutable struct Hurdle{Z<:RegressionModel,P<:RegressionModel} <: TwoPartModel{Z,P}
  mzero::Z                # model for zeros
  mpos::P                 # model for positive counts
  fittedzero::Bool        # whether the model for zeros was fitted
  fittedpos::Bool         # whether the model for positives was fitted
end

# does nothing for hurdle
excessy!(ypos, ::Type{Hurdle}) = ypos

# minimum ypos element
minypos(::Type{Hurdle}) = 1.0

"Returns an (offsetzero, offsetpos) tuple of offset vector"
function setoffsets(y::AbstractVector, ixpos::BitVector, offset::AbstractVector, offsetzero::AbstractVector, offsetpos::AbstractVector)
    # set offsets
    if length(offset) != 0
      if length(offsetpos) == 0
        offsetpos = offset
      end
      if length(offsetzero) == 0
        offsetzero = offset
      end
    end

    if length(offsetpos) == length(y)
      offsetpos = offsetpos[ixpos]
    end

    offsetzero, offsetpos
end

"Returns positives indicators for y"
function getIy(y::AbstractVector{T}) where {T}
  # find positive y entries
  n = length(y)
  ixpos = BitVector(undef, n)
  Iy = Vector{T}(undef, n)
  @inbounds for i=1:n
    if y[i] == zero(T)
      Iy[i] = zero(T)
      ixpos[i] = false
    else
      Iy[i] = one(T)
      ixpos[i] = true
    end
  end

  ixpos, Iy
end

"Drops observations with infinite offset"
function finiteoffsetobs(part::Symbol, X, y, offset, wts, showwarnings)
  if any(isinf,offset)
      ixfinite = isfinite.(offset)
      showwarnings && @warn("omitting $(length(offset)-sum(ixfinite)) observations with infinite offset from $part model")
      offset = offset[ixfinite]
      X = X[ixfinite, :]
      y = y[ixfinite]
      wts = wts[ixfinite]
  end

  X, y, offset, wts
end

mildexception(e) = false
mildexception(e::PosDefException) = true
mildexception(e::DomainError) = true
mildexception(e::ConvergenceException) = true
mildexception(e::ErrorException) = (occursin("step-halving", e.msg) || occursin("failure to converge", e.msg) || occursin("failed to converge", e.msg))

function getlogger(showwarnings::Bool)
  if showwarnings
    current_logger()
  else
    MinLevelLogger(current_logger(), Logging.Error)
  end
end

# fitblank is wrapped in try catch only for M<:RegularizationPath and when the penalty_factor is specified
fitblank(::Type{M}, args...; penalty_factor::Union{Vector,Nothing}=nothing, kwargs...) where M<:RegularizationPath =
  fitblank(penalty_factor, M, args...; kwargs...)

function fitblank(penalty_factor::Vector, ::Type{M}, args...; kwargs...) where M<:RegularizationPath
  try
    fit(M, args...; λ=[0.0], dofit=false, penalty_factor=penalty_factor, kwargs...)
  catch e
    @warn "fitblank($M) caught a $e and will ignore penalty_factor for this part"
    fit(M, args...; λ=[0.0], dofit=false, penalty_factor=nothing, kwargs...)
  end
end

fitblank(penalty_factor::Nothing, ::Type{M}, args...; kwargs...) where M<:RegularizationPath =
  fit(M, args...; λ=[0.0], dofit=false, kwargs...)
fitblank(::Type{M}, args...; kwargs...) where M<:RegressionModel =
  fit(M, args...; dofit=false, kwargs...)

"Fits the model for zeros Iy ~ X"
function fitzero(::Type{M},
  X::AbstractMatrix{T}, Iy::V,
  dzero::UnivariateDistribution,
  lzero::Link,
  dofit::Bool,
  wts::V,
  offsetzero::AbstractVector,
  verbose::Bool,
  showwarnings::Bool,
  fitargs...) where {M<:RegressionModel,T<:FP,V<:FPVector}

  X, Iy, offsetzero, wts = finiteoffsetobs(:zeros, X, Iy, offsetzero, wts, showwarnings)

  local mzero
  fittedzero = false

  with_logger(getlogger(showwarnings)) do
    if var(Iy) > zero(T)
      try
        mzero = fit(M, X, Iy, dzero, lzero; dofit=dofit, wts=wts, offset=offsetzero, verbose=verbose, fitargs...)
        fittedzero = dofit
      catch e
        @warn("failed to fit zero counts model, possibly not enough variation in I(y). countmap(Iy)=$(countmap(Iy))")
        if mildexception(e)
          fittedzero = false
        else
          @warn("X'=$(X')")
          @warn("Iy=$Iy)")
          rethrow(e)
        end
      end
    else
      if verbose
        if all(iszero,Iy)
          @warn("I(y) is all zeros. There is nothing to explain.")
        else
          @warn("I(y) is all ones. Data may be fully described by a poisson model.")
        end
      end
    end

    if !fittedzero
      # create blank zeros model without fitting
      mzero = fitblank(M, X, Iy, dzero, lzero; wts=wts, offset=offsetzero, verbose=verbose, fitargs...)
    end
  end

  mzero, fittedzero
end

"Fits the model for positives ypos ~ Xpos"
function fitpos(::Type{TPM},::Type{M},
  Xpos::AbstractMatrix{T}, ypos::V,
  dpos::UnivariateDistribution,
  lpos::Link,
  dofit::Bool,
  wtspos::V,
  offsetpos::AbstractVector,
  verbose::Bool,
  showwarnings::Bool,
  fitargs...) where {TPM<:TwoPartModel, M<:RegressionModel,T<:FP,V<:FPVector}

  excessy!(ypos, TPM)
  Xpos, ypos, offsetpos, wtspos = finiteoffsetobs(:positives, Xpos, ypos, offsetpos, wtspos, showwarnings)

  local mpos
  fittedpos = false

  with_logger(getlogger(showwarnings)) do
    if any(x->x>minypos(TPM), ypos)
      try
        mpos = fit(M, Xpos, ypos, dpos, lpos; dofit=dofit, wts=wtspos, offset=offsetpos, verbose=verbose, fitargs...)
        fittedpos = dofit
      catch e
        @warn("failed to fit truncated counts model to positive subsample, possibly not enough variation in ypos. countmap(ypos)=$(sort(countmap(ypos)))")
        if mildexception(e)
          fittedpos = false
        else
          @warn("Xpos'=$(Xpos')")
          @warn("ypos=$ypos)")
          rethrow(e)
        end
      end
    else
      if length(ypos) == 0
        error("y is all zeros! There is nothing to explain.")
      else
        @warn("ypos has no elements larger than $(minypos(TPM))! Data may be fully described by a probability model.")
      end
    end

    if !fittedpos
      # create blank positives model without fitting
      mpos = fitblank(M, Xpos, ypos, dpos, lpos; wts=wtspos, offset=offsetpos, verbose=verbose, fitargs...)
    end
  end

  mpos, fittedpos
end

"""
    fit(Hurdle,M,X,y; Xpos=Xpos, <keyword arguments>)

Fit a Hurdle (Mullahy, 1986) of count vector y on X with potentially another
covariates matrix Xpos used to model positive counts.

# Example with GLM:
```julia
  m = fit(Hurdle,GeneralizedLinearModel,X,y; Xpos=Xpos)
  yhat = predict(m, X; Xpos=Xpos)
```

# Example with Lasso regularization:
```julia
  m = fit(Hurdle,GammaLassoPath,X,y; Xpos=Xpos)
  yhat = predict(m, X; Xpos=Xpos, select=MinAICc())
```

# Arguments
- `M::RegressionModel`
- `counts` n-by-d matrix of counts (usually sparse)
- `dzero::UnivariateDistribution = Binomial()` distribution for zeros model
- `dpos::UnivariateDistribution = PositivePoisson()` distribution for positives model
- `lzero::Link=canonicallink(dzero)` link function for zeros model
- `lpos::Link=canonicallink(dpos)` link function for positives model

# Keywords
- `Xpos::Union{AbstractMatrix{T},Nothing} = nothing` covariates matrix for positives
  model or nothing to use X for both parts
- `dofit::Bool = true` fit the model or just construct its shell
- `wts::V = ones(y)` observation weights
- `offsetzero::AbstractVector = similar(y, 0)` offsets for zeros model
- `offsetpos::AbstractVector = similar(y, 0)` offsets for positives model
- `offset::AbstractVector = similar(y, 0)` offsets for both model parts
- `verbose::Bool=true`
- `showwarnings::Bool=false`
- `fitargs...` additional keyword arguments passed along to fit(M,...)
"""
function StatsBase.fit(::Type{Hurdle},::Type{M},
  X::AbstractMatrix{T}, y::V,
  dzero::UnivariateDistribution = Binomial(),
  dpos::UnivariateDistribution = PositivePoisson(),
  lzero::Link = canonicallink(dzero),
  lpos::Link = canonicallink(dpos);
  Xpos::Union{AbstractMatrix{T},Nothing} = nothing,
  dofit::Bool = true,
  wts::V = fill(one(eltype(y)),size(y)),
  offsetzero::AbstractVector = similar(y, 0),
  offsetpos::AbstractVector = similar(y, 0),
  offset::AbstractVector = similar(y, 0),
  verbose::Bool = false,
  showwarnings::Bool = false,
  fitargs...) where {M<:RegressionModel,T<:FP,V<:FPVector}

  ixpos, Iy = getIy(y)

  offsetzero, offsetpos = setoffsets(y, ixpos, offset, offsetzero, offsetpos)

  mzero, fittedzero = fitzero(M, X, Iy, dzero, lzero, dofit, wts, offsetzero, verbose, showwarnings, fitargs...)

  # Xpos optional argument allows to specify a data matrix only for positive counts
  if Xpos == nothing
    # use X for Xpos too
    Xpos = X[ixpos,:]
  elseif size(Xpos,1) == length(y)
    # Xpos has same dimensions as X, take only positive y ones
    Xpos = Xpos[ixpos,:]
  end

  mpos, fittedpos = fitpos(Hurdle, M, Xpos, y[ixpos], dpos, lpos, dofit, wts[ixpos], offsetpos, verbose, showwarnings, fitargs...)

  Hurdle(mzero,mpos,fittedzero,fittedpos)
end

"""
    fit(Hurdle,M,f,df; fpos=Xpos, <keyword arguments>)

Takes dataframe and two formulas, one for each model part. Otherwise same arguments
as [`fit(::Hurdle)`](@ref)

# Example
```julia
  fit(Hurdle,GeneralizedLinearModel,@formula(y ~ x1*x2), df; fpos=@formula(y ~ x1*x2+x3))
```
"""
function StatsBase.fit(::Type{Hurdle},::Type{M},
                      f::FormulaTerm,
                      data,
                      dzero::UnivariateDistribution = Binomial(),
                      dpos::UnivariateDistribution = PositivePoisson(),
                      lzero::Link = canonicallink(dzero),
                      lpos::Link = canonicallink(dpos);
                      fpos::FormulaTerm = f,
                      dofit::Bool = true,
                      wts = fill(1.0,size(data,1)),
                      offsetzero = Float64[],
                      offsetpos = Float64[],
                      offset = Float64[],
                      verbose::Bool = false,
                      showwarnings::Bool = false,
                      contrasts::Dict{Symbol,<:Any} = Dict{Symbol,Any}(),
                      fitargs...) where {M<:RegressionModel}

  Tables.istable(data) || throw(ArgumentError("expected data in a Table, got $(typeof(data))"))

  cols = columntable(data)

  mfzero = ModelFrame(f, cols, model=M, contrasts=contrasts)
  mmzero = ModelMatrix(mfzero)
  y = response(mfzero)

  ixpos, Iy = getIy(y)

  offsetzero, offsetpos = setoffsets(y, ixpos, offset, offsetzero, offsetpos)

  # fit zero model to entire sample
  # TODO: should be wrapped in TableRegressionModel but can't figure out right now why it complains it is not defined
  # mzero = TableRegressionModel(fit(M, mmzero.m, Iy, dzero, lzero; wts=wts, offset=offsetzero, fitargs...), mfzero, mmzero)
  mzero, fittedzero = fitzero(M, mmzero.m, Iy, dzero, lzero, dofit, wts, offsetzero, verbose, showwarnings, fitargs...)

  mfpos = (f===fpos) ? mfzero : ModelFrame(fpos, cols, model=M, contrasts=contrasts)
  mmpos = (f===fpos) ? mmzero : ModelMatrix(mfpos)
  mmpos.m = mmpos.m[ixpos,:]

  mpos, fittedpos = fitpos(Hurdle, M, mmpos.m, y[ixpos], dpos, lpos, dofit, wts[ixpos], offsetpos, verbose, showwarnings, fitargs...)

  Hurdle(mzero,mpos,fittedzero,fittedpos)
end

function Base.show(io::IO, tpm::TwoPartModel)
  name = typeof(tpm).name
  println(io, "$name regression\n")

  if typeof(tpm.mpos) <: RegularizationPath
    println(io, "Positive part regularization path ($(distfun(tpm.mpos)) with $(linkfun(tpm.mpos)) link):")
    println(io, tpm.fittedpos ? tpm.mpos : "Not Fitted")

    println(io, "Zero part regularization path ($(distfun(tpm.mzero)) with $(linkfun(tpm.mzero)) link):")
    println(io, tpm.fittedzero ? tpm.mzero : "Not Fitted")
  else
    println(io, "Positive part coefficients ($(Distribution(tpm.mpos)) with $(Link(tpm.mpos)) link):")
    println(io, tpm.fittedpos ? tpm.mpos : "Not Fitted")

    println(io, "Zero part coefficients ($(Distribution(tpm.mzero)) with $(Link(tpm.mzero)) link):")
    println(io, tpm.fittedzero ? tpm.mzero : "Not Fitted")
  end
end

"""
    coef(m::Hurdle; select=MinAICc())

Returns a selected segment of the coefficient matrices of the fitted the TwoPartModel.

# Example:
```julia
  m = fit(Hurdle,GeneralizedLinearModel,X,y; Xpos=Xpos)
  coefspos, coefszero = coef(m)
```

# Keywords
- `kwargs...` are passed along to two coef() calls on the two model parts.
"""
function StatsBase.coef(tpm::TwoPartModel; kwargs...)
  czero = coef(tpm.mzero; kwargs...)
  cpos = coef(tpm.mpos; kwargs...)
  cpos,czero
end

# the following coef definitions override those in Lasso for Hurdles
function StatsBase.coef(tpm::TwoPartModel, select::SegSelect; kwargs...)
  czero = coef(tpm.mzero, select; kwargs...)
  cpos = coef(tpm.mpos, select; kwargs...)
  cpos,czero
end

function StatsBase.coef(tpm::TwoPartModel, select::S; kwargs...) where {S<:CVSegSelect}
  error("""
    Specifying an instance of a `CVSegSelect` is not supported because there is
    more than one path and its generator has a fixed number of observation indices.
    Instead, consider passing a `MinCVKfold{$S}(k)`.
    """)
end

"Selects the RegularizationPath segment according to `CVSegSelect` with `k`-fold cross-validation"
struct MinCVKfold{S<:CVSegSelect} <: CVSegSelect
  k::Int  # number of CV folds
end

"Selects the RegularizationPath segment coefficients according to `S` with `k`-fold cross-validation"
function StatsBase.coef(path::RegularizationPath, select::MinCVKfold{S};
  kwargs...) where {S<:CVSegSelect}

  selector = S(path, select.k)
  coef(path, selector; kwargs...)
end

"Selects the RegularizationPath segment coefficients according to `S` with `k`-fold cross-validation"
function StatsBase.coef(tpm::TwoPartModel, select::MinCVKfold{S};
  kwargs...) where {S<:CVSegSelect}

  selectzero = S(tpm.mzero, select.k)
  selectpos = S(tpm.mpos, select.k)
  czero = coef(tpm.mzero, selectzero; kwargs...)
  cpos = coef(tpm.mpos, selectpos; kwargs...)
  cpos,czero
end

# predicted counts vector given expected inclusion (μzero) and repetition (μpos)
μtpm(m::Hurdle, μzero, μpos) = μzero .* μpos

"""
    predict(m,X; Xpos=Xpos, <keyword arguments>)

Predict using a fitted TwoPartModel given new X (and potentially Xpos).

# Example with GLM:
```julia
  m = fit(Hurdle,GeneralizedLinearModel,X,y; Xpos=Xpos)
  yhat = predict(m, X; Xpos=Xpos)
```

# Example with Lasso regularization:
```julia
  m = fit(Hurdle,GammaLassoPath,X,y; Xpos=Xpos)
  yhat = predict(m, X; Xpos=Xpos, select=MinAICc())
```

# Arguments
- `m::Hurdle` fitted Hurdle model
- `X` n-by-p matrix of covariates of same dimensions used to fit m.

# Keywords
- `Xpos::Union{AbstractMatrix{T},Nothing} = nothing` covariates matrix for positives
  model or nothing to use X for both parts
- `kwargs...` additional keyword arguments passed along to predict() for each
  of the two model parts.
"""
function StatsBase.predict(tpm::TwoPartModel, X::AbstractMatrix{T};
  Xpos::AbstractMatrix{T} = X,
  offsetzero::AbstractVector = Array{T}(undef, 0),
  offsetpos::AbstractVector = Array{T}(undef, 0),
  offset::AbstractVector=Array{T}(undef, 0),
  kwargs...) where {T<:AbstractFloat}

  # set offsets
  if length(offset) != 0
    if length(offsetpos) == 0
      offsetpos = offset
    end
    if length(offsetzero) == 0
      offsetzero = offset
    end
  end

  μzero = predict(tpm.mzero, X; offset=offsetzero, kwargs...)
  μpos = predict(tpm.mpos, Xpos; offset=offsetpos, kwargs...)

  @assert size(μzero) == size(μpos) "Predicted values from zero and positives models have different dimensions, $(size(μzero)) != $(size(μpos))\nCan result from fitting a RegularizationPath with autoλ and select=AllSeg() is specified."

  μtpm(tpm, μzero, μpos)
end


text:
using DataFrames, Unitful, UnitfulAstro, Statistics

export compare_df, strip_blanks!, gpc

gaussian(x, mean=0., sigma=1.) =
    (1 / sqrt(2pi) / sigma) * exp(-((x - mean) / sigma)^2 / 2)

showv(df::DataFrameRow) =
    show(DataFrame(field=names(df), value=[values(df)...]), allrows=true, allcols=true)

function askpass(msg="")
    @assert Base.Sys.isunix()
    io = Base.getpass(msg)
    out = String(read(io))
    Base.shred!(io)
    return out
end

function compare_df(a::DataFrame, b::DataFrame)
    if ncol(a) != ncol(b)
        @warn "Number of cols is different: " ncol(a) ncol(b)
    end
    if nrow(a) != nrow(b)
        @warn "Number of rows is different: " nrow(a) nrow(b)
    end
    na = names(a)
    nb = names(b)

    function sd(v1, v2, i1, i2)
        if v1[i1] < v2[i2]
            return -1
        elseif v1[i1] == v2[i2]
            return 0
        end
        return 1
    end
    jj = sortmerge(na, nb, sd=sd)
    if  (length(jj[1]) < length(na))  ||
        (length(jj[2]) < length(nb))
        @warn "Only " * string(length(jj[1])) * " columns have the same name"
    end

    if nrow(a) != nrow(b)
        @warn "Number of rows is different"
    else
        for i in 1:length(jj[1])
            colname = string(na[jj[1][i]])
            @info colname
            da = a[:, jj[1][i]]
            db = b[:, jj[2][i]]

            k = findall(ismissing.(da) .!= ismissing.(db))
            if length(k) != 0
                @warn string(length(k)) * " rows have different missing condition on column $colname"
            else
                if count(.!ismissing.(da)) > 0
                    da = disallowmissing(collect(skipmissing(da)))
                    db = disallowmissing(collect(skipmissing(db)))

                    (length(da) == 0)  &&  continue
                    if eltype(da) <: AbstractFloat
                        k = findall(da .!== db)  # handle NaNs
                    else
                        k = findall(da .!= db)  # handle NaNs
                    end
                    if length(k) != 0
                        @warn string(length(k)) * " rows have different values on column $colname"
                    end
                end
            end
        end
    end
end

function strip_blanks!(df::DataFrame)
    for i in 1:ncol(df)
        if nonmissingtype(eltype(df[:, i])) == String
            df[:, i] .= string.(strip.(df[:, i]))
        end
    end
end

# Reference: http://www.astro.wisc.edu/~dolan/constants.html
Base.@kwdef struct gpc
    # Dimensionless or common to both CGS and MKS
    jd_mjd    = 2400000.5              # Offset to change from Julian date to MJD
    NA        = 6.0221367    * 1e23    # Avagadro's number
    α         = 7.29735308   * 1e-3    # Fine structure constant
    day       = 86400.                 # days to sec fator
    month     = 2592000.               # months to sec fator
    year      = 31536000.              # year to sec fator

    # CGS
    c         = 2.99792458  * 1e10     * u"cm" / u"s"                          # Vacuum speed of light
    h         = 6.6260755   * 1e-27    * u"erg" * u"s"                         # Planck's constant
    ħ         = 1.05457266  * 1e-27    * u"erg" * u"s"                         # Reduced Planck's constant
    G         = 6.67259     * 1e-8     * u"cm"^3 * u"g"^-1 * u"s"^-2           # Gravitational constant
    e         = 4.8032068   * 1e-10                                            # Electron charge (esu)
    me        = 9.1093897   * 1e-28    * u"g"                                  # Mass of electron
    mp        = 1.6726231   * 1e-24    * u"g"                                  # Mass of proton
    mn        = 1.6749286   * 1e-24    * u"g"                                  # Mass of neutron
    mH        = 1.6733      * 1e-24    * u"g"                                  # Mass of hydrogen
    amu       = 1.6605402   * 1e-24    * u"g"                                  # Atomic mass unit
    k         = 1.380658    * 1e-16    * u"erg" / u"K"                         # Boltzmann constant
    eV        = 1.6021772   * 1e-12    * u"erg" / u"eV"                        # Electron volt to erg factor
    a         = 7.5646      * 1e-15    * u"erg" * u"cm"^-3 * u"K"^-4           # Radiation density constant
    σ         = 5.67051     * 1e-5     * u"erg" * u"cm"^-2 * u"K"^-4 * u"s"^-1 # Stefan-Boltzmann constant
    R_inf     = 1.097373    * 1e5      * u"cm"^-1                              # R_infinity
    au        = 1.496       * 1e13     * u"cm" / u"AU"                         # Astronomical unit to cm
    pc        = 3.086       * 1e18     * u"cm" / u"pc"                         # Parsec to cm factor
    ly        = 9.463       * 1e17     * u"cm" / u"ly"                         # Light year to cm factor
    ld        = 2.5902      * 1e15                                             # Light day to cm factor (cm ly^-1)
    Msun      = 1.99        * 1e33     * u"g"                                  # Solar mass
    Rsun      = 6.96        * 1e10     * u"cm"                                 # Solar radius
    Lsun      = 3.9         * 1e33     * u"erg" / u"s"                         # Solar luminosity
    Tsun      = 5.780       * 1e3      * u"K"                                  # Solar temperature
    Mearth    = 5.974       * 1e27     * u"g"                                  # Earth mass
    Rearth    = 6372.8      * 1e5      * u"cm"                                 # Earth mean radius
    g         = 981.52                 * u"cm" * u"s"^-2                       # Acceleration at Earth surface
    thom      = 0.66524616  * 1e-24    * u"cm"^2                               # Thomson cross section
    jansky    = 1e-23                  * u"erg" * u"cm"^-2 * u"s"^-1 * u"Hz"^-1# Flux density (keV cm^-2 s^-1 keV^/1)
    wien      = 2.82 * k / h
    deg       = pi/180                 * u"rad" * u"°"^-1
    arcsec    = pi/180/3600                                                    # rad arcsec^-1
    mas       = pi/180/3600/1000                                               # rad milliarcsec^-1
    edd       = 1.26e38                * u"erg" * u"s"^-1 * u"Msun"^-1         # Eddington luminosity (erg s^-1 M_sun^-1)
    r_cm      = 2 * pi^2 * me * e^4 / (h^3) / c                                # Rydberg constant (cm^-1)
    r2_cm     = r_cm / (1 + me/mp)                                             # Rydberg constant (cm^-1, reduced mass)
    A         = 1.e-8                  * u"cm" / Unitful.angstrom
end

function smooth(y, n)
    out = y .* 1.
    @assert mod(n, 2) == 1
    @assert n >= 3
    h = div(n-1, 2)
    for i in 1+h:length(y)-h
        out[i] = mean(y[i-h:i+h])
    end
    return out
end

function ppvalunc(v, e)
    f = 10^ceil(-log10(e))
    pe = round(e * f) / f
    if isinteger(log10(pe))  &&  pe <= e
        # if the only sig. digit is 1
        f *= 10
        pe = round(e * f) / f
    end
    pv = round(v * f) / f
    # println(v, " ", e, " ", f, " ", pv, " ", pe)
    return pv, pe
end


text:
module ScanImageIO
using Images,FileIO, ImageMetadata, ImageAxes, ImageMagick, SharedArrays, Distributed

function scanImage2016Reader(files::Array{String,1};binFile=nothing)

    @info "Reading metadata"
    (pixelsPerLine,linesPerFrame,nSlices,nFrames,resolutionXY,resolutionZ,samplingTime,realSlices) = read_metadata(files[1],read_pos=false)

    nFrames = [nFrames]

    for f in files[2:end]
        tempMetadata = read_metadata(f,read_pos=false)
        @assert tempMetadata[[1:3;5:7]] == (pixelsPerLine,linesPerFrame,nSlices,resolutionXY,resolutionZ,samplingTime) "Images do not have the same size"
        nFrames = vcat(nFrames,tempMetadata[4])
    end

    ## Finding the longest recording and get the frame positions from this one
    long_run = findmax(nFrames)[2]
    im_pos = read_metadata(files[long_run],read_pos=true)[9]

    nTotalFrames = sum(nFrames)
    if binFile == nothing
        @info "Creating shared Array"
        out = SharedArray{Int16}((linesPerFrame,pixelsPerLine,realSlices,nTotalFrames))
    else
        binFile = abspath(binFile)
        @info "Creating bin file $(binFile)"
        out = SharedArray{Int16}(binFile,(linesPerFrame,pixelsPerLine,realSlices,nTotalFrames))
    end

    im_length = pixelsPerLine * linesPerFrame * realSlices
    startPoints = vcat(1,cumsum(nFrames*im_length)[1:end-1].-1)

    frame_length = pixelsPerLine * linesPerFrame

    pos_idx = [vcat([i:(i+realSlices-1) for i in 1:nSlices:nF]...) for nF in nFrames*nSlices]
    im_pos = [im_pos[pI] for pI in pos_idx]

    @sync begin
        for p in procs(out)
            @async remotecall_wait(writeim_toshared_native_chunk, p, out,files,startPoints,im_pos,frame_length)
        end
    end

    ## Returning the SharedArray and a dict of metadata
    (out,Dict("framesPerTrial" => nFrames,"resolutionXY"=>resolutionXY,"resolutionZ"=>resolutionZ,"samplingTime"=>samplingTime))

end

function read_metadata(f::String;read_pos=true)
    ## WIP. Just reads the main metadata for now. TOD : MROI situations etc...
    fi = open(f)

    (ifd_header_pos,nv_frame_data_length,roi_group_data_length) = read_SI_header(fi)

    extracomment = read_SI_meta(fi,nv_frame_data_length)

    if read_pos
        ## Collecting where the frames are and their length
        img_pos = []
        while ifd_header_pos != 0
            img_pos_temp,ifd_header_pos = read_ifd_header(fi,ifd_header_pos)
            img_pos = vcat(img_pos,img_pos_temp)
        end
    end

    close(fi)

    improps = extract_image_properties(extracomment)

    if read_pos
        return (improps...,img_pos)
    else
        return improps
    end
end

function read_ifd_header(fi,ifd_header_pos)
    seek(fi,ifd_header_pos)
    n_ifd = read(fi,Int)

    ## Skip to the image location field
    skip(fi,20*6)

    @assert read(fi,Int16) == 273 "Failed to find strip offset field"
    skip(fi,10)
    img_pos = read(fi,Int)

    ## Finding the next frame
    skip(fi,20*11)
    ifd_header_pos = read(fi,Int)
    return (img_pos,ifd_header_pos)
end

function read_SI_header(fi)
    seek(fi,8)
    ifd_header_pos = read(fi,Int)

    @assert read(fi,Int32) == 117637889 "Not a scanimage file"
    @info "ScanImage TIFF version $(read(fi,Int32))"

    nv_frame_data_length = read(fi,Int32)

    roi_group_data_length = read(fi,Int32)
    return (ifd_header_pos,nv_frame_data_length,roi_group_data_length)
end

function read_SI_meta(fi,header_length)
    parse_SI_meta(String(read(fi,header_length)))
end

function parse_SI_meta(extracomment)
    extracomment = split(extracomment,"\n")[1:end-1]
    extracomment = [split(replace(replace(str,"'"=>"\""),"SI."=>""),"= ") for str in extracomment]
    extracomment = Dict(strip(des[1]) => des[2] for des in extracomment)
end

function extract_image_properties(extracomment)
    framesPerVolume = Meta.parse(extracomment["hFastZ.numFramesPerVolume"])
    if typeof(framesPerVolume) !== Expr
        nSlices = Meta.parse(extracomment["hFastZ.numFramesPerVolume"])
        nFrames = Meta.parse(extracomment["hFastZ.numVolumes"])
        samplingTime = (1/Meta.parse(extracomment["hRoiManager.scanVolumeRate"]))
    else
        nFrames = Meta.parse(extracomment["hStackManager.framesPerSlice"])
        nSlices = 1
        samplingTime = Meta.parse(extracomment["hRoiManager.scanFramePeriod"])
    end
    pixelsPerLine = Meta.parse(extracomment["hRoiManager.pixelsPerLine"])
    linesPerFrame = Meta.parse(extracomment["hRoiManager.linesPerFrame"])
    realNSlices = nSlices - Meta.parse(extracomment["hFastZ.numDiscardFlybackFrames"])

    resolutionXY = 2*eval(Meta.parse(extracomment["hRoiManager.imagingFovUm"]))[2,1]/Meta.parse(extracomment["hRoiManager.pixelsPerLine"])
    resolutionZ = Meta.parse(extracomment["hStackManager.stackZStepSize"])

    return (pixelsPerLine,linesPerFrame,nSlices,nFrames,resolutionXY,resolutionZ,samplingTime,realNSlices)
end

function myrange(q::SharedArray,bigrange)
    idx = indexpids(q)
    if idx == 0 # This worker is not assigned a piece
        return 1:0, 1:0
    end
    nchunks = length(procs(q))
    splits = [round(Int, s) for s in range(bigrange[1],stop=bigrange[2],length=nchunks+1)]
    splits[idx]+1:splits[idx+1]
end

function writeim_toshared_native_chunk(out,files,startPoints,im_pos,imlength)
   writeim_toshared_native_chunk(out,files,startPoints,im_pos,imlength,myrange(out,[0,length(files)]))

end

function writeim_toshared_native_chunk(out,files,startPoints,im_pos,imlength,filerange)
    for i in filerange
        writeim_toshared_native(out,files[i],startPoints[i],im_pos[i],imlength)
    end
end

function writeim_toshared_native(out,file,startPoint,im_pos,imlength)

    imsz = (size(out,2),size(out,1))

    try
        fi = open(file)
        sP = startPoint
        for fr in 1:length(im_pos)
            seek(fi,im_pos[fr])
            out[sP:(sP+imlength-1)] = permutedims(reshape([read(fi,Int16) for j in 1:imlength],imsz),(2,1))
            sP+=imlength
        end
        @info "Loaded file $(file)"
        close(fi)
        return out
    catch
        @info "Need to read the image positions for file $(file)"
        metad = read_metadata(file,read_pos=true)
        realSlices = metad[8]
        im_pos = metad[9]
        nSlices = metad[3]
        nFrames =metad[4]
        pos_idx = vcat([i:(i+realSlices-1) for i in 1:nSlices:nFrames]...)
        im_pos = im_pos[pos_idx]
        fi = open(file)
        for fr in 1:length(im_pos)
            seek(fi,im_pos[fr])
            out[startPoint:(startPoint+imlength-1)] = permutedims(reshape([read(fi,Int16) for j in 1:imlength],imsz),(2,1))
            startPoint+=imlength
        end
        @info "Loaded file $(file)"
        close(fi)
        return out
    end

end

## Reads a single file.
function scanImage2016Reader(file::String;binFile=nothing)
    scanImage2016Reader([file];binFile=binFile)
end

function scanImage5Reader(f;objMag=40,shared=false,mapped=false,binFile="imgFile")

    img = load(File(format"TIFF",f))
    extraprops = magickinfo(f,["comment"])

    ## Correcting the fact that ImageMagick reads the image as unsigned integers (and return the integer value)
    img = reinterpret(Int16,img)

    comment = extraprops["comment"]
    comment = split(comment,"\n")[1:end-1]
    comment = [split(replace(replace(str,"'","\""),"scanimage.SI5.","")," = ") for str in comment]
    comment = Dict(strip(des[1]) => des[2] for des in comment)

    ## Extract the dimensions/number of frames from the SI comment
    framesAcq = size(img,3)#parse(comment["frameNumberAcquisition"])
    nSlices = nSlices =  parse(comment["stackNumSlices"])
    nFrames = div(framesAcq,nSlices)
    if nSlices>1
        samplingTime = samplingTime = parse(comment["fastZPeriod"])
    else
        samplingTime = samplingTime = parse(comment["scanFramePeriod"])
    end

    img = img[:,:,1:(nFrames*nSlices)]

    img = reshape(img,(size(img)[1],size(img)[2],nSlices,nFrames))
    img = img[:,:,1:(nSlices-parse(comment["fastZNumDiscardFrames"])),:]
    #img = Gray.(img)

    resolutionXY = ((tan(10*2*pi/360)*45000/objMag)/(parse(comment["pixelsPerLine"])*parse(comment["zoomFactor"])))
    resolutionZ = (parse(comment["stackZStepSize"])/10)

    if shared
        img = SharedArray(img)
    elseif mapped
        img = SharedArray(binFile,img)
    end

    img = ImageMeta(AxisArray(img,Axis{:x}(range(0,resolutionXY,size(img)[1])),
                              Axis{:y}(range(0,resolutionXY,size(img)[2])),
                              Axis{:z}(range(0,resolutionZ,size(img)[3])),
                              Axis{:time}(range(0,samplingTime,size(img)[4]))),comment=comment)

    img

end

export scanImage2016Reader, scanImage5Reader

end # module


text:
export
    run_callback

# run callback and return whether simlation should terminate
run_callback(callback::Any, rec::EntityQueueRecord{S,D,I}, roadway::R, models::Dict{I,M}, tick::Int) where {S,D,I,R,M<:DriverModel} = error("run_callback not implemented for callback $(typeof(callback))")

function _run_callbacks(callbacks::C, rec::EntityQueueRecord{S,D,I}, roadway::R, models::Dict{I,M}, tick::Int) where {S,D,I,R,M<:DriverModel,C<:Tuple{Vararg{Any}}}
    isdone = false
    for callback in callbacks
        isdone |= run_callback(callback, rec, roadway, models, tick)
    end
    return isdone
end
function simulate!(
    ::Type{A},
    rec::EntityQueueRecord{S,D,I},
    scene::EntityFrame{S,D,I},
    roadway::R,
    models::Dict{I,M},
    nticks::Int,
    callbacks::C,
    ) where {S,D,I,A,R,M<:DriverModel,C<:Tuple{Vararg{Any}}}

    empty!(rec)
    update!(rec, scene)

    # potential early out right off the bat
    if _run_callbacks(callbacks, rec, roadway, models, 0)
        return rec
    end

    actions = Array{A}(undef, length(scene))
    for tick in 1 : nticks
        get_actions!(actions, scene, roadway, models)
        tick!(scene, roadway, actions, get_timestep(rec))
        update!(rec, scene)
        if _run_callbacks(callbacks, rec, roadway, models, tick)
            break
        end
    end

    return rec
end
function simulate!(
    rec::EntityQueueRecord{S,D,I},
    scene::EntityFrame{S,D,I},
    roadway::R,
    models::Dict{I,M},
    nticks::Int,
    callbacks::C,
    ) where {S,D,I,R,M<:DriverModel,C<:Tuple{Vararg{Any}}}

    return simulate!(Any, rec, scene, roadway, models, nticks, callbacks)
end

max_stars_repo_name: intelligent-control-lab/Autoenv

max_issues_repo_name: intelligent-control-lab/Autoenv

max_forks_repo_name: intelligent-control-lab/Autoenv