text:
function rebin2D(df::DataFrame, _prec::Float64 = 0.001, _dx = 0.001)
    df_rebinned = DataFrame(E1 = Float64[],
                            minE = Float64[], maxE = Float64[],
                            minG = Float64[], maxG = Float64[],
                            a = Float64[], b = Float64[])

    for e in unique(df.E1)
        df_rebinned = vcat(df_rebinned, rebin1D(df[df.E1 .== e, 2], df[df.E1 .== e, 3], e, _prec, _dx))
    end

    return df_rebinned
end

function rebin1D(_x::Vector{Float64}, _y::Vector{Float64}, _E1::Float64, _prec::Float64 = 0.001, _dx = 0.001)
    e0=1 #initialize index of the first point in fit
    ef=2 #initialize index of the second point in fit

    a,b = get_line_params(0.0, 0.0, _x[e0], _y[e0])                 # first approximation is made from origin to first point in spectrum

    df_rebinned = DataFrame(E1 = _E1,                               # initialized DataFrame with first approximation
                            minE = 0.0, maxE = _x[e0],
                            minG = 0.0, maxG = _y[e0],
                            a = a, b = b )

    if length(_x) == 1                                              # if there is only one point in the vector, first point is approximated from (0,0) to (x1,y1)
        a,b = get_line_params(_x[e0], _y[e0], _x[e0] + _dx, 0.0)    # second line is approximated from (x1,y1) to (x1+dx, 0) creating a triangle with area equal to
                                                                    # rectangle with sides dx and y1 (dx is the step in x-direction...we assume x-array is evenly spaced)
        push!(df_rebinned, [_E1,
                            _x[e0], _x[e0] + _dx,
                            _y[e0], 0.0     ,
                            a     , b       ] )
        return df_rebinned
    end

    a,b = get_line_params(_x[e0], _y[e0], _x[ef], _y[ef])           # first fit points 1,2
    Δy  = abs( (get_line_point(_x[ef], a, b) - _y[ef]) / _y[ef] )   # Δy of points 1,2

    for _ in 1:length(_x)

        while Δy <= _prec && ef < length(_x)
            Δy = abs((get_line_point(_x[ef], a, b) - _y[ef])/_y[ef])
            ef += 1
        end

        if ef != length(_x)
            a, b = get_line_params(_x[e0], _y[e0], _x[ef-1], _y[ef-1])

            push!(df_rebinned, [_E1,
                                _x[e0], _x[ef-1],                   # refit line with previous point
                                _y[e0], _y[ef-1],                   # unless it's last point of the series
                                a     , b       ] )
        else
            a, b = get_line_params(_x[e0], _y[e0], _x[ef], _y[ef])
            push!(df_rebinned, [_E1,
                                _x[e0], _x[ef],
                                _y[e0], _y[ef],                     # for last point line stops at the endpoint
                                a     , b       ] )
        end

        if ef >= length(_x)                                         # if the next point is more than the size of data, check ends
            break
        end

        e0  = ef -1                                                 # last point from fit becomes new point in next fit
        a,b = get_line_params(_x[e0], _y[e0], _x[ef], _y[ef])       # find new line from the next set of points
        Δy  = abs((get_line_point(_x[ef], a, b) - _y[ef])/_y[ef])   # find new Δy

    end

    a, b = get_line_params(_x[ef], _y[ef], _x[ef]+_dx, 0.0)       # finally to add very last approximation by going from (xf, yf) to (xf+dx, 0.0)
    push!(df_rebinned, [_E1,
                        _x[ef], _x[ef]+_dx,
                        _y[ef], 0.0         ,                       # for last point line stops at the endpoint
                        a     , b           ] )

    return df_rebinned
end

function normalize2D!(df_rebinned::DataFrame, thickness = 0.001) # CHANGE THIS, NO NEED TO PREALLOCATE df_normed. JUST CHANGE THE INPUT DF and RETURN IT!
    volume_rebinned = get_total_volume(df_rebinned, thickness)

    df_normed = select(df_rebinned,
                    :E1            => :E1,
                    :minE          => :minE,
                    :maxE          => :maxE,
                    :minG=> ByRow(g->g/volume_rebinned) => :minG,
                    :maxG=> ByRow(g->g/volume_rebinned) => :maxG,
                    :a   => ByRow(a->a/volume_rebinned) => :a,
                    :b   => ByRow(b->b/volume_rebinned) => :b,
                    )

    return df_normed  #return select(df_rebinned,
                    #  :E1            => :E1,
                    #  :minE          => :minE,
                    #  :maxE          => :maxE,
                    #  :minG=> ByRow(g->g/volume_rebinned) => :minG,
                    #  :maxG=> ByRow(g->g/volume_rebinned) => :maxG,
                    #  :a   => ByRow(a->a/volume_rebinned) => :a,
                    #  :b   => ByRow(b->b/volume_rebinned) => :b,
                    #  )
end

"""
#### function ```get_cdf(df::DataFrame, thickness::Real = 0.001)```
<br>

    Description of ```get_cdf```
    ------------------------------
Calculate CDF for pdf = a + bx. CDF is cummulative density function defined as integral from -inf to x of pdf dx.

"""
function get_cdf(df::DataFrame, thickness::Real = 0.001)
    n = nrow(df)
    cdf = Vector{Float64}(undef, n)
    cdf[1] = get_integral_linear(df[1,2], df[1,3], df[1,6], df[1,7])*thickness

    for i in 2:n
        cdf[i] = cdf[i-1] + get_integral_linear(df[i,2], df[i,3], df[i,6], df[i,7])*thickness
    end
    return cdf
end

function get_cdf!(df::DataFrame, thickness::Real = 0.001)
    n = nrow(df)
    cdf = Vector{Float64}(undef, n)
    cdf[1] = get_integral_linear(df[1,2], df[1,3], df[1,6], df[1,7])*thickness

    for i in 2:n
        cdf[i] = cdf[i-1] + get_integral_linear(df[i,2], df[i,3], df[i,6], df[i,7])*thickness
    end

    @transform! df :cdf = cdf
end

ext: jl


text:
workspace()
using HoTRG_lemon.SimulatorModule
using HoTRG_lemon.SpinModule
using HoTRG_lemon.LatticeModule
using TensorOperations
using TensorMatrices_lemon

dimM = 4
trotterparameter = 0.01
trotteriteration = 40

isinginfo = SpinInfo("quantum_ising_2",1.0e-13,0.0)
trotterinfo = TrotterInfo(trotterparameter, trotteriteration)
ising = QuantumIsingModel(isinginfo, trotterinfo)
lattice = buildLattice("quantum_2d_fractal", ising)
simulator = Quantum2dFractalSimulator(lattice,  4, 0)
println(getInititeration(simulator))
# simulator = buildSimulator(lattice,dimM)

# fieldrange = linspace(0.1, 3.5, 20)
# simulatorQuantum(fieldrange, simulator; filename = "data_dim" * string(dimM)  * ".txt")


text:
sess = Gallium.FakeMemorySession([(UInt64(0),reinterpret(UInt8,[UInt64(0)]))],
  Gallium.X86_64.X86_64Arch(), nothing)
@test Gallium.load(sess, RemotePtr{UInt64}(0)) == UInt64(0)
@test_throws ErrorException Gallium.load(sess, RemotePtr{UInt64}(8))


text:
module CellMaps

using Test
using Gridap
using Gridap.Helpers
using Gridap.CachedArrays
using Gridap.CellValues: _test_iter_cell_value
using Gridap.CellValues: _test_index_cell_value
using Gridap.CellValues: _eq

export CellMap
export IterCellMap
export IndexCellMap
import Gridap: evaluate
export test_iter_cell_map
export test_index_cell_map
export test_index_cell_map_with_index_arg
export test_iter_cell_map_with_index_result

import Base: iterate
import Base: length
import Base: size
import Base: getindex
import Base: show

"""
Abstract object that traverses a set of cells and at every cell returns a
`Map{S,M,T,N}`
"""
const IterCellMap{S,M,T,N,R<:Map{S,M,T,N}} = IterCellValue{R}

"""
Abstract array indexed by cells that returns a `Map{S,M,T,N}`
"""
const IndexCellMap{S,M,T,N,C,R<:Map{S,M,T,N}} = IndexCellValue{R,C}

"""
Abstract object that for a given cell index returns a `Map{S,M,T,N}`
"""
const CellMap{S,M,T,N} = Union{IterCellMap{S,M,T,N},IndexCellMap{S,M,T,N}}

# evaluation

"""
Return the cellwise maps of a `CellMap` on a cellwise set of points
"""
function evaluate(cm::CellMap{S,M,T,N},ca::CellArray{<:S,M}) where {S,M,T,N}
  IterCellMapValue(cm,ca)
end

function evaluate(cm::IndexCellMap{S,M,T,N},ca::IndexCellArray{<:S,M}) where {S,M,T,N}
  IndexCellMapValue(cm,ca)
end

struct IterCellMapValue{T,N,V,A} <: IterCellArray{T,N,CachedArray{T,N,Array{T,N}}}
  cm::V
  ca::A
end

function IterCellMapValue(cm::CellMap{S,M,T,N},ca::CellArray{<:S,M}) where {S,M,T,N}
  V = typeof(cm)
  A = typeof(ca)
  IterCellMapValue{T,N,V,A}(cm,ca)
end

length(v::IterCellMapValue) = length(v.ca)

@inline function iterate(v::IterCellMapValue{T,N}) where {T,N}
  cmnext = iterate(v.cm)
  canext = iterate(v.ca)
  r = CachedArray(T,N)
  _iterate(cmnext,canext,r)
end

@inline function iterate(v::IterCellMapValue,state)
  cmstate, castate, r = state
  cmnext = iterate(v.cm,cmstate)
  canext = iterate(v.ca,castate)
  _iterate(cmnext,canext,r)
end

@inline function _iterate(cmnext,canext,r)
  if cmnext === nothing; return nothing end
  if canext === nothing; return nothing end
  m, cmstate = cmnext
  a, castate = canext
  s = return_size(m,size(a))
  setsize!(r,s)
  evaluate!(m,a,r)
  state = (cmstate,castate,r)
  (r, state)
end

struct IndexCellMapValue{T,N,V,A} <: IndexCellArray{T,N,CachedArray{T,N,Array{T,N}},1}
  cm::V
  ca::A
  r::CachedArray{T,N,Array{T,N}}
end

function IndexCellMapValue(cm::IndexCellMap{S,M,T,N},ca::IndexCellArray{<:S,M}) where {S,M,T,N}
  V = typeof(cm)
  A = typeof(ca)
  r = CachedArray(T,N)
  IndexCellMapValue{T,N,V,A}(cm,ca,r)
end

length(v::IndexCellMapValue) = length(v.ca)

size(v::IndexCellMapValue) = (length(v),)

function getindex(v::IndexCellMapValue,i::Integer)
  m = v.cm[i]
  a = v.ca[i]
  s = return_size(m,size(a))
  setsize!(v.r,s)
  evaluate!(m,a,v.r)
  v.r
end

# Testers

function test_iter_cell_map(
  m::CellMap{S,M,T,N},
  a::CellArray{<:S,M},
  b::AbstractArray{<:AbstractArray{T,N}}) where {S,M,T,N}

  @test length(m) == length(b)
  @test length(m) == length(a)

  c = evaluate(m,a)
  _test_iter_cell_value(c,b)

  for (mi,ai,bi) in zip(m,a,b)
    @assert isa(mi,Map{S,M,T,N})
    @assert _eq(evaluate(mi,ai),bi)
    @assert typeof(mi) == eltype(m)
  end

end

function test_index_cell_map(
  m::IndexCellMap{S,M,T,N},
  a::CellArray{<:S,M},
  b::AbstractArray{<:AbstractArray{T,N}}) where {S,M,T,N}

  @test size(m) == size(b)

  @test IndexStyle(m) == IndexStyle(b)

  test_iter_cell_map(m,a,b)

  for (i,ai) in enumerate(a)
    mi = m[i]
    bi = b[i]
    @assert isa(mi,Map{S,M,T,N})
    @assert _eq(evaluate(mi,ai),bi)
    @assert typeof(mi) == eltype(m)
  end

end

function test_index_cell_map_with_index_arg(
  m::IndexCellMap{S,M,T,N},
  a::IndexCellArray{<:S,M},
  b::AbstractArray{<:AbstractArray{T,N}}) where {S,M,T,N}

  test_index_cell_map(m,a,b)

  c = evaluate(m,a)
  @test isa(c,IndexCellArray)

  @test length(c) == length(m)
  @test length(c) == length(a)

  for i in 1:length(c)
    mi = m[i]
    bi = b[i]
    ci = c[i]
    @assert _eq(bi,ci)
  end

  for i in 1:length(m)
    mi = m[i]
    ai = a[i]
    bi = b[i]
    @assert isa(mi,Map{S,M,T,N})
    @assert _eq(evaluate(mi,ai),bi)
    @assert typeof(mi) == eltype(m)
  end

end

function test_iter_cell_map_with_index_result(
  m::CellMap{S,M,T,N},
  a::CellArray{<:S,M},
  b::AbstractArray{<:AbstractArray{T,N}}) where {S,M,T,N}

  test_iter_cell_map(m,a,b)

  c = evaluate(m,a)

  _test_index_cell_value(c,b)

end

function show(io::IO,self::CellMap)
  s = "$(nameof(typeof(self))) object"
  print(io,s)
end

function show(io::IO,::MIME"text/plain",self::CellMap)
  show(io,self)
  print(io,":")
  print(io,"\n length: $(length(self))")
end

end # module CellMaps


text:
using IndirectCoStateInit
using Documenter

DocMeta.setdocmeta!(IndirectCoStateInit, :DocTestSetup, :(using IndirectCoStateInit); recursive=true)

makedocs(;
    modules=[IndirectCoStateInit],
    authors="Grant Hecht",
    repo="https://github.com/GrantHecht/IndirectCoStateInit.jl/blob/{commit}{path}#{line}",
    sitename="IndirectCoStateInit.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://GrantHecht.github.io/IndirectCoStateInit.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/GrantHecht/IndirectCoStateInit.jl",
    devbranch="main",
)


text:
using DelayDiffEq, OrdinaryDiffEq, RecursiveArrayTools, Test

function f_lotka(du,u,h,p,t)
    du[1] = 0.5*u[1] - p[1]*u[1]*u[2]
    du[2] = p[2]*h(p,t-0.5)[1]*u[2] - 0.5*u[2]
end
h = (p,t) -> [0.5, 0.5]
u0 = [0.5, 0.5]
tspan = (0., 10.)
p = [0.5, 1.]

prob = DDEProblem(f_lotka,u0,h,tspan,p,
                   constant_lags = [0.5])
sol = solve(prob, MethodOfSteps(Tsit5()))

t = collect(range(0, stop=10, length=30))
randomized = VectorOfArray([(sol(t[i]) + .03randn(2)) for i in 1:length(t)])
data = convert(Array, randomized)

using DiffEqParamEstim, NLopt

function f_lotka2(du,u,h,p,t)
    du[1] = 0.5*u[1] - p[1]*u[1]*u[2]
    du[2] = h(p,t-0.5)[1]*u[2] - 0.5*u[2]
end
h = (p,t) -> [0.5, 0.5]
u0 = [0.5, 0.5]
p = [0.5]

prob_opt = DDEProblem(f_lotka2, u0, h, tspan, p, constant_lags = [0.5])
cost_function = build_loss_objective(prob_opt, MethodOfSteps(Tsit5()),
            L2Loss(t,data),maxiter=10000,abstol=1e-8,reltol=1e-8,verbose=false)

opt = Opt(:GN_ESCH, 1)
min_objective!(opt, cost_function.cost_function2)
lower_bounds!(opt, [0.0])
upper_bounds!(opt, [1.0])
maxeval!(opt, 10000)

(minf, minx, ret) = NLopt.optimize(opt,[0.2])
@test minx[1] ≈ 0.5 atol=5e-3


text:
using Preferences
function get_provider()
    default_provider = "Polyester"

    # Load the preference
    provider = @load_preference("provider", default_provider)

    # Ensure the provider matches one of the ones we support
    if provider ∉ ("Threads", "Polyester")
        @error("Invalid Threads setting \"$(provider)\"; valid settings include [\"Threads\", \"Polyester\"], defaulting to \"Threads\"")
        provider = default_provider
    end
    return provider
end

# Read in preferences, see if any users have requested a particular backend
const threads_provider = get_provider()

"""
    set_provider!(provider; export_prefs::Bool = false)
Convenience wrapper for setting the Threads provider.  Valid values include `"Threads"`, `"Polyester"`.
Also supports `Preferences` sentinel values `nothing` and `missing`; see the docstring for
`Preferences.set_preferences!()` for more information on what these values mean.
"""
function set_provider!(provider; export_prefs::Bool = false)
    if provider !== nothing && provider !== missing && provider ∉ ("Threads", "Polyester")
        throw(ArgumentError("Invalid provider value '$(provider)'; valid settings include [\"Threads\", \"Polyester\"]"))
    end
    set_preferences!(@__MODULE__, "provider" => provider; export_prefs, force = true)
    if provider != threads_provider
        # Re-fetch to get default values in the event that `nothing` or `missing` was passed in.
        provider = get_provider()
        @info("Threads provider changed; restart Julia for this change to take effect", provider)
    end
end


text:
#--------------------------------------------------------------# plot of coefficients
@recipe function f(o::OnlineStat{XYInput})
    β = coef(o)
    nonzero = collect(β .== 0)
    x = 1:length(β)
    try
        if o.intercept  # if intercept, make indices start at 0
            x -= 1
        end
    end
    seriestype --> :scatter
    legend --> length(unique(nonzero)) > 1
    group --> nonzero
    label --> ["Nonzero" "Zero"]
    ylabel --> "Value"
    xlims --> (minimum(x) - 1, maximum(x) + 1)
    xlabel --> "Index of Coefficient Vector"
    x, β
end


text:
#Struct to hold information on wavelet-grid (v-component).
struct VWavelets{T <: Real, N <: Integer}
          nxvw :: N                                       # Number of grid points in x in UWavelets (equal to UGrid)
          nyvw :: N                                       # Number of grid points in y in UWavelets (equal to UGrid)
        mask :: Array{Bool,2}                           # Model domain on the U grid
           n :: Base.RefValue{N}                        # Number of grid points in domain
        crop :: Base.RefValue{Diagonal{T,Array{T,1}}}   # Crop matrix: diagonal matrix with mask entries on diag
        samp :: Base.RefValue{SparseMatrixCSC{T,N}}     # Sampling matrix: take full domain to model domain
      spread :: Base.RefValue{SparseMatrixCSC{T,N}}     # Spread matrix: take model domain to full domain
      levels :: N                                       # Number of wavelet levels
        idwt :: KronType{T,N}                           # Wavelet matrix cross produce
    wavelets :: Array{T,2}                              # Wavelet matrix
end

"""
    VWavelets(;
            nxvw,
            nyvw,
            mask = trues(nxvw,nyvw),
            levels)

Construct a WAVI.jl VWavelets structure with size (nxvw,nyvw)
UWavelets stores wavelet matrix and associated fields for velocities in y direction

Keyword arguments
=================
    - 'nxvw': (required) Number of grid cells in x-direction in VWavelets
    - 'nyvw': (required) Number of grid cells in y-direction in VWavelets
    - 'mask': Mask specifying the model domain with respect to UWavelet grid (VGrid)
    - levels: (required) Number of levels in the preconditioner
"""
function VWavelets(;
                    nxvw,
                    nyvw,
                    mask = trues(nxvw,nyvw),
                    levels)

    (size(mask) == (nxvw,nyvw)) || throw(DimensionMismatch("Sizes of inputs to UWavelets must all be equal to nxvw x nyvw (i.e. $nxvw x $nyvw)"))

    #compute non-inputs
    n = Ref(count(mask));  @assert n[] == count(mask)
    crop = Ref(Diagonal(float(mask[:]))); @assert crop[] == Diagonal(float(mask[:]));
    samp  = Ref(crop[][mask[:],:]); @assert samp[] == crop[][mask[:],:]
    spread = Ref(sparse(samp[]')); @assert spread[] == sparse(samp[]')
    idwt =  wavelet_matrix(nyvw,levels,"reverse" ) ⊗ wavelet_matrix(nxvw,levels,"reverse")
    wavelets = zeros(nxvw,nyvw); @assert size(wavelets)==(nxvw,nyvw)

    #make sure boolean type rather than bitarray
    mask = convert(Array{Bool,2}, mask)

    return VWavelets(
                    nxvw,
                    nyvw,
                    mask,
                    n,
                    crop,
                    samp,
                    spread,
                    levels,
                    idwt,
                    wavelets)
end


text:
using fci
using JLD2

@load "_testdata_h4_triplet_alphaconfigs.jld2"
norbs = 4
nalpha = 3
ndets = factorial(norbs)÷(factorial(nalpha)*factorial(norbs-nalpha))
configs = alpha_configs
y_matrix = ya

@testset "lookup table" begin
    index_table_test = zeros(Int, configs[1].norbs, configs[1].norbs, ndets)
    orbs = [1:configs[1].norbs;]
    for I in 1:ndets
        vir = filter!(x->!(x in configs[I].config), [1:configs[I].norbs;])
        for p in configs[I].config
            for q in vir
                new_config, sorted_config, sign_s = fci.excit_config(deepcopy(configs[I].config), [p,q])
                idx = fci.get_index(new_config, y_matrix, configs[I].norbs)
                index_table_test[p,q,configs[I].label]=sign_s*idx
            end
        end
    end

    @test isapprox(index_table_a, index_table_test, atol=0.05)

end

max_stars_repo_name: nbraunsc/fci

max_issues_repo_name: nbraunsc/fci

max_forks_repo_name: nbraunsc/fci


text:
# ---
# title: 138. Copy List with Random Pointer
# id: problem138
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table, Linked List
# link: <https://leetcode.com/problems/copy-list-with-random-pointer/description/>
# hidden: true
# ---
#
# A linked list is given such that each node contains an additional random
# pointer which could point to any node in the list or null.
#
# Return a [**deep
# copy**](https://en.wikipedia.org/wiki/Object_copying#Deep_copy) of the list.
#
# The Linked List is represented in the input/output as a list of `n` nodes.
# Each node is represented as a pair of `[val, random_index]` where:
#
#   * `val`: an integer representing `Node.val`
#   * `random_index`: the index of the node (range from `0` to `n-1`) where random pointer points to, or `null` if it does not point to any node.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/12/18/e1.png)
#
#
#
#     Input: head = [[7,null],[13,0],[11,4],[10,2],[1,0]]
#     Output: [[7,null],[13,0],[11,4],[10,2],[1,0]]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2019/12/18/e2.png)
#
#
#
#     Input: head = [[1,1],[2,1]]
#     Output: [[1,1],[2,1]]
#
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2019/12/18/e3.png)**
#
#
#
#     Input: head = [[3,null],[3,0],[3,null]]
#     Output: [[3,null],[3,0],[3,null]]
#
#
# **Example 4:**
#
#
#
#     Input: head = []
#     Output: []
#     Explanation: Given linked list is empty (null pointer), so return null.
#
#
#
#
# **Constraints:**
#
#   * `-10000 <= Node.val <= 10000`
#   * `Node.random` is null or pointing to a node in the linked list.
#   * The number of nodes will not exceed 1000.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function itensor(I::UniformScaling, is...)
  return ITensor(I, is...)
end

function ITensor(I::UniformScaling, is...)
  return ITensor(I(isqrt(dim(is))), is...)
end

# Using ITensors.jl definitions
function _ITensor(
  which_op::AbstractString, sites::Tuple, params::NamedTuple, s::Vector{<:Index}
)
  return op(which_op, s, sites; params...)
end

function _ITensor(
  which_op::Union{AbstractMatrix,UniformScaling},
  sites::Tuple,
  params::NamedTuple,
  s::Vector{<:Index},
)
  sₙ = s[collect(sites)]
  return itensor(which_op, sₙ', dag(sₙ))
end

function hassamesites(o)
  if length(o) ∈ (0, 1)
    return true
  end
  return reduce(issetequal, Ops.sites.(o))
end

function ITensor(o::Op, s::Vector{<:Index})
  return _ITensor(Tuple(o)..., s)
end

# Extend the operator `o` to the sites `n`
# by filling the rest of the sites with `op`.
function insert_ids(o, n)
  insert_n = @ignore_derivatives setdiff(n, Ops.sites(o))
  for i in insert_n
    o *= Op(I, i)
  end
  return o
end

# TODO: Merge these two.
function insert_ids(o::∏)
  n = Ops.sites(o)
  return ∏([insert_ids(oₙ, n) for oₙ in o])
end
function insert_ids(o::∑)
  n = Ops.sites(o)
  return ∑([insert_ids(oₙ, n) for oₙ in o])
end

# TODO: Does this work for fermions?
function ITensor(o::∏, s::Vector{<:Index})
  o_id = insert_ids(o)
  ∏layers = ∏([∏{ITensor}(oₙ, s) for oₙ in o_id])
  res = ITensor(Op(I, Ops.sites(o)), s)
  for layer in reverse(∏layers)
    res = layer(res)
  end
  return res
end

function ITensor(o::∑{ITensor})
  res = ITensor()
  for oₙ in o
    res += oₙ
  end
  return res
end

function ITensor(o::∑, s::Vector{<:Index})
  o_id = insert_ids(o)
  ∑layers = ∑([ITensor(oₙ, s) for oₙ in o_id])
  return ITensor(∑layers)
  ## if hassamesites(o)
  ##   return o.f([ITensor(arg, s) for arg in o])
  ## end
  ## return error("Trying to make an ITensor from operator expression $o. When making an ITensor from a sum of operators, the operators need to have the same sites.")
end

function ITensor(o::α, s::Vector{<:Index})
  return coefficient(o) * ITensor(Ops.op(o), s)
end

function ITensor(o::Applied{typeof(exp)}, s::Vector{<:Index})
  return o.f(ITensor(Ops.op(o), s))
end

itensor_adjoint(T::ITensor) = swapprime(dag(T), 0 => 1)

function ITensor(o::Applied{typeof(adjoint)}, s::Vector{<:Index})
  return itensor_adjoint(ITensor(Ops.op(o), s))
end

function ITensor(o::Applied, s::Vector{<:Index})
  return error("Trying to make ITensor from expression $(o), not yet implemented.")
end

∏{ITensor}(o::∏, s::Vector{<:Index}) = ∏([ITensor(oₙ, s) for oₙ in o])
∏{ITensor}(o::Union{Op,Applied}, s::Vector{<:Index}) = ∏{ITensor}(∏([o]), s)
∑{ITensor}(o::∑, s::Vector{<:Index}) = ∑([ITensor(oₙ, s) for oₙ in o])
∑{ITensor}(o::Union{Op,Applied}, s::Vector{<:Index}) = ∑{ITensor}(∑([o]), s)

(o::∏{ITensor})(x; kwargs...) = apply(o, x; kwargs...)
# Apply it in reverse to follow the linear algebra convention:
# (O₁ * O₂)|x⟩ = O₁ * (O₂|x⟩)
apply(o::∏{ITensor}, x; kwargs...) = apply([oₙ for oₙ in reverse(o)], x; kwargs...)


text:
using Test
using BoundaryIntegralMethod
using BoundaryIntegralMethod.DirectSweep

# Create a circular boundary
# Internal n
ni=2.0
# External n
ne=1.0
# Radius
r = 1.0
b = BoundaryIntegralMethod.Circle(r, ni, ne)

kTrue = 13.708675947990717 - 0.0006855557536092232im
k0= 13.7 -.0001im
n=100
(ks, xs) = directSweep(b, n, k0)
convergence = 1e-6
i = 1
(kf, xf) = residualInverseIteration(25, convergence, b, n, ks[i], xs[i,:])
@test abs(kTrue-kf) ≈ 0 atol=10e-10


text:
# Autogenerated wrapper script for SurfelsLib_jll for x86_64-linux-gnu-cxx11
export libsurfels

using GLEW_jll
using GLU_jll
JLLWrappers.@generate_wrapper_header("SurfelsLib")
JLLWrappers.@declare_library_product(libsurfels, "libsurfels.so")
function __init__()
    JLLWrappers.@generate_init_header(GLEW_jll, GLU_jll)
    JLLWrappers.@init_library_product(
        libsurfels,
        "lib/libsurfels.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module SMM

using NLsolve
using Distributions
using QuadGK
using CompStat
using Macros
using CustomTypes

export graph!, simulate!, fill_RegUncSS!, average, test_eq, loss!

@def reg_unc begin

    F_1 = (1-Δ)*F
    F_2 = (1+Δ)*F

    U_2 = ((1+ϵ/r)*b + η*γ*(θ_2 + ϵ*0.5*(θ_2+θ_1)/r)/(1-η))/(r+ϵ)
    U_1 = ((1+ϵ/r)*b + η*γ*(θ_1 + ϵ*0.5*(θ_2+θ_1)/r)/(1-η))/(r+ϵ)

    zf_2 = (r+δ+λ)*((r+ϵ)*U_2 - δ*0.5*ϵ*(U_2+U_1)/(r+δ))/(r+δ+ϵ) - λ*Ez/(r+δ)
    zf_1 = (r+δ+λ)*((r+ϵ)*U_1 - δ*0.5*ϵ*(U_2+U_1)/(r+δ))/(r+δ+ϵ) - λ*Ez/(r+δ)

    q_2 = m*θ_2^(-σ)
    q_1 = m*θ_1^(-σ)

    zc_1 = (r+ϵ)*U_1+λ*F_1-0.5*ϵ*(U_1+U_2)-λ*I0(zp_1)/(r+s+λ)
    zc_2 = (r+ϵ)*U_2+λ*F_2-0.5*ϵ*(U_1+U_2)-λ*I0(zp_1)/(r+s+λ)-λ*(I0(zp_2)-I0(zp_1))/(r+s+λ+ϵ/2) + 0.5*ϵ*zc_1/(r+s+λ+ϵ/2)
    zc_2 *= (r+s+λ+ϵ/2)/(r+s+λ+ϵ)

    zs_2 = (zc_2/(r+s+λ) - zf_2/(r+δ+λ))/(1/(r+s+λ) - 1/(r+δ+λ))
    zs_1 = (zc_1/(r+s+λ) - zf_1/(r+δ+λ))/(1/(r+s+λ) - 1/(r+δ+λ))

end

"""
    graph!(t, sim_np, sim_nf, sim_Wp, shocks, states, ss, param, D; k=10)

"""
function graph!(t, sim_np, sim_nf, sim_Wp, shocks, states, ss, param, D; k=10)
    @load_struct param η r s λ α m b γ σ F δ ϵ
    @spec_funs

    #Starting from one of the steady-states
    T = shocks[1]
    T_next = shocks[2]
    i0 = states[1]
    @load_index_cat ss i0 old zp zs zc npᶜ np⁰ np nf Wp Wpᶜ Wp⁰ μᵖ μᶠ ξ q

	n = size(shocks,1)
    push!(t,T)
    push!(sim_np, np_old)
    push!(sim_nf, nf_old)
	push!(sim_Wp, Wp_old)

    if n > 2
        for i in 2:n-1
            T = shocks[i]
            T_next = shocks[i+1]
            sh = states[i]

            # println((T, T_next, sh, zp_old, zs_old, zc_old, npᶜ_old, np⁰_old, nf_old, Wpᶜ_old, Wp⁰_old))

            @load_index ss sh θ zp zc zs ξ μᵖ μᶠ npᶜ np⁰ np nf u Wpᶜ Wp⁰ Wp F
            F_bar = 0.5*sum(ss.F)
            @dual_to_dual
            @dual_to_dual_wages

            #Updating old values
            zp_old = zp
            zs_old = zs
            zc_old = zc
            npᶜ_old = np_c(T_next-T)
            np⁰_old = np_0(T_next-T)
            nf_old = nᶠ(T_next-T)
            Wpᶜ_old = Wp_c(T_next-T)
            Wp⁰_old = Wp_0(T_next-T)

			time_interval = range(T,T_next,length=k)
			append!(t,time_interval)
			append!(sim_np,nᵖ.(time_interval .- T))
			append!(sim_nf,nᶠ.(time_interval .- T))
			append!(sim_Wp, Wᵖ.(time_interval .- T))

        end
    end
end

"""

    Eq!(R, x, Δ, param)

Fills vector `R` with the residuals of the system characterizing the model with regulatory uncertainty given variables `x`, firing cost wedge `Δ` and parameters `param`.
"""
function Eq!(R, x, Δ, param)

	@load_struct param η r s λ α m b γ σ F δ ϵ

    @spec_funs

    θ_2 = exp(x[1])
    θ_1 = exp(x[2])
    zp_2 = x[3]
    zp_1 = x[4]

    @reg_unc

    R[1] = zp_1-(r+ϵ)*U_1+(r+s+ϵ)*F_1+0.5*ϵ*(U_1+U_2-F_1-F_2)+λ*I0(zp_1)/(r+s+λ)+0.5*ϵ*(zp_1-zp_2)/(r+s+λ+ϵ/2)
    R[2] = zp_2-(r+ϵ)*U_2+(r+s+ϵ)*F_2+0.5*ϵ*(U_1+U_2-F_1-F_2)+λ*I0(zp_1)/(r+s+λ)+λ*(I0(zp_2)-I0(zp_1))/(r+s+λ+ϵ/2)
    R[3] = γ - q_2*(1-η)*( I0(max(zs_2,zc_2))/(r+s+λ) + (I0(zf_2)-I0(max(zs_2,zf_2)))/(r+δ+λ) )
    R[4] = γ - q_1*(1-η)*( I0(max(zs_1,zc_1))/(r+s+λ) + (I0(zf_1)-I0(max(zs_1,zf_1)))/(r+δ+λ) )

end

"""
    fill_RegUncSS!(ss, Δ, param, x0, show_trace)

Fills the `RegUncSS` structure `ss` with the proper steady-state values for both values of steady-state associated with `Δ`.
"""
function fill_RegUncSS!(ss, Δ, param, x0, show_trace)
    R = zeros(4)
    sol = nlsolve((R,x)->Eq!(R,x,Δ,param), x0, show_trace=show_trace)
    θ_2 = exp(sol.zero[1])
    θ_1 = exp(sol.zero[2])
    zp_2 = sol.zero[3]
    zp_1 = sol.zero[4]
    @load_struct param η r s λ α m b γ σ F δ ϵ

    @spec_funs

    @reg_unc

    ss.θ = [θ_1, θ_2]
    ss.zp = [zp_1, zp_2]
    ss.zc = [zc_1, zc_2]
    ss.zf = [zf_1, zf_2]
    ss.zs = [zs_1, zs_2]
    ss.F = [F_1, F_2]

    @load_struct ss θ zp zc zf zs F
    F_bar = 0.5*sum(F)

    p = @. m*θ^(1-σ)
    q = @. p/θ
    ξ = @. s+λ*G(zp)
    μᵖ = @. p*(1-G(max(zs,zc)))
    μᶠ = @. p*(G(max(zs,zf))-G(zf))
    npᶜ = @. λ*(1-G(zp))*μᵖ*δ/((s+λ)*(μᵖ*δ + ξ*δ + μᶠ*ξ))
    np⁰ = @. ξ*μᵖ*δ/((s+λ)*(μᵖ*δ + ξ*δ + μᶠ*ξ))
    np = @. npᶜ + np⁰
    nf = @. μᶠ*ξ/(μᵖ*δ + ξ*δ + μᶠ*ξ)
    u = @. 1-np-nf
    v = @. u*θ
    z_bar_p = @. (λ*np*I(zp)+p*u*I(max(zs,zc)))/((s+λ)*np)
    z_bar_f = @. (λ*nf*Ez+p*u*(I(zf)-I(max(zs,zf))))/((δ+λ)*nf)
    z_bar = @. (np*z_bar_p+nf*z_bar_f)/(np+nf)
    W = @. z_bar_p*np + z_bar_f*nf + b*u - γ*v - λ*G(zp)*np*F
    Eᶜ = @. η*(I(zp)/(1-G(zp)) + (r+s+ϵ)*F - ϵ*F_bar + γ*θ) + (1-η)*b
    E⁰ = @. η*(I(max(zs,zc))/(1-G(max(zs,zc))) - λ*F + γ*θ) + (1-η)*b
    Wpᶜ = @. λ*(1-G(zp))*Eᶜ*np/(s+λ)
    Wp⁰ = @. μᵖ*E⁰*u/(s+λ)
    Wp = @. Wpᶜ + Wp⁰

    ss.flag = ss.flag && converged(sol) && !(any(isnan.(Wp⁰))) && !(any(zs .< zc[2])) && !(any(np⁰ .< 0.01)) && !(any(npᶜ .< 0.01)) && !(any(nf .< 0.01))

    @save_struct ss zp zc zs zf npᶜ np⁰ np nf θ ξ μᵖ μᶠ z_bar_p z_bar_f z_bar W u Wpᶜ Wp⁰ Wp q

    return
end

macro test_regunc_eq(a,regunc,dual)
	expr = Expr(:block)
	for name in intersect(fieldnames(RegUncSS),fieldnames(DualSS))
		push!(expr.args,:(push!($a,$regunc.$name .- $dual.$name)))
	end
	esc(expr)
end

"""

	test_eq(param)

Verifies for parameters `param` that steady-state values of the dual model without uncertainty and of the regulatory uncertainty model with Δ = 0 and ϵ = 0 coincide.
"""
function test_eq(param)
	param.ϵ = 0.
    Δ = 0.
	regunc = RegUncSS()
	init!(regunc, 2)
	fill_RegUncSS!(regunc, Δ, param, [log(2.2), log(2.2), 0.93, 0.93], false)
	dual = DualSS()
	fill_DualSS!(dual, param, [0.9, 1.3])
	t = []
	@test_regunc_eq t regunc dual
	@assert maximum([(t...)...]) < 1e-7
    return
end

"""

    simulate!(shocks, states, D)

Simulates a path of shocks on firing costs of duration `D` and stores the resulting shocks and states in `shocks` and `states`.
"""
function simulate!(shocks, states, param, D)
    @load_struct param ϵ
    T = 0.
    c = 1
    u = rand()
    s = (u<0.5)*1+(u>=0.5)*2
    while (T < D)
        push!(shocks,T)
        T += -log(rand())/(ϵ/2)
        c += 1
        push!(states,s)
        s = 3-s
    end
    push!(states,states[end])
    push!(shocks,D)
end

"""

    average(shocks, states, ss, param, D)

Computes the empirical means of targets given a shock instants `shocks`, states `states`, steady-state values `ss`, parameters `param` and duration `D`.
"""
function average(shocks, states, ss, param, D)
    @load_struct param η r s λ α m b γ σ F δ ϵ
    @spec_funs

    #Starting from one of the steady-states
    T = shocks[1]
    T_next = shocks[2]
    i0 = states[1]
    @load_index_cat ss i0 old zp zs zc npᶜ np⁰ np nf Wpᶜ Wp⁰ Wp μᵖ μᶠ ξ q

    np_m = np_old*(T_next-T)
    nf_m = nf_old*(T_next-T)
	Fswp_m = F*(T_next-T)/(Wp_old/np_old)
	μᶠsμ_m = μᶠ_old*(T_next-T)/(μᵖ_old+μᶠ_old)
	ssξ_m = s*(T_next-T)/ξ_old
	q_m = q_old*(T_next-T)

    n = size(shocks,1)
    if n > 2
        for i in 2:n-1
            T = shocks[i]
            T_next = shocks[i+1]
            sh = states[i]

            # println((T, T_next, sh, zp_old, zs_old, zc_old, npᶜ_old, np⁰_old, nf_old, Wpᶜ_old, Wp⁰_old))

            @load_index ss sh θ zp zc zs ξ μᵖ μᶠ npᶜ np⁰ np nf u Wpᶜ Wp⁰ Wp F
            F_bar = 0.5*sum(ss.F)
            @dual_to_dual
            @dual_to_dual_wages

            #Means on [T,T_next]
            np_m += (T_next-T)*np + ρ₁ᵖ*(exp(α₁*(T_next-T))-1)/α₁ + ρ₂ᵖ*(exp(α₂*(T_next-T))-1)/α₂
            nf_m += (T_next-T)*nf + ρ₁ᶠ*(exp(α₁*(T_next-T))-1)/α₁ + ρ₂ᶠ*(exp(α₂*(T_next-T))-1)/α₂
            ssξ_m += (T_next-T)*s/ξ
            μᶠsμ_m += (T_next-T)*μᶠ/(μᶠ+μᵖ)
            q_m += (T_next-T)*m*θ^(-σ)

            Δ_wpm1,_ = quadgk(t->F/(Wᵖ(t)/nᵖ(t)), 0, T_next-T)
            Fswp_m += Δ_wpm1

            #Updating old values
            zp_old = zp
            zs_old = zs
            zc_old = zc
            npᶜ_old = np_c(T_next-T)
            np⁰_old = np_0(T_next-T)
            nf_old = nᶠ(T_next-T)
            Wpᶜ_old = Wp_c(T_next-T)
            Wp⁰_old = Wp_0(T_next-T)

        end
    end
    return [Fswp_m, μᶠsμ_m, np_m, nf_m, ssξ_m, q_m] ./ D

end

function moments!(ss, shocks, states, Δ, param, D, show_trace; x0=[log(2.2), log(2.2), 0.93, 0.93])
    mom = Inf*ones(6)
    fill_RegUncSS!(ss, Δ, param, x0, show_trace)
    if ss.flag
        mom .= average(shocks, states, ss, param, D)
    end

	return mom

end

function loss!(x, param, ss, shocks, states, Δ, D, targets; show_trace=false)
    @load_vec_into_struct x param s α m b γ F
    ss.flag = !any(x .< 0.)
    err = Inf
	if ss.flag
        m = moments!(ss, shocks, states, Δ, param, D, show_trace)
		e = (m' - targets)./targets
		# err = (e*W*e')[1,1]
		err = maximum(abs.(e))
	end
    return err
end

# function EqDual!(R,x,param)

# 	@load_vec param s α m b γ F

# 	Ez = exp(α^2/2)
# 	G(x) = cdf(LogNormal(0,α),x)
# 	g(x) = pdf(LogNormal(0,α),x)
# 	I(x) = exp(α^2/2)*cdf(Normal(0,1),((α^2-log(max(x,0)))/α))
# 	I0(x) = I(x) - x*(1-G(x))

# 	zp = x[1]
# 	θ = x[2]

# 	zc = zp+(r+s+λ)*F
# 	q = m*θ^(-σ)
# 	rU = b + η*γ*θ/(1-η)
# 	zf = rU - λ*(Ez-rU)/(r+δ)
# 	zs = (zp/(r+s+λ) + F - zf/(r+δ+λ))/(1/(r+s+λ) - 1/(r+δ+λ))

# 	R[1] = zp+(r+s)*F+λ*I0(zp)/(r+s+λ)-rU
# 	R[2] = γ - q*(1-η)*( I0(max(zs,zc))/(r+s+λ) + (I0(zf)-I0(max(zs,zf)))/(r+δ+λ) )
# end

# @everywhere function compute_n_dual(sol)

# 	@load_vec param s α m b γ F

# 	Ez = exp(α^2/2)
# 	G(x) = cdf(LogNormal(0,α),x)
# 	g(x) = pdf(LogNormal(0,α),x)
# 	I(x) = exp(α^2/2)*cdf(Normal(0,1),((α^2-log(max(x,0)))/α))
# 	I0(x) = I(x) - x*(1-G(x))

# 	zp = sol[1]
# 	θ = sol[2]

# 	zc = zp+(r+s+λ)*F
# 	q = m*θ^(-σ)
# 	rU = b + η*γ*θ/(1-η)
# 	zf = rU - λ*(Ez-rU)/(r+δ)
# 	zs = (zp/(r+s+λ) + F - zf/(r+δ+λ))/(1/(r+s+λ) - 1/(r+δ+λ))
# 	ξ = s+λ*G(zp)
# 	p = m*θ^(1-σ)
# 	μᵖ = p*(1-G(max(zs,zc)))
# 	μᶠ = p*(G(max(zs,zf))-G(zf))
# 	np = μᵖ*δ/(μᵖ*δ + ξ*δ + μᶠ*ξ)
# 	nf = μᶠ*ξ/(μᵖ*δ + ξ*δ + μᶠ*ξ)

# 	return (np,nf)

# end

end


text:
struct MeasureNout{T} <: NaiveNASflux.AbstractMutableComp
    wrapped::T
    actdim::Int
    outsize::Ref{Int}
end
MeasureNout(l; actdim=actdim(l), outsize=0) = MeasureNout(l, actdim, Ref(outsize))
NaiveNASflux.wrapped(mn::MeasureNout) = mn.wrapped
NaiveNASflux.layer(mn::MeasureNout) = layer(mn.wrapped)

NaiveNASlib.nout(mn::MeasureNout) = mn.outsize[]
NaiveNASlib.nin(mn::MeasureNout) = nin(mn.wrapped)

function NaiveNASlib.Δsize!(mn::MeasureNout, ins::AbstractVector, outs::AbstractVector; kwargs...)
    mn.outsize[] = length(outs)
    Δsize!(mn.wrapped, ins, outs; kwargs...)
end

function (mn::MeasureNout)(x...)
    out = mn.wrapped(x...)
    measuresize!(mn, out)
    return out
end
measuresize!(mn::MeasureNout, x::AbstractArray) = mn.outsize[] = size(x, mn.actdim)
function measuresize!(::MeasureNout, x) end # Need methods for AbstractProbe?

calc_outsize(m::NaiveNASflux.AbstractMutableComp, v) = calc_outsize(wrapped(m), v)
function calc_outsize(mn::MeasureNout, v)
    if mn.outsize[] == 0
        mn.outsize[] = calc_outsize(mn.wrapped, v)
    end
    return mn.outsize[]
end

function try_infer_sizes!(g, insizes...)
    all(v -> nout(v) > 0, vertices(g)) && return
    try
        insizes_nobatch = map(zip(insizes, layertype.(inputs(g)))) do (insize, lt)
            # not sure this can happen...
            lt isa FluxRecurrent && length(insize) == 3 && return insize[1:end-2]
            insize[1:end-1]
        end
        if length(insizes) === length(inputs(g)) && all(inshape -> !isempty(inshape) && all(s -> s isa Number && s > 0, inshape), insizes_nobatch)
            # This will make any MeasureNout to become aware of the size
            Flux.outputsize(g,insizes_nobatch...; padbatch=true)
            update_size_meta!(g)
        else
            @warn  "No valid input sizes provided. Shape inference could not be done. Either provide Integer insizes manually or use load(...; infer_shapes=false) to disable. If disabled, graph mutation might not work."
        end
        # If insizes are not provided there was an old method called fix_zerosizes! which might be refactored to work with NaiveNASlib 2.0 and MeasureNouts
        # it was not 100% to succeed either and required a CompGraphBuilder
    catch e
        throw(CompositeException([ErrorException("Size inference failed! Use load(...; infer_shapes=false) to disable! If disabled, graph mutation might not work."), e]))
    end
end

update_size_meta!(g::CompGraph) = foreach(update_size_meta!, vertices(g))
update_size_meta!(v::AbstractVertex) = update_size_meta!(base(v))
update_size_meta!(::InputVertex) = nothing
update_size_meta!(v::CompVertex) = update_size_meta!(v.computation)
update_size_meta!(m::NaiveNASflux.AbstractMutableComp, args...) = update_size_meta!(wrapped(m), args...)

function update_size_meta!(m::ActivationContribution, args...)
    outsize = update_size_meta!(wrapped(m), args...)
    if (outsize !== missing) && (m.contribution === missing || length(m.contribution) != outsize)
        m.contribution = zeros(eltype(m.contribution), outsize)
    end
    return outsize
end

function update_size_meta!(m::LazyMutable, args...)
    outsize = update_size_meta!(wrapped(m), args...)
    if (outsize !== missing) && (m.outputs === missing || length(m.outputs) != outsize)
        m.outputs = 1:outsize
    end
    return outsize
end

update_size_meta!(m::MeasureNout, args...) = update_size_meta!(wrapped(m), nout(m))
update_size_meta!(f, outsize::Integer) = outsize
update_size_meta!(f) = missing


text:
using Distributions, HypothesisTests

p = [0.08, 0.12, 0.2, 0.2, 0.15, 0.25]
O = [3, 2, 9, 11, 8, 27]
M = length(O)
n = sum(O)
E = n*p

testStatistic = sum((O-E).^2 ./E)
pVal = ccdf(Chisq(M-1), testStatistic)

println("Manually calculated test statistic: ", testStatistic)
println("Manually calculated p-value: ", pVal,"\n")

println(ChisqTest(O,p))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
include("icp_trajectory_gen.jl")
include("foot_trajectories.jl")
include("notebooks.jl")


text:
"""
    find(registry, [; col = filter, ...])

Find entries in a registry that match a set of filters given as keyword
arguments. The keyword name must correspond to a field in `registry`,
while the value is a filter that is checked against each field value.

See [how to use registries](/documents/docs/using.md) for examples.
"""
function find(registry::Registry; filters...)
    mask = [matchesentry(getfield(registry, :fields), row; filters...)
            for row in getfield(registry, :data)]
    return view(registry, mask)
end

Base.filter(registry::Registry; kwargs...) = find(registry; kwargs...)

function Base.findfirst(registry::Registry; filters...)
    for row in getfield(registry, :data)
        if matchesentry(getfield(registry, :fields), row; filters...)
            return RegistryEntry(row, registry)
        end
    end
    nothing
end

function matchesentry(fields, row; filters...)
    return all(matchesentry(fields[col], row[col], f) for (col, f) in filters)
end

# The default for any function is to match against a value
matchesentry(f::Field, value, query) = isnothing(f.filterfn) ? value == query : f.filterfn(value, query)
matchesentry(field::Field, value::Missing, ::typeof(ismissing)) = true
matchesentry(field::Field, value::Missing, _) = false

# Passing a function as a query uses it as predicate
matchesentry(::Field, value, query::Function) = query(value)
matchesentry(::Field{<:Function}, value::Function, query::Function) = query(value)

# For strings, do fuzzy matching and allow regex
matchesentry(::Field, str::String, query::String) = occursin(query, str)
matchesentry(::Field, str::String, query::Regex) = match(query, str) !== nothing

# TODO: efficient `Base.findfirst`

@testset "find" begin
    function testregistry()
        registry = Registry((
                id = Field(String, name="ID"),
                name = Field(String, name="name"),
            ),
        name = "Test registry")
        push!(registry, (id = "id1", name = "name1"))
        push!(registry, (id = "id2", name = "name2"))
        return registry
    end
    @test_nowarn testregistry()
    registry = testregistry()

    @test length(find(registry, id = "id1")) == 1
    @test length(find(registry, id = "1")) == 1
    @test length(find(registry, id = r"1")) == 1
    @test length(find(registry, id = "id")) == 2
    @test length(find(registry, id = "id", name = "name")) == 2
    @test length(find(registry, id = x -> startswith(x, "id"), name = "name")) == 2

    @test NamedTuple(findfirst(registry, id = "id")) == (id = "id1", name = "name1")
    @test isnothing(findfirst(registry, id = "404"))

    @testset "filtering" begin
        @testset "Fallback: equality" begin
            @test matchesentry(Field(Symbol), :x, :x)
            @test !matchesentry(Field(Symbol), :x, :y)
        end

        @testset "Predicates" begin
            @test matchesentry(Field(Symbol), :x, ==(:x))
            @test matchesentry(Field(Function), sum, ==(sum))
        end

        @testset "Missing values" begin
            @test !matchesentry(Field(Symbol, optional = true), missing, :x)
            @test matchesentry(Field(Symbol, optional = true), missing, ismissing)
        end

        @testset "String" begin
            @test matchesentry(Field(String), "hello", "hel")
            @test !matchesentry(Field(String), "hello", "yo")
            @test matchesentry(Field(String), "hello", r"^h")
            @test !matchesentry(Field(String), "hello", r"^h$l")
        end

        @testset "Row" begin
            fields = (f1 = Field(String), f2 = Field(Symbol, optional = true))
            @test matchesentry(fields, (f1 = "", f2 = :x))
            @test matchesentry(fields, (f1 = "", f2 = :x), f2 = :x, f1 = "")
            @test !matchesentry(fields, (f1 = "", f2 = :s), f2 = :x, f1 = "")
            @test !matchesentry(fields, (f1 = "", f2 = :x), f2 = :x, f1 = "s")
        end
    end

end


text:
#=

This is a script that we use to run speed tests on Chimera graphs.
Run like julia st_chimera.jl 123456 10
Which will run for approximately 10 hours on graphs with 123456 vertices

=#

n = parse(ARGS[1])
hours = parse(ARGS[2])

fn = "st_chimera_$(n)h$(hours).jld"
println(fn)

using Laplacians
using MATLAB

lapdir = Pkg.dir("Laplacians")

include("$(lapdir)/compare/matlabSafe.jl")
include("$(lapdir)/compare/compare_solvers_TL.jl")

ac_deg = function(a; verbose=false, args...)
    approxchol_lap(a; params=ApproxCholParams(:deg), verbose=verbose, args...)
end
ac_wdeg = function(a; verbose=false, args...)
    approxchol_lap(a; params=ApproxCholParams(:wdeg), verbose=verbose, args...)
end

test_ac = SolverTest(ac_deg, "ac")
test_acw = SolverTest(ac_wdeg, "ac_fast")
test_amg = SolverTest(AMGLapSolver, "pyamg")

# removed chol because killing it can cause a crash
#test3 = SolverTest(chol_lap, "chol")
#test4 = SolverTest(cgLapSolver, "cg")
tests = [test_ac test_acw test_amg]

# warm up the code
a = chimera(10000,1)
b = randn(10000)
b = b - mean(b)
for solver in tests
    f = solver.solver(a,verbose=true)
    x = f(b)
end

dic = Dict()

using JLD

i = 0
t0 = time()

while time() - t0 < 60*60*hours

    i += 1
    println("-----------")
    println(i)

    a = chimera(n,i)
    aw = randWeight(a)

    tn = "chimera($n,$i)"
    b = randn(n)
    b = b - mean(b)
    b = b / norm(b)
    x = testVMatlabLap(tests, dic, a, b, verbose=true, tol=1e-8, testName=tn, test_icc=false)

    tn = "wtedChimera($n,$i)"
    x = testVMatlabLap(tests, dic, aw, b, verbose=true, tol=1e-8, testName=tn, test_icc=false)

    save(fn,"dic",dic)

end


text:
using SweepContractor

n=8; d=2

TN=TensorNetwork(undef,n)
for i=1:n
    TN[i]=Tensor(
        setdiff(1:n,i),
        randn(d*ones(Int,n-1)...),
        randn(),
        randn()
    )
end

println("example_rand")
for χ∈d.^(4:9)
    sweep=sweep_contract(TN,χ,2*χ;connected=false)
    println("χ=$χ:\t",ldexp(sweep...))
end


text:
using Dates: DateTime
using Dates: DatePeriod, Day, Hour, Millisecond, Minute, Month, Second, Week, Year
using LaxZonedDateTimes
using LaxZonedDateTimes: isrepresentable, NonExistent
using Test
using TimeZones
using TimeZones: Transition, timezone

const winnipeg = TimeZone("America/Winnipeg")

@testset "LaxZonedDateTimes" begin
    include("accessors.jl")
    include("conversions.jl")
    include("intervals.jl")
    include("laxzoneddatetimes.jl")
    include("ranges.jl")
    include("rounding.jl")
end


text:
import IntervalArithmetic: ∩, ∪, entireinterval, ∈

export ∩, ∪, entireinterval, ∈

for f in (:∩, :∪, :∈)
    @eval $f(a::NumberInterval, b::NumberInterval) =
        NumberInterval($f(Interval(a), Interval(b)))
    @eval $f(a::Interval{T}, b::NumberInterval{T}) where T =
        NumberInterval($f(a, Interval(b)))
    @eval $f(a::NumberInterval, b::Interval) = $f(b, a)
end
∈(a::Number, b::NumberInterval) = ∈(a, Interval(b))

for f in (:entireinterval, )
    @eval $f(a::NumberInterval) = NumberInterval($f(Interval(a)))
end


text:
#This implementation is very similar to original implementation in https://github.com/stanfordnlp/mac-network
using JSON,Knet,Images,HDF5,Plots
#import AutoGrad: cat1d
using Printf,Random
include("model.jl")
savemodel(filename,m,mrun,o) = Knet.save(filename,"m",m,"mrun",mrun,"o",o)

function loadmodel(filename;onlywrun=false)
    d = Knet.load(filename)
    if onlywrun
        mrun=d["mrun"];o=d["o"]
        m=nothing
    else
        m=d["m"];o=d["o"];mrun=d["mrun"];
    end
    return m,mrun,o;
end

function getQdata(dhome,set)
    JSON.parsefile(dhome*set*".json")
end

function invert(vocab)
       int2tok = Array{String}(undef,length(vocab))
       for (k,v) in vocab; int2tok[v] = k; end
       return int2tok
end

function getDicts(dhome,dicfile)
    dic  = JSON.parsefile(dhome*dicfile*".json")
    qvoc = dic["word_dic"]
    avoc = dic["answer_dic"]
    i2w  = invert(qvoc)
    i2a  = invert(avoc)
    return qvoc,avoc,i2w,i2a
end

function loadFeatures(dhome,set;h5=false)
    if h5
        return h5open(dhome*set*".hdf5","r")["data"]
    else
        feats = reinterpret(Float32,read(open(dhome*set*".bin")))
        return reshape(feats,(14,14,1024,div(length(feats),200704)))
    end
end

function miniBatch(data;shfl=true,srtd=false,B=32)
    L = length(data)
    shfl && shuffle!(data)
    srtd && sort!(data;by=x->length(x[2]))
    batchs = [];
    for i=1:B:L
        b         = min(L-i+1,B)
        questions = Any[]
        answers   = zeros(Int,b)
        images    = Any[]
        families  = zeros(Int,b)

        for j=1:b
            crw = data[i+j-1]
            push!(questions,reverse(Array{Int}(crw[2])))
            push!(images,parse(Int,crw[1][end-9:end-4])+1)
            answers[j]  = crw[3]
            families[j] = crw[4]
        end

        lngths     = length.(questions);
        srtindices = sortperm(lngths;rev=true)

        lngths     = lngths[srtindices]
        Tmax       = lngths[1]
        questions  = questions[srtindices]
        answers    = answers[srtindices]
        images     = images[srtindices]
        families   = families[srtindices]

        qs = Int[];
        batchSizes = Int[];
        pads = falses(b,Tmax)

        for k=1:b
           pads[k,lngths[k]+1:Tmax].=true
        end

        if sum(pads)==0
           pads=nothing
        end

        while true
            batch = 0
            for j=1:b
                if length(questions[j]) > 0
                    batch += 1
                    push!(qs,pop!(questions[j]))
                end
            end
            if batch != 0
                push!(batchSizes,batch)
            else
                break;
            end
        end
        push!(batchs,(images,qs,answers,batchSizes,pads,families))
    end
    return batchs
end

function loadTrainingData(dhome="data/";h5=false)
    !h5 && println("Loading pretrained features for train&val sets.
                It requires minimum 70GB RAM!!!")
    trnfeats = loadFeatures(dhome,"train";h5=h5)
    valfeats = loadFeatures(dhome,"val";h5=h5)
    println("Loading questions ...")
    trnqstns = getQdata(dhome,"train")
    valqstns = getQdata(dhome,"val")
    println("Loading dictionaries ... ")
    qvoc,avoc,i2w,i2a = getDicts(dhome,"dic")
    return (trnfeats,valfeats),(trnqstns,valqstns),(qvoc,avoc,i2w,i2a)
end

function loadDemoData(dhome="data/demo/")
    println("Loading demo features ...")
    feats = loadFeatures(dhome,"demo")
    println("Loading demo questions ...")
    qstns = getQdata(dhome,"demo")
    println("Loading dictionaries ...")
    dics = getDicts(dhome,"dic")
    return feats,qstns,dics
end

function modelrun(M,data,feats,o,Mrun=nothing;train=false)
    getter(id) = view(feats,:,:,:,id)
    cnt=total=0.0; L=length(data);
    Mparams   = params(M)
    Rparams   = Mrun !== nothing ? params(Mrun) : nothing
    # results   = similar(Array{Float32},200704*48) #uncomment for INPLACE
    println("Timer Starts");
    for i in 1:L
        ids,questions,answers,batchSizes,pad,families = data[i]
        B    = batchSizes[1]
        xB   = atype(ones(Float32,1,B))
        #x    = inplace_batcher(results,feats,ids) #uncomment for INPLACE
        x    = batcher1(feats,ids) #comment for INPLACE
        xS   = atype(x)
        #xS   = atype(reshape(cat1d(map(getter,ids)...),14,14,1024,B))
        xP   = pad==nothing ? nothing : atype(pad*Float32(1e22))
        if train
            J = @diff M(questions,batchSizes,xS,xB,xP;answers=answers,p=o[:p],selfattn=o[:selfattn],gating=o[:gating])
            cnt += value(J)*B; total += B;
            for w in Mparams
                update!(w.value,grad(J,w),w.opt)
            end
            if Mrun != nothing
                for (wr,wi) in zip(Rparams,Mparams);
                    axpy!(1.0f0-o[:ema],wi.value-wr.value,wr.value);
                end
            end
        else
            preds  = M(questions,batchSizes,xS,xB,xP;p=o[:p],selfattn=o[:selfattn],gating=o[:gating])
            cnt   += sum(preds.==answers)
            total += B
        end
        i % 1000 == 0 && println(@sprintf("%.2f Accuracy|Loss", train ? cnt/total : 100cnt/total))
    end
    train && savemodel(o[:prefix]*".jld2",M,Mrun,o);
end

function train!(M,Mrun,sets,feats,o)
    @info "Training Starts...."
    setoptim!(M,o)
    for i=1:o[:epochs]
        println("Epoch $(i) starts...")
        modelrun(M,sets[1],feats[1],o,Mrun;train=true)
        modelrun(Mrun,sets[2],feats[2],o;train=false)
    end
    return M,Mrun;
end

function train(sets,feats,o)
     if o[:mfile]==nothing
         M    = MACNetwork(o);
         Mrun = deepcopy(M)
     else
         M,Mrun,o = loadmodel(o[:mfile])
     end
     train!(M,Mrun,sets,feats,o)
     return M,Mrun;
end

function train(dhome="data/",o=nothing)
     if o==nothing
         o=Dict(:h5=>false,:mfile=>nothing,:epochs=>10,
                :lr=>0.0001,:p=>12,:ema=>0.999f0,:batchsize=>32,
                :selfattn=>false,:gating=>false,:d=>512,
                :shuffle=>true,:sorted=>false,:prefix=>string(now())[1:10],
                :vocab_size=>90,:embed_size=>300, :dhome=>"data/", :loadresnet=>false)
     end
     feats,qdata,dics = loadTrainingData(dhome;h5=o[:h5])
     sets = []
     for q in qdata; push!(sets,miniBatch(q;shfle=o[:shuffle],srtd=o[:sorted])); end
     qdata = nothing; #gc();
     M,Mrun = train(sets,feats,o)
     return M,Mrun,sets,feats,dics;
end

function validate(Mrun,valset,valfeats,o)
     modelrun(Mrun,valset,valfeats,o;train=false)
end

function validate(mfile,valset,valfeats,o)
     _,Mrun,_ = loadmodel(mfile)
     modelrun(Mrun,valset,valfeats;train=false)
     return Mrun
end

function validate(mfile,dhome,o)
     _,Mrun,_,o   = loadmodel(mfile)
     valfeats     = loadFeatures(dhome,"val")
     qdata        = getQdata(dhome,"val")
     dics         = getDicts(dhome,"dic")
     valset       = miniBatch(qdata;shfle=o[:shuffle],srtd=o[:sorted])
     modelrun(Mrun,valset,valfeats,o;train=false)
     return Mrun,valset,valfeats
end

function singlerun(Mrun,feat,question;p=12,selfattn=false,gating=false)
    results        = Dict{String,Any}("cnt"=>1)
    batchSizes     = ones(Int,length(question))
    xB             = atype(ones(Float32,1,1))
    outputs = Mrun(question,batchSizes,feat,xB,nothing;tap=results,p=p,selfattn=selfattn,gating=gating,allsteps=true)
    prediction = argmax(results["y"])
    return results,prediction,outputs
end

function visualize(img,results;p=12)
    s_y,s_x = size(img)./14
    for k=1:p
        α = results["w_attn_$(k)"][:]
        wrds    = i2w[question]
        p = bar(α;xticks=(collect(1:length(wrds)),wrds),xrotation=90,bar_width = 1,
            xtickfont = font(8, "Courier"),yticks=0:.1:(maximum(α)+.1),
            legend=false,size=(600,100+400*(maximum(α))),aspect_ratio=10)
        savefig(p,"plots/$(k).png")
        display(RGB.(load("plots/$(k).png")))
        println("Image Attention Map: ")
        flush(stdout)
        hsvimg = HSV.(img);
        attn = results["KB_attn_$(k)"]
        for i=1:14,j=1:14
            rngy          = floor(Int,(i-1)*s_y+1):floor(Int,min(i*s_y,320))
            rngx          = floor(Int,(j-1)*s_x+1):floor(Int,min(j*s_x,480))
            hsvimg[rngy,rngx]  = scalepixel.(hsvimg[rngy,rngx],attn[LinearIndices((1:14,1:14))[i,j]])
        end
        display(hsvimg)
    end
end

function scalepixel(pixel,scaler)
     return HSV(pixel.h,pixel.s,min(1.0,pixel.v+5*scaler))
end

max_stars_repo_path: Courses/World of machine learning/visualize/src/main.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
module MPISimulations

using Pkg.TOML, StaticArrays, LinearAlgebra, FastGaussQuadrature, DocStringExtensions

include("TOML.jl")

include("Quadrature.jl")
include("TracerDistributions.jl")
include("MagneticFields.jl")
include("MagneticMoments.jl")
include("Magnetizations.jl")
include("Receivers.jl")
include("ExcitationFields.jl")

end


text:
function reflection_matrix(delta, mx, my, MinvT, bar::Barrier, grd::Grid)

    # define function we want to map onto grd
    function f(x, y, grd)
        x, y = reflection(x, y, delta, bar)
        row_dat, row_js = interpolation_row(x, y, mx, my, MinvT, grd)
        return hcat(row_dat, row_js)
    end

    M = grid_map(f, (mx*my,2), grd)
    dat = M[:,1]
    js = Int64[M[:,2]...]
    is = vcat([k*ones(Int64, mx*my) for k=1:grd.Nk]...)

    R2d = sparse(is, js, dat, grd.Nk, grd._Nx*grd._Ny) * transpose(grd.Proj)
    R3d = kron(sparse(I, grd.Nz, grd.Nz), R2d)

    return R3d
end

function boundary_operators(mx, my, MinvT, deltas::Vector{Float64}, bars::Vector{Barrier}, qs::Vector{Float64}, grd::Grid)

    Nb = length(bars)

    PEN = SparseMatrixCSC[sparse(I, grd.Nz*grd.Nk, grd.Nz*grd.Nk)]
    REF = SparseMatrixCSC[]
    DCHLT = SparseMatrixCSC[]
    NMANN = SparseMatrixCSC[]

    for i=1:Nb
        P = Diagonal(f_to_grid((x,y)-> smoothstep(x, y, deltas[i], bars[i]), grd))
        for j=1:length(PEN)
            PEN[j] = P * PEN[j]
        end
        append!(PEN, [(I - P)/qs[i]])
    end

    for i=1:Nb
        R = reflection_matrix(deltas[i], mx, my, MinvT, bars[i], grd)
        @info "Calculated reflection matrix $i"
        append!(REF, [R])
        append!(DCHLT, [qs[i]^2 * 0.5*PEN[i+1]*(I + R)])
        append!(NMANN, [0.5*PEN[i+1]*(I - R)])
    end

    return PEN, REF, DCHLT, NMANN
end

function grid_map(f::Function, outsize::Tuple{Int64, Int64}, grd::Grid)

    n, m = outsize
    output = SharedArray{Float64}((n*grd.Nk,m))

    # @distributed unfortunately only works if nprocs > 1
    if nprocs() > 1
        @distributed for i=1:grd.Nk
            k = 1 + n*(i-1)
            output[k:k+n-1,:] = f(grd.points[:,i]..., grd)
        end
    else
        for i=1:grd.Nk
            k = 1 + n*(i-1)
            output[k:k+n-1,:] = f(grd.points[:,i]..., grd)
        end
    end

    return output[:,:]
end


text:
initialize(c::IterativeMapCarrier) = rand(c.density)

function generate!(c::IterativeMapCarrier, x::Vector{Float64})
  x[1] = initialize(c)

  for i = 2:c.len
    x[i] = c.iterate(x[i-1])
  end

  return x
end

generate(c::IterativeMapCarrier) = generate!(c, Array(Float64, c.len))

mean(c::IterativeMapCarrier) = mean(c.density)

var(c::IterativeMapCarrier) = var(c.density)

immutable LogisticCarrier <: IterativeMapCarrier
  density::Distribution
  len::Int
  iterate::Function
end

LogisticCarrier(len::Int) = LogisticCarrier(Beta(0.5, 0.5), len, logistic)
LogisticCarrier() = LogisticCarrier(Beta(0.5, 0.5), 0, logistic)
LogisticCarrier(; density::Distribution=Beta(0.5, 0.5), len::Int=5, iterate::Function=logistic) =
  LogisticCarrier(density, len, iterate)

immutable FSPWL2BCarrier <: FSPWLCarrier # 2-branch fully-stretching piece-wise linear maps
  l::Float64 # minimum of domain
  u::Float64 # maximum of domain
  nc::Float64 # non-centrality parameter
  density::Distribution
  len::Int
  iterate::Function

  function FSPWL2BCarrier(l::Float64, u::Float64, nc::Float64, density::Distribution, len::Int, iterate::Function)
    @assert l<nc && nc<u "Non-centrality parameter must be in the interior of domain of Bernoulli map."
    new(l, u, nc, density, len, iterate)
  end
end

fspwl2b_types = (:bernoulli, :nbernoulli, :tent, :valley)

function FSPWL2BCarrier(l::Float64, u::Float64, nc::Float64, len::Int, mtype::Symbol)
  @assert in(mtype, fspwl2b_types) "$mtype is not a supported type of FSPWL2BCarrier map."

  if mtype == :bernoulli
    return FSPWL2BCarrier(l, u, nc, Uniform(l, u), len, (x::Float64)->bernoulli(x, l, u, nc))
  elseif mtype == :nbernoulli
    return FSPWL2BCarrier(l, u, nc, Uniform(l, u), len, (x::Float64)->nbernoulli(x, l, u, nc))
  elseif mtype == :tent
    return FSPWL2BCarrier(l, u, nc, Uniform(l, u), len, (x::Float64)->tent(x, l, u, nc))
  elseif mtype == :valley
    return FSPWL2BCarrier(l, u, nc, Uniform(l, u), len, (x::Float64)->valley(x, l, u, nc))
  end
end

FSPWL2BCarrier(; l::Float64=0.0, u::Float64=1.0, nc::Float64=0.5, density::Distribution=Uniform(0.0, 1.0), len::Int=5,
  iterate::Function=(x::Float64)->bernoulli(x, l, u, nc)) = FSPWL2BCarrier(l, u, nc, density, len, iterate)

BernoulliCarrier(l::Float64, u::Float64, nc::Float64, len::Int) = FSPWL2BCarrier(l, u, nc, len, :bernoulli)
BernoulliCarrier(l::Float64, u::Float64, len::Int) = BernoulliCarrier(l, u, 0.5*(l-u), len)
BernoulliCarrier(len::Int) = BernoulliCarrier(0.0, 1.0, 0.5, len)
BernoulliCarrier() = BernoulliCarrier(0.0, 1.0, 0.5, 0)

NBernoulliCarrier(l::Float64, u::Float64, nc::Float64, len::Int) = FSPWL2BCarrier(l, u, nc, len, :nbernoulli)
NBernoulliCarrier(l::Float64, u::Float64, len::Int) = NBernoulliCarrier(l, u, 0.5*(l-u), len)
NBernoulliCarrier(len::Int) = NBernoulliCarrier(0.0, 1.0, 0.5, len)
NBernoulliCarrier() = NBernoulliCarrier(0.0, 1.0, 0.5, 0)

TentCarrier(l::Float64, u::Float64, nc::Float64, len::Int) = FSPWL2BCarrier(l, u, nc, len, :tent)
TentCarrier(l::Float64, u::Float64, len::Int) = TentCarrier(l, u, 0.5*(l-u), len)
TentCarrier(len::Int) = TentCarrier(0.0, 1.0, 0.5, len)
TentCarrier() = TentCarrier(0.0, 1.0, 0.5, 0)

ValleyCarrier(l::Float64, u::Float64, nc::Float64, len::Int) = FSPWL2BCarrier(l, u, nc, len, :valley)
ValleyCarrier(l::Float64, u::Float64, len::Int) = ValleyCarrier(l, u, 0.5*(l-u), len)
ValleyCarrier(len::Int) = ValleyCarrier(0.0, 1.0, 0.5, len)
ValleyCarrier() = ValleyCarrier(0.0, 1.0, 0.5, 0)

immutable CircularCarrier <: IterativeMapCarrier
  nc::Float64 # non-centrality parameter
  density::Distribution
  len::Int
  iterate::Function

  function CircularCarrier(nc::Float64, density::Distribution, len::Int, iterate::Function)
    @assert -1.0<nc && nc<1.0 "Non-centrality parameter must be in the interior of domain of circular map."
    new(nc, density, len, iterate)
  end
end

CircularCarrier(nc::Float64, len::Int) = CircularCarrier(nc, VDist(nc), len, (x::Float64)->circular(x, nc))
CircularCarrier(len::Int) = CircularCarrier(0.42, VDist(0.42), len, (x::Float64)->circular(x, 0.42))
CircularCarrier() = CircularCarrier(0.42, VDist(0.42), 0, (x::Float64)->circular(x, 0.42))
CircularCarrier(; nc::Float64=0.42, density::Distribution=VDist(0.42), len::Int=5,
  iterate::Function=(x::Float64)->circular(x, 0.42)) = CircularCarrier(nc, density, len, iterate)


text:
using Pkg

function has_cov(pkg_dir)
    for (root, dirs, files) in walkdir(pkg_dir)
        if any(endswith(".cov"), files)
            return true
        end
    end
    return false
end

depot_packages = filter(isdir, readdir(joinpath(Base.DEPOT_PATH[1], "packages"), join = true))

found_dict = Dict{String,Bool}()
for dir in depot_packages
    pkgname = splitpath(dir)[end]
    try
        found_dict[pkgname] = has_cov(dir)
    catch ex
        @error "Could not search coverage for $dir" exception = (ex, catch_backtrace())
    end
end

deps = copy(Pkg.Types.Context().env.project.deps)
for dep in deps
    pkgname = first(dep)
    if !haskey(found_dict, pkgname) || found_dict[pkgname]
        delete!(deps, first(dep))
    end
end

if !isempty(deps)
    @error "Some direct dependencies were not used by the test code" deps

    @warn """This information should be used with caution.
    The check is only _at best_ as good as the coverage of the package tests, or provided test script.
    Also consider that the setup of the CI machine may impact coverage, with platform-guarded code
    usage hiding real code usage on other platforms.
    """
    exit(1)
end

max_stars_repo_name: julia-actions/julia-treeshake

max_issues_repo_name: julia-actions/julia-treeshake

max_forks_repo_name: julia-actions/julia-treeshake


text:
# Small script to read MATLAB data into Julia

using MAT

function load_EMPS_data()
	"Load from EMPS.mat files"

	# Open file
	file = matopen("./data/DATA_EMPS.mat")

	# Extract input-output signals
	vir = read(file, "vir")
	qm = read(file, "qm")

	# gc
	close(file)

	return vir, qm
end

max_stars_repo_name: wmkouw/NSI-EMPS

max_issues_repo_name: wmkouw/NSI-EMPS

max_forks_repo_name: wmkouw/NSI-EMPS


text:
# The following two test is based on a 4 layer network.
# We can consider this network is in the form of y=f(x), where x and y are 1d variables.
# The graph of this function looks like \/\/\/\/
# To get a precise rechable set, the verification algorithm needs to fully split
# the input set to reduce all over-approximation.

w_nnet = read_nnet(net_path * "spiky_nnet.nnet", last_layer_activation = NeuralVerification.Id())

ϵ = 1e-1

x_min = 1.0
x_max = 100.0
y_min = 42.51485148514851
y_max = 100.0

problem_holds = Problem(w_nnet,
                           Hyperrectangle(low = [x_min], high = [x_max]),
                           Hyperrectangle(low = [y_min - ϵ], high = [y_max + ϵ]));

problem_violated = Problem(w_nnet,
                           Hyperrectangle(low = [x_min], high = [x_max]),
                           Hyperrectangle(low = [y_min + ϵ], high = [y_max - ϵ]));

# NOTE: 'max_iter' of the solver must be large enough to fully split.
@time @testset "Fully split, ReluVal" begin
    solver = ReluVal(max_iter = 1000)
    @test solve(solver, problem_holds).status == :holds
    @test solve(solver, problem_violated).status == :violated
end

@time @testset "Fully split, Neurify" begin
    solver = Neurify(max_iter = 1000)
    @test solve(solver, problem_holds).status == :holds
    @test solve(solver, problem_violated).status == :violated
end


text:
"""
    module PickleIO

Functions to read and write Python pickle files.
"""
module PickleIO

export readpkl

# We use the following as a reference:
# https://formats.kaitai.io/python_pickle/

using ..PyObjects

include("opcode.jl")
include("stack.jl")
include("memo.jl")
include("unpickle.jl")

const HIGHEST_PROTOCOL = 5
const DEFAULT_PROTOCOL = 4

"""
    readpkl(io_or_filename; simplify=false)

Reads a .pkl file.
"""
function readpkl(io::IO; simplify=false)
    ans = unpickle(io)
    if simplify !== nothing
        ans = PyObjects.simplify(ans)
        ans = PyObjects.walk(simplify, ans)
    end
    return ans
end
readpkl(fn::AbstractString; kw...) = open(io->readpkl(io; kw...), fn)

end # module


text:
module GNSSSimulator

    using
        DocStringExtensions,
        Rotations,
        Unitful,
        CoordinateTransformations,
        StaticArrays,
        GNSSSignals,
        LinearAlgebra,
        Parameters,
        PhasedArray,
        Random,
        StructArrays,
        LoopVectorization

    import Base.transpose
    import Unitful: Hz, rad, s, m, dB, °, dBHz, upreferred
    import PhasedArray.get_steer_vec

    const EARTH_RADIUS = 6_360_000m
    const SPEED_OF_LIGHT = 299_792_458m/s

export
    NoisyStaticAttitude,
    DynamicAttitude,
    RandomWalkAttitude,
    NoisyDynamicAttitude,
    LinearDynamicAttitude,
    NoisyLinearDynamicAttitude,
    DynamicDOA,
    LinearDynamicDOA,
    RandomWalkDOA,
    DynamicExistence,
    ConstantDopplerSatellite,
    ConstantDopplerStructuralInterference,
    SyntheticSatellite,
    AsymptoticGainPhaseMismCrosstalk,
    RandomWalkGainPhaseMismCrosstalk,
    CWJammer,
    NoiseJammer,
    DynamicJNR,
    Noise,
    Receiver,
    get_attitude,
    get_doa,
    get_existence,
    get_gain_phase_mism_crosstalk,
    get_amplitude,
    get_carrier_doppler,
    get_carrier_phase,
    get_code_doppler,
    get_code_phase,
    get_gnss_system,
    get_carrier_to_noise_density_ratio,
    get_noise_density,
    get_noise_std,
    get_prn,
    get_id,
    get_jammer_to_noise_ratio,
    get_sampling_frequency,
    get_intermediate_frequency,
    get_measurement!,
    get_measurement,
    DynamicCN0

    propagate(val::Number, Δt, rng) = val
    const cart2sph = SphericalFromCartesian()
    const sph2cart = CartesianFromSpherical()

    function get_sampled_value(time, sample_freq, values)
        index = floor(Int, time * sample_freq) + 1
        index < size(values, 1) ? values[index] : values[end]
    end

    include("covdraw.jl")
    include("amplitude.jl")
    include("carrier.jl")
    include("attitude.jl")
    include("doa.jl")
    include("existence.jl")
    include("emitter.jl")
    include("satellite.jl")
    include("structural_interference.jl")
    include("synthetic_satellite.jl")
    include("jammer.jl")
    include("gain_phase_mism_crosstalk.jl")
    include("receiver.jl")
    include("measurement.jl")
end


text:
using LinearAlgebra, Statistics, BenchmarkTools, Interpolations, Parameters, Plots, QuantEcon, Roots, Optim, Random, IterTools, BasicInterpolators

function K!(Kg, g, grid, β, ∂u∂c, f, f′, shocks)
# This function requires the container of the output value as argument Kg
    # Construct radial basis functions
    #g_func = RBFInterpolator(grid, vec(g), 1)
    g_func = LinearInterpolation(vec(grid), vec(g), extrapolation_bc=Line())
    # solve for updated consumption value
    for i in 1:length(grid[:,1])
        y=grid[i,:]
        function h(c) # implies c is similar object to grid, this is not compatible with our root finding method currently
            vals = ∂u∂c.(g_func.(sum(f.(y - c), dims=2)  * shocks)) .* f′(y - c) .* shocks
            return ∂u∂c(c) - β * mean(vals)
        end
        Kg[i] = find_zero(h, (1e-10, y[1] .- 1e-10)) # would need to make this N (or at least 2-dimensional for a test)
    end
    return Kg
end

# The following function does NOT require the container of the output value as argument
K(g, grid, β, ∂u∂c, f, f′, shocks) = K!(similar(g), g, grid, β, ∂u∂c, f, f′, shocks)

# returns a cartesian product matrix that is compatible with the RBF interpolation
function cartProd(gridMin,gridMax,gridSize, numSec)
    # creates a vector of all the ranges we need
    ranges = fill(LinRange(gridMin,gridMax,gridSize),numSec)
    # initialise empty array
    grid = Array{Float64}(undef,gridSize^numSec,numSec)
    for (index, p) in enumerate(product(ranges...))
        # coordinates in each dimension are across rows
        grid[index, :] = collect(p)
    end
    return grid
end

function verify_true_policy(m, shocks, c_star)
    # compute (Kc_star)
    @unpack grid, β, ∂u∂c, f, f′ = m
    c_star_new = K(c_star, grid, β, ∂u∂c, f, f′, shocks)

    #compare
    return (c_star-c_star_new)./c_star
end

function check_convergence(m, shocks, c_star, g_init; n_iter = 15)
    @unpack grid, β, ∂u∂c, f, f′ = m
    g = g_init;
    for i in 1:n_iter
        new_g = K(g, grid, β, ∂u∂c, f, f′, shocks)
        g = new_g
    end
    #compare
    return (c_star-g)./c_star
end

Model=@with_kw (α = 0.65,                            # Productivity parameter
                β = 0.95,                            # Discount factor
                γ = 1.0,                             # Risk aversion
                μ = 0.0,                             # First parameter in lognorm(μ, σ)
                s = 0.1,                             # Second parameter in lognorm(μ, σ)
                numSec = 1,
                gridMin = 1e-6,                     # Smallest grid point
                gridMax = 4.0,                      # Largest grid point
                gridSize = 20,                     # Number of grid points
                grid = cartProd(gridMin, gridMax, gridSize, numSec), # Grid
                u = (c, γ = γ) -> isoelastic(c, γ),  # utility function
                ∂u∂c = c -> c^(-γ),                  # u′
                f = k -> k^α,                        # production function
                f′ = k -> α * k^(α - 1),             # f′
                )

m = Model();

Random.seed!(42) # for reproducible results.
shock_size = 250 # number of shock draws in Monte Carlo integral
shocks = collect(exp.(m.μ .+ m.s * randn(shock_size))); # generate shocks

c_star = sum((1 - m.α * m.β) .* m.grid, dims=2) # true policy (c_star)
verify_true_policy(m, shocks, c_star)

isoelastic(c, γ) = isone(γ) ? log(c) : (c^(1 - γ) - 1) / (1 - γ)

check_convergence(m, shocks, c_star, m.grid, n_iter = 15)

max_stars_repo_name: uq-aibe/aim

max_issues_repo_name: uq-aibe/aim

max_forks_repo_name: uq-aibe/aim


text:
export integral, integrate

const RTOL_DEFAULT = 1e-4

integral(args...; kwargs...) = integrate(args...; kwargs...).val

function integrate(vc::VoxelwiseConstant{K,AN}) where {K,AN}
    CI = CartesianIndices(size(vc))
    val = zero(eltype(AN))
    for ci = CI
        ∆lvxl_ci = map((∆lvxlₖ,iₖ) -> ∆lvxlₖ[iₖ], vc.∆lvxl, ci.I)
        val += vc.val[ci] * prod(∆lvxl_ci)
    end

    return (val=val, err=0.0)
end

function integrate(vc::VoxelwiseConstant,
                   dom;  # dom(x) returns true if x is inside domain; false otherwise
                   kwargs...)
    val, err = hcubature(x -> vc(x) * dom(x), bounds(vc)...; rtol=RTOL_DEFAULT, kwargs...)
    return (val=val, err=err)
end


text:
# Schrodinger equation
function SchrodingerEvo(H₀::Operator)
    # Constant Hamiltonian term
    L₀ = -1im*data(H₀)
    return Liouvillian(dims(H₀),L₀)
end
function SchrodingerEvo(H₀::Operator, Hₙ::Tuple{Vararg{Tuple}})
    dimsmatch(H₀,first.(Hₙ))
    # Constant Hamiltonian term
    L₀ = -1im*data(H₀)
    # Time-dependent Hamiltonian terms
    Lₙ,fₙ,pₙ = unpack_operators(-1im,data,Hₙ)
    return Liouvillian(dims(H₀),L₀,Lₙ,fₙ,pₙ)
end
SchrodingerEvo(H₀::Operator,Hₙ::Vararg{Tuple}) = SchrodingerEvo(H₀,Hₙ)
SchrodingerEvo(H::Tuple{Operator,Vararg{Tuple}}) = SchrodingerEvo(first(H),tail(H))
SchrodingerEvo(::Any) = throw(ArgumentError("invalid Hamiltonian specification"))

# von Neumann equation
# TODO

# Lindblad form master equation
function LindbladEvo(H₀::Operator, Cₘ::Tuple{Vararg{Operator}})
    dimsmatch(H₀,Cₘ)
    Id = data(qeye(dims(H₀)))
    # Constant Hamiltonian term
    L₀ = -1im*(Id⊗data(H₀) - transpose(data(H₀))⊗Id)
    # Add constant collapse operator terms
    L₀ += sum_collapse(Cₘ,Id,1)
    return Liouvillian(dims(H₀),L₀)
end
function LindbladEvo(H₀::Operator, Hₙ::Tuple{Vararg{Tuple}}, Cₘ::Tuple{Vararg{Operator}})
    dimsmatch(H₀,first.(Hₙ)); dimsmatch(H₀,Cₘ)
    Id = data(qeye(dims(H₀)))
    # Constant Hamiltonian term
    L₀ = -1im*(Id⊗data(H₀) - transpose(data(H₀))⊗Id)
    # Add constant collapse operator terms
    L₀ += sum_collapse(Cₘ,Id,1)
    # Time-dependent Hamiltonian terms
    Lₙ,fₙ,pₙ = unpack_operators(-1im,data,Hₙ)
    return Liouvillian(dims(H₀),L₀,@.((Id,)⊗Lₙ-transpose(Lₙ)⊗(Id,)),fₙ,pₙ)
end
LindbladEvo(H₀::Operator,Cₘ::Vararg{Operator}) = LindbladEvo(H₀,Cₘ)
LindbladEvo(H::Tuple{Operator,Vararg{Tuple}},Cₘ::Vararg{Operator}) = LindbladEvo(first(H),tail(H),Cₘ)
LindbladEvo(H::Tuple{Operator,Vararg{Tuple}},Cₘ::Tuple{Vararg{Operator}}) = LindbladEvo(first(H),tail(H),Cₘ)
LindbladEvo(::Any) = throw(ArgumentError("invalid Hamiltonian specification"))

# Schrodinger evo propagator
SchrodingerProp(H₀::Operator, tspan) = SchrodingerProp(H₀,float(tspan[2]-tspan[1]))
function SchrodingerProp(H₀::Operator, Δt::Real)
    # Constant Hamiltonian term
    U = expim(Hermitian(-full(H₀).*Δt))
    return Propagator(U,Δt,dims(H₀))
end
SchrodingerProp(H₀::Operator, Hₙ::Tuple, tspan, steps::Integer) = SchrodingerProp(H₀,(Hₙ,),tspan,steps)
function SchrodingerProp(H₀::Operator, Hₙ::Tuple{Vararg{Tuple}}, tspan, steps::Integer)
    dimsmatch(H₀,first.(Hₙ))
    F = real(eltype(H₀))
    # Sampling times and spacing dt
    t₁, t₂ = tspan; dt = (t₂-t₁)/steps
    # Unpack constant and time dep operators
    H0, Hn = full(H₀), unpack_operators(1,full,Hₙ)
    # Multiply sampled propagators together to generate total evolution
    U = Matrix{Complex{F}}(I, size(H0))
    H = Hermitian(zeros(compute_H_type(H0,Hn),size(H0)...))
    A = similar(U); B = similar(U); C = similar(H.data); D = similar(U)
    Λ = Vector{F}(undef,size(H,1))
    for i = 1:steps
        step_hamiltonian!(H.data,H0,Hn,(t₁,dt,i)) # calc H for this time step
        expim!(A,H,Λ,C,D) # A = exp(-1im*H*dt)
        mul!(B,A,U) # U = A*U
        U,B = B,U # swappitty swap for the next step
    end
    return Propagator(U,float(t₂-t₁),dims(H₀))
end
SchrodingerProp(H) = throw(ArgumentError("invalid Propagator specification"))

# Lindblad evo propagator
LindbladProp(H₀::Operator, Cₘ::Operator, tspan) = LindbladProp(H₀,(Cₘ,),tspan)
LindbladProp(H₀::Operator, Cₘ::Tuple{Vararg{Operator}}, tspan) = LindbladProp(H₀,Cₘ,tspan[2]-tspan[1])
function LindbladProp(H₀::Operator, Cₘ::Tuple{Vararg{Operator}}, Δt::Real)
    dimsmatch(H₀,Cₘ)
    F = real(eltype(H₀)); d = prod(dims(H₀))
    Id = Matrix{F}(I, d, d)
    # Constant Hamiltonian term
    L₀Δt = -1im*Δt.*(Id⊗full(H₀) .- transpose(full(H₀))⊗Id)
    # Add constant collapse operator terms
    L₀Δt .+= sum_collapse(Cₘ,Id,Δt)
    # Build constant propagator
    U = exp!(L₀Δt)
    return Propagator(U,Δt,dims(H₀))
end
LindbladProp(H₀::Operator, Hₙ::Tuple, Cₘ::Operator, tspan, steps::Integer) = LindbladProp(H₀,(Hₙ,),(Cₘ,),tspan,steps)
LindbladProp(H₀::Operator, Hₙ::Tuple, Cₘ::Tuple{Vararg{Operator}}, tspan, steps::Integer) = LindbladProp(H₀,(Hₙ,),Cₘ,tspan,steps)
LindbladProp(H₀::Operator, Hₙ::Tuple{Vararg{Tuple}}, Cₘ::Operator, tspan, steps::Integer) = LindbladProp(H₀,Hₙ,(Cₘ,),tspan,steps)
function LindbladProp(H₀::Operator, Hₙ::Tuple{Vararg{Tuple}}, Cₘ::Tuple{Vararg{Operator}}, tspan, steps::Integer)
    dimsmatch(H₀,first.(Hₙ)); dimsmatch(H₀,Cₘ)
    F = real(eltype(H₀)); d = prod(dims(H₀))
    Id = Matrix{F}(I, d, d)
    # Sampling times and spacing dt
    t₁, t₂ = tspan; dt = (t₂-t₁)/steps
    # Unpack constant and time dep operators
    H0, Hn = full(H₀), unpack_operators(1,full,Hₙ)
    # Build constant collapse propagator part
    U₀ = exp!(sum_collapse(Cₘ,Id,dt))
    # Multiply sampled propagators together to generate total evolution
    U = Matrix{Complex{F}}(I, size(U₀))
    H = Hermitian(zeros(compute_H_type(H0,Hn),size(H0)...))
    A = Matrix{Complex{F}}(undef,size(H0)...); B = similar(U); C = similar(H.data); D = similar(A)
    Λ = Vector{F}(undef,size(H,1))
    for i = 1:steps
        step_hamiltonian!(H.data,H0,Hn,(t₁,dt,i))
        expim!(A,H,Λ,C,D) # A = exp(-1im*H*dt)
        invA = Compat.LinearAlgebra.inv!(luf(A))
        mul!(B,U₀,U) # use the Lie product formula here for better results
        @static if VERSION < v"0.7.0-"
            At_mul_B!(U,invA⊗Id,B)
        else
            mul!(U,transpose(invA⊗Id),B)
        end
        I_kron_mul!(B,A,U)
        U,B = B,U # swappitty swap for the next step
    end
    return Propagator(U,float(t₂-t₁),dims(H₀))
end
LindbladProp(H) = throw(ArgumentError("invalid Propagator specification"))

# Utils
function unpack_operators(x::Number,f::Function,t::Tuple{Tuple{Operator,Function,Array},Vararg{Tuple}})
    a, b = first(t), unpack_operators(x,f,tail(t))
    return (x*f(a[1]), b[1]...), (a[2], b[2]...), (a[3], b[3]...)
end
function unpack_operators(x::Number,f::Function,t::Tuple{Tuple{Operator,Function},Vararg{Tuple}})
    a, b = first(t), unpack_operators(x,f,tail(t))
    return (x*f(a[1]), b[1]...), (a[2], b[2]...), ([], b[3]...)
end
unpack_operators(x::Number,f::Function,t::Tuple{}) = (), (), ()

function sum_collapse(Cₘ,Id,dt)
    # sum the Lindblad superoperators corresponding to the action of
    # D[Cᵢ](ρ) = CᵢρCᵢ† - 1/2 * (Cᵢ†Cᵢρ + ρCᵢ†Cᵢ)
    mapreduce(+,Cₘ) do Cᵢ
        C = full(Cᵢ); CdC = C'*C
        dt*(conj(C)⊗C - 0.5*(Id⊗CdC + transpose(CdC)⊗Id))
    end
end

function step_hamiltonian!(H,H0,Hn,tspec)
    Hₙ,fₙ,pₙ = Hn
    t₁, dt, i = tspec; tᵢ = t₁ + (i-1)*dt
    H .= -dt.*H0 # constant part
    for j = 1:length(Hₙ) # sample function at 3 points
        Hᵢ,fᵢ,pᵢ = Hₙ[j],fₙ[j],pₙ[j]
        H .+= -dt*(fᵢ(tᵢ,pᵢ) + fᵢ(tᵢ+0.5dt,pᵢ) + fᵢ(tᵢ+dt,pᵢ))/3 .* Hᵢ
    end
end

function compute_H_type(H0,Hn)
    T1 = promote_eltype(H0,Hn[1]...)
    T2 = promote_type(_typeof_H_funs(Hn[2],Hn[3])...)
    return promote_type(T1,T2)
end

@inline function _typeof_H_funs(funs::NTuple{N,Function},ps::NTuple{N,Vector}) where N
    this_T = typeof(first(funs)(0,first(ps)))
    return (this_T, _typeof_H_funs(tail(funs),tail(ps))...)
end

@inline _typeof_H_funs(funs::Tuple{},ps::Tuple{}) = ()


text:
immutable kNNClassifier
    t::BruteTree
    y::Vector
end

function knn(X::Matrix,
             y::Vector;
             metric::Metric = Euclidean())
    return kNNClassifier(BruteTree(X, metric), y)
end

# TODO: Don't construct copy of model.y just to extract majority vote
function StatsBase.predict(model::kNNClassifier,
                           x::Vector,
                           k::Integer = 1)
    inds, dists = NearestNeighbors.knn(model.t, x, k)
    return majority_vote(model.y[inds])
end

function StatsBase.predict!(predictions::Vector,
                            model::kNNClassifier,
                            X::Matrix,
                            k::Integer = 1)
    n = size(X, 2)
    # @assert eltype(predictions) == eltype(model.y)
    # @assert length(predictions) == n
    for i in 1:n
        predictions[i] = predict(model, X[:, i], k)
    end
    return predictions
end

function StatsBase.predict(model::kNNClassifier,
                           X::Matrix,
                           k::Integer = 1)
    predictions = Array{eltype(model.y)}(size(X, 2))
    predict!(predictions, model, X, k)
    return predictions
end


text:
function alpha_bounds(umin,umax)
    #= Upper and lower bounds for function e_1(u) = (u-sin(u))/u
       given upper/lower bounds on u
    Inputs:
        umin (array of floats): Lower bound on u for each cell
        umax (array of floats): Upper bound on u for each cell
    Outputs:
        (array of floats): Lower and upper bounds of error function
    =#
    num = length(umin)
    e1_umin = zeros(num); e1_umax = zeros(num)
    min_map = umin.!=0.0
    max_map = umax.!=0.0
    e1_umin[min_map] = 1 .-sin.(umin[min_map])./umin[min_map]
    e1_umax[max_map] = 1 .-sin.(umax[max_map])./umax[max_map]

    min_map = e1_umax .< e1_umin

    e1_min = copy(e1_umin)
    e1_max = copy(e1_umax)
    e1_max[min_map] = e1_umin[min_map]
    e1_min[min_map] = e1_umax[min_map]
    sign_map = (umin.<0.0) .& (umax.>0.0)
    e1_min[sign_map] .*=0.0
    return hcat(e1_min,e1_max)
end

function beta_bounds(umin,umax)
    #= Upper and lower bounds for function e_2(u) = (1-cos(u))/u
       given upper/lower bounds on u
    Inputs:
        umin (array of floats): Lower bound on u for each cell
        umax (array of floats): Upper bound on u for each cell
    Outputs:
        (array of floats): Lower and upper bounds of error function
    =#
    num = length(umin)
    e2_min = zeros(num); e2_max = zeros(num)
    min_map = umin.!=0.0
    max_map = umax.!=0.0

    e2_min[min_map] = (1 .-cos.(umin[min_map]))./umin[min_map]
    e2_max[max_map] = (1 .-cos.(umax[max_map]))./umax[max_map]
    return hcat(e2_min,e2_max)
end

function getNextSetBounds!(xBounds,yBounds,psiBounds,u0min,u0max,u1min,u1max;v0=200,v1=200)
    #= Get an over-approximated region of space reachable from each cell
    Inputs:
        xBounds (array of floats): lower/upper bounds on x for each cell
        yBounds (array of floats): lower/upper bounds on y for each cell
        psiBounds (array of floats): lower/upper bounds on psi for each cell
        u0min (array of floats): lower bounds on u0 (ownship turn rate) for each cell
        u0max (array of floats): Upper bounds on u0 (ownship turn rate) for each cell
        u1min (array of floats): lower bounds on u1 (intruder turn rate) for each cell
        u1max (array of floats): Upper bounds on u1 (intruder turn rate) for each cell
        v0 (float): Ownship speed (ft/s)
        v1 (float): Intruder speed (ft/s)
    Outputs:
        None (xBounds, yBounds, and psiBounds are changed to be the boundaries of reachable space for each cell
    =#
    dpsi = Float32(0.5)*(psiBounds[:,2].-psiBounds[:,1])
    v1=Float32(v1)
    v0=Float32(v0)

    xBounds.+= v0.*alpha_bounds(u0min,u0max) .+ v1.*cos.((psiBounds[:,1].+psiBounds[:,2])./Float32(2.0)).-v0
    yBounds.-=v0*beta_bounds(u0min,u0max)[:,end:-1:1] .- v1.*sin.((psiBounds[:,1].+psiBounds[:,2])./Float32(2.0))

    v1cpsi = v1.*cos.((psiBounds[:,1].+psiBounds[:,2])./Float32(2.0))
    v1cpsi_pos = v1cpsi.>0.0
    v1cpsi_neg = .!v1cpsi_pos
    xBounds[v1cpsi_pos,:].-= v1cpsi[v1cpsi_pos].*alpha_bounds(u1min[v1cpsi_pos],u1max[v1cpsi_pos])[:,end:-1:1]
    xBounds[v1cpsi_neg,:].-= v1cpsi[v1cpsi_neg].*alpha_bounds(u1min[v1cpsi_neg],u1max[v1cpsi_neg])
    yBounds[v1cpsi_pos,:].+= v1cpsi[v1cpsi_pos].*beta_bounds(u1min[v1cpsi_pos],u1max[v1cpsi_pos])
    yBounds[v1cpsi_neg,:].+= v1cpsi[v1cpsi_neg].*beta_bounds(u1min[v1cpsi_neg],u1max[v1cpsi_neg])[:,end:-1:1]

    xBounds[v1cpsi_pos,1].+= v1cpsi[v1cpsi_pos].*(cos.(dpsi[v1cpsi_pos]).-1)
    yBounds[v1cpsi_pos,2].+= v1cpsi[v1cpsi_pos].*sin.(dpsi[v1cpsi_pos])
    yBounds[v1cpsi_pos,1].-= v1cpsi[v1cpsi_pos].*sin.(dpsi[v1cpsi_pos])

    xBounds[v1cpsi_neg,2].+= v1cpsi[v1cpsi_neg].*(cos.(dpsi[v1cpsi_neg]).-1)
    yBounds[v1cpsi_neg,2].-= v1cpsi[v1cpsi_neg].*sin.(dpsi[v1cpsi_neg])
    yBounds[v1cpsi_neg,1].+= v1cpsi[v1cpsi_neg].*sin.(dpsi[v1cpsi_neg])
    v1cpsi = nothing; v1cpsi_pos=nothing; v1cpsi_neg=nothing; GC.gc()

    v1spsi = v1.*sin.((psiBounds[:,1].+psiBounds[:,2])./Float32(2.0))
    v1spsi_pos = v1spsi.>0.0
    v1spsi_neg = .!v1spsi_pos
    xBounds[v1spsi_pos,2].+= v1spsi[v1spsi_pos].*sin.(dpsi[v1spsi_pos])
    xBounds[v1spsi_pos,1].-= v1spsi[v1spsi_pos].*sin.(dpsi[v1spsi_pos])
    xBounds[v1spsi_pos,:].-= v1spsi[v1spsi_pos].*beta_bounds(u1min[v1spsi_pos],u1max[v1spsi_pos])[:,end:-1:1]
    yBounds[v1spsi_pos,1].+= v1spsi[v1spsi_pos].*(cos.(dpsi[v1spsi_pos]).-1)
    yBounds[v1spsi_pos,:].-= v1spsi[v1spsi_pos].*alpha_bounds(u1min[v1spsi_pos],u1max[v1spsi_pos])[:,end:-1:1]

    xBounds[v1spsi_neg,2].-= v1spsi[v1spsi_neg].*sin.(dpsi[v1spsi_neg])
    xBounds[v1spsi_neg,1].+= v1spsi[v1spsi_neg].*sin.(dpsi[v1spsi_neg])
    xBounds[v1spsi_neg,:].-= v1spsi[v1spsi_neg].*beta_bounds(u1min[v1spsi_neg],u1max[v1spsi_neg])
    yBounds[v1spsi_neg,2].+= v1spsi[v1spsi_neg].*(cos.(dpsi[v1spsi_neg]).-1)
    yBounds[v1spsi_neg,:].-= v1spsi[v1spsi_neg].*alpha_bounds(u1min[v1spsi_neg],u1max[v1spsi_neg])
    v1spsi=nothing; v1spsi_pos=nothing; v1spsi_cos=nothing; GC.gc()

    psiBounds[:,2].+=u1max.-u0min
    psiBounds[:,1].+=u1min.-u0max

    cos_u0_min = cos.(max.(abs.(u0max),abs.(u0min)))
    sin_u0_min = sin.(u0min)
    sin_u0_max = sin.(u0max)

    xBounds_rotMin = min.(xBounds[:,1],cos_u0_min.*xBounds[:,1])+min.(sin_u0_min.*yBounds[:,1],sin_u0_max.*yBounds[:,1],sin_u0_min.*yBounds[:,2],sin_u0_max.*yBounds[:,2])
    xBounds_rotMax = max.(xBounds[:,2],cos_u0_min.*xBounds[:,2])+max.(sin_u0_min.*yBounds[:,1],sin_u0_max.*yBounds[:,1],sin_u0_min.*yBounds[:,2],sin_u0_max.*yBounds[:,2])

    yBounds_rotMin = min.(yBounds[:,1],cos_u0_min.*yBounds[:,1])-max.(sin_u0_min.*xBounds[:,1],sin_u0_max.*xBounds[:,1],sin_u0_min.*xBounds[:,2],sin_u0_max.*xBounds[:,2])
    yBounds_rotMax = max.(yBounds[:,2],cos_u0_min.*yBounds[:,2])-min.(sin_u0_min.*xBounds[:,1],sin_u0_max.*xBounds[:,1],sin_u0_min.*xBounds[:,2],sin_u0_max.*xBounds[:,2])

    xBounds[:,:] = hcat(xBounds_rotMin,xBounds_rotMax)
    yBounds[:,:] = hcat(yBounds_rotMin,yBounds_rotMax)
    cos_u0_min=nothing; sin_u0_min=nothing; sin_u0_max=nothing; xBounds_rotMin = nothing; xBounds_rotMax=nothing; yBounds_rotMin=nothing; yBounds_rotMax=nothing; GC.gc()

end

function getCutpoints(XS)
    xStepsAll = XS[2:end]-XS[1:end-1]
    xSteps = Array{Float32,1}()
    xNums = Array{Int32,1}()
    xCuts = [XS[1]]
    currentStep = -1
    for (i,step) in enumerate(xStepsAll)
        if step != currentStep
            currentStep = step
            xSteps = vcat(xSteps,[step])
            xNums = vcat(xNums,[i])
            if length(xNums)>1
                xCuts = vcat(xCuts,[xCuts[end]+(xNums[end]-xNums[end-1])*xSteps[end-1]])
            end
        end
    end
    return (xSteps, xCuts, xNums)
end

function getIndexBounds(xBoundsNext,yBoundsNext,psiBoundsNext)
    #= Compute which cells overlap with region of state space for each region
    Inputs:
        xBoundsNext (array of floats): Lower/upper bounds on x for each region
        yBoundsNext (array of floats): Lower/upper bounds on y for each region
        psiBoundsNext (array of floats): Lower/upper bounds on psi for each region
    Outputs:
        xInds (array of ints): First/last index of x-array that overlap with the region
        yInds (array of ints): First/last index of y-array that overlap with the region
        psiInds (array of ints): First/last index of psi-array that overlap with the region
    =#

    # Only use cell if it is inside, not just touching the boundary. This small adjustment
    # prevents cells from being added that just touch the boundary of the region
    psiBoundsNext[:,1].+=1e-6
    psiBoundsNext[:,2].-=1e-6
    xBoundsNext[:,1].+=1e-6
    xBoundsNext[:,2].-=1e-6
    yBoundsNext[:,1].+=1e-6
    yBoundsNext[:,2].-=1e-6

    # The X/Y arrays of the grid are not uniform, but because it is piecewise linearly spaced, we
    # can compute the indices more efficiently than a brute force search
    xInds = ones(Int32,size(xBoundsNext)[1],2)
    xMap = xInds.>0
    xSteps, xCuts, xNums = getCutpoints(XS)
    for i=1:length(xSteps)
        xMap = xMap .& (xBoundsNext.>xCuts[i])
        xInds[xMap] = Int32.(div.(xBoundsNext[xMap].+(xSteps[i]*xNums[i]-xCuts[i]),xSteps[i]))
    end

    yInds = ones(Int32,size(yBoundsNext)[1],2)
    yMap = yInds.>0
    ySteps, yCuts, yNums = getCutpoints(YS)
    for i=1:length(ySteps)
        yMap = yMap .& (yBoundsNext.>yCuts[i])
        yInds[yMap] = Int32.(div.(yBoundsNext[yMap].+(ySteps[i]*yNums[i]-yCuts[i]),ySteps[i]))
    end

    # Psi is split uniformly, so mapping to indices is much quicker
    psiInds = ones(Int32,size(psiBoundsNext)[1],2)
    psiMap = psiInds.>0
    psiSteps, psiCuts, psiNums = getCutpoints(PSIS_DEG)
    for i=1:length(psiSteps)
        psiMap = psiMap .& (psiBoundsNext.>psiCuts[i])
        psiInds[psiMap] = Int32.(div.(psiBoundsNext[psiMap].*180.0./pi.+(psiSteps[i]*psiNums[i]-psiCuts[i]),psiSteps[i]))
    end

    xInds[xInds.>NUMX] .= NUMX
    yInds[yInds.>NUMY] .= NUMY
    psiInds[psiInds.>NUMP] .= NUMP

    psiInds[psiInds[:,1].>psiInds[:,2],1] .-=NUMP

    return xInds,yInds,psiInds
end

function getReachDynamics(xBoundInds,yBoundInds,psiBoundInds)
    #= Convert the first/last indices for each dimension to a sparse
       representation of the reachable dynamics. Output r has length the same as the number of cells.
       The cells reachable from cell i are given by c[r[i]:r[i+1]-1]

    This function is slow due to for-loops

    Inputs:
        xBoundInds (array of ints): First/last index of x-array that overlap with the region
        yBoundInds (array of ints): First/last index of y-array that overlap with the region
        psiBoundInds (array of ints): First/last index of psi-array that overlap with the region
    Outputs:
        r (array of ints): Pointers to the c array to determine where to begin reading
        c (array of ints): Array of cell indices
    =#
    dX = Int32(maximum(xBoundInds[:,2]-xBoundInds[:,1]))
    dY = Int32(maximum(yBoundInds[:,2]-yBoundInds[:,1]))
    dPsi = Int32(maximum(psiBoundInds[:,2]-psiBoundInds[:,1]))

    r = Vector{Int32}()
    c = Vector{Int32}()
    inds = Int32(1):NUMREGIONS

    mapX = BitArray(undef,NUMREGIONS)
    mapX .= true
    mapY = BitArray(undef,NUMREGIONS)
    mapPsi = BitArray(undef,NUMREGIONS)

    for i = Int32(0):dX
        println(@sprintf("%d of %d",i+1,dX+1))
        mapX .&= (xBoundInds[:,1].+i .<= xBoundInds[:,2])
        mapY.=mapX
        @time for j = Int32(0):dY
            mapY .&= (yBoundInds[:,1].+j .<= yBoundInds[:,2])
            mapPsi.=mapY
            for k = Int32(0):dPsi
                mapPsi .&= (psiBoundInds[:,1].+k .<= psiBoundInds[:,2])
                append!(c,inds[mapPsi])
                append!(r,indicesToIndex_Vector(psiBoundInds[mapPsi,1].+k,yBoundInds[mapPsi,1].+j,xBoundInds[mapPsi,1].+i))
            end
        end
    end
    sortIdx = sortperm(c)
    return r[sortIdx],c[sortIdx]
end

max_stars_repo_name: cherrywoods/HorizontalCAS

max_issues_repo_name: ahmed-irfan/HorizontalCAS

max_forks_repo_name: ahmed-irfan/HorizontalCAS


text:
mutable struct joint_angles_t <: LCMType
    utime::Int64
    robot_name::String
    num_joints::Int32
    joint_name::Vector{String}
    joint_position::Vector{Float64}
end

@lcmtypesetup(joint_angles_t,
    joint_name => (num_joints,),
    joint_position => (num_joints,)
)


text:
export OutOfSampleBootstrap, recursive_bootstrap!

## Simple dummy object type that will allow us to dispatch on the name
## of the bootstrap.
type OutOfSampleBootstrap{T}
end
OutOfSampleBootstrap(s::Symbol) = OutOfSampleBootstrap{s}()

## Main implementations of the different oos bootstraps. Different
## methods are invoked by dispatching on the OutOfSampleBootstrap
## abstract type (second argument).

## Naive approach used for bootstrap reality check:
function recursive_bootstrap!(b::AbstractVector,
                              ::OutOfSampleBootstrap{:naive},
                              block::Block,
                              f::AbstractVector)
    P = length(f)
    bootindex = Array(Int, P)
    fsum = sum(f)
    for i in 1:length(b)
        fsum_i = 0.0
        bootstrap_index!(bootindex, block)
        for t_i in bootindex
            fsum_i += f[t_i]
        end
        b[i] = (fsum_i - fsum) / sqrt(P)
    end
end

## Corradi-Swanson (2007)
function recursive_bootstrap!(b::AbstractVector,
                              boot::OutOfSampleBootstrap{:cs07_ols},
                              block::Block,
                              y::AbstractVector,
                              x::AbstractMatrix,
                              βhat::AbstractMatrix,
                              L::Function)
    P = size(βhat, 2)
    n = length(y)
    R = n - P
    bootindex = Array(Int, n)
    correction = cs07_correction(boot, y, x, βhat, L)
    for i in 1:length(b)
        fsum_boot = 0.0
        bootstrap_index!(bootindex, block)
        for j in 1:P
            mAdj_j = 2 * mean((y[1:(n-1)] - x[1:(n-1),:] * βhat[:,j]) .* x[1:(n-1),:], 1)'
            # Ugly and slow; replace with the right BLAS calls laters
            βboot_ij = ((x[bootindex[1:(R+j-1)],:]'*x[bootindex[1:(R+j-1)],:]) \
                        (x[bootindex[1:(R+j-1)],:] '* y[bootindex[1:(R+j-1)]] - mAdj_j))
            fsum_boot += L(y[bootindex[R+j]] - scalar(x[bootindex[R+j],:] * βboot_ij))
        end
        b[i] = (fsum_boot - correction) / sqrt(P)
    end
end

function cs07_correction(::OutOfSampleBootstrap{:cs07_ols},
                         y::AbstractVector,
                         x::AbstractMatrix,
                         βhat::AbstractMatrix,
                         L::Function)
    P = size(βhat, 2)
    n = length(y)
    fsum = 0.0
    for j in 1:P
        for s in 2:n
            fsum += L(y[s] - scalar(x[s,:] * βhat[:,j]))
        end
    end
    return fsum / (n-1)
end

## The new one proposed by my paper
function recursive_bootstrap!(b::AbstractVector,
                              boot::OutOfSampleBootstrap{:mine_ols},
                              block::Block,
                              y::AbstractVector,
                              x::AbstractMatrix,
                              R::Integer,
                              L::Function)
    n = length(y)
    P = n - R
    bootindex = Array(Int, n)
    correction = sum(L(y - x * (x \ y)))
    βhat_i = Array(Float64, size(x, 2), P)
    errors_i = Array(Float64, P)

    for i in 1:length(b)
        bootstrap_index!(bootindex, block)
        ## This is a terrible approach; I should pass the index as a
        ## separate argument and reorder in recursive_ols!
        recursive_ols!(βhat_i, errors_i, y, x, bootindex)
        b[i] = (sum(L(errors_i)) - correction) / sqrt(P)
    end
end

## Several variations that wrap the main versions for
## convenience.

## Allow statistic type to be specified as a symbol
recursive_bootstrap!(b, boot::Symbol, block, y, x, R, L) =
    recursive_bootstrap!(b, OutOfSampleBootstrap(boot), block, y, x, R, L)

## Default to MSE if loss function not specified
function recursive_bootstrap!(b, boot, block, y, x, R)
    L(x) = x .* x
    recursive_bootstrap!(b, boot, block, y, x, R, L)
end

## Alternative invokations for the naive oos bootstrap
function recursive_bootstrap!(b, boot::OutOfSampleBootstrap{:naive}, block, y, x, R::Integer, L::Function)
    P = length(y) - R
    βhat = Array(Float64, size(x, 2), P)
    errors = Array(Float64, P)
    recursive_ols!(βhat, errors, y, x)
    recursive_bootstrap!(b, boot, block, L(errors))
end

## Alternative invokation for the Corradi-Swanson oos bootstrap
function recursive_bootstrap!(b, boot::OutOfSampleBootstrap{:cs07_ols}, block, y, x, R::Integer, L::Function)
    βhat = recursive_ols(y, x, R)
    recursive_bootstrap!(b, boot, block, y, x, βhat, L)
end


text:
@show isfile("input.txt")
@show isfile("/input.txt")
@show isdir("docs")
@show isdir("/docs")
@show isfile("")
@show isfile("`Abdu'l-Bahá.txt")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@safetestset "Path traversal" begin

  @safetestset "Returns 401 unauthorised" begin
    using Genie
    using HTTP

    isdir(Genie.config.server_document_root) || mkdir(Genie.config.server_document_root)

    port = rand(8000:10_000)
    server = up(port)

    req = HTTP.request("GET", "http://localhost:$port////etc/hosts"; status_exception = false)
    @test req.status == (Sys.iswindows() ? 404 : 401)

    req = HTTP.request("GET", "http://localhost:$port/../../src/mimetypes.jl"; status_exception = false)
    @test req.status == 401

    Genie.AppServer.down!()
    server = nothing
  end

end


text:
using MultipleScattering

p1 = Particle([-2.,2.])
p2 = Particle([-2.,-2.])
particles = [p1,p2]

w_arr = collect(0.1:0.01:1.)
simulation = FrequencySimulation(particles, w_arr)

using Plots
pyplot()
plot(simulation)

savefig("plot_simulation.png")

plot(simulation,0.8)
savefig("plot_field.png")