text:
include("/Users/dbassen/work/KwateeServer-v1.0/src/DataFile.jl")
include("/Users/dbassen/work/KwateeServer-v1.0/src/SolveBalances.jl")
using Debug

# =============================================================================== #
# Solve.jl
# Function library to solve the EMT model for different experimental/test cases
# Author: J.Varner
#
# Version: 1.0
# Date: 12/17/15
# =============================================================================== #
# Run model to steady state and then give pertubation

# Simulate the addition of TGFB -
function SolveModelAddTGFB(TSTART,TSTOP,Ts,data_dictionary)

  # Do we have a data_dictionary (if not load a new one)
  if (isempty(data_dictionary))
    data_dictionary = DataFile(TSTART,TSTOP,Ts)
  end

  # Setup simulation time -
  TSIM = collect(TSTART:Ts:TSTOP);

  # Run the model to steady-state *before* we add the TGFB -
  XSS = EstimateSteadyState(data_dictionary);

  # ok, we established a steady state -
  # reset the initial conditions -
  initial_condition_array = vec(XSS);

  # Add TGFB1/2 -
  initial_condition_array[56] = 10.0;
  data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;

  # Run the simulation -
  (TSIM,X) = SolveBalances(TSTART,TSTOP,Ts,data_dictionary);

  # return tuple -
  return (TSIM,X)

end

# Run the model to steady state -
function EstimateSteadyState(data_dictionary)

  # Do we have a data_dictionary (if not load a new one)
  if (isempty(data_dictionary))
    data_dictionary = DataFile(TSTART,TSTOP,Ts)
  end

  # Setup loop -
  EPSILON = 1e-3;
  TSTART = 0.0;
  Ts = 1.0;
  TSTOP = 100;
  did_reach_steady_state = false
  while (!did_reach_steady_state)

    # solve the balances -
    (TSIM,X1) = SolveBalances(TSTART,TSTOP,Ts,data_dictionary);

    # Take a few additional steps -
    TNEXT_START = TSTOP+Ts;
    TNEXT_STOP = TNEXT_START+1.0;
    Ts = 0.1;

    # solve the balances again 0
    initial_condition_array = vec(X1[end,:])
    data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
    (TSIM,X2) = SolveBalances(TNEXT_START,TNEXT_STOP,Ts,data_dictionary);

    # Find the difference -
    DIFF = norm((X2[end,:] - X1[end,:]));

    # Should we stop -or- go around again?
    if (DIFF<EPSILON)
      did_reach_steady_state = true;
      return (X2[end,:]);
    else

      # No, we did *not* reach steady state ....
      TSTART = TSTOP+Ts
      TSTOP = 1000 + TSTART;
      Ts = 100.0;

      initial_condition_array = vec(X2[end,:])
      data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
    end
  end

  # return
  return XSS;

end


text:
#=Copyright (c) 2017 Gabriel Goh

Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
the Software without restriction, including without limitation the rights to
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
the Software, and to permit persons to whom the Software is furnished to do so,
subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.=#

using Plots

function getiterk(Q,b,α)
  Σ,U = eig(Q)
  c = U'b
  lambdak(k,λ) = (1 - (1 - α*λ).^k)./λ
  return k -> U*(lambdak(k,Σ).*c)
end

Box = zeros(16,16)
Cavity = ones(16,16)
Cavity[8:9,:] = 0
R = [Box Cavity Box Cavity Box Cavity Box]
A = ones(size(R,1)+2, size(R,2) + 2)
A[2:end-1,2:end-1] = R

n = size(A,1) - 2
m = size(A,2) - 2

Ind = zeros(Integer, size(A))
IndInv = Tuple{Int64,Int64}[]
count = 1
for i = 1:size(A,1), j = 1:size(A,2)
  if A[i,j] == 0
    Ind[i,j] = count
    push!(IndInv, (i,j))
    count = count + 1
  end
end

fp = open("Uval.json","w")
JSON.print(fp, convert(Array{Float16,2},U))
close(fp)

fp = open("Sigma.json","w")
JSON.print(fp, Σ)
close(fp)
# ==================================================
# Generate Laplacian
# ==================================================

Q = zeros(length(IndInv), length(IndInv))
for i = 2:(size(A,1)-1), j = 2:(size(A,2)-1)
  k = Ind[i,j]
  if A[i,j] == 0
    A[i,j+1] == 0 ? Q[Ind[i,j+1],k] = 1 : nothing
    A[i+1,j] == 0 ? Q[Ind[i+1,j],k] = 1 : nothing
    A[i,j-1] == 0 ? Q[Ind[i,j-1],k] = 1 : nothing
    A[i-1,j] == 0 ? Q[Ind[i-1,j],k] = 1 : nothing
  end
end

Q = (Q + Q')/2
L = diagm(sum(Q,1)[:]) - Q
L[449,449] = L[449,449] + 0.1
L[561,561] = L[561,561] + 0.1

Σ,U = eig(L)

function arrange(x)
  A = zeros(n,m)
  for ind = 1:length(x)
    i,j = IndInv[ind]
    A[i-1,j-1] = x[ind]
  end
  return A
end

plot(size = (880,140)); heatmap!(arrange(U[:,10]), map = "seismic")

imshow(arrange(U[:,100]), interpolation = "None", cmap = "spectral"); show()

pyplot(leg=false, ticks=nothing)

@gif for i in 2:500
  plot(size = (880,140)); heatmap!(reshape(U[:,224*3+i],n,m))
end

e1 = zeros(length(IndInv))
e1[449] = 1
e1[561] = 1

α = 1/max(Σ...)

iterk = getiterk(L,e1,α)

imshow(arrange(iterk(10000)))
show()

# α = 1/max(eig(L)[1]...)
# iters = 5000
# X = zeros(n*m,iters)
# x = 10*zeros(n*m)
# for i = 1:iters
#   X[:,i] = x
#   for i = 1:100
#     x = x - (1*α)*(L*x - e1)
#   end
#   println(norm(L*x - e1))
# end

# plot(size = (880,140)); heatmap!(reshape(X[:,100],n,m))
# pyplot(leg=false, ticks=nothing)

# @gif for i in 10:10:2000
#   println(i)
#   xi = reshape(X[:,i],n,m)
#   plot(size = (880,140)); heatmap!(xi)
# end

# lambdak(k,λ) = α*(1 - (1 - α*λ).^k)./(α*λ)
# c = U'e1

# @gif for i i guein 10:500:300000
#   println(i)
#   z = U*(lambdak(i,Σ).*(c)
#   plot(size = (880,140)); heatmap!(reshape(z, n, m))
# end

# A = randn(10,10); A = A*A'

# b = rand(10)
# x = zeros(10)
# for i = 1:1000
#   x = x - 0.01*(A*x - b)
# end
# x - getiterk(A,b,0.01)(1000)

ext: jl

max_stars_repo_name: sbeleidy/post--momentum

max_issues_repo_name: sbeleidy/post--momentum

max_forks_repo_name: sbeleidy/post--momentum


text:
#===========================================
Tests on the RTN generation routine function
===========================================#

#=
Here we calculate the autocorrelation functions of the two methods of generating RTN histories. The functions need a huge time in input (40000. is ok) in order to have ergodicity, and they return the autocorrelation array [<RTN(0)RTN(0)>-<RTN(0)><RTN(0)>,<RTN(1)RTN(0)>-<RTN(1)><RTN(0)>,...] up to the lag corresponding to an evolution time of tc=4./mu.
=#
using Test
using Statistics
using StatsBase

function autocorrelation(noisefunction, time, mu, dt)
    t = range(0., stop=time, step=dt)

    signal = noisefunction(t, mu)[:]
    tc = 4. / mu
    dt = Float64(t.step)

    lags = [x for x = 0:(Int(fld(tc, dt))-1)]
    autoC = autocov(signal, lags)

    tvec = range(0., stop=tc, length=Int(fld(tc,dt)))
    return (tvec, autoC)
end

function spectrum(t, signal)
    n = length(signal)
    p = fft(signal)
    sampFreq = 1/(t[2] - t[1])
    nUniquePts = ceil(Int, (n+1)/2)
    p = p[1:nUniquePts]
    p = abs2.(p) / n^2
    # the fourier transform of the tone returned by the fft function contains both magnitude and phase information and is given in a complex representation (i.e. returns complex numbers). By taking the absolute value of the fourier transform we get the information about the magnitude of the frequency components.

    # odd nfft excludes Nyquist point
    if n % 2 > 0
        p[2 : length(p)] = p[2:length(p)]*2 # we've got odd number of   points fft
    else
        p[2: (length(p) -1)] = p[2: (length(p) -1)]*2 # we've got even number of points fft
    end

    freqArray = (0:(nUniquePts-1)) * (sampFreq / n)
    return (freqArray, p)
end

function test_autocorrelation(noisefunction, mu; atol=1e-1)
    (t, ac) = autocorrelation(noisefunction, 20000. /mu, mu, .05 / mu)
    @test ac ≈ exp.(-2 * mu * t) atol = atol
end

# Test that the mean is zero
@testset "RTN mean" begin
    @testset "γ = $γ" for γ in [0.01, 0.1 , 1., 10.]
        @test mean(generateRTN(0. : .05 / γ : 1000 / γ, γ, 100)) ≈ 0 atol=5e-3
    end
end

# Test that the autocorrelation function is the decaying exponential
@testset "RTN autocorrelation" begin
    @testset "γ = $γ" for γ in [0.01, 0.1 , 1., 10.]
        @testset "RTNgenTime" for i in 1:100
            test_autocorrelation(generateRTN, γ, atol=1e-1)
        end
    end
end

max_stars_repo_name: matteoacrossi/QuantumSpatialSearch

max_issues_repo_name: matteoacrossi/QuantumSpatialSearch

max_forks_repo_name: matteoacrossi/QuantumSpatialSearch


text:
@generated function apodization!(p::NFFTPlan{D,0,T}, f::AbstractArray{U,D}, g::StridedArray{Complex{T},D}) where {D,T,U}
    quote
        @nexprs $D d -> offset_d = round(Int, p.n[d] - p.N[d]/2) - 1

        @nloops $D l f d->(gidx_d = rem(l_d+offset_d, p.n[d]) + 1) begin
            v = @nref $D f l
            @nexprs $D d -> v *= p.windowHatInvLUT[d][l_d]
            (@nref $D g gidx) = v
        end
    end
end

@generated function apodization_adjoint!(p::NFFTPlan{D,0,T}, g::AbstractArray{Complex{T},D}, f::StridedArray{U,D}) where {D,T,U}
    quote
        @nexprs $D d -> offset_d = round(Int, p.n[d] - p.N[d]/2) - 1

        @nloops $D l f begin
            v = @nref $D g d -> rem(l_d+offset_d, p.n[d]) + 1
            @nexprs $D d -> v *= p.windowHatInvLUT[d][l_d]
            (@nref $D f l) = v
        end
    end
end

function convolve!(p::NFFTPlan{D,0,T}, g::AbstractArray{Complex{T},D}, fHat::StridedVector{U}) where {D,T,U}
  if isempty(p.B)
    convolve_LUT!(p, g, fHat)
  else
    convolve_sparse_matrix!(p, g, fHat)
  end
end

function convolve_LUT!(p::NFFTPlan{D,0,T}, g::AbstractArray{Complex{T},D}, fHat::StridedVector{U}) where {D,T,U}
    scale = 1.0 / p.m * (p.K-1)

    #Threads.@threads
    for k in 1:p.M
        fHat[k] = _convolve_LUT(p, g, scale, k)
    end
end

@generated function _convolve_LUT(p::NFFTPlan{D,0,T}, g::AbstractArray{Complex{T},D}, scale, k) where {D,T}
    quote
        @nexprs $D d -> xscale_d = p.x[d,k] * p.n[d]
        @nexprs $D d -> c_d = floor(Int, xscale_d)

        fHat = zero(T)

        @inbounds @nloops $D l d -> (c_d-p.m):(c_d+p.m) d->begin
            # preexpr
            gidx_d = rem(l_d+p.n[d], p.n[d]) + 1
            idx = abs( (xscale_d - l_d)*scale ) + 1
            idxL = floor(idx)
            idxInt = Int(idxL)
            tmpWin_d = p.windowLUT[d][idxInt] + ( idx-idxL ) * (p.windowLUT[d][idxInt+1] - p.windowLUT[d][idxInt])
        end begin
            # bodyexpr
            v = @nref $D g gidx
            @nexprs $D d -> v *= tmpWin_d
            fHat += v
        end

        return fHat
    end
end

function convolve_sparse_matrix!(p::NFFTPlan{D,0,T}, g::AbstractArray{Complex{T},D}, fHat::StridedVector{U}) where {D,T,U}
  fill!(fHat, zero(T))

  mul!(fHat, transpose(p.B), vec(g))
end

function convolve_adjoint!(p::NFFTPlan{D,0,T}, fHat::AbstractVector{U}, g::StridedArray{Complex{T},D}) where {D,T,U}
  if isempty(p.B)
    convolve_adjoint_LUT!(p, fHat, g)
  else
    convolve_adjoint_sparse_matrix!(p, fHat, g)
  end
end

@generated function convolve_adjoint_LUT!(p::NFFTPlan{D,0,T}, fHat::AbstractVector{U}, g::StridedArray{Complex{T},D}) where {D,T,U}
    quote
        fill!(g, zero(T))
        scale = 1.0 / p.m * (p.K-1)

        @inbounds @simd for k in 1:p.M
            @nexprs $D d -> xscale_d = p.x[d,k] * p.n[d]
            @nexprs $D d -> c_d = floor(Int, xscale_d)

            @nloops $D l d -> (c_d-p.m):(c_d+p.m) d->begin
                # preexpr
                gidx_d = rem(l_d+p.n[d], p.n[d]) + 1
                idx = abs( (xscale_d - l_d)*scale ) + 1
                idxL = floor(idx)
                idxInt = Int(idxL)
                tmpWin_d = p.windowLUT[d][idxInt] + ( idx-idxL ) * (p.windowLUT[d][idxInt+1] - p.windowLUT[d][idxInt])
            end begin
                # bodyexpr
                v = fHat[k]
                @nexprs $D d -> v *= tmpWin_d
                (@nref $D g gidx) += v
            end
        end
    end
end

function convolve_adjoint_sparse_matrix!(p::NFFTPlan{D,0,T},
                        fHat::AbstractVector{U}, g::StridedArray{Complex{T},D}) where {D,T,U}
  fill!(g, zero(T))

  mul!(vec(g), p.B, fHat)
end


text:
using SeuratRDS
using Documenter

makedocs(;
    modules=[SeuratRDS],
    authors="Matt Karikomi <mattkarikomi@gmail.com> and contributors",
    repo="https://github.com/mkarikom/SeuratRDS.jl/blob/{commit}{path}#L{line}",
    sitename="SeuratRDS.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://mkarikom.github.io/SeuratRDS.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/mkarikom/SeuratRDS.jl",
)


text:
# function to generate a Kitaev colorcode
function generateBondColorcodeKitaev(
        lattice :: L
    ) :: Dict{LB,Vector{Int64}} where {LS,D,LB,N,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N},U,L<:AbstractLattice{S,B,U}}

    # put the default kitaev labels with default colors into a dictonary
    return Dict{LB,Vector{Int64}}(
        getDefaultLabelX(LB) => [255,0,0],
        getDefaultLabelY(LB) => [0,255,0],
        getDefaultLabelZ(LB) => [0,0,255]
    )
end

# function to generate a Kitaev colorcode
function generateBondColorcodeRandom(
        lattice :: L
    ) :: Dict{LB,Vector{Int64}} where {LS,D,LB,N,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N},U,L<:AbstractLattice{S,B,U}}

    # put the default kitaev labels with default colors into a dictonary
    d = Dict{LB,Vector{Int64}}()
    # insert all labels
    for lb in label.(bonds(lattice))
        d[lb] = round.(Int64,[rand()*255, rand()*255, rand()*255])
    end
    # return the colorcode dict
    return d
end


text: using HTTP2.Session


text:
import SatSolver: AstNode, VariableNode, NotNode, AndNode, OrNode, create_ast,
                  ast2cnf

@testset "NNF formula parser" begin
  function test_case(input, reference)
    buffer = IOBuffer(input)
    ast = create_ast(buffer)
    @test ast == reference
  end

  test_case("a", VariableNode("a"))
  test_case("(not a)", NotNode("a"))
  test_case("(and a b)", AndNode(VariableNode("a"), VariableNode("b")))
  test_case("(and (not a) b)", AndNode(NotNode("a"), VariableNode("b")))
  test_case("(and (or (not a) b) (and b c))",
            AndNode(OrNode(NotNode("a"), VariableNode("b")),
                    AndNode(VariableNode("b"), VariableNode("c"))))
end

@testset "Ast2Cnf" begin
  function test_case(ast::T,
                     clauses::Vector{Clause},
                     olr::Bool=false) where {T <: AstNode}
    res = ast2cnf(ast, only_left_to_right=olr)
    fla = Formula(clauses)
    @test Formula(res.clauses) == fla
  end

  test_case(VariableNode("a"), [Clause(1)])
  test_case(NotNode("a"), [Clause(-1)])
  test_case(AndNode(VariableNode("a"), NotNode("b")),
                   [Clause(-3, 1), Clause(-3, -2), Clause(-1, 2, 3)])
  test_case(AndNode(VariableNode("a"), NotNode("b")),
                   [Clause(-3, 1), Clause(-3, -2)],
                   true)
  test_case(OrNode(VariableNode("c"), AndNode(VariableNode("a"), VariableNode("b"))),
            [Clause(-4, 2), Clause(-4, 3), Clause(-2, -3, 4), Clause(-5, 1, 4),
             Clause(-1, 5), Clause(-4, 5)])
end


text:
using Iterators
using RLESUtils, RunUtils, IFTTTUtils

const COMP = gethostname()
const DATA = "nvn_dascfilt"
const CONFIG = "normal"

V_seeds = collect(1:20)
V_efs = [0.1, 0.25, 0.5, 0.6]
V_wnews = [0.1, 0.2, 0.4, 0.6, 0.8, 1.0]
V_wprior = [0.01, 0.05, 0.1, 0.2, 0.4]
V_maxsteps = [20, 40, 60]

template_ce(seed, elite_frac, w_new, w_prior, maxsteps) =
"""
using GrammarExpts,ACASX_CE,RLESUtils,IFTTTUtils
config=configure(ACASX_CE,"$CONFIG","$DATA")
config[:seed]=$seed
config[:elite_frac]=$(elite_frac)
config[:w_new]=$(w_new)
config[:w_prior]=$(w_prior)
config[:maxsteps]=$maxsteps
acasx_ce(;outdir=Pkg.dir("GrammarExpts/results/ACASX_CE_ParamStudy/seed$(seed)ef$(elite_frac)wnew$(w_new)wprior$(w_prior)maxsteps$maxsteps"), config...)
"""

notifydone() = sendifttt(;value1="$COMP done")

A = JuliaSource[]
for args in product(V_seeds, V_efs, V_wnews, V_wprior, V_maxsteps)
    push!(A, JuliaSource(template_ce(args...)))
end

#include("script.jl")
#fs=julia_process(A, np); notifydone()


text:
function aggregate(obj::AggregatableStatistic, arg0::Iterable)
    return jcall(obj, "aggregate", void, (Iterable,), arg0)
end

function aggregate(obj::AggregatableStatistic, arg0::Object)
    return jcall(obj, "aggregate", void, (Object,), arg0)
end

function aggregate(obj::AggregatableStatistic, arg0::Vector{Object})
    return jcall(obj, "aggregate", void, (Vector{Object},), arg0)
end


text:
"""
    autocorrelation(x)

Calculate autocorrelation of a vector for all positive lags.
"""
autocorrelation(x) = xcorr(x,x)[length(x):end]

# Numbers of bin
# Freeman-Diaconis rule
function n_fd(x::AbstractVector)
    first, last = extrema(x)
    iqr = quantile(x,.75) - quantile(x,.25)
    ceil(Integer, (last-first) / (2iqr * length(x)^(-1/3)) )
end

# Sturges
n_sturges(x::AbstractVector) = ceil(Integer, 1+log2(length(x)))

makebins(x, n) = linspace(extrema(x)..., n+1)

# Own simplified implementation of hist2d for square edges
function hist2!(h, x, edges)
    n = length(edges)-1
    for r=1:n, c=1:n
        h[r,c] = 0;
    end
    for i = 1:size(x,1)
        r = searchsortedfirst(edges, x[i,1]) - 1
        c = searchsortedfirst(edges, x[i,2]) - 1
        if 1 <= r <= n && 1 <= c <= n
            h[r,c] += 1
        end
    end
    h
end

"""
    ami(x, delay[, nbins])

Calculate the average mutual information of a signal for given delays.

The delay can either be an integer (in which case the function returns an integer),
or a series of numbers (returning a vector) specified as an array of ascending integers,
a range, or a tuple with the limits of the range (minimum, maximum).

The number of bins used to estimate the marginal and joint entropies of the
original and delayed signal can be given as a fixed integer, or as a
string that specifies a binning criterion. Currently available criteria are
Sturges (`"Sturges"`, the default value) or Freeman-Diaconis (`"FD"`).
"""
function ami(x, delay::Integer, nbins::Integer)
    n = length(x)-delay
    log2n = log2(n)
    edges = makebins(x, nbins)
    pxy = zeros(nbins,nbins)
    hist2!(pxy, x[(1:n).+[0 delay]], edges)
    px = sum(pxy, 2)
    py = sum(pxy, 1)
    ret = 0.
    for ix=1:nbins, iy=1:nbins
        if pxy[ix,iy] != 0
            pxpy = px[ix]*py[iy]
            ret += pxy[ix,iy]/n*(log2(pxy[ix,iy]/pxpy)+log2n)
        end
    end
    # Maximum entropy for a given nbins is log2(nbins)
    ret/log2(nbins)
end

function ami(x, delay::Union{Array, Range}, nbins::Integer)
    d2 = maximum(delay)
    miv = zeros(d2+1)
    n = length(x)-d2
    log2n = log2(n)
    # Pre-allocated arrays
    pxy = zeros(nbins,nbins)
    px = zeros(nbins,1)
    py = zeros(1,nbins)
    for d in delay
        edges = makebins(x[1:n+d], nbins)
        hist2!(pxy, x[(1:n).+[0 d]], edges)
        sum!(px, pxy)
        sum!(py, pxy)
        for ix=1:nbins, iy=1:nbins
            if pxy[ix,iy] != 0
                pxpy = px[ix]*py[iy]
                miv[d+1] += pxy[ix,iy]/n*(log2(pxy[ix,iy]/pxpy)+log2n)
            end
        end
    end
    miv[1+delay]/log2(nbins)
end

ami(x, delay::Tuple{Integer,Integer}, nbins) = ami(x, colon(delay...), nbins)

function ami(x, delay, nbins="Sturges")
    # Define number of bins
    dict_binfuns = Dict("Sturges" => n_sturges,
        "FD" => n_fd)
    !haskey(dict_binfuns,nbins) && error("incorrect definition of bins.")
    binfun = dict_binfuns[nbins]
    nbins = binfun(x)
    ami(x, delay, nbins)
end

"""
    gmi(x, delay, radius)

Calculate the generalized mutual information of a signal for given delays, based
on Rényi entropies.

The delay can either be an integer (in which case the function returns an integer),
or a series of numbers (returning a vector) specified as an array of ascending integers,
a range, or a tuple with the limits of the range (minimum, maximum).

The radius is a fixed value in the absolute scale of the signal,
used for the calculation of recurrence matrices to estimate
Rényi entropies (see `?recurrencematrix` for details).
Maximum norm and no scaling of distances are assumed.
"""
function gmi(x, delay::Integer, radius::Real)
    rmfull = recurrencematrix(x,radius,scale=1)
    n = length(x) - delay
    # Entropy of the non-delayed series
    rm1 = rmfull[1:n, 1:n]
    h2 = -log2(nnz(rm1)/n^2)
    # Joint entropy
    rmj = rm1 .* rmfull[delay+(1:n), delay+(1:n)]
    h2tau = -log2(nnz(rmj)/n^2)
    # Maximum entropy for a given radius corresponds to uniform distribution
    maxh2 = log2((maximum(x)-minimum(x))/radius)
    (2h2 - h2tau)/maxh2
end

function gmi(x, delay::Union{Array, Range}, radius::Real)
    d2 = maximum(delay)
    rmfull = recurrencematrix(x,radius,scale=1)
    h2tau = zeros(d2+1)
    n = length(x) - d2
    # Entropy of the non-delayed series
    rm1 = rmfull[1:n, 1:n]
    h2 = -log2(nnz(rm1)/n^2)
    # Joint entropy
    rmj = spzeros(n, n)
    for d in delay
        rmj = rm1 .* rmfull[d+(1:n), d+(1:n)]
        h2tau[d+1] = -log2(nnz(rmj)/n^2)
    end
    # Maximum entropy for a given radius corresponds to uniform distribution
    maxh2 = log2((maximum(x)-minimum(x))/radius)
    (2h2 - h2tau[1+delay])/maxh2
end

gmi(x, delay::Tuple{Integer,Integer}, radius) = gmi(x, colon(delay...), radius)

# gmi(x::AbstractVector, delay) = gmi(x, delay, radius_mrr(x, .01))


text:
if !isdefined(Main, :CustomStencil)
    include("../module/CustomStencil.jl")
    using Main.CustomStencil
end
include("misc.jl")
include("cpu_stencils.jl")

## Stencil Definition
coefs = Float64.([1; 0.5; 0.25; 0.125; 0.06])
star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1, 4, c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
            D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))
st_def = CreateStencilDefinition(star_stencil, coefs)
st_inst = NewStencilInstance(st_def)

## Input Data size Definition
radius = 4
nx = 8
ny = 8
nz = 8
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

t_steps = 4
## GPU Implementation
original_data = CreateData(dx,dy,dz)
data = PadData(radius, original_data)
dout = ApplyStencil(st_inst, data, t_steps)
gpu_out = Array(dout)


text:
mutable struct Boss
    name::String
    hp::Int
end

Boss_dict = Dict(
            "normal" => Dict(
                1 => Boss("ゴブリングレート",600),
                2 => Boss("ライライ",800),
                3 => Boss("ニードルクリーパー",1000),
                4 => Boss("サイクロプス",1200),
                5 => Boss("レサトパルト",1500)
            ),
            "very_hard" => Dict(
                1 => Boss("ゴブリングレート",700),
                2 => Boss("ライライ",900),
                3 => Boss("ニードルクリーパー",1200),
                4 => Boss("サイクロプス",1400),
                5 => Boss("レサトパルト",1700)
            )
)

max_stars_repo_name: hibibol/clan-battle

max_issues_repo_name: hibibol/clan-battle

max_forks_repo_name: hibibol/clan-battle


text:
"""
Package: Forecast

Store results from the function `forecast`

# Arguments
   model:  Model object containing information about the fitted model.
   x:      Original time series.
   alpha:  The confidence levels associated with the prediction intervals.
   mean:   Point forecasts.
   lower:  Lower limits for prediction intervals.
   upper:  Upper limits for prediction intervals.
   se:     Standard Error.
"""
mutable struct FORECAST
    model
    alpha::Tuple
    mean::DataFrame
    upper::DataFrame
    lower::DataFrame
    se::DataFrame
    call::String
end

function Base.show(io::IO, fc::FORECAST)
    ts = eltype(fc.mean[:,1]) in [Date, DateTime, Time,
                                  Year, Month, Quarter, Week, Day,
                                  Hour, Second, Millisecond,
                                  Microsecond, Nanosecond]
    printstyled("Forecast Information\n",bold=true,color=:underline)
    print("\n",fc.call, "\n")
    printstyled("\nMean Forecasting\n",bold=true,color=:underline)
    pretty_table(fc.mean, noheader = false, nosubheader = true, show_row_number=false)
    printstyled("\nPrediction Intervals alpha at: "*string(fc.alpha)*"\n",bold=true,color=:underline)
    printstyled("\nUpper:\n",color=:underline)
    pretty_table(fc.upper, noheader = false, nosubheader = true, show_row_number=false,
                 vlines = [0,1,div(size(fc.upper,2),2)+1,size(fc.upper,2)])
    printstyled("\nLower:\n",color=:underline)
    pretty_table(fc.lower, noheader = false, nosubheader = true, show_row_number=false,
                 vlines = [0,1,div(size(fc.lower,2),2)+1,size(fc.lower,2)])
end

"""
Rename data in FORECAST object
"""
function setnames!(fc::FORECAST,new_names::AbstractVector{String})
    n_mean = names(fc.mean)
    n_mean[2:end] = new_names
    rename!(fc.mean,n_mean)

    nn_upper = []
    for n in new_names
        push!(nn_upper,"upper1_"*n)
        push!(nn_upper,"upper2_"*n)
    end
    n_upper = names(fc.upper)
    n_upper[2:end] = nn_upper
    rename!(fc.upper,n_upper)

    nn_lower = []
    for n in new_names
        push!(nn_lower,"lower1_"*n)
        push!(nn_lower,"lower2_"*n)
    end
    n_lower = names(fc.lower)
    n_lower[2:end] = nn_lower
    rename!(fc.lower,n_lower)

    if fc.model.x isa DataFrame
        n_model = names(fc.model.x)
        n_model[2:end] = new_names
        rename!(fc.model.x,n_model)
        fc.model.varnames=new_names
    end
end

"""
Transform a FORECAST object value with given function
"""
function transform(fc::FORECAST, f::Function, vari::Union{Integer,UnitRange} = 1:size(fc.mean,2)-1)

    vari = collect(vari)

    # Renaming
    tfc = deepcopy(fc)
    fnames = tfc.model.varnames
    fnames[vari] .= "$(string(f))_" .* fnames[vari]
    setnames!(tfc, fnames)

    # Transform
    vari2 =  vcat(vari, vari .+ size(fc.mean,2) .- 1)
    tfc.model.x[:,vari .+ 1] .= f.(Array(tfc.model.x[:,vari .+ 1]))
    tfc.mean[:,vari .+ 1] .= f.(Array(tfc.mean[:,vari .+ 1]))
    tfc.lower[:,vari2 .+ 1] .= f.(Array(tfc.lower[:,vari2 .+ 1]))
    tfc.upper[:,vari2 .+ 1] .= f.(Array(tfc.upper[:,vari2 .+ 1]))

    tfc.call = fc.call * "\nData transformed with function: $(f)"

    return(tfc)
end


text:
######################################################################################################
# Copy from https://github.com/xuedong/LinBAI.jl  by Xuedong Shang                                   #
# Please note that XY-Adaptive can only run for one fixed δ, not multiple δs!                        #
######################################################################################################
using JLD2;
using Distributed;
using Printf;
using IterTools;
using Distributions
@everywhere include("runit_bai.jl");
@everywhere include("../utilities/thresholds.jl");
include("../utilities/experiment_helpers.jl");

# Setup: replicate the benchmark setting (ii) in (Jedra and Proutiere, 2020)
#        which uses the example considered by (Soare et al., 2014) but with ω = 0.1
dist = Gaussian();
dim = 6;
# best arm
μ = zeros(dim); µ[1] = 2.0;
arms = Vector{Float64}[];
for k = 1:dim
    v = zeros(dim); v[k] = 1.0; push!(arms, v);
end
# the last arm which is hard to distinguish with the 1st arm
ω = 0.1; v = zeros(dim); v[1] = cos(ω); v[2] = sin(ω); push!(arms, v)

pep = LinearBestArm(dist, arms);
# methods to be compared
srs = [FWSampling(), LazyTrackAndStop(), ConvexGame(CTracking), LearnerK(CTracking), XYAdaptive(), RoundRobin()]; # (1) run only for δs = (0.01,) for plot figure
#srs = [FWSampling(), LazyTrackAndStop(), ConvexGame(CTracking), LearnerK(CTracking), RoundRobin()]; # (2) run all except XYAdaptive for δs = (0.1,0.01,0.001,0.0001)
#srs = [XYAdaptive()]; # (3) separately run for different δs

"""
Note that XY-Adaptive can only run for one δ.
If you'd like to run multiple δs, e.g., δs=(0.1, 0.01, 0.0001, 0.00001),
then please remove XYAdaptive() from the srs
"""
δs = (0.01,); # (1)
#δs = (0.1,0.01,0.001,0.0001); # (2)
#δs = (0.1,); # (3) separately run for each δ
βs = GK16.(δs);
repeats = 1000;
seed = 1234;

println("ω=$ω, dim=$dim, repeat $repeats times");
# compute
@time data = pmap(
    ((sr, i),) -> runit(seed + i, sr, μ, pep, βs, δs),
    Iterators.product(srs, 1:repeats),
);
dump_stats(pep, μ, δs, βs, srs, data, repeats);
# save
@save isempty(ARGS) ? "BAI_0.01.dat" : ARGS[1] dist μ pep srs data δs βs repeats seed
# visualise by loading viz.jl


text:
## GCMGridSpec function with default gridName argument:

GCMGridSpec() = GCMGridSpec("LLC90")

## GCMGridSpec function with gridName argument:

"""
    GCMGridSpec(gridName)

Set global variables in the module scope for grDir, nFaces, grTopo, ioSize,
facesSize, ioPrec using hard-coded values for LLC90, CS32, LL360 (for now).
"""
function GCMGridSpec(gridName)

global grDir, nFaces, grTopo, ioSize, facesSize, ioPrec;

if gridName=="LLC90";
    grDir="GRID_LLC90/";
    nFaces=5;
    grTopo="llc";
    ioSize=[90 1170];
    facesSize=[(90, 270), (90, 270), (90, 90), (270, 90), (270, 90)]
    ioPrec=Float64;
elseif gridName=="CS32";
    grDir="GRID_CS32/";
    nFaces=6;
    grTopo="cs";
    ioSize=[32 192];
    facesSize=[(32, 32), (32, 32), (32, 32), (32, 32), (32, 32), (32, 32)]
    ioPrec=Float32;
elseif gridName=="LL360";
    grDir="GRID_LL360/";
    nFaces=1;
    grTopo="ll";
    ioSize=[360 160];
    facesSize=[(360, 160)]
    ioPrec=Float32;
else;
    error("unknown gridName case");
end;

return "GCMGridSpec: passed"

end

## GCMGridLoad function

"""
    GCMGridLoad()

Loads grid variables from files located in grDir set by GCMGridSpec.

Grid variables are XC, XG, YC, YG, RAC, RAZ, DXC, DXG, DYC, DYG, hFacC,
hFacS, hFacW, Depth based on the MITgcm naming convention.
"""
function GCMGridLoad()

    global XC, XG, YC, YG, RAC, RAZ, DXC, DXG, DYC, DYG, hFacC, hFacS, hFacW, Depth;

    list0=("XC","XG","YC","YG","RAC","RAZ","DXC","DXG","DYC","DYG","hFacC","hFacS","hFacW","Depth");
    for ii=1:length(list0);
        tmp1=read_bin(grDir*list0[ii]*".data",ioPrec);
        tmp2=Symbol(list0[ii]);
        @eval (($tmp2) = ($tmp1))
    end

    return "GCMGridLoad: passed"

end

"""
    GCMGridOnes(grTp,nF,nP)

Define all-1 grid variables instead of using GCMGridSpec & GCMGridLoad.
"""
function GCMGridOnes(grTp,nF,nP)

    global grDir, nFaces, grTopo, ioSize, facesSize, ioPrec;
    global XC, XG, YC, YG, RAC, RAZ, DXC, DXG, DYC, DYG, hFacC, hFacS, hFacW, Depth;

    grDir=""
    grTopo=grTp
    nFaces=nF
    ioSize=[nP nP*nF]
    facesSize=Array{NTuple{2, Int},1}(undef,nFaces)
    facesSize[:].=[(nP,nP)]
    ioPrec=Float32

    list0=("XC","XG","YC","YG","RAC","RAZ","DXC","DXG","DYC","DYG","hFacC","hFacS","hFacW","Depth");
    for ii=1:length(list0);
        tmp1=fill(1.,nP,nP*nF);
        tmp1=convert2gcmfaces(tmp1);
        tmp2=Symbol(list0[ii]);
        @eval (($tmp2) = ($tmp1))
    end

    return "GCMGridOnes: passed"

end


text:
"""
    $(FUNCTIONNAME)(FT=Float64)

Creates a `LeafDistribution` following a planophile (mostly horizontal) distribution using a Beta distribution

Standard values for θₗ (26.76) and s (18.51) from Bonan https://doi.org/10.1017/9781107339217.003, Table 2.1
"""
function planophile_leaves(FT=Float64)
    θ = FT(26.76)
    s = FT(18.51)
    α, β =  βparameters(deg2rad(θ), deg2rad(s));
    return LeafDistribution(Beta(FT(α), FT(β)), FT(2/π));
end

function planophile_leaves2(FT=Float64)
    return LeafDistribution(Beta(FT(1), FT(2.768)), FT(2/π));
end

"""
    $(FUNCTIONNAME)(FT=Float64)

Creates a `LeafDistribution` following a uniform distribution (all angles equally likely)
"""
function uniform_leaves(FT=Float64)
    return LeafDistribution(Uniform(FT(0), FT(1)),FT(2/π));
end

"""
    $(FUNCTIONNAME)(FT=Float64)

Creates a `LeafDistribution` following a plagiophile (mostly 45degrees) distribution using a Beta distribution

Standard values for θₗ (45.0) and s (16.27) from Bonan https://doi.org/10.1017/9781107339217.003, Table 2.1
"""
function plagiophile_leaves(FT=Float64)
    θ = FT(45)
    s = FT(16.27)
    α, β =  βparameters(deg2rad(θ), deg2rad(s));
    return LeafDistribution(Beta(FT(α), FT(β)), FT(2/π));
end

"""
    $(FUNCTIONNAME)(FT=Float64)

Creates a `LeafDistribution` following a erectophile (mostly vertical) distribution using a Beta distribution

Standard values for θₗ (63.24) and s (18.5) from Bonan https://doi.org/10.1017/9781107339217.003, Table 2.1
"""
function erectophile_leaves(FT=Float64)
    θ = FT(63.24)
    s = FT(18.5)
    α, β =  βparameters(deg2rad(θ), deg2rad(s));
    return LeafDistribution(Beta(FT(α), FT(β)), FT(2/π));
end

"""
    $(FUNCTIONNAME)(FT=Float64)

Creates a `LeafDistribution` following a spherical (distributed as facets on a sphere) distribution using a Beta distribution

Standard values for θₗ (57.3) and s (21.55) from Bonan https://doi.org/10.1017/9781107339217.003, Table 2.1
"""
function spherical_leaves(FT=Float64)
    θ = FT(57.3)
    s = FT(21.55)
    α, β =  βparameters(deg2rad(θ), deg2rad(s));
    return LeafDistribution(Beta(FT(α), FT(β)), FT(2/π));
end

function flat_leaves(FT=Float64)
    θ = FT(57.3)
    s = FT(21.55)
    α, β =  βparameters(deg2rad(θ), deg2rad(s));
    return LeafDistribution(Beta(FT(α), FT(β)), FT(2/π));
end


text:
#    Copyright 2015 Raytheon BBN Technologies
#
#     Licensed under the Apache License, Version 2.0 (the "License");
#     you may not use this file except in compliance with the License.
#     You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#     Unless required by applicable law or agreed to in writing, software
#     distributed under the License is distributed on an "AS IS" BASIS,
#     WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#     See the License for the specific language governing permissions and
#     limitations under the License.

using Convex

"""
 ϕ represents the isomorphism between complex and real matrices.
 e.g., see [Invariant semidefinite programs](http://arxiv.org/abs/1007.2905)
 by Bachoc et al.  If two arguments are given, they are taken to be the real and
 imaginary parts of a complex matrix.
"""
function ϕ(r,i)
    if size(r) != size(i)
        error("ϕ requires both arguments to be of the same size.")
    end
    [r i; -i r]
end

function ϕ(c)
    r = real(c)
    i = imag(c)
    ϕ(r,i)
end

"""
Extract the real part of a complex matrix represented as a real matrix
"""
function ϕr(m)
    m[1:end/2,1:end/2]
end

"""
Extract the real part of a complex matrix represented as a real matrix.
"""
function ϕi(m)
    m[1:end/2,end/2+1:end]
end

"""
Extract the real part of a complex matrix represented as a real matrix.
"""
function ϕinv(m)
   ϕr(m), ϕi(m)
end

"""
Compute the trace of the real representation of a complex matrix.
"""
function retrϕ(m)
    trace(ϕr(m))
end

function ket(i,d)
    v = spzeros(Float64,d,1)
    v[i+1] = 1.0
    return v
end

function bra(i,d)
    return ket(i,d)'
end

"""
Generates linear map E_ such that E_(ρ) → 1 ⊗ ρ
"""
function E_(id_dim, ρ_dim)
    M = spzeros(Float64,id_dim^2*ρ_dim^2,ρ_dim^2)
    for m in 0:ρ_dim-1
        for n in 0:ρ_dim-1
            for k in 0:id_dim-1
                M += kron(ket(k,id_dim),ket(m,ρ_dim),ket(k,id_dim),ket(n,ρ_dim))*kron(bra(m,ρ_dim),bra(n,ρ_dim))
            end
        end
    end
    return M
end

"""
Permutes the elements of a matrix so that it transforms the column
major representation of a linear map L into a matrix C that is
positive iff L is completely positive, Hermitian iff L maps
(vectorized) Hermitian matrices to Hermitian matrices. In other words,
it corresponds to the Choi-Jamiolkoski isomorphism.
"""
function involution(m)
    dsq = size(m,1) # we assume the matrix is square
    d   = Int(sqrt(dsq))
    return reshape(permutedims(reshape(m,(d,d,d,d)),[2,4,1,3]),(dsq,dsq))
end

let # wat13b
    global dnorm
    local prev_dx, M

    prev_dx = -1

    """
    Computes the diamond norm of a linear superoperator `L` (i.e., a
    linear transformation of operators). The superoperator must be
    represented in column major form. In other words, it must be given
    by a matrix that, when multiplying a vectorized (column major)
    operator, it should result in the vectorized (column major)
    representation of the result of the transformation.
    """
    function dnorm(L)
        J = involution(L)

        dx = size(J,1) |> sqrt |> x -> round(Int,x)
        dy = dx

        if prev_dx != dx
            M = E_(dy,dx)
            prev_dx = dx
        end

        Jr = real(J)
        Ji = imag(J)

        Xr  = Variable(dy*dx, dy*dx)
        Xi  = Variable(dy*dx, dy*dx)
        ρ0r = Variable(dx, dx)
        ρ0i = Variable(dx, dx)
        ρ1r = Variable(dx, dx)
        ρ1i = Variable(dx, dx)

        prob = maximize( retrϕ( ϕ(Jr,Ji)'*ϕ(Xr,Xi) ) )

        prob.constraints += trace(ρ0r) == 1
        prob.constraints += trace(ρ0i) == 0
        prob.constraints += trace(ρ1r) == 1
        prob.constraints += trace(ρ1i) == 0

        Mρ0r = reshape(M * vec(ρ0r), dy*dx, dy*dx)
        Mρ0i = reshape(M * vec(ρ0i), dy*dx, dy*dx)
        Mρ1r = reshape(M * vec(ρ1r), dy*dx, dy*dx)
        Mρ1i = reshape(M * vec(ρ1i), dy*dx, dy*dx)

        prob.constraints += isposdef( ϕ(ρ0r,ρ0i) )

        prob.constraints += isposdef( ϕ(ρ1r,ρ1i) )

        prob.constraints += isposdef( ϕ( [ Mρ0r Xr ; Xr' Mρ1r ], [ Mρ0i Xi ; -Xi' Mρ1i ] ) )

        solve!(prob)

        if prob.status != :Optimal
            println("DNORM error.")
            println("Input: $(L)")
            println("Input's Choi spectrum: $(eigvals(liou2choi(L)))")
            error("Could not compute the diamond norm.")
        end

        return prob.optval
    end
end


text:
# PSET 3: Estimating Mixture Models #
# Danielle Nemschoff #

# Note: code adapted from Thibaut's R code version and Florian Oswald #

# ------------------------------------------------------------------------------------------------------------- #
# Set Up ------------------------------------------------------------------------------------------------------ #
# ------------------------------------------------------------------------------------------------------------- #
# models of wage formaation with unobserved heterogeneity in firms and workers #

# firms have types l: 1,...,L
# workers have types by k: 1,...,K

# worker i is of type k and workers for firm l in a certain period
# wages are drawn from the normal distribution N(μkl, σkl)

# one firm type and 2 worker types

# add time dimension

using(RCall)
using(Distributions)
using(StatsFuns)
using(Random)
using(Plots)
using(DataFrames)

# ------------------------------------------------------------------------------------------------------------- #
# Data Generator Function ------------------------------------------------------------------------------------- #
# ------------------------------------------------------------------------------------------------------------- #
# sets true mixture params
# creates a distribution mixture mmodoel
# draws n random realizations

function sdata(n; t=40, doplot = false)
    Random.seed!(6)

    # true values
    μ = [2.0,5.0]
    σ = [0.5,0.7]
    pk = [0.3,0.7]

    m = MixtureModel([Normal(μ[i], σ[i]) for i in 1:2], pk)
    if doplot
        plot(
            plot(m,linewidth=2),
            plot(m,linewidth=2, fill=(0,:red,0.5), components = false, title="Mixture"),dpi = 300
            )
        #savefig("mixtures.png")
    end
    y = rand(m,n)

    return Dict(:y => y, :μ => μ, :σ => σ, :pk => pk)
end

# ------------------------------------------------------------------------------------------------------------- #
# Estimation Function ----------------------------------------------------------------------------------------- #
# ------------------------------------------------------------------------------------------------------------- #
# 1. Take vector y
# 2. set (the same) wrong starting values
# 3. run the EM algorithm for iters iterations to find true values of
#    proportion weights pk, means \mu and variances \sigma for each component.

function bm_jl(y::Vector{Float64}; t= 40, iters=100)
    y = fake_data[:y]
    # poor starting values
    μ = [4.0,6.0]
    σ = [1.0,1.0]
    pk = [0.5,0.5]

    N = length(y)
    K = length(μ)

    # initialize objects
    L = zeros(N,K)
    p = similar(L)

    for it in 1:iters

        dists = [Normal(μ[ik], σ[ik] ) for ik in 1:K]

        # evaluate likelihood for each type
        for i in 1:N
            for k in 1:K
                L[i,k] = log(pk[k]) + logpdf.(dists[k], y[i])
            end
        end
        L
        # get posterior of each type
        p[:,:] = exp.(L .- logsumexp(L))

        # with p in hand, update
        pk[:] .= vec(sum(p,dims=1) ./ N)
        μ[:] .= vec(sum(p .* y, dims = 1) ./ sum(p, dims = 1))
        σ[:] .= vec(sqrt.(sum(p .* (y .- μ').^2, dims = 1) ./ sum(p, dims = 1)))
        pk
    end
    return Dict(:pk => pk, :μ => μ, :σ => σ)
end

# ------------------------------------------------------------------------------------------------------------- #
# Simulate ---------------------------------------------------------------------------------------------------- #
# ------------------------------------------------------------------------------------------------------------- #

fake_data = sdata(10, t=40)

bm_jl(fake_data[:y], t = 40, iters=100)
using GaussianMixtures
function bm_jl_GMM(y::Vector{Float64};iters=100)
    gmm = GMM(2,1)  # initialize an empty GMM object
    # stick in our starting values
    gmm.μ[:,1] .= [4.0;6.0]
    gmm.Σ[:,1] .= [1.0;1.0]
    gmm.w[:,1] .= [0.5;0.5]

    # run em!
    em!(gmm,y[:,:],nIter = iters)
    return gmm
end

res = bm_jl_GMM(fake_data[:y])
res
res.μ
res.w

# ------------------------------------------------------------------------------------------------------------- #
# Implement on PSID ------------------------------------------------------------------------------------------- #
# ------------------------------------------------------------------------------------------------------------- #

#data = data.table(read.dta13("~/Dropbox/Documents/Teaching/ECON-24030/lectures-laborsupply/homeworks/data/AER_2012_1549_data/output/data4estimation.dta"))


text:
# Main design principles:
# TODO: manifest design

abstract type Layer end
abstract type Activation <: Layer end

abstract type Loss end


text:
const DIMC = 500
const DIMR = 500
const DIMD = 1
const MAXCOLOR = 255

const ROTATIONINDICES = Dict(:x => [2 2; 2 3; 3 2; 3 3],
                             :y => [1 1; 3 1; 1 3; 3 3],
                             :z => [1 1; 1 2; 2 1; 2 2])

const CIRCSTEPS = 360
const SPHSTEPS  = 100
const TORUSTEPS = 100
const BEZSTEPS  = 100
const HERMSTEPS = 100

const CURVES = Dict(:bezier  => [-1  3 -3  1;  3 -6  3  0; -3  3  0  0;  1  0  0  0],
                    :hermite => [ 2 -2  1  1; -3  3 -2 -1;  0  0  1  0;  1  0  0  0])

const CAMBIENT = [ 50,  50,  50]

const LIGHTS = [Dict(:location => [ 0.50,  0.75,  1.00], :color => [  0, 255, 255]),
                Dict(:location => [-0.50,  0.75,  1.00], :color => [255,   0,   0])]

const VIEWVEC = [ 0.00,  0.00,  1.00]

const REFLECTION = Dict(:ambient  => [0.10, 0.10, 0.10],
                        :diffuse  => [0.50, 0.50, 0.50],
                        :specular => [0.50, 0.50, 0.50])
const SPEC_EXP = 16

const ANIM_DELAY = 10

max_stars_repo_name: AnAverageHuman/jlgraphics

max_issues_repo_name: AnAverageHuman/jlgraphics

max_forks_repo_name: AnAverageHuman/jlgraphics


text:
using Documenter
using HDF5
using H5Zblosc
using H5Zbzip2
using H5Zlz4
using H5Zzstd
using MPI  # needed to generate docs for parallel HDF5 API

# Used in index.md to filter the autodocs list
not_low_level_api(m::Method) = !endswith(String(m.file), "src/api.jl")
not_low_level_api(f::Function) = all(not_low_level_api, methods(f))
not_low_level_api(o) = true

# defined on separate page
not_low_level_api(::Type{<:HDF5.Properties}) = false

# Manually-defined low-level API (in source file src/api_helpers.jl)
not_low_level_api(::typeof(HDF5.API.h5p_get_class_name)) = false
not_low_level_api(::typeof(HDF5.API.h5t_get_member_name)) = false
not_low_level_api(::typeof(HDF5.API.h5t_get_tag)) = false

makedocs(;
    sitename="HDF5.jl",
    modules=[HDF5, H5Zblosc, H5Zbzip2, H5Zlz4, H5Zzstd],
    authors="Mustafa Mohamad <mus-m@outlook.com> and contributors",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaIO.github.io/HDF5.jl",
        assets=String[],
        sidebar_sitename=false,
    ),
    pages=[
        "Home" => "index.md",
        "Interface" => [
            "dataset.md",
            "properties.md",
            "filters.md"
        ],
        "Low-level library bindings" => "api_bindings.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaIO/HDF5.jl.git",
    push_preview=true,
)


text:
module Direct

using jInv.Mesh
using jInv.Utils
using jInv.InverseSolve
using ShapeReconstructionPaLS.Utils
using ShapeReconstructionPaLS.ParamLevelSet
using MAT
using SparseArrays
using Distributed
import jInv.ForwardShare.getData
import jInv.ForwardShare.getSensTMatVec
import jInv.ForwardShare.getSensMatVec
import jInv.ForwardShare.ForwardProbType

## The rotation planes are:
# first rotation (theta) is for planes x1 and x2
# second rotation (phi) is for planes x1 and x3
# Sampling is done by summing planes of x1 and x2 for all values of x3.

export DirectParam, getDirectParam
mutable struct DirectParam <: ForwardProbType
    Mesh      				:: RegularMesh
	ndipsAll				:: Int64
	workerSubIdxs
    #theta_phi_rad			:: Array{Float64,2}
	#b						:: Array{Float64,2}
	samplingBinning			:: Int64
	method					:: String
	Jacobian				:: SparseMatrixCSC{Float32,Int32}
	#S						:: SparseMatrixCSC{Float32,Int32}
end

function getDirectParamInternal(Mesh::RegularMesh,ndipsAll::Int64,workerSubIdxs,samplingBinning::Int64,method)
	return DirectParam(Mesh, ndipsAll, workerSubIdxs, samplingBinning, method,spzeros(Float32,Int32,0,0));
end

function getDirectParam(Mesh::RegularMesh,samplingBinning::Int64,numWorkers::Int64,method = MATBased)
	## This function does use the parallel mechanism of jInv (i.e., returns a RemoteChannel), even if numWorkers=1.
	if numWorkers > nworkers()
		numWorkers = nworkers();
	end
	#SourcesSubInd = Array{Array{Int64,1}}(numWorkers);
	SourcesSubInd = Array{Array{Int64}}(undef,numWorkers);
	ActualWorkers = workers();
	if numWorkers < nworkers()
		ActualWorkers = ActualWorkers[1:numWorkers];
	end
	#pFor   = Array{RemoteChannel}(numWorkers)
	pFor   = Array{RemoteChannel}(undef,numWorkers);
	i = 1; nextidx() = (idx=i; i+=1; idx)
	idx=i;
	ndips  = 1;
	# send out jobs
	@sync begin
		for p = ActualWorkers
			@async begin
				while true
					idx = nextidx();
					if idx > numWorkers
						break
					end
					I_p = getIndicesOfKthWorker(numWorkers,idx,ndips);
					# println("Sending ",collect(I_p)," to worker ", p);
					# find src and rec on mesh
					#SEE IF CORRECT!!!!!!!!!!!!!!!!
					tmp = initRemoteChannel(getDirectParamInternal,p, Mesh,ndips,I_p,samplingBinning,method);
					pFor[1]  = initRemoteChannel(getDirectParamInternal,p, Mesh,ndips,I_p,samplingBinning,method);
					wait(pFor[1]);

				end
			end
		end
	end
	return pFor # Array of Remote Refs
end

import jInv.Utils.clear!
function clear!(pFor::DirectParam)
	pFor.SampleMat = spzeros(0);
	clear!(pFor.Mesh);
	return pFor;
end

#include("generateDipMatrix.jl");
#include("prepareSyntheticDipData.jl");
include("getData.jl");
include("sensitivityFuncs.jl");
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
module ImageAxes

using Base: @pure, tail
using Reexport, Colors, SimpleTraits, MappedArrays

# maybe return to "@reexport AxisArrays" if AxisArrays is fixed

import AxisArrays
using AxisArrays: AxisArray, Axis, axisnames, axisvalues, axisdim, atindex, atvalue, collapse
export AxisArray, Axis, axisnames, axisvalues, axisdim, atindex, atvalue, collapse

@reexport using ImageCore

export # types
    HasTimeAxis,
    IndexAny,
    IndexIncremental,
    StreamingContainer,
    TimeAxis,
    StreamIndexStyle,
    # functions
    colordim,
    arraydata,
    getindex!,
    istimeaxis,
    timeaxis,
    timedim

"""
    TimeAxis{Ax}

A trait (from SimpleTraits) indicating whether axis `Ax` corresponds
to time. This decision is based on the symbol-name given to `Ax`. For
example, the following declares that all `Axis{:time}` objects
correspond to time:

    @traitimpl TimeAxis{Axis{:time}}

This definition has already been made in ImageAxes, but you can add
new names as well.
"""
@traitdef TimeAxis{X}

@traitimpl TimeAxis{Axis{:time}}

# Note: any axis not marked as a TimeAxis is assumed to correspond to
# space. It might be useful to allow users to add other possibilities,
# but this is not currently possible with SimpleTraits.

"""
    timeaxis(A)

Return the time axis, if present, of the array `A`, and `nothing` otherwise.
"""
@inline timeaxis(A::AxisArray) = _timeaxis(A.axes...)
timeaxis(A::AbstractArray) = nothing
timeaxis(A::AbstractMappedArray) = timeaxis(parent(A))
@traitfn _timeaxis(ax::Ax, axes...) where {Ax<:Axis; !TimeAxis{Ax}} = _timeaxis(axes...)
@traitfn _timeaxis(ax::Ax, axes...) where {Ax<:Axis;  TimeAxis{Ax}} = ax
_timeaxis() = nothing

"""
    istimeaxis(ax)

Test whether the axis `ax` corresponds to time.
"""
istimeaxis(ax::Axis) = istimeaxis(typeof(ax))
@traitfn istimeaxis(::Type{Ax}) where {Ax<:Axis; !TimeAxis{Ax}} = false
@traitfn istimeaxis(::Type{Ax}) where {Ax<:Axis;  TimeAxis{Ax}} = true

@traitdef HasTimeAxis{X}
"""
    HasTimeAxis{AA}

A trait for testing whether type `AA` has a time axis. Time axes must
be declared before use.

# Examples

```julia
using ImageAxes, SimpleTraits

# Declare that all axes named `:time` are time axes
@traitimpl TimeAxis{Axis{:time}}

# Define functions that dispatch on AxisArrays that may or may not have time axes
@traitfn got_time{AA<:AxisArray;  HasTimeAxis{AA}}(img::AA) = "yep, I've got time"
@traitfn got_time{AA<:AxisArray; !HasTimeAxis{AA}}(img::AA) = "no, I'm too busy"

julia> A = AxisArray(1:5, Axis{:time}(1:5));

julia> got_time(A)
"yep, I've got time"

julia> A = AxisArray(1:5, Axis{:x}(1:5));

julia> got_time(A)
"no, I'm too busy"
```
"""
HasTimeAxis

axtype(::Type{AxisArray{T,N,D,Ax}}) where {T,N,D,Ax} = Ax
axtype(A::AxisArray) = axtype(typeof(A))

Base.@pure function SimpleTraits.trait(t::Type{HasTimeAxis{AA}}) where AA<:AxisArray
    axscan = map(S->istimeaxis(S), axtype(AA).parameters)
    any(axscan) ? HasTimeAxis{AA} : Not{HasTimeAxis{AA}}
end

# Specializations to preserve the AxisArray wrapper
function ImageCore.permuteddimsview(A::AxisArray, perm)
    axs = AxisArrays.axes(A)
    AxisArray(permuteddimsview(A.data, perm), axs[[perm...]]...)
end
function ImageCore.channelview(A::AxisArray)
    Ac = channelview(A.data)
    _channelview(A, Ac)
end
# without extra dimension:
_channelview(A::AxisArray{C,N}, Ac::AbstractArray{T,N}) where {C,T,N} = AxisArray(Ac, AxisArrays.axes(A)...)
# with extra dimension: (bug: the type parameters shouldn't be necessary, but julia 0.5 dispatches incorrectly without them)
_channelview(A::AxisArray{C,M}, Ac::AbstractArray{T,N}) where {C,T,M,N} = AxisArray(Ac, Axis{:color}(axes(Ac,1)), AxisArrays.axes(A)...)

### Image properties based on traits ###

"""
    timedim(img) -> d::Int

Return the dimension of the array used for encoding time, or 0 if not
using an axis for this purpose.

Note: if you want to recover information about the time axis, it is
generally better to use `timeaxis`.
"""
timedim(img::AxisArray{T,N}) where {T,N} = _timedim(filter_time_axis(AxisArrays.axes(img), ntuple(identity, Val(N))))
_timedim(dim::Tuple{Int}) = dim[1]
_timedim(::Tuple{}) = 0

ImageCore.nimages(img::AxisArray) = _nimages(timeaxis(img))
_nimages(::Nothing) = 1
_nimages(ax::Axis) = length(ax)

function colordim(img::AxisArray)
    d = _colordim(1, AxisArrays.axes(img))
    d > ndims(img) ? 0 : d
end
_colordim(d, ax::Tuple{Ax,Vararg{Any}}) where {Ax<:Axis{:color}} = d
_colordim(d, ax::Tuple{Any,Vararg{Any}}) = _colordim(d+1, tail(ax))
_colordim(d, ax::Tuple{}) = d+1

ImageCore.pixelspacing(img::AxisArray) = map(step, filter_space_axes(AxisArrays.axes(img), axisvalues(img)))

ImageCore.spacedirections(img::AxisArray) = ImageCore._spacedirections(pixelspacing(img))

ImageCore.coords_spatial(img::AxisArray{T,N}) where {T,N} = filter_space_axes(AxisArrays.axes(img), ntuple(identity, Val(N)))

ImageCore.spatialorder(img::AxisArray) = filter_space_axes(AxisArrays.axes(img), axisnames(img))

ImageCore.size_spatial(img::AxisArray)    = filter_space_axes(AxisArrays.axes(img), size(img))
ImageCore.indices_spatial(img::AxisArray) = filter_space_axes(AxisArrays.axes(img), axes(img))

arraydata(img::AxisArray) = img.data
@deprecate data arraydata

### Utilities for writing "simple algorithms" safely ###

# Check that the time dimension, if present, is last
@traitfn function ImageCore.assert_timedim_last(img::AA) where {AA<:AxisArray; HasTimeAxis{AA}}
    ax = AxisArrays.axes(img)[end]
    istimeaxis(ax) || error("time dimension is not last")
    nothing
end
@traitfn ImageCore.assert_timedim_last(img::AA) where {AA<:AxisArray; !HasTimeAxis{AA}} = nothing

### Convert ###
function Base.convert(::Type{Array{C,n}},
                      img::AxisArray{C,n}) where {C<:Colorant,n}
    copyto!(Array{C}(undef, size(img)), img)
end
function Base.convert(::Type{Array{Cdest,n}},
                      img::AxisArray{Csrc,n}) where {Cdest<:Colorant,n,Csrc<:Colorant}
    copyto!(Array{ccolor(Cdest, Csrc)}(undef, size(img)), img)
end

### StreamingContainer ###

checknames(axnames, ::Type{P}) where {P} = checknames(axnames, axisnames(P))
@noinline function checknames(axnames, parentnames::Tuple{Symbol,Vararg{Symbol}})
    mapreduce(x->in(x, parentnames), &, axnames) || throw(DimensionMismatch("names $axnames are not included among $parentnames"))
    nothing
end

"""
    A = StreamingContainer{T}(f!, parent, streamingaxes::Axis...)

An array-like object possessing one or more axes for which changing "slices" may
be expensive or subject to restrictions. A canonical example would be
an AVI stream, where addressing pixels within the same frame is fast
but jumping between frames might be slow.

Here's a simple example of dividing by the mean of each slice of an image before returning values.

    A = AxisArrays.AxisArray(reshape(1:36, 3, 3, 4))

    function f!(buffer, slice)
        meanslice = mean(slice)
        buffer .= slice./meanslice
    end

    B = StreamingContainer{Float64}(f!, A, AxisArrays.axes(A)[3])

    julia> A[:,:,1]
    3×3 AxisArray{Int64,2,Array{Int64,2},Tuple{Axis{:row,Base.OneTo{Int64}},Axis{:col,Base.OneTo{Int64}}}}:
     1  4  7
     2  5  8
     3  6  9

    julia> B[:,:,1]
    3×3 Array{Float64,2}:
     0.2  0.8  1.4
     0.4  1.0  1.6
     0.6  1.2  1.8

The user-provided `f!` function should take arguments:

    f!(buffer, slice)

Where `buffer` will be an empty array that can hold a slice of your series, and `slice` will hold the current input slice.

It's worth noting that `StreamingContainer` is *not* a subtype of
`AbstractArray`, but that much of the array interface (`eltype`,
`ndims`, `axes`, `size`, `getindex`, and `IndexStyle`) is
supported. A StreamingContainer `A` can be built from an AxisArray,
but it may also be constructed from other "parent" objects, even
non-arrays, as long as they support the same functions. In either
case, the parent should also support the standard AxisArray functions
`axes`, `axisnames`, `axisvalues`, and `axisdim`; this support will be
extended to the `StreamingContainer`.

Additionally, a StreamingContainer `A` supports

    getindex!(dest, A, axt::Axis{:time}, ...)

to obtain slices along the streamed axes (here it is assumed that
`:time` is a streamed axis of `A`). You can implement this directly
(dispatching on the parameters of `A`), or (if the parent is an
`AbstractArray`) rely on the fallback

    A.getindex!(dest, view(parent, axs...))

where `A.getindex! = f!` as passed as an argument at construction. `dest` should
have dimensionality `ndims(parent)-length(streamingaxes)`.

Optionally, define [`StreamIndexStyle(typeof(parent),typeof(f!))`](@ref).
"""
struct StreamingContainer{T,N,streamingaxisnames,P,GetIndex}
    getindex!::GetIndex
    parent::P
end
function StreamingContainer{T}(f!::Function, parent, axs::Axis...) where T
    N = ndims(parent)
    axnames = axisnames(axs...)
    checknames(axnames, typeof(parent))
    StreamingContainer{T,N,axnames,typeof(parent),typeof(f!)}(f!, parent)
end

Base.parent(S::StreamingContainer) = S.parent
Base.axes(S::StreamingContainer) = axes(S.parent)
Base.size(S::StreamingContainer)    = size(S.parent)
Base.axes(S::StreamingContainer, d) = axes(S.parent, d)
Base.size(S::StreamingContainer, d)    = size(S.parent, d)

AxisArrays.axes(S::StreamingContainer) = AxisArrays.axes(parent(S))
AxisArrays.axisnames(S::StreamingContainer)  = axisnames(AxisArrays.axes(S)...)
AxisArrays.axisvalues(S::StreamingContainer) = axisvalues(AxisArrays.axes(S)...)
function AxisArrays.axisdim(S::StreamingContainer, ::Type{Axis{name}}) where name
    isa(name, Int) && return name <= ndims(S) ? name : error("axis $name greater than array dimensionality $(ndims(S))")
    names = axisnames(S)
    idx = findfirst(isequal(name), names)
    idx isa Nothing && error("axis $name not found in array axes $names")
    idx
end
AxisArrays.axisdim(S::StreamingContainer, ax::Axis) = axisdim(S, typeof(ax))
AxisArrays.axisdim(S::StreamingContainer, ::Type{Axis{name,T}}) where {name,T} = axisdim(S, Axis{name})

ImageCore.nimages(S::StreamingContainer) = _nimages(timeaxis(S))
ImageCore.coords_spatial(S::StreamingContainer{T,N}) where {T,N} =
    filter_space_axes(AxisArrays.axes(S), ntuple(identity, Val(N)))
ImageCore.spatialorder(S::StreamingContainer) = filter_space_axes(AxisArrays.axes(S), axisnames(S))
ImageCore.size_spatial(img::StreamingContainer)    = filter_space_axes(AxisArrays.axes(img), size(img))
ImageCore.indices_spatial(img::StreamingContainer) = filter_space_axes(AxisArrays.axes(img), axes(img))
function ImageCore.assert_timedim_last(S::StreamingContainer)
    istimeaxis(AxisArrays.axes(S)[end]) || error("time dimension is not last")
    nothing
end

Base.eltype(::Type{StreamingContainer{T,N,names,P,GetIndex}}) where {T,N,names,P,GetIndex} = T
Base.eltype(S::StreamingContainer) = eltype(typeof(S))
Base.ndims(::Type{StreamingContainer{T,N,names,P,GetIndex}}) where {T,N,names,P,GetIndex} = N
Base.ndims(S::StreamingContainer) = ndims(typeof(S))
Base.length(S::StreamingContainer) = prod(size(S))

streamingaxisnames(::Type{StreamingContainer{T,N,names,P,GetIndex}}) where {T,N,names,P,GetIndex} =
    names
streamingaxisnames(S::StreamingContainer) = streamingaxisnames(typeof(S))

timeaxis(S::StreamingContainer) = _timeaxis(AxisArrays.axes(S)...)

@inline function getindex!(dest, S::StreamingContainer, axs::Axis...)
    all(ax->isstreamedaxis(ax,S), axs) || throw(ArgumentError("$axs do not coincide with the streaming axes $(streamingaxisnames(S))"))
    _getindex!(dest, S.getindex!, parent(S), axs...)
end

@inline function getindex!(dest, S::StreamingContainer, I...)
    axs = getslicedindices(S, I)
    all(x->isa(x, Colon), filter_notstreamed(I, S)) || throw(ArgumentError("positional indices must use `:` for any non-streaming axes"))
    _getindex!(dest, S.getindex!, parent(S), axs...)
end

# _getindex! just makes it easy to specialize for particular parent or function types
@inline function _getindex!(dest, f!, P::AbstractArray, axs::Axis...)
    v = view(P, axs...)
    f!(dest, v)
end

function isstreamedaxis(ax::Axis{name},
                        S::StreamingContainer{T,N,saxnames}) where {name,T,N,saxnames}
    in(name, saxnames)
end

sliceindices(S::StreamingContainer) = filter_notstreamed(axes(S), S)
sliceaxes(S::StreamingContainer) = filter_notstreamed(AxisArrays.axes(S), S)
getslicedindices(S::StreamingContainer, I) = filter_streamed(map((ax, i) -> ax(i), AxisArrays.axes(S), I), S)

filter_streamed(inds, S)    = _filter_streamed(inds, AxisArrays.axes(S), S)
filter_notstreamed(inds, S) = _filter_notstreamed(inds, AxisArrays.axes(S), S)
filter_streamed(inds::Tuple{Axis,Vararg{Axis}}, S)    = _filter_streamed(inds, inds, S)
filter_notstreamed(inds::Tuple{Axis,Vararg{Axis}}, S) = _filter_notstreamed(inds, inds, S)

@generated function _filter_streamed(a, axs::NTuple{N,Axis}, S::StreamingContainer) where N
    inds = findall(x->x in streamingaxisnames(S), axisnames(axs.parameters...))
    Expr(:tuple, Expr[:(a[$i]) for i in inds]...)
end
@generated function _filter_notstreamed(a, axs::NTuple{N,Axis}, S::StreamingContainer) where N
    inds = findall(x->x in streamingaxisnames(S), axisnames(axs.parameters...))
    inds = setdiff(1:N, inds)
    Expr(:tuple, Expr[:(a[$i]) for i in inds]...)
end

@inline function Base.getindex(S::StreamingContainer{T,N}, inds::Vararg{Union{Colon,Base.ViewIndex},N}) where {T,N}
    tmp = similar(Array{T}, sliceindices(S))
    getindex!(tmp, S, getslicedindices(S, inds)...)
    tmp[filter_notstreamed(inds, S)...]
end
@inline function Base.getindex(S::StreamingContainer, ind1::Axis, inds_rest::Axis...)
    axs = sliceaxes(S)
    tmp = AxisArray(Array{eltype(S)}(undef, map(length, axs)), axs)
    inds = (ind1, inds_rest...)
    getindex!(tmp, S, _filter_streamed(inds, inds, S)...)
    getindex_rest(tmp, _filter_notstreamed(inds, inds, S))
end
getindex_rest(tmp, ::Tuple{}) = tmp
getindex_rest(tmp, inds) = tmp[inds...]

"""
    style = StreamIndexStyle(A)

A trait that indicates the degree of support for indexing the streaming axes
of `A`. Choices are [`IndexAny()`](@ref) and
[`IndexIncremental()`](@ref) (for arrays that only permit advancing
the time axis, e.g., a video stream from a webcam). The default value
is `IndexAny()`.

This should be specialized for the type rather than the instance. For
a StreamingContainer `S`, you can define this trait via

```julia
StreamIndexStyle(::Type{P}, ::typeof(f!)) = IndexIncremental()
```

where `P = typeof(parent(S))`.
"""
abstract type StreamIndexStyle end

"""
    IndexAny()

Indicates that an axis supports full random-access indexing.
"""
struct IndexAny <: StreamIndexStyle end
"""
    IndexIncremental()

Indicates that an axis supports only incremental indexing, i.e., from `i` to `i+1`.
This is commonly used for the temporal axis with media streams.
"""
struct IndexIncremental <: StreamIndexStyle end

StreamIndexStyle(::Type{A}) where {A<:AbstractArray} = IndexAny()
StreamIndexStyle(A::AbstractArray) = StreamIndexStyle(typeof(A))

StreamIndexStyle(::Type{StreamingContainer{T,N,axnames,P,GetIndex}}) where {T,N,axnames,P,GetIndex} = StreamIndexStyle(P, GetIndex)
StreamIndexStyle(::Type{P},::Type{GetIndex}) where {P,GetIndex} = IndexAny()

StreamIndexStyle(S::StreamingContainer) = StreamIndexStyle(typeof(S))

### Low level utilities ###

filter_space_axes(axes::NTuple{N,Axis}, items::NTuple{N,Any}) where {N} =
    _filter_space_axes(axes, items)
@inline @traitfn _filter_space_axes(axes::Tuple{Ax,Vararg{Any}}, items) where {Ax<:Axis;  TimeAxis{Ax}} =
    _filter_space_axes(tail(axes), tail(items))
@inline @traitfn _filter_space_axes(axes::Tuple{Ax,Vararg{Any}}, items) where {Ax<:Axis; !TimeAxis{Ax}} =
    (items[1], _filter_space_axes(tail(axes), tail(items))...)
_filter_space_axes(::Tuple{}, ::Tuple{}) = ()
@inline _filter_space_axes(axes::Tuple{Ax,Vararg{Any}}, items) where {Ax<:Axis{:color}} =
    _filter_space_axes(tail(axes), tail(items))

filter_time_axis(axes::NTuple{N,Axis}, items::NTuple{N}) where {N} =
    _filter_time_axis(axes, items)
@inline @traitfn _filter_time_axis(axes::Tuple{Ax,Vararg{Any}}, items) where {Ax<:Axis; !TimeAxis{Ax}} =
    _filter_time_axis(tail(axes), tail(items))
@inline @traitfn _filter_time_axis(axes::Tuple{Ax,Vararg{Any}}, items) where {Ax<:Axis;  TimeAxis{Ax}} =
    (items[1], _filter_time_axis(tail(axes), tail(items))...)
_filter_time_axis(::Tuple{}, ::Tuple{}) = ()

# summary: print color types & fixed-point types compactly
function AxisArrays._summary(io, A::AxisArray{T,N}) where {T<:Union{Fractional,Colorant},N}
    print(io, "$N-dimensional AxisArray{")
    if T<:Colorant
        ColorTypes.colorant_string_with_eltype(io, T)
    else
        ColorTypes.showcoloranttype(io, T)
    end
    println(io, ",$N,...} with axes:")
end

end # module


text:
module InteractiveFormatter

using JuliaFormatter: format_text
using REPL: LineEdit, LineEditREPL, REPL
using REPL.LineEdit: InputAreaState, edit_insert, input_string, refresh_line

function format_prompt(s, _...)
    code = try
        format_text(String(take!(copy(LineEdit.buffer(s)))))
    catch err
        println(stderr)
        @error "Error in `format_prompt`" exception = (err, catch_backtrace())
        println(stderr)
        refresh_line(s)
        return
    end
    take!(LineEdit.buffer(s))  # empty prompt
    edit_insert(s, code)
    refresh_line(s)
end

"""
    afterreplinit(f)

Like `atreplinit` but triggers `f` even after REPL is initialized when
it is called.
"""
function afterreplinit(f)
    # See: https://github.com/JuliaLang/Pkg.jl/blob/v1.0.2/src/Pkg.jl#L338
    function wrapper(repl)
        if isinteractive() && repl isa REPL.LineEditREPL
            f(repl)
        end
    end
    if isdefined(Base, :active_repl)
        wrapper(Base.active_repl)
    else
        atreplinit() do repl
            @async begin
                wait_repl_interface(repl)
                wrapper(repl)
            end
        end
    end
end

function wait_repl_interface(repl)
    for _ in 1:20
        try
            repl.interface.modes[1].keymap_dict
            return
        catch
        end
        sleep(0.05)
    end
end

function init_repl(repl)
    formatter_keymap = Dict{Any,Any}("^V" => format_prompt)

    main_mode = repl.interface.modes[1]
    main_mode.keymap_dict = LineEdit.keymap_merge(
        main_mode.keymap_dict,
        LineEdit.keymap([formatter_keymap]),
    )
end
# See: https://github.com/JuliaInterop/RCall.jl/blob/master/src/RPrompt.jl

function __init__()
    afterreplinit(init_repl)
end

end # module


text:
module CovidAgent
using Agents
mutable struct GenericAgent <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int} #region
    type::Int # 1 - nsec+1 HH, nsec + 2 - 2*(nsec+1) Firm, man / service / Food/public
    age:: Int # 1 young, 2 old
    cor::Int #infection status: 1 susceptible, 2 infected, 3 recovered, 4 infected with mutant
    vaccinated::Bool # true, if agent has been vaccinated
    cortime:: Int # time of infection
    infect:: Int # number of people infected by this agent
    emp:: Int #employer ID, 0 is unemployed
    shorttime:: Bool # on short time
    shorttimecounter:: Int # Number of weeks on short time
    homeoffice:: Bool # homeoffice
    worker:: Dict{Int,CovidAgent.GenericAgent} # list of firm employees
    workershorttime:: Dict{Int,CovidAgent.GenericAgent} # list of firm employees on short time
    income:: Float64 # average income
    saving::Float64 # savings of HH and firms
    consbud::Array{Float64,1} # consumption budget for HH for three cons types
    prevconsbud::Array{Float64,1} # consumption budget from previous period
    shopday::Array{Int,1} # array of shopping days in the week
    shoprepeat:: Array{Bool,1} # indicator repetition of shopping trip possible
    product::Float64 # productiviy of firm
    stock::Float64 # stock of product of firm
    demand::Float64 #firm's demand in previous week
    demandexp::Float64 # firm’s demand estimate
    rev::Float64 # firm revenues
    pastrev::Array{Float64,1} # array of last 4 weekly revenes
    markup::Float64 # firm markup
    fixedcosts::Float64 # firms fixed costs
end

end

# function for roulette wheel selection based on fitness vector vec
function conselect(vect::Array{Float64,1})
    vecn = size(vect)[1]
    vecsum = sum(vect)
    vecr = rand()*vecsum
    for i=1:vecn
        if sum(vect[1:i]) >= vecr
            return i
            break
        end
    end
end

function getagent(id::Int)
    return covidmodel[id]
end

# one step in the dynamics
function model_step!(covidmodel)
    global hh, firms, lochh, unemplist, shorttimelist, cas, daycounter, noinf, noinfmut, unemp, consumption, weeklyconsumption, shoppers, nsec
    global regcons, regemp, actshoplist, pubacc, empcount, tau,trigger, proftax, totaldiv, totalfixedcosts
    global virus, unempcount, shorttimecount, oldcount, R0count,gdp, genhomeoffice, noinftraj
    global bankruptcypossible, inactivefirms, sumbailouts, totalaccounts, divperhh, totalsav

    global contact_count_traj, contact_work_traj, contact_social_traj, contact_shop_traj
    global contact_count, contact_work, contact_social, contact_shop, contact_num, nosuscep

    global updatefixedcosts, bankruptcount, curinfcount, curinfregcount

    # reset consumption counter
    consumption = zeros(nsec, k1, k2)
    #reset total savings
    totalsav = 0

    #determine mortaliy
    noinf = 0
    noinfmut = 0
    aginf = 0
    for i in union(hh[1],hh[2])
        if (getagent(i).cor ==2) || (getagent(i).cor ==4)
            noinf = noinf +1
        end
        if (getagent(i).cor ==4)
            noinfmut = noinfmut +1
        end
        if getagent(i).cor > 1
            aginf +=1
        end
    end
    push!(totinftraj,aginf)  # number of infected and recoverd
    mort = (min(hcap,icufrac*noinf) /(icufrac*noinf))*mortl  + (1 - min(hcap,icufrac*noinf) / (icufrac*noinf))*morth

    # determine first which of the young hh is travelling
    #lochhytemp = deepcopy(lochh[1,:,:]) # temp list of young at locations
    postemp = [(0,0) for i=1:nhh] # temp list of young positions

    lochhytemp = [Dict{Int,CovidAgent.GenericAgent}() for a = 1:2, i=1:k1, j=1:k2]
    for a=1:2
        for i=1:k1
            for j=1:k2
                for k in keys(lochh[a,i,j])
                    lochhytemp[a,i,j][k] = lochh[a,i,j][k]
                end
            end
        end
    end

    for i in hh[1]  # only young households travel
        agent = getagent(i)
        if rand() < ptravel
            posh = agent.pos
            posn = [rand(1:k1),rand(1:k2)]
            delete!(lochhytemp[posh[1],posh[2]],i)
            lochhytemp[posn[1],posn[2]][i] = agent
            postemp[i] = (posn[1],posn[2]) # save new position
        else
            postemp[i] =agent.pos # stick to default position
        end
    end

    for i in hh[2] # old agents do not travel
        postemp[i] =getagent(i).pos # stick to default position
    end

    # firm decisions once a week
    if daycounter == 0
        totaldiv = 0. # sums dividends of all firms
        totalfixedcosts = 0. # sums fixed costs of all firms
        proftax = 0. # total profit tax
        gdp = 0.
        inactivefirms = 0
        sumbailouts = 0.0
        badbankacc = 0.0
        actshoplist = [[] for tt=1:nsec, i=1:k1, j2 = 1:k2]

        totprice = zeros(nsec)
        pricecount = zeros(nsec)
        permlist = shuffle(firms)
        for id in permlist
            firm = getagent(nhh+id)
            if firm.type-(nsec+1) <= nsec # not a public sector office
                numworker = length(firm.worker)

                # firms collect profits, pay dividends, wages, taxes, update demand estimate
                profit = firm.rev - wage[firm.type-(nsec+1)]*numworker - firm.fixedcosts
                totalfixedcosts += firm.fixedcosts

                # determine dividend
                for revt = 1:3
                    firm.pastrev[revt] = firm.pastrev[revt+1]
                end
                firm.pastrev[4] = firm.rev
                divthr = savtar[firm.type-(nsec+1)]*mean(firm.pastrev) # dividend threhshold
                if firm.saving + (1-tau)*max(0,profit) > divthr
                    div = (1-tau)*max(0,profit)
                else
                    div = divrate*(1-tau)*max(0,profit)
                end
                totaldiv += div  # add to total dividends
                proftax += tau*max(0,profit)
                firm.saving += profit-tau*max(0,profit) - div

                # pricing
                if weeklyconsumption[firm.type-(nsec+1), firm.pos[1], firm.pos[2]] > 0
                    marketshare = firm.demand / weeklyconsumption[firm.type-(nsec+1), firm.pos[1], firm.pos[2]]
                    firm.markup = etmin[firm.type-(nsec+1)] + marketshare*(etmax[firm.type-(nsec+1)]-etmin[firm.type-(nsec+1)])
                else
                    firm.markup = etmin[firm.type-(nsec+1)]
                end

                firm.demandexp = (1 - rhodem)*firm.demandexp + rhodem*firm.demand # update demand expectation
                firm.demand = 0. # reset weekly demand
                firm.rev = 0. # reset weekly revenue

                # production, hiring/firing
                if firm.stock > 0
                    # mormal plan
                    prodplan = max(0,(firm.demandexp*(1 + buff[firm.type-(nsec+1)]) - (1-de[firm.type-(nsec+1)])*firm.stock))
                else
                    # stronger expansion if stock is empty
                    prodplan = max(0,(firm.demandexp*(1 + zerostockboost*buff[firm.type-(nsec+1)]) - (1-de[firm.type-(nsec+1)])*firm.stock))
                end
                labdem = ceil(prodplan / firm.product) - numworker

                if firm.saving < 0 && bailoutprogram
                    sumbailouts += - firm.saving
                    firm.saving = 0
                end

                if firm.saving < 0 && bankruptcypossible
                    prodplan = 0
                    labdem = -numworker
                end

                for (id2, stworker) in firm.workershorttime
                    stworker.shorttimecounter += 1

                    # move worker from short term to unemployed if:
                    # - short time program has ended
                    # - worker excedded maximum time in short time
                    # - firm is bankrupt
                    if (!shorttimeprogram || stworker.shorttimecounter > shorttimeperiod || (firm.saving < 0 && bankruptcypossible))
                        delete!(firm.workershorttime,id2)
                        delete!(shorttimelist[firm.type-(nsec+1),firm.pos[1],firm.pos[2]],id2)
                        unemplist[firm.type-(nsec+1),stworker.pos[1],stworker.pos[2]][id2] = stworker
                        stworker.shorttime = false
                        stworker.shorttimecounter = 0
                        stworker.emp = 0
                    end
                end

                if labdem < 0 # random firing
                    firelist = shuffle(collect(keys(firm.worker)))[1:Int(-labdem)]
                    for id2 in firelist
                        fireagent = getagent(id2)
                        if shorttimeprogram && rand() < shorttimeprob && !(firm.saving < 0 && bankruptcypossible)
                            delete!(firm.worker,id2)
                            firm.workershorttime[id2] = fireagent
                            shorttimelist[firm.type-(nsec+1),fireagent.pos[1],fireagent.pos[2]][id2] = fireagent # put on regional short time list
                            fireagent.shorttime = true
                            fireagent.shorttimecounter = 0
                        else
                            delete!(firm.worker,id2)
                            unemplist[firm.type-(nsec+1),fireagent.pos[1],fireagent.pos[2]][id2] = fireagent # put on regional unemploymentlist
                            fireagent.emp = 0
                        end
                    end
                    numworker += labdem
                elseif labdem > 0 # hiring
                    numhire = min(labdem)
                    numhirefromshorttime = min(numhire, length(firm.workershorttime))
                    numhirefrommarket = min(numhire - numhirefromshorttime, length(unemplist[firm.type-(nsec+1),firm.pos[1],firm.pos[2]])) # check if firm is rationed on labor market

                    if numhirefromshorttime > 0
                        hirelist = shuffle(collect(keys(firm.workershorttime)))[1:Int(numhirefromshorttime)] # search among unemployd from firms short time list
                        for id2 in hirelist
                            hireagent = getagent(id2)
                            firm.worker[id2] = hireagent
                            delete!(firm.workershorttime,id2) # remove from firms short time list
                            delete!(shorttimelist[firm.type-(nsec+1),firm.pos[1],firm.pos[2]],id2) # remove from regional unemploymentlist
                            hireagent.shorttime = false
                            hireagent.shorttimecounter = 0
                        end
                        numworker += numhirefromshorttime
                    end

                    if numhirefrommarket > 0
                        hirelist = shuffle(collect(keys(unemplist[firm.type-(nsec+1),firm.pos[1],firm.pos[2]])))[1:Int(numhirefrommarket)] # search among unemployd in region
                        for id2 in hirelist
                            hireagent = getagent(id2)
                            firm.worker[id2] = hireagent
                            delete!(unemplist[firm.type-(nsec+1),firm.pos[1],firm.pos[2]],id2) # remove from regional unemploymentlist
                            hireagent.emp = id
                        end
                        numworker += numhirefrommarket
                    end
                end

                firm.stock = (1-de[firm.type-(nsec+1)])*firm.stock + numworker*firm.product # update product stock with production
                gdp += (numworker*wage[firm.type-(nsec+1)]+firm.fixedcosts)*(1+firm.markup)/(1 - buff[firm.type-(nsec+1)]*de[firm.type-(nsec+1)])
                if firm.stock > 0
                    push!(actshoplist[firm.type-(nsec+1),firm.pos[1],firm.pos[2]],id) # if sotcvk positive add to actve shops
                end
            else # public sector office
                gdp += length(firm.worker)*wage[nsec+1]*(1+fcvcratio[nsec+1])
            end
            if updatefixedcosts  # update fixed costs
                firm.fixedcosts = fcvcratio[firm.type-(nsec+1)] * wage[firm.type-(nsec+1)]*length(firm.worker) # fixed costs here prop to vc
            else
                if firm.fixedcosts < fcvcratio[firm.type-(nsec+1)] * wage[firm.type-(nsec+1)]*length(firm.worker) # adjust fc only if firm expands
                    firm.fixedcosts = fcvcratio[firm.type-(nsec+1)] * wage[firm.type-(nsec+1)]*length(firm.worker)
                end
            end
            if firm.type < 2*(nsec+1) # temp delete
                totprice[firm.type-(nsec+1)] += (wage[firm.type-(nsec+1)] + firm.fixedcosts / (length(firm.worker)+10^(-5)))*(1+firm.markup)/ (firm.product*(1 - buff[firm.type-(nsec+1)]*de[firm.type-(nsec+1)]))
                pricecount[firm.type-(nsec+1)] += 1
            end
            if firm.saving < 0 && bankruptcypossible
                badbankacc -= firm.saving
                setdiff!(firms,(firm.id-nhh))
                setdiff!(actshoplist[firm.type-(nsec+1),firm.pos[1],firm.pos[2]],(firm.id-nhh))
                setdiff!(locf[firm.type-(nsec+1),firm.pos[1],firm.pos[2]],(firm.id-nhh))
                bankruptcount[firm.type-(nsec+1)] += 1
                kill_agent!(firm,covidmodel)
            end
        end # end firm loop
        #println(totprice ./ pricecount)

        weeklyconsumption = zeros(nsec, k1, k2)

        # government pays unemployment benefits, psensions and public sector, adjusts tax rate
        pubsur = proftax + tau*sum(empcount.*wage) + tau*(sum(empcount)+sum(unempcount)+sum(shorttimecount)+oldcount)*divperhh- (1-tau)*(sum(unempcount.*unempbenefit) + sum(shorttimecount.*shorttimewage) + oldcount*pension) - empcount[nsec+1]*wage[nsec+1]
        pubacc += pubsur

        # Cost of bailouts to the public account
        pubacc -= (sumbailouts + badbankacc)
        # --- End of tax period ---

        tauhat = max(0,(sum(unempcount.*unempbenefit) + sum(shorttimecount.*shorttimewage) + oldcount*pension + empcount[nsec+1]*wage[nsec+1] - debtfrac*pubacc) /
        ((proftax / tau) + sum(empcount.*wage) + sum(unempcount.*unempbenefit) + sum(shorttimecount.*shorttimewage) + oldcount*pension +
        (sum(empcount)+sum(unempcount)+sum(shorttimecount)+oldcount)*divperhh)) #taxrate needed to balance budget

        tau = (1-tauadj)*tau + tauadj*tauhat

        # reset employment, unemployment and shorttime counter
        empcount = zeros(nsec+1)
        unempcount = zeros(nsec)
        shorttimecount = zeros(nsec)
        oldcount = 0

        # calculate div per hh
        divperhh = (totaldiv+totalfixedcosts) / size(union(hh[1],hh[2]))[1]  # equasl share of total dividends for each hh
    end # end daycounter == 0

    #reset consumption counter
    consumption = zeros(nsec,k1,k2)

    #resent local shoppers list
    shoppers = [Dict{Int,CovidAgent.GenericAgent}() for tt=1:nsec, i=1:k1, j=1:k2]

    curinfcount = zeros(2)
    curinfregcount = zeros(2,k1,k2)
    # Household loop
    permlist = shuffle(union(hh[1],hh[2]))
    for ii in permlist
        agent = getagent(ii)
        # potential death of agent

        if (agent.cor == 2) || (agent.cor == 4) # if infected
            if rand() < mort[agent.age]
                setdiff!(hh[agent.age],agent.id)
                delete!(lochh[agent.age,agent.pos[1],agent.pos[2]],agent.id)
                delete!(lochh[agent.age,agent.pos[1],agent.pos[2]],agent.id)
                cas[agent.pos[1],agent.pos[2]][agent.age] += 1
                push!(R0count,agent.infect)
                if agent.age == 1
                    delete!(lochhytemp[postemp[agent.id][1],postemp[agent.id][2]],agent.id)
                    if agent.emp ==0
                        delete!(unemplist[agent.type,agent.pos[1],agent.pos[2]],agent.id)
                    else
                        employer = getagent(nhh+agent.emp)
                        if agent.shorttime
                            delete!(employer.workershorttime,agent.id)
                            delete!(shorttimelist[agent.type,agent.pos[1],agent.pos[2]],agent.id)
                        else
                            delete!(employer.worker,agent.id)
                        end
                    end
                end
                # transfer savings of HH to random other HH
                permlisttemp = shuffle(union(hh[1],hh[2]))
                heiragent = getagent(permlisttemp[1])
                heiragent.saving += agent.saving
                kill_agent!(agent,covidmodel)
                continue
            end
            agent.cortime = agent.cortime-1 # if not dying reduce recovery time by one
            if agent.cortime ==0
                agent.cor = 3
                push!(R0count,agent.infect)
            end
        end

    # hh determine consumption budgets and count unemp once a week
        if daycounter ==0
            agent.saving += sum(agent.consbud)   # non-spent consumption budget is added to savings
            actualincome = 0
            if agent.emp > 0 # employed

                if agent.shorttime
                    shorttimecount[agent.type] += 1
                    actualincome = (1-tau)*(shorttimewage[agent.type] + divperhh)
                    agent.income = (1-rhoinc)*agent.income +  rhoinc*(1-tau)*(shorttimewage[agent.type] + divperhh) # update income expectations
                else
                    empcount[agent.type] += 1
                    actualincome = (1-tau)*(wage[agent.type] + divperhh)
                    agent.income = (1-rhoinc)*agent.income +  rhoinc*(1-tau)*(wage[agent.type] + divperhh) # update income expectations

                end
            elseif agent.age==1 # unemployed young
                unempcount[agent.type] += 1
                agent.income = (1-rhoinc)*agent.income +  rhoinc*(1-tau)*(unempbenefit[agent.type] + divperhh)
                actualincome = (1-tau)*(unempbenefit[agent.type] + divperhh)
            else # old agent
                oldcount += 1
                agent.income = (1-rhoinc)*agent.income +  rhoinc*(1-tau)*(pension +  divperhh)
                actualincome = (1-tau)*(pension + divperhh)
            end

            totcons = max(0,agent.income + kap*(agent.saving - phi*agent.income))
            agent.consbud = fraccons*totcons

            if essentialsec > 0
                if agent.consbud[essentialsec] < (1-essentialadj) * agent.prevconsbud[essentialsec]
                    if totcons < (1-essentialadj) * agent.prevconsbud[essentialsec]
                        for i = 1:size(agent.consbud, 1)
                            agent.consbud[i] = 0
                        end
                        agent.consbud[essentialsec] = totcons
                    else
                        remainingcons = totcons - (1-essentialadj) * agent.prevconsbud[essentialsec]
                        tempfraccons = deepcopy(fraccons)
                        tempfraccons[essentialsec] = 0
                        sumfrac = sum(tempfraccons)
                        for i = 1:size(tempfraccons, 1)
                            tempfraccons[i] = tempfraccons[i] / sumfrac
                        end
                        agent.consbud = tempfraccons*remainingcons
                        agent.consbud[essentialsec] = (1-essentialadj) * agent.prevconsbud[essentialsec]
                    end
                end
            end

            agent.prevconsbud = deepcopy(agent.consbud)
            agent.saving += actualincome - totcons
            agent.shopday = rand(0:6,nsec) # determine shopping days for three sectors
            agent.shoprepeat = [true, true, true] # consumer repeats shopping once if rationed
            # calculate total savings
            totalsav += agent.saving
        end # end daycounter = 0

        for tt = 1:nsec
            if agent.shopday[tt] == daycounter && rand() < pshop[agent.age][tt] # agent goes shopping
                numshoph = min(size(actshoplist[tt,postemp[agent.id][1],postemp[agent.id][2]])[1],numshops[tt])
                shoppers[tt,postemp[agent.id][1],postemp[agent.id][2]][agent.id] = agent # if agent shops in region where she currently is
                if numshoph > 0 # only if there are potential suppliers
                    storelist = shuffle(actshoplist[tt,postemp[agent.id][1],postemp[agent.id][2]])[1:numshoph] # list of supplier agent observes
                    pricelist = Float64[]
                    for storeid in storelist
                        store = getagent(nhh + storeid)
                        push!(pricelist, exp( - gac*log((wage[store.type-(nsec+1)] + store.fixedcosts / (length(store.worker)+10^(-5)))*(1+store.markup)/ (store.product*(1 - buff[store.type-(nsec+1)]*de[store.type-(nsec+1)]) ))))
                    end
                    store = getagent(nhh + storelist[conselect(pricelist)]) # selected supplier
                    price = (wage[store.type-(nsec+1)] + store.fixedcosts / (length(store.worker)+10^(-5)))*(1+store.markup)/ (store.product*(1 - buff[store.type-(nsec+1)]*de[store.type-(nsec+1)]))
                    quant = min(agent.consbud[tt] / price, store.stock) # minimum of demand and stock
                    agent.consbud[tt] -= quant*price # update agent's consumer budget
                    store.stock -= quant # update stock of supplier
                    if store.stock < 0.0000001
                        setdiff!(actshoplist[store.type-(nsec+1),store.pos[1],store.pos[2]],(store.id-nhh)) # if stock is empty make store inactive
                    end
                    consumption[tt,postemp[agent.id][1],postemp[agent.id][2]] += quant # update total consumption
                    weeklyconsumption[tt,postemp[agent.id][1],postemp[agent.id][2]] += quant
                    store.rev += quant*price # update firm revenues
                    store.demand += quant # update demand
                end
                if agent.consbud[tt] > 0.000001 && agent.shoprepeat[tt]
                    agent.shopday[tt] += 1 # try again next day if rationed
                    agent.shoprepeat[tt] = false # only one repetition possible
                end
            end
        end
    end # end HH loop

    totalaccounts = pubacc
    for i in firms
        totalaccounts += getagent(nhh+i).saving
    end
    for i in union(hh[1], hh[2])
        totalaccounts += getagent(i).saving
    end
    # calculate unemployment
    unemp = 1 - (sum(empcount)+sum(shorttimecount))/size(hh[1])[1]

    #update total consumption
    push!(regcons, consumption)

    # start counting contatcs
    contact_count = 0
    contact_work = 0
    contact_social = 0
    contact_shop = 0
    contact_num = 0
    nosuscep = 0 # count susceptible

### infection loop from perspective of infected
    if virus
        for ii in union(hh[1],hh[2])
            agent = getagent(ii)
            if (agent.cor == 2) || (agent.cor == 4)
                global curinfcount[agent.age] += 1
                global curinfregcount[agent.age,agent.pos[1], agent.pos[2]] += 1
            end
            # possible infection
            if ((agent.cor==2)|| (agent.cor == 4))  && (agent.cortime < trec - corlatent) && (agent.cortime > trec - corlatent- corinf)
                contact_num += 1
                # first infection at the workplace
                if agent.emp > 0 && !agent.shorttime
                    if !agent.homeoffice || !genhomeoffice || rand() > phomeoffice # no home-office for agent
                        employer = getagent(nhh+agent.emp)
                        idlist = employer.worker
                        nummeet = rand(0:workmeet[employer.type-(nsec+1)])
                        for i in shuffle(collect(keys(idlist)))[1:min(nummeet,length(idlist))]
                            partner = getagent(i)
                            contact_work = contact_work + 1 # add contatcs
                            if !partner.vaccinated && (rand() < pinf[Int(agent.cor/2)]) && partner.cor == 1 && !(partner.homeoffice && genhomeoffice && rand() < phomeoffice) && !partner.shorttime
                                partner.cor=agent.cor
                                partner.cortime = trec
                                agent.infect +=1
                            end
                        end
                    end
                end

                # second infection through consumption
                for tt=1:nsec
                    if length(shoppers[tt,postemp[agent.id][1],postemp[agent.id][2]]) > 2
                        if sum(collect(keys(shoppers[tt,postemp[agent.id][1],postemp[agent.id][2]])) .== agent.id) > 0 # agent shopped in sector tt at that day
                            nummeet = rand(0:shopmeet[tt])
                            shopmeeth = Int(ceil(nummeet * length(shoppers[tt,postemp[agent.id][1],postemp[agent.id][2]]) * 7 / (length(lochh[1,postemp[agent.id][1],postemp[agent.id][2]]) + length(lochh[2,postemp[agent.id][1],postemp[agent.id][2]])))) # number of contacts when shopping, equal to nummeet if 1/7 of local pop is shopping
                            shopcontact = rand(collect(keys(delete!(shoppers[tt,postemp[agent.id][1],postemp[agent.id][2]],agent.id))),shopmeeth) #random selection from list of shoppers in same region
                            for i in shopcontact
                                partner = getagent(i)
                                contact_shop = contact_shop + 1 # add contacts
                                if !partner.vaccinated && (rand() < pinf[Int(agent.cor/2)]) &&  (partner.cor ==1)
                                    partner.cor=agent.cor
                                    partner.cortime = trec
                                    agent.infect +=1
                                end
                            end
                            shoppers[tt,postemp[agent.id][1],postemp[agent.id][2]][agent.id] = agent
                        end
                    end
                end

                # third infection through other social contacts
                if agent.age ==1 # young agent
                        # meet young
                        helplist = collect(keys(delete!(lochhytemp[postemp[agent.id][1],postemp[agent.id][2]],agent.id))) # list of potential young contacts in the region
                        if socialmaxyy >= 1
                            numcontact = rand(0:socialmaxyy) # number of contacts
                        else
                            if rand() < mean([0,socialmaxyy])
                                numcontact = 1
                            else
                                numcontact = 0
                            end
                        end
                        if numcontact > 0
                            helpcontact = shuffle(helplist)[1:min(numcontact,size(helplist)[1])] #list of contacts
                            for i in helpcontact
                                partner = getagent(i)
                                contact_social = contact_social + 1 # add contacts
                                if !partner.vaccinated && (rand() < pinf[Int(agent.cor/2)]) &&  (partner.cor ==1)
                                    partner.cor=agent.cor
                                    partner.cortime = trec
                                    agent.infect +=1
                                end
                            end
                        end
                        lochhytemp[postemp[agent.id][1],postemp[agent.id][2]][agent.id] = agent

                        #meet old
                        helplist = collect(keys(lochh[2,postemp[agent.id][1],postemp[agent.id][2]])) # list of potential old in the region
                        if socialmaxyo >= 1
                            numcontact = rand(0:socialmaxyo) # number of contacts
                        else
                            if rand() < mean([0,socialmaxyo])
                                numcontact = 1
                            else
                                numcontact = 0
                            end
                        end
                        if numcontact > 0
                            helpcontact = shuffle(helplist)[1:min(numcontact,size(helplist)[1])] #list of contacts
                            for i in helpcontact
                                partner = getagent(i)
                                contact_social = contact_social + 1 # add contacts
                                if !partner.vaccinated && (rand() < pinf[Int(agent.cor/2)]) &&  (partner.cor ==1)
                                    partner.cor=agent.cor
                                    partner.cortime = trec
                                    agent.infect +=1
                                end
                            end
                        end
                else # old agent
                    # meet young
                    helplist = collect(keys(lochhytemp[agent.pos[1],agent.pos[2]])) # list of potential young contacts in the region
                    if socialmaxoy >= 1
                        numcontact = rand(0:socialmaxoy) # number of contacts
                    else
                        if rand() < mean([0,socialmaxoy])
                            numcontact = 1
                        else
                            numcontact = 0
                        end
                    end
                    if numcontact > 0
                        helpcontact = shuffle(helplist)[1:min(numcontact,size(helplist)[1])] #list of contacts
                        for i in helpcontact
                            partner = getagent(i)
                            contact_social = contact_social + 1 # add contacts
                            if !partner.vaccinated && (rand() < pinf[Int(agent.cor/2)]) &&  (partner.cor ==1)
                                partner.cor=agent.cor
                                partner.cortime = trec
                                agent.infect +=1
                            end
                        end
                    end
                    #meet old
                    helplist = collect(keys(delete!(lochh[2,agent.pos[1],agent.pos[2]],agent.id))) # list of potential old in the region
                    if socialmaxoo >= 1
                        numcontact = rand(0:socialmaxoo) # number of contacts
                    else
                        if rand() < mean([0,socialmaxoo])
                            numcontact = 1
                        else
                            numcontact = 0
                        end
                    end
                    if numcontact > 0
                        helpcontact = shuffle(helplist)[1:min(numcontact,size(helplist)[1])] #list of contacts
                        for i in helpcontact
                            partner = getagent(i)
                            contact_social = contact_social + 1 # add contacts
                            if !partner.vaccinated && (rand() < pinf[Int(agent.cor/2)]) &&  (partner.cor ==1)
                                partner.cor=agent.cor
                                partner.cortime = trec
                                agent.infect +=1
                            end
                        end
                    end
                    lochh[2,agent.pos[1],agent.pos[2]][agent.id] = agent
                end
            end # if cor = 1

        end # summing over all young hh

        # collect contact data
        contact_count = contact_work + contact_social + contact_shop
        push!(contact_count_traj,contact_count/(contact_num + 10^(-5)))
        push!(contact_work_traj,contact_work/(contact_num + 10^(-5)))
        push!(contact_social_traj,contact_social/(contact_num + 10^(-5)))
        push!(contact_shop_traj,contact_shop/(contact_num + 10^(-5)))

    end # virus from infected perspective

    if daycounter == 0
        daycounter = 7
    end
    # reduce daycounter
    daycounter -=1
    return
end


text:
function parseMsg(msg::AbstractString)
    pairs = split(msg, ";")
    d = Dict{AbstractString, Any}()
    for i = 1:length(pairs)

        if strip(pairs[i]) == ""
            continue
        end

        key, val = split(pairs[i], ":")
        key = String(key)
        val = String(val)
        d[key] = val
    end
    return d
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# makie
using CairoMakie
let
    x = 1:10
    fig = Figure(resolution = (600, 400), backgroundcolor = :honeydew,
          font = "CMU Serif")
    ax = Axis(fig[1, 1], xlabel = L"x", ylabel = L"x^2",
            xlabelsize = 22, ylabelsize = 22)
    lines!(x, x .^ 2; linewidth = 2, linestyle = ".-", label = L"x^2",
      color = :black)
    limits!(0, 10, 0, 100)
    axislegend(; merge = true, framecolor = :mediumpurple, position = :lt)
    current_figure()
    display(fig)
end;

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations