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

"""
    split(object, fraction, [normal])

Split spatial `object` into two parts where the first
part has a `fraction` of the elements. Optionally, the
split is performed perpendicular to a `normal` direction.
"""
function split(object, fraction::Real, normal=nothing)
  if isnothing(normal)
    partition(object, FractionPartitioner(fraction))
  else
    partition(object, BisectFractionPartitioner(normal, fraction))
  end
end


text:
module X

using CompEcon
include("../../ModelAbstraction.jl")
using .ModelAbstraction

immutable BFZStateSpace <: AbstractStateSpace{2}
    k::EndogenousState{Float64}
    x::ExogenousState{Float64}
    g::ExogenousState{Float64}
    grid::Matrix{Float64}
    grid_transpose::Matrix{Float64}
    basis::Basis{2}
    basis_struct::BasisStructure{Direct}
end

immutable BFZModel
    agent::EpsteinZinAgent
    producer::CESProducer
    exog::ConstantVolatility
    state_space::BFZStateSpace
end

function BFZcv(;ρ::Real=-1, α::Real=-9, β::Real=0.99,
              η::Real=1/3, ν::Float64=0.1, δ::Real=0.025,
              lgbar::Real=0.004, A::Real=0., B::Real=1., vbar::Real=0.015^2,
              nk::Int=20, nx::Int=8, nϵ1::Int=3)
    # (ρ, α, β, η, ν, δ, lgbar, A, B, vbar, nk, nx, nϵ1) = (-1., -9., .99, .33, .1, .025, .004, 0., 1., .015^2, 20, 8, 9)
    # Create the Consumer
    BFZagent = EpsteinZinAgent(ρ, α, β)

    # Create means of Production
    BFZprod = CESProducer(η, ν, δ)

    # Create Exogenous Process
    exog = ConstantVolatility(A, B, vbar, nϵ1)

    #
    # Create the State Space
    #
    # First Create Basis
    basis = Basis(SplineParams(collect(linspace(15.0, 30.0, nk)), 0, 3),
                  SplineParams(collect(linspace(-.05, .05, nx)), 0, 1))
    basis_struc = BasisStructure(basis, Direct())

    # Get the implied grids
    grid, (kgrid, xgrid) = nodes(basis)
    grid_transpose = grid'
    ggrid = exp(lgbar + xgrid)

    xpgrid = Array(Float64, nϵ1, nx)
    gpgrid = Array(Float64, nϵ1, nx)
    for i=1:nx
        xpgrid[:, i] = step(exog, xgrid[i], exog.ϵ)
        gpgrid[:, i] = exp(lgbar + xpgrid[:, i])
    end

    # Create the states
    kstate = EndogenousState(kgrid)
    xstate = ExogenousState(xgrid, xpgrid)
    gstate = ExogenousState(ggrid, gpgrid)

    # Finally create space
    statespace = BFZStateSpace(kstate, xstate, gstate, grid, grid_transpose,
                               basis, basis_struc)

    return BFZModel(BFZagent, BFZprod, exog, statespace)
end

end

ext: jl

max_stars_repo_name: NYUEcon/GrowthModels

max_issues_repo_name: NYUEcon/GrowthModels

max_forks_repo_name: NYUEcon/GrowthModels


text:
# packages -
using LinearAlgebra
using GLPK

# my codes -
include("Stoichiometric.jl")
include("Parser.jl")
include("Flux.jl")


text:
using Documenter, PoincareInvariants

makedocs(
    sitename = "PoincareInvariants.jl",
    format = Documenter.HTML(),
    pages = ["Home" => "index.md",
             "Poincaré Invariants"    => "poincare_invariants.md",
             "1st Poincaré Invariant" => "poincare_invariant_1st.md",
             "2nd Poincaré Invariant" => "poincare_invariant_2nd.md",
             ]
)

deploydocs(
    repo   = "github.com/DDMGNI/PoincareInvariants.jl.git"
)


text:
import ModelParameters.Flatten

"""
    ParameterSource

Abstract supertypes for parameter source wrappers such as `Aux` and `Grid`.
These allow flexibly in that parameters can be retreived from various data sources.
"""
abstract type ParameterSource end

"""
    get(data::AbstractSimData, key::ParameterSource, I...)
    get(data::AbstractSimData, key::ParameterSource, I::Union{Tuple,CartesianIndex})

Allows parameters to be taken from a single value or a [`ParameterSource`](@ref)
such as another [`Grid`](@ref), an [`Aux`](@ref) array, or a [`Delay`](@ref).
"""
@propagate_inbounds Base.get(data::AbstractSimData, val, I...) = val
@propagate_inbounds Base.get(data::AbstractSimData, key::ParameterSource, I...) = get(data, key, I)
@propagate_inbounds Base.get(data::AbstractSimData, key::ParameterSource, I::CartesianIndex) = get(data, key, Tuple(I))

"""
    Aux <: ParameterSource

    Aux{K}()
    Aux(K::Symbol)

Use auxilary array with key `K` as a parameter source.

Implemented in rules with:

```julia
get(data, rule.myparam, I)
```

When an `Aux` param is specified at rule construction with:

```julia
rule = SomeRule(; myparam=Aux{:myaux})
output = ArrayOutput(init; aux=(myaux=myauxarray,))
```

If the array is a DimensionalData.jl `DimArray` with a `Ti` (time)
dimension, the correct interval will be selected automatically,
precalculated for each timestep so it has no significant overhead.

Currently this is cycled by default. Note that cycling may be incorrect
when the simulation timestep (e.g. `Week`) does not fit
equally into the length of the time dimension (e.g. `Year`).
This will reuire a `Cyclic` index mode in DimensionalData.jl in future
to correct this problem.
"""
struct Aux{K} <: ParameterSource end
Aux(key::Symbol) = Aux{key}()

_unwrap(::Aux{X}) where X = X
_unwrap(::Type{<:Aux{X}}) where X = X

@propagate_inbounds Base.get(data::AbstractSimData, key::Aux, I::Tuple) = _getaux(data, key, I)

@inline aux(nt::NamedTuple, ::Aux{Key}) where Key = nt[Key]

# If there is no time dimension we return the same data for every timestep
@propagate_inbounds _getaux(data::AbstractSimData, key::Union{Aux,Symbol}, I::Tuple) = _getaux(aux(data, key), data, key, I)
@propagate_inbounds _getaux(A::AbstractMatrix, data::AbstractSimData, key::Union{Aux,Symbol}, I::Tuple) = A[I...]
function _getaux(A::AbstractDimArray{<:Any,3}, data::AbstractSimData, key::Union{Aux,Symbol}, I::Tuple)
    A[I..., auxframe(data, key)]
end

function boundscheck_aux(data::AbstractSimData, auxkey::Aux{Key}) where Key
    a = aux(data)
    (a isa NamedTuple && haskey(a, Key)) || _auxmissingerror(Key, a)
    gsize = gridsize(data)
    asize = size(a[Key], 1), size(a[Key], 2)
    if asize != gsize
        _auxsizeerror(Key, asize, gsize)
    end
    return true
end

@noinline _auxmissingerror(key, aux) = error("Aux data $key is not present in aux")
@noinline _auxsizeerror(key, asize, gsize) =
    error("Aux data $key is size $asize does not match grid size $gsize")

_calc_auxframe(data) = _calc_auxframe(aux(data), data)
_calc_auxframe(aux::NamedTuple, data) = map(A -> _calc_auxframe(A, data), aux)
function _calc_auxframe(A::AbstractDimArray, data)
    hasdim(A, Ti) || return nothing
    curtime = currenttime(data)
    firstauxtime = first(dims(A, TimeDim))
    auxstep = step(dims(A, TimeDim))
    # Use julias range objects to calculate the distance between the
    # current time and the start of the aux
    i = if curtime >= firstauxtime
        length(firstauxtime:auxstep:curtime)
    else
        1 - length(firstauxtime-timestep(data):-auxstep:curtime)
    end
    _cyclic_index(i, size(A, 3))
end
_calc_auxframe(aux, data) = nothing

function _cyclic_index(i::Integer, len::Integer)
    return if i > len
        rem(i + len - 1, len) + 1
    elseif i <= 0
        i + (i ÷ len -1) * -len
    else
        i
    end
end

"""
    Grid <: ParameterSource

    Grid{K}()
    Grid(K::Symbol)

Use grid with key `K` as a parameter source.

Implemented in rules with:

```julia
get(data, rule.myparam, I)
```

And specified at rule construction with:

```julia
SomeRule(; myparam=Grid(:somegrid))
```
"""
struct Grid{K} <: ParameterSource end
Grid(key::Symbol) = Grid{key}()

_unwrap(::Grid{X}) where X = X
_unwrap(::Type{<:Grid{X}}) where X = X

@propagate_inbounds Base.get(data::AbstractSimData, key::Grid{K}, I::Tuple) where K = data[K][I...]

"""
    AbstractDelay <: ParameterSource

Abstract supertype for [`ParameterSource`](@ref)s that use data from a grid
with a time delay.

$EXPERIMENTAL
"""
abstract type AbstractDelay{K} <: ParameterSource end

@inline frame(delay::AbstractDelay) = delay.frame

@propagate_inbounds function Base.get(data::AbstractSimData, delay::AbstractDelay{K}, I::Tuple) where K
    _getdelay(frames(data)[frame(delay, data)], delay, I)
end

@propagate_inbounds function _getdelay(frame::NamedTuple, delay::AbstractDelay{K}, I::Tuple) where K
    _getdelay(frame[K], delay, I)
end
@propagate_inbounds function _getdelay(frame::AbstractArray, ::AbstractDelay{K}, I::Tuple) where K
    frame[I...]
end

"""
    Delay <: AbstractDelay

    Delay{K}(steps)

`Delay` allows using a [`Grid`](@ref) from previous timesteps as a parameter source as
a field in any `Rule` that uses `get` to retrieve it's parameters.

It must be coupled with an output that stores all frames, so that `@assert
DynamicGrids.isstored(output) == true`.  With [`GraphicOutput`](@ref)s this may be
acheived by using the keyword argument `store=true` when constructing the output object.

# Type Parameters

- `K::Symbol`: matching the name of a grid in `init`.

# Arguments

- `steps`: As a user supplied parameter, this is a multiple of the step size of the output
    `tspan`. This is automatically replaced with an integer for each step. Used within the
    code in a rule, it must be an `Int` number of frames, for performance.

# Example

```julia
SomeRule(;
    someparam=Delay(:grid_a, Month(3))
    otherparam=1.075
)
`` `

$EXPERIMENTAL
"""
struct Delay{K,S} <: AbstractDelay{K}
    steps::S
end
Delay{K}(steps::S) where {K,S} = Delay{K,S}(steps)

ConstructionBase.constructorof(::Delay{K}) where K = Delay{K}
steps(delay::Delay) = delay.steps

# _to_frame
# Replace the delay step size with an integer for fast indexing
# checking that the delay matches the simulation step size.
# Delays at the start just use the init frame, for now.
function _to_frame(delay::Delay{K}, data) where K
    nsteps = steps(delay) / step(tspan(data))
    isteps = trunc(Int, nsteps)
    nsteps == isteps || _delaysteperror(delay, step(tspan(data)))
    frame = max(currentframe(data) - isteps, 1)
    Frame{K}(frame)
end

"""
    Lag <: AbstractDelay

    Lag{K}(frames::Int)

`Lag` allows using a [`Grid`](@ref) from a specific previous frame from within a rule,
using `get`. It is similar to [`Delay`](@ref), but an integer amount of steps should be
used, instead of a quantity related to the simulation `tspan`. Used within rule code,
the lower bound will not be checked. Do this manually, or use [`Frame`](@ref) instead.

# Type Parameter

- `K::Symbol`: matching the name of a grid in `init`.

# Argument

- `frames::Int`: number of frames to lag by, 1 or larger.

# Example

```julia
SomeRule(;
    someparam=Delay(:grid_a, Month(3))
    otherparam=1.075
)
`` `

$EXPERIMENTAL
"""
struct Lag{K} <: AbstractDelay{K}
    nframes::Int
end

function _to_frame(lag::Lag{K}, data) where K
    frame = max(currentframe(data) - lag.nframes, 1)
    Frame{K}(frame)
end

@inline frame(lag::Lag, data) = max(1, currentframe(data) - lag.nframes)

"""
    Frame <: AbstractDelay

    Frame{K}(frame)

`Frame` allows using a [`Grid`](@ref) from a specific previous timestep from within
a rule, using `get`. It should only be used within rule code, not as a parameter.

# Type Parameter

- `K::Symbol`: matching the name of a grid in `init`.

# Argument

- `frame::Int`: the exact frame number to use.

$EXPERIMENTAL
"""
struct Frame{K} <: AbstractDelay{K}
    frame::Int
    Frame{K}(frame::Int) where K = new{K}(frame)
end

ConstructionBase.constructorof(::Frame{K}) where K = Frame{K}

@inline frame(delay::Frame) = delay.frame
@inline frame(delay::Frame, data) = frame(delay)

# Delay utils

const DELAY_IGNORE = Union{Function,SArray,AbstractDict,Number}

@inline function hasdelay(rules::Tuple)
    # Check for Delay as parameter or used in rule code
    length(_getdelays(rules)) > 0 || any(map(needsdelay, rules))
end

needsdelay(rule::Rule) = false

# _setdelays
# Update any Delay anywhere in the rules Tuple.
function _setdelays(rules::Tuple, data)
    delays = _getdelays(rules)
    if length(delays) > 0
        newdelays = map(d -> _to_frame(d, data), delays)
        _setdelays(rules, newdelays)
    else
        rules
    end
end

_getdelays(rules::Tuple) = Flatten.flatten(rules, AbstractDelay, DELAY_IGNORE)
_setdelays(rules::Tuple, delays::Tuple) =
    Flatten.reconstruct(rules, delays, AbstractDelay, DELAY_IGNORE)

_notstorederror() =
    throw(ArgumentError("Output does not store frames, which is needed for a Delay. Use ArrayOutput or any GraphicOutput with `store=true` keyword"))
_delaysteperror(delay::Delay{K}, simstep) where K =
    throw(ArgumentError("Delay $K size $(steps(delay)) is not a multiple of simulations step $simstep"))


text:
using InterpretableValidation
using Test
using Distributions
using Random

## IID tests
x = collect(1.:100.)
uniform  = IID(x, Uniform(2,3))
@test uniform isa IID{Uniform{Float64}}
normal = IID(x, Normal(0.5, 1.))
@test normal isa IID{Normal{Float64}}
categorical = IID(x, Categorical(5))
@test categorical isa IID{Categorical{Float64, Array{Float64,1}}}
bernoulli = IID(x, Bernoulli(0.1))
@test bernoulli isa IID{Bernoulli{Float64}}

@test N_pts(uniform) == 100

@test logpdf(uniform, rand(100)) == -Inf
@test logpdf(uniform, 2.5*ones(100)) == 0
@test logpdf(normal, ones(100)) ≈ 100*logpdf(Normal(0.5, 1.), 1.)
@test logpdf(categorical, ones(Int, 100)) ≈ 100*log(0.2)
@test logpdf(bernoulli, fill(true, 100)) ≈ 100*log(0.1)

v = rand(Random.GLOBAL_RNG, uniform)
@test all((v .> 2.) .& (v .< 3.))
v = rand(Random.GLOBAL_RNG, uniform, 2.5*ones(100), 2.7*ones(100))
@test all((v .> 2.5) .& (v .< 2.7))

v = rand(Random.GLOBAL_RNG, normal)
@test all((v .> -4.) .& (v .< 4.))
v = rand(Random.GLOBAL_RNG, normal, 2.5*ones(100), 2.7*ones(100))
@test all((v .> 2.5) .& (v .< 2.7))

v = rand(Random.GLOBAL_RNG, categorical)
@test all([v[i] in [1,2,3,4,5] for i=1:100])

v = rand(Random.GLOBAL_RNG, categorical, fill([1], 100))
@test all(v.==1)

v = rand(Random.GLOBAL_RNG, bernoulli)
@test sum(v) < 20

v = rand(Random.GLOBAL_RNG, bernoulli, fill([1], 100))
@test all(v)

## GaussianProcess Tests
gp = GaussianProcess(m = (x) -> 0, k = squared_exp_kernel(l=100), x = x)

@test gp.m(gp.x) == 0
@test gp.k(gp.x[1], gp.x[1]) == 1
@test gp.k(gp.x[1], gp.x[2]) < 1
@test gp.x == x
@test gp.σ2 == 1e-6

@test logpdf(gp, rand(100)) < 0

v = rand(Random.GLOBAL_RNG, gp)
@test length(v) == 100
l = collect(range(-1.5, 0, length=100))
u = collect(range(0, 1.5, length=100))
v = rand(Random.GLOBAL_RNG, gp, l, u)
ϵ = 0.18
@test all((v .>= l .- ϵ)) && all((v .<= u .+ ϵ))


text:
using Gadfly
using ElectronDisplay
using Electron
using VegaDatasets
using Test

Electron.prep_test_env()

struct DummyDisplayable{M}
    data
end

Base.show(io::IO, m::M, x::DummyDisplayable{M}) where {M} = true
Base.show(io::IO, m::M, x::DummyDisplayable{M}) where {M} = print(io, x.data)

include("construct_vega_specs.jl")

@testset "ElectronDisplay" begin

p1 = plot(y=[1,2,3])

f = display(p1)

@test f isa Electron.Window

ElectronDisplay.CONFIG.single_window = true

f = display(p1)

@test f === ElectronDisplay._getglobalwindow()

p2 = plot(y=[1,2,3])
f2 = display(p2)

@test f2 === f   # Window is reused

eldt = ElectronDisplay.ElectronDisplayType()

@test electrondisplay(vl2) isa Electron.Window
@test electrondisplay(vl3) isa Electron.Window
@test electrondisplay(vg3) isa Electron.Window
@test electrondisplay(vg4) isa Electron.Window
@test electrondisplay(vg5) isa Electron.Window

@test_logs(
    (:warn, r"The size of JSON representation.*exceeds.*max_json_bytes"),
    electrondisplay(vl3png; max_json_bytes=-1)::Electron.Window
)

mdo = DummyDisplayable{MIME"text/markdown"}("""foo""")
@test electrondisplay(mdo) isa Electron.Window

pngo = DummyDisplayable{MIME"image/png"}("""fakedata""")
@test electrondisplay(pngo) isa Electron.Window

dro = DummyDisplayable{MIME"application/vnd.dataresource+json"}("""{"schema":{"fields":[{"name": "Miles_per_Gallon","type": "number"}]},"data":[{"Miles_per_Gallon":18}]}""")
@test electrondisplay(dro) isa Electron.Window

@test displayable(eldt, "text/html") == true
@test displayable(eldt, "text/markdown") == true
@test displayable(eldt, "image/png") == true
@test displayable(eldt, "image/svg+xml") == true
@test displayable(eldt, "application/vnd.vegalite.v2+json") == true
@test displayable(eldt, "application/vnd.vegalite.v3+json") == true
@test displayable(eldt, "application/vnd.vegalite.v4+json") == true
@test displayable(eldt, "application/vnd.vega.v3+json") == true
@test displayable(eldt, "application/vnd.vega.v4+json") == true
@test displayable(eldt, "application/vnd.vega.v5+json") == true
@test displayable(eldt, "application/vnd.plotly.v1+json") == true
@test displayable(eldt, "application/vnd.dataresource+json") == true

@testset "smoke test: single_window=$single_window focus=$focus " for
        single_window in [false, true],
        focus in [false, true]
    ElectronDisplay.CONFIG.single_window = single_window
    ElectronDisplay.CONFIG.focus = focus
    @test electrondisplay(dataset("cars")) isa Electron.Window
    @test electrondisplay(@doc reduce) isa Electron.Window

    config = (single_window=single_window, focus=focus)
    @test electrondisplay(dataset("cars"); config...) isa Electron.Window
    @test electrondisplay(@doc reduce; config...) isa Electron.Window
end

end


text:
using DataStructures
# author: Younis Shah

# Gates: H, CNOT, S, Pauli gates

b₁ = [0, 1]
b₀ = [1, 0]

# Hadamard gate
# H = [1\√2  1\√2]
#     [1\√2 -1\√2]
H = round(1/sqrt(2) * [1 1; 1 -1], 3)

# Pauli X gate
σₓ = [0 1;
      1 0]

# Pauli Y gate
Y = [0 -im;
     im 0]

# Pauli Z gate
Z = [1 0;
     0 -1]

# Phase gate
# S= [1 0]
#    [0 ι]
S = [1 0;
     0 im]

# Controlled NOT gate
CNOT = [1 0 0 0;
        0 1 0 0;
        0 0 0 1;
        0 0 1 0]

# 2x2 Identity matrix
Ι = eye(Int64, 2)

gate_symbol_list = (:X, :Y, :Z, :H, :S)

gate_list = ("I", "H", "X", "Y", "Z", "CNOT", "S")

# |0⟩ = [1]
#       [0]
pure_state_0 = [1, 0]
#|1⟩ = [0]
#      [1]
pure_state_1 = [0, 1]

# Define a new Qubit type
# Rest -TODO
type Qubit
  index::Int64
  ket::Array
end

# Define a QuantumRegister type
# Rest - TODO
type QuantumRegister
  register::Array{Qubit}
end

# Define a new Stage output type
type StageOutput
  index::Int64
  output::Union(Array, SparseMatrixCSC)
end

# Define a new type to store all the Stage outputs
type StageOutputDeque
  stages::Deque{StageOutput}
end

# Define a new Gate type using which gates are represnted as matrices
type Gate
  qubit::Int64
  name::String
  repr::Array
end

# Define a new Gates type: the number of gates and qubits
# Represents total gates per stage
type Gates
  gates::Array{Gate}
  nqubits::Int64
end

max_stars_repo_name: younisshah/cherokee

max_issues_repo_name: younisshah/cherokee

max_forks_repo_name: younisshah/cherokee


text:
#  using NetABM, LightGraphs
using NetABM
using LightGraphs
using StatsBase
using Plots
using Statistics
using Bootstrap
using LaTeXStrings

num_agents = 1000
#  g = erdos_renyi(num_agents,0.4)
g = barabasi_albert(num_agents,9)
agents = [Agent(i) for i = 1:num_agents]
init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
set_coop_agents!(agents;p_cop=0.1)
map(x -> assign_contacts!(g,x), agents);
get_coop!(agents)

S = Array{Float64}(undef,0)
I = Array{Float64}(undef,0)
R = Array{Float64}(undef,0)

for i in 1:100
    next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.7, R=true)
    update_state!(agents)
    push!(S,([ag.state for ag in agents if ag.state=="S"] |> length)/num_agents)
    push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
    push!(R,([ag.state for ag in agents if ag.state=="R"] |> length)/num_agents)
end

plot(1:length(S),S,label="S")
plot!(1:length(S),I,label="I")
plot!(1:length(S),R,label="R")

#  What about levels of cooperation?
#

n_boot = num_agents
some_data = randn(100)
bs1 = bootstrap(std, some_data, BasicSampling(n_boot))
cil = 0.95
bci1 = confint(bs1, BasicConfInt(cil))
II = [zeros(100) for i in 1:100]

g = barabasi_albert(5000,10)
global p = plot()
begin
    Threads.@threads for j in [0.0,0.3,0.6,0.9]
        agents = [Agent(i) for i = 1:5000]
        init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
        set_coop_agents!(agents;p_cop=j)
        map(x -> assign_contacts!(g,x), agents);
        get_coop!(agents)
        I = Array{Float64}(undef,0)
        Threads.@threads for i in 1:100
            next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.7, R=true)
            update_state!(agents)
            push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/5000)
        end
        p = plot!(1:length(I),I,label=string(j))
    end
end
p

#  Changing the average degree of the network

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 10
global cil = 0.95
global j
        #  Threads.@threads
global p = plot()
begin
    for j in [2, 6, 10, 20]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,j)
            agents = [Agent(i) for i = 1:num_agents]
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
            set_coop_agents!(agents;p_cop=0.9)
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents)
            I = Array{Float64}(undef,0)
            for i in 1:100
                next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.5, R=true)
                update_state!(agents)
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
            end
            push!(Iall,I)
            #  p = plot!(1:length(I),I,label=string(j))
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.35,label=string(j))
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

#  step = 1
#  for step in eachindex(Iall[1])
#      bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
#      bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
#      bci1 = confint(bs1, BasicConfInt(cil))
#      bci2 = confint(bs2, BasicConfInt(cil))
#      push!(Imean,bci1[1][1])
#      push!(Ilow,bci2[1][2])
#      push!(Ihigh,bci2[1][3])
#  end
#
#  Imean
#  Ilow
#  Ihigh
#  plot(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.5)
#  plot(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.5)

#  Changing the cooperation according to neighboors

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 10
global cil = 0.95
global j
        #  Threads.@threads
global p = plot()
begin
    for j in [0.1, 0.4, 0.7]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents]
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
            set_coop_agents!(agents;p_cop=0.4)
            set_adapt_agents!(agents;p_cop=j)
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents)
            I = Array{Float64}(undef,0)
            for i in 1:100
                next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.8, R=true)
                update_state!(agents)
                update_coop!(agents,0.3)
                get_coop!(agents)
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
            end
            push!(Iall,I)
            #  p = plot!(1:length(I),I,label=string(j))
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.35,label=string(j))
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

#  Changing the cooperation according to previous infections

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 10
global cil = 0.95
global j = 0.1
        #  Threads.@threads
global p = plot()
begin
    for j in [0.1, 0.4, 0.7, 1]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents]
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
            set_coop_agents!(agents;p_cop=0.4)
            set_adapt_agents!(agents;p_cop=j)
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents)
            I = Array{Float64}(undef,0)
            for i in 1:100
                next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.8, R=false)
                update_state!(agents)
                update_coop_infections!(agents,i;lrt=false)
                get_coop!(agents)
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
            end
            push!(Iall,I)
            #  p = plot!(1:length(I),I,label=string(j))
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.35,label=string(j))
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

p=plot!(xlims=(0,100),ylims=(0,1))
p=plot!(xlabel = L"t_i", ylabel = L"I/N")
savefig(p,"../figs/sis_adapting_infections.png")

#  Changing the cooperation according to previous infections including lrt

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 10
global cil = 0.95
global j
        #  Threads.@threads
global p = plot()
begin
    for j in [0.1, 0.4, 0.7, 1]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents]
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
            set_coop_agents!(agents;p_cop=0.4)
            set_adapt_agents!(agents;p_cop=j)
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents)
            I = Array{Float64}(undef,0)
            for i in 1:100
                next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.8, R=false)
                update_state!(agents)
                update_coop_infections!(agents,i;lrt=true)
                get_coop!(agents)
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
            end
            push!(Iall,I)
            #  p = plot!(1:length(I),I,label=string(j))
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.35,label=string(j))
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

p=plot!(xlims=(0,100),ylims=(0,1))
p=plot!(xlabel = L"t_i", ylabel = L"I/N")
savefig(p,"../figs/sir_adapting_infections_incl_lrt.png")

#  Changing the cooperation according to neighboors given distance

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 10
global cil = 0.95
global j
        #  Threads.@threads
global p = plot()
begin
    for j in [1,2,3]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents]
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
            set_coop_agents!(agents;p_cop=0.4)
            set_adapt_agents!(agents;p_cop=0.4)
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents)
            I = Array{Float64}(undef,0)
            for i in 1:100
                next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.8, R=true)
                update_state!(agents)
                update_coop_distance!(agents,g,j,0.3;lrt=true)
                get_coop!(agents)
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
            end
            push!(Iall,I)
            #  p = plot!(1:length(I),I,label=string(j))
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.35,label=string(j))
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

p=plot!(xlims=(0,100),ylims=(0,1))
p=plot!(xlabel = L"t_i", ylabel = L"I/N")
savefig(p,"../figs/sir_adapting_distance.png")

#  Changing the cooperation according to neighboors infections given distance

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 10
global cil = 0.95
global j
        #  Threads.@threads
global p = plot()
begin
    for j in [0,1]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents]
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
            set_coop_agents!(agents;p_cop=0.4)
            set_adapt_agents!(agents;p_cop=0.4)
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents)
            I = Array{Float64}(undef,0)
            for i in 1:100
                next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.8, R=false)
                update_state!(agents)
                update_coop_distance_inf!(agents,g,j,i;lrt=true)
                get_coop!(agents)
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
            end
            push!(Iall,I)
            #  p = plot!(1:length(I),I,label=string(j))
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.35,label=string(j))
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

p=plot!(xlims=(0,100),ylims=(0,1))
p=plot!(xlabel = L"t_i", ylabel = L"I/N")
savefig(p,"figs/sis_adapting_distance_infections_neighborhood.png")

#  Changing the cooperation according to first neighboors infection

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 10
global cil = 0.95
global j
        #  Threads.@threads
global p = plot()
begin
    for j in [0.1, 0.4, 0.7, 1]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents]
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
            set_coop_agents!(agents;p_cop=0.4)
            set_adapt_agents!(agents;p_cop=j)
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents)
            I = Array{Float64}(undef,0)
            for i in 1:100
                next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.8, R=false)
                update_state!(agents)
                update_coop_distance_inf!(agents,g,1,i;lrt=false)
                get_coop!(agents)
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
            end
            push!(Iall,I)
            #  p = plot!(1:length(I),I,label=string(j))
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.35,label=string(j))
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

p=plot!(xlims=(0,100),ylims=(0,1))
p=plot!(xlabel = L"t_i", ylabel = L"I/N")
savefig(p,"figs/sis_adapting_first_neighboors.png")

#  Changing the cooperation using infection counter
#  first infection.

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 10
global cil = 0.95
#  global d = 0
#  global probs = [1]
global j
global adoptprob = (1,[0.5,0.3])
        #  Threads.@threads
global p = plot(xlims=(0,100),ylims=(0,1),xlabel = L"t_i", ylabel = L"I/N")
begin
    for j in [0.1, 0.4, 0.7, 1]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents];
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1]);
            set_coop_agents!(agents;p_cop=0.4);
            set_adapt_agents!(agents;p_cop=j);
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents);
            I = Array{Float64}(undef,0);
            for i in 1:100
                next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.8, R=false);
                update_state!(agents);
                update_coop_given_distance!(agents,g,adoptprob[1],1,adoptprob[2];lrt=false);
                get_coop!(agents);
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents);
            end
            push!(Iall,I)
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.35,label=string(j))
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

savefig(p,"figs/sis_adapting_d1.png")
p

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 10
global cil = 0.95
#  global d = 0
#  global probs = [1]
global j
global adoptprob = (1,[0.5,0.3])
        #  Threads.@threads
global p = plot(xlims=(0,100),ylims=(0,1),xlabel = L"t_i", ylabel = L"I/N")
begin
    for j in [0.1, 0.4, 0.7, 1]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents];
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1]);
            set_coop_agents!(agents;p_cop=0.4);
            set_adapt_agents!(agents;p_cop=j);
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents);
            I = Array{Float64}(undef,0);
            for i in 1:100
                next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.8, R=true);
                update_state!(agents);
                update_coop_given_distance!(agents,g,adoptprob[1],1,adoptprob[2];lrt=false);
                get_coop!(agents);
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents);
            end
            push!(Iall,I)
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.35,label=string(j))
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

savefig(p,"figs/sir_adapting_d1.png")
p

###################################################
###################################################
###################################################
###################################################
# Adaptive responses with risk attitudes
# first let's just increase the effectiveness
# given the times they got infected themselves
# this would be an equivalent to increase biosecurity

num_agents = 1000
global n_boot = 100
global cil = 0.95
global j
global adoptstep = (0,[0.1])
global p = plot(xlims=(0,100),ylims=(0,1),xlabel = L"t_i", ylabel = L"I/N")
begin
    for j in [0.1, 0.4, 0.7]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents];
            #  the_selected = rand(1:num_agents)
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1]);
            #  init_demographics!(agents[setdiff(1:num_agents,the_selected)];states=["S","I"],initial=[1.,0.]);
            #  init_demographics!([agents[the_selected]];states=["S","I"],initial=[0.,1.]);
            set_coop_agents!(agents;p_cop=0.5);
            set_adapt_agents!(agents;p_cop=j);
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents);
            I = Array{Float64}(undef,0);
            for i in 1:100
                next_state!(agents;fun=SI_attitude!,inf_prob=0.1, rec_prob=0.1, R=true);
                update_state!(agents);
                update_effect_given_distance!(agents,g,adoptstep[1],1,adoptstep[2]);
                get_coop!(agents);
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents);
            end
            push!(Iall,I)
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

savefig(p,"figs/sir_attitudes_1I.png")
p

###################################################
###################################################
###################################################
###################################################
# This is the same as the previous but looking at the
# mean ρ of each group. Didn't want to mess with the
# already working code.

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 100
global cil = 0.95
#  global d = 0
#  global probs = [1]
global j
#  global adoptprob = (1,[0.5,0.3])
global adoptstep = (0,[0.1])
global p = plot(xlims=(0,100),ylims=(0,1),xlabel = L"t_i", ylabel = L"I/N")
        global rtall = Array{Float64}(undef,0);
        global lrtall = Array{Float64}(undef,0);
        global raall = Array{Float64}(undef,0);
        global lraall = Array{Float64}(undef,0);
begin
    for j in [0.4]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents];
            init_demographics!(agents;states=["S","I"],initial=[0.7,0.3]);
            set_coop_agents!(agents;p_cop=0.5);
            set_adapt_agents!(agents;p_cop=j);
            map(x -> assign_contacts!(g,x), agents);
            #  get_coop!(agents);
            rt = Array{Float64}(undef,0);
            ra = Array{Float64}(undef,0);
            lrt = Array{Float64}(undef,0);
            lra = Array{Float64}(undef,0);
            for i in 1:100
                next_state!(agents;fun=SI_attitude!,inf_prob=0.1, rec_prob=0.1, R=false);
                update_state!(agents);
                update_effect_given_distance!(agents,g,adoptstep[1],1,adoptstep[2]);
                #  get_coop!(agents);
                #  agents[1].adapter = true
                #  agents[3].counter = 1
                #  agents[1].coop_effect
                #  update_single_effect_distance!(agents,g,3,0,1,adoptstep[2])
                #  agents[1].coop_effect
                push!(rt,mean([ag.coop_effect for ag in agents if ag.attitude == "rt" && !ag.adapter]));
                push!(lrt,mean([ag.coop_effect for ag in agents if ag.attitude == "rt" && ag.adapter]));
                push!(ra,mean([ag.coop_effect for ag in agents if ag.attitude == "ra" && !ag.adapter]));
                push!(lra,mean([ag.coop_effect for ag in agents if ag.attitude == "ra" && ag.adapter]));
            end
            push!(rtall,I)
            push!(lrtall,I)
            push!(raall,I)
            push!(lraall,I)
        end
        end
        #  for step in eachindex(Iall[1])
        #      bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
        #      bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
        #      bci1 = confint(bs1, BasicConfInt(cil))
        #      bci2 = confint(bs2, BasicConfInt(cil))
        #      push!(Imean,bci1[1][1])
        #      push!(Ilow,bci2[1][2])
        #      push!(Ihigh,bci2[1][3])
        #  end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
        #  p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

savefig(p,"figs/sis_attitudes_1I.png")
p

####################################
# Memory test
#
using NetABM
using LightGraphs
using StatsBase
using Plots
using Statistics
using Bootstrap
using LaTeXStrings

global num_agents = 1000
global n_boot = 100
global cil = 0.95
global j
global adoptstep = (0,[0.1])
global p = plot(xlims=(0,100),ylims=(0,1),xlabel = L"t_i", ylabel = L"I/N")

@time let
    for j in [0.0, 0.1, 0.4, 0.7, 1.0]
    #  for j in [0.7]
        Imean = Array{Float64}(undef,0)
        Ilow = Array{Float64}(undef,0)
        Ihigh = Array{Float64}(undef,0)
        Iall = Array{Array}(undef,0)
        for a in 1:100
            g = barabasi_albert(num_agents,3)
            agents = [Agent(i) for i = 1:num_agents];
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1]);
            set_coop_agents!(agents;p_cop=0.5);
            set_adapt_agents!(agents;p_cop=j);
            map(x -> assign_contacts!(g,x), agents);
            #  get_coop!(agents);
            I = Array{Float64}(undef,0);
            for i in 1:100
                next_state!(agents;fun=SI_attitude!,inf_prob=0.1, rec_prob=0.03, R=false);
                update_state!(agents);
                update_effect_given_distance!(agents,g,adoptstep[1],1,adoptstep[2]);
                #  get_coop!(agents);
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents);
                #  GC.gc()
            end
            push!(Iall,I)
            GC.gc()
        end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
            #  GC.gc()
        end
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
    GC.gc()
end
p
savefig(p,"figs/sis_attitudes_01.png")

global rtall = Array{Array}(undef,0);
global lrtall = Array{Array}(undef,0);
global raall = Array{Array}(undef,0);
global lraall = Array{Array}(undef,0);
global p = plot(xlims=(0,100),ylims=(0,1),xlabel = L"t_i", ylabel = L"I/N")
@time let
    #  for j in [0.0, 0.1, 0.4, 0.7, 1.0]
    for j in [0.99]
        Imean = Array{Float64}(undef,0)
        Ilow = Array{Float64}(undef,0)
        Ihigh = Array{Float64}(undef,0)
        Iall = Array{Array}(undef,0)
        for a in 1:100
            g = barabasi_albert(num_agents,3)
            agents = [Agent(i) for i = 1:num_agents];
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1]);
            set_coop_agents!(agents;p_cop=0.5);
            set_adapt_agents!(agents;p_cop=j);
            map(x -> assign_contacts!(g,x), agents);
            #  get_coop!(agents);
            I = Array{Float64}(undef,0);
            rt = Array{Float64}(undef,0);
            ra = Array{Float64}(undef,0);
            lrt = Array{Float64}(undef,0);
            lra = Array{Float64}(undef,0);
            for i in 1:100
                next_state!(agents;fun=SI_attitude!,inf_prob=0.1, rec_prob=0.03, R=false);
                update_state!(agents);
                update_effect_given_distance!(agents,g,adoptstep[1],1,adoptstep[2]);
                #  get_coop!(agents);
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents);
                push!(rt,mean([ag.coop_effect for ag in agents if ag.attitude == "rt" && !ag.adapter]));
                push!(lrt,mean([ag.coop_effect for ag in agents if ag.attitude == "rt" && ag.adapter]));
                push!(ra,mean([ag.coop_effect for ag in agents if ag.attitude == "ra" && !ag.adapter]));
                push!(lra,mean([ag.coop_effect for ag in agents if ag.attitude == "ra" && ag.adapter]));
            end
            push!(Iall,I)
            push!(rtall,rt)
            push!(lrtall,lrt)
            push!(raall,ra)
            push!(lraall,lra)
            GC.gc()
        end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
            #  GC.gc()
        end
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
    GC.gc()
end

Imean = Array{Float64}(undef,0)
Ilow = Array{Float64}(undef,0)
Ihigh = Array{Float64}(undef,0)
Iall = Array{Array}(undef,0)
for step in eachindex(lrtall[1])
    bs1 = bootstrap(mean,getindex.(lraall,step),BasicSampling(n_boot))
    bs2 = bootstrap(std,getindex.(lraall,step),BasicSampling(n_boot))
    bci1 = confint(bs1, BasicConfInt(cil))
    bci2 = confint(bs2, BasicConfInt(cil))
    push!(Imean,bci1[1][1])
    push!(Ilow,bci2[1][2])
    push!(Ihigh,bci2[1][3])
end
p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label="lra")
p

p
#  savefig(p,"figs/sis_attitudes_01.png")


text:
#using Pkg
#Pkg.add("Luxor")
#Pkg.add("Color")

using Luxor, Colors
Drawing(1000,1000,"turtle.png")
origin()
background("midnightblue")

turtle = Turtle()
Pencolor(turtle, "cyan")
Penwidth(turtle, 1.5)
n = 5
for i in 1:4000
    global n
    Forward(turtle, n)
    Turn(turtle, 89.5)
    HueShift(turtle)
    n += 0.75
end
fontsize(20)
#Message(turtle, "finished")
finish()

using Luxor, Colors
Drawing(1000,1000,"graph.png")
origin()
background("white")

turtle = Turtle()
Pencolour(turtle, "black")
Penwidth(turtle, 2)
Forward(turtle, 100)
Reposition(t::turtle, 0,0)
Turn(turtle, 90)
Forward(turtle, 100)
Reposition(t::turtle, 0,0)
finish()


text:
using MLJLinearModels
using BenchmarkTools, Random, LinearAlgebra
DO_COMPARISONS = false; include("../testutils.jl")

n, p = 50_000, 500
((X, y, θ), (X1, y1, θ1)) = generate_continuous(n, p;  seed=512, sparse=0.5)

# ======================== #
# ELEMENTARIES
#
# > apply_X!        ✅ Sept 13, 2019
# > apply_Xt!       ✅ Sept 13, 2019
# > get_residuals!  ✅ Sept 13, 2019
# ======================== #

# ------------
# No intercept

Xθ = similar(y)
Xtθ = similar(θ)
r = zeros(n)

# Sept 13, 2019 :: 7.83 ms (0 allocations: 0 bytes)
# May 12, 2020 :: 8.091 ms (8 allocations: 128 bytes)
#              :: 7.612 ms (8 allocations: 128 bytes)
@btime R.apply_X!($Xθ, $X, $θ);

# Sept 13, 2019 :: 7.94 ms (0 allocations: 0 bytes) [pretty much only apply_X!]
# May 12, 2020 :: 8.305 ms (12 allocations: 256 bytes)
#              :: 7.826 ms (12 allocations: 256 bytes)
@btime R.get_residuals!($r, $X, $θ, $y);

# Sept 13, 2019 :: 7.9 ms (0 allocations: 0 bytes)
# May 12, 2020 :: 8.047 ms (8 allocations: 128 bytes)
#              :: 7.616 ms (8 allocations: 128 bytes)
@btime R.apply_Xt!($Xtθ, $X, $y);

# --------------
# With intercept

Xθ = similar(y)
Xtθ = similar(θ1)

# Sept 13, 2019 :: 7.949 ms (1 allocation: 48 bytes) -- alloc for the view
# May 12, 2020 :: 8.190 ms (12 allocations: 272 bytes)
#              :: 7.810 ms (12 allocations: 272 bytes)
@btime R.apply_X!($Xθ, $X, $θ1);

# Sept 13, 2019 :: 7.955 ms (1 allocation: 48 bytes) -- mostly apply_X!
# May 12, 2020 :: 8.375 ms (16 allocations: 400 bytes)
#              :: 7.945 ms (16 allocations: 400 bytes)
@btime R.get_residuals!($r, $X, $θ1, $y);

# Sept 13, 2019 :: 7.883 ms (1 allocation: 48 bytes)
# May 12, 2020 :: 8.320 ms (9 allocations: 176 bytes)
#              :: 7.884 ms (9 allocations: 176 bytes)
@btime R.apply_Xt!($Xtθ, $X, $y);


text:
using Test, MPI
include(joinpath("..", "..", "testhelpers.jl"))

@testset "ODE Solvers" begin
    runmpi(joinpath(@__DIR__, "ode_tests_basic.jl"))
    runmpi(joinpath(@__DIR__, "callbacks.jl"))
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
export iauAticq
"""
Quick CIRS RA,Dec to ICRS astrometric place, given the star-
independent astrometry parameters.

Use of this function is appropriate when efficiency is important and
where many star positions are all to be transformed for one date.
The star-independent astrometry parameters can be obtained by
calling one of the functions iauApci[13], iauApcg[13], iauApco[13]
or iauApcs[13].

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   ri,di  double     CIRS RA,Dec (radians)
   astrom iauASTROM* star-independent astrometry parameters:
      pmt    double       PM time interval (SSB, Julian years)
      eb     double[3]    SSB to observer (vector, au)
      eh     double[3]    Sun to observer (unit vector)
      em     double       distance from Sun to observer (au)
      v      double[3]    barycentric observer velocity (vector, c)
      bm1    double       sqrt(1-|v|^2): reciprocal of Lorenz factor
      bpn    double[3][3] bias-precession-nutation matrix
      along  double       longitude + s' (radians)
      xpl    double       polar motion xp wrt local meridian (radians)
      ypl    double       polar motion yp wrt local meridian (radians)
      sphi   double       sine of geodetic latitude
      cphi   double       cosine of geodetic latitude
      diurab double       magnitude of diurnal aberration vector
      eral   double       "local" Earth rotation angle (radians)
      refa   double       refraction constant A (radians)
      refb   double       refraction constant B (radians)

Returned:
   rc,dc  double     ICRS astrometric RA,Dec (radians)

Notes:

   1. Only the Sun is taken into account in the light deflection
   correction.

   2. Iterative techniques are used for the aberration and light
   deflection corrections so that the functions iauAtic13 (or
   iauAticq) and iauAtci13 (or iauAtciq) are accurate inverses;
   even at the edge of the Sun's disk the discrepancy is only about
   1 nanoarcsecond.

Called:
   iauS2c       spherical coordinates to unit vector
   iauTrxp      product of transpose of r-matrix and p-vector
   iauZp        zero p-vector
   iauAb        stellar aberration
   iauLdsun     light deflection by the Sun
   iauC2s       p-vector to spherical
   iauAnp       normalize angle into range +/- pi

This revision:   2013 October 9

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
# void iauAticq(double ri, double di, iauASTROM *astrom,
#               double *rc, double *dc)

function iauAticq(ri::Real, di::Real, astrom::iauASTROM)
   # Allocate return values
   ref_astrom = Ref{iauASTROM}(astrom)
   ref_rc     = Ref{Float64}(0.0)
   ref_dc     = Ref{Float64}(0.0)

   status = ccall((:iauAticq, libsofa_c), Cint,
            (Cdouble, Cdouble, Ref{iauASTROM},
            Ref{Cdouble}, Ref{Cdouble}),
            convert(Float64, ri), convert(Float64, di),
            ref_astrom, ref_rc, ref_dc)

    return ref_rc[], ref_dc[]
end


text:
function findk(f, A, B, C, delta, maxk = 170)
    # Compute the total degree of the truncated bivariate Taylor expansion
    # for approximating f{A, B}(C) with accuracy delta, for triangular A and B
    #
    #--------------------------INPUT-------------------------------------------
    #
    # f			handle function for evaluating f and its derivatives: f(x, y, i, j)
    # A,B		mxm and nxn triangular matrix coefficients (likely the Schur forms)
    # C 		mxn right hand side
    # delta     upper bound for the magnitude of the residual
    # maxk		(optional) maximum value for k, default = 170
    #
    #--------------------------------------------------------------------------
    #

    (m, n) = size(C);
    lA = diag(A); lB = diag(B);
    l = sum(lA) / m; mu = sum(lB) / n;
    nrmC = norm(C);
    nrmA = norm(diag(A) .- l, Inf); # opnorm(A - l * one(A), 1);
    nrmB = norm(diag(B) .- mu, Inf); # opnorm(B - mu * one(B), 1);

    k = 0;

    ff = norm([ f(x,y,0,0) for x = lA, y = lB ], Inf)

    while k < maxk
        k = k + 1;
        y = findk_res(f, l, mu, k, nrmA, nrmB);
        if y < delta * ff
            mx = -Inf;
            for x = lA  # TODO it might be done more efficiently if f(x, y, i, j) allows vectorization
                for y = lB
                    tmp = findk_res(f, x, y, k, nrmA, nrmB);
                    if tmp > mx
                        mx = tmp;
                    end
                end
            end
            if mx < delta * ff
                break
            end
        end
    end

    return k
end

#------------------Auxiliary function for findk --------------------------------
function findk_res(f, x, y, k, nrmA, nrmB)
    yy = 0.0;

    fj = 1.0;
    fk = convert(Float64, factorial(big(k+1)));

    for j = 0 : k + 1
        # yy = yy + abs(f(x, y, k + 1 - j, j)) / factorial(k + 1 - j) * nrmA^(k + 1 - j) / factorial(j) * nrmB^j;
        yy = yy + abs(f(x, y, k + 1 - j, j)) / fk * nrmA^(k + 1 - j) / fj * nrmB^j;

        # Update factorials for the next round
        fj = (j+1) * fj;
        fk = fk / (k+1+j);
    end

    return yy;
end

function taylor_eval(f, A, B, C, delta)
    # Returns [X, digits]
    # Evaluate the truncated bivariate taylor expansion of accuracy delta of f{A, B}(C)
    # for triangular coefficients A and B
    #
    #--------------------------INPUT-------------------------------------------
    #
    # f			handle function for evaluating f and its derivatives: f(x, y, i, j)
    # A,B		mxm and nxn triangular matrix coefficients (likely the Schur forms)
    # C 		mxn right hand side
    # delta     upper bound for the magnitude of the residual
    #
    #--------------------------------------------------------------------------
    #

    (m, n) = size(C);
    l = sum(diag(A)) / m; mu = sum(diag(B)) / n;
    k = findk(f, A, B, C, delta);
    digits = k;
    A = A - l * one(A);
    B = B - mu * one(B);

    coef = zeros(ComplexF64, k, k)
    for i = 1 : k
        for j = 1 : k
            if i + j <= k + 1
                coef[i, j] = f(l, mu, i - 1, j - 1)
            end
        end
    end

    # XX at the generic step is A^i / i! * C
    XX = C
    X = zero(C)

    for i = 1 : k

        YY = XX
        for j = 1 : k + 1 - i
            X = X + coef[i,j] * YY
            YY = YY * B / j
        end

        XX = A * XX / i
    end

    return X, digits;
end

max_stars_repo_name: numpi/BivMatFun

max_issues_repo_name: numpi/BivMatFun

max_forks_repo_name: numpi/BivMatFun


text:
# Copyright(c) Rafael Dias
# MIT Licensed

module PropriedadesSecoes
    mutable struct Coordenada
        x::Float32
        y::Float32
    end

    function comprimentoEntrePontos(coordenadas::Array{Coordenada})::Float32
        quantidade = length(coordenadas)
        atual = 1
        proximo = 2
        comprimento = 0f0

        for contador in 1:quantidade - 1
            comprimento += sqrt(((coordenadas[proximo].x - coordenadas[atual].x)^2) + ((coordenadas[proximo].y - coordenadas[atual].y)^2))

            if proximo < quantidade
                proximo += 1
            else
                proximo = 1
            end
            if contador == quantidade
                atual = 1
            else
                atual += 1
            end
        end

        return comprimento
    end

    function areaPoligono(coordenadas::Array{Coordenada})::Float32
        quantidade = length(coordenadas)
        atual = 1
        proximo = 2
        area = 0f0

        for contador in 1:quantidade
            area += (coordenadas[proximo].x * coordenadas[atual].y) - (coordenadas[atual].x * coordenadas[proximo].y)

            if proximo < quantidade
                proximo += 1
            else
                proximo = 1
            end
            if contador == quantidade
                atual = 1
            else
                atual += 1
            end
        end

        return area / 2f0
    end

    function perimetroPoligono(coordenadas::Array{Coordenada})::Float32
        quantidade = length(coordenadas)
        atual = 1
        proximo = 2
        perimetro = 0f0

        for contador in 1:quantidade
            perimetro += sqrt(((coordenadas[proximo].x - coordenadas[atual].x)^2) + ((coordenadas[proximo].y - coordenadas[atual].y)^2))

            if proximo < quantidade
                proximo += 1
            else
                proximo = 1
            end
            if contador == quantidade
                atual = 1
            else
                atual += 1
            end
        end

        return perimetro
    end

    function centroGravidadePoligonoX(coordenadas::Array{Coordenada})::Float32
        quantidade = length(coordenadas)
        atual = 1
        proximo = 2
        xgc = 0f0

        for contador in 1:quantidade
            xgc += (coordenadas[proximo].x + coordenadas[atual].x) * ((coordenadas[proximo].x * coordenadas[atual].y) - (coordenadas[atual].x * coordenadas[proximo].y))

            if proximo < quantidade
                proximo += 1
            else
                proximo = 1
            end
            if contador == quantidade
                atual = 1
            else
                atual += 1
            end
        end

        return xgc * (1 / (6 * areaPoligono(coordenadas)))
    end

    function centroGravidadePoligonoY(coordenadas::Array{Coordenada})::Float32
        quantidade = length(coordenadas)
        atual = 1
        proximo = 2
        ygc = 0f0

        for contador in 1:quantidade
            ygc += (coordenadas[proximo].y + coordenadas[atual].y) * ((coordenadas[proximo].x * coordenadas[atual].y) - (coordenadas[atual].x * coordenadas[proximo].y))

            if proximo < quantidade
                proximo += 1
            else
                proximo = 1
            end
            if contador == quantidade
                atual = 1
            else
                atual += 1
            end
        end

        return ygc * (1 / (6 * areaPoligono(coordenadas)))
    end
end

max_stars_repo_name: rsd-consultoria/scripts-engenharia

max_issues_repo_name: rsd-consultoria/scripts-engenharia

max_forks_repo_name: rsd-consultoria/scripts-engenharia


text:
function estimate(entity::HomographyMatrix, method::DirectLinearTransform, 𝒟::Tuple{AbstractArray, Vararg{AbstractArray}})
    ℳ, ℳʹ =  𝒟
    N = length(ℳ)
    if (N != length(ℳʹ))
          throw(ArgumentError("There should be an equal number of points for each view."))
    end
    𝒪, 𝐓 = transform(HomogeneousCoordinates(),CanonicalToHartley(),ℳ)
    𝒪ʹ, 𝐓ʹ  = transform(HomogeneousCoordinates(),CanonicalToHartley(),ℳʹ)
    𝐀 = moments(HomographyMatrix(), (𝒪, 𝒪ʹ))
    λ, h = smallest_eigenpair(Symmetric(𝐀))
    𝐇 = reshape(h,(3,3))
    𝐇 = SMatrix{3,3,Float64,9}(𝐇 / norm(𝐇))
    # Transform estimate back to the original (unnormalised) coordinate system.
    inv(𝐓ʹ)*𝐇*𝐓
end


text:
function get_kernel_argdiffs(argdiffs::Tuple)
    n_args = length(argdiffs)
    kernel_argdiffs = Dict{Int,Vector}()
    for (i, diff) in enumerate(argdiffs)
        if isa(diff, VectorDiff)
            for (key, element_diff) in diff.updated
                if !haskey(kernel_argdiffs, key)
                    kernel_argdiff = Vector{Any}(undef, n_args)
                    fill!(kernel_argdiff, NoChange())
                    kernel_argdiffs[key] = kernel_argdiff
                end
                kernel_argdiffs[key][i] = diff.updated[key]
            end
        end
    end
    kernel_argdiffs
end

function process_all_retained!(gen_fn::Map{T,U}, args::Tuple, argdiffs::Tuple,
                               choices_or_selection, prev_length::Int, new_length::Int,
                               retained_and_targeted::Set{Int}, state) where {T,U}
    kernel_no_change_argdiffs = map((_) -> NoChange(), args)
    kernel_unknown_change_argdiffs = map((_) -> UnknownChange(), args)

    if all(diff == NoChange() for diff in argdiffs)

        # only visit retained applications that were targeted
        for key in retained_and_targeted
            @assert key <= min(new_length, prev_length)
            process_retained!(gen_fn, args, choices_or_selection, key, kernel_no_change_argdiffs, state)
        end

    elseif any(diff == UnknownChange() for diff in argdiffs)

        # visit every retained application
        for key in 1:min(prev_length, new_length)
            @assert key <= min(new_length, prev_length)
            process_retained!(gen_fn, args, choices_or_selection, key, kernel_unknown_change_argdiffs, state)
        end

    else

        key_to_kernel_argdiffs = get_kernel_argdiffs(argdiffs)

        # visit every retained applications that either has an argdiff or was targeted
        for key in union(keys(key_to_kernel_argdiffs), retained_and_targeted)
            @assert key <= min(new_length, prev_length)
            if haskey(key_to_kernel_argdiffs, key)
                kernel_argdiffs = tuple(key_to_kernel_argdiffs[key]...)
            else
                kernel_argdiffs = kernel_no_change_argdiffs
            end
            process_retained!(gen_fn, args, choices_or_selection, key, kernel_argdiffs, state)
        end

    end
end

"""
Process all new applications.
"""
function process_all_new!(gen_fn::Map{T,U}, args::Tuple, choices_or_selection,
                          prev_len::Int, new_len::Int, state) where {T,U}
    for key=prev_len+1:new_len
        process_new!(gen_fn, args, choices_or_selection, key, state)
    end
end


text:
#=
Copyright (c) 2016, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module HPAT

ENV["CGEN_MPI_COMPILE"]=1

using CompilerTools.DebugMsg
DebugMsg.init()

using CompilerTools
using CompilerTools.AstWalker
using CompilerTools.LambdaHandling
using CompilerTools.Helper
using CompilerTools.OptFramework
using CompilerTools.OptFramework.OptPass
using ParallelAccelerator
using ParallelAccelerator.Driver
using ParallelAccelerator.CGen
import ParallelAccelerator.API.runStencil

import MPI

export hpat, hpat_checkpoint, hpat_debug, @acc, @noacc, runStencil
using Base.typed_hcat
export typed_hcat

ParallelAccelerator.CGen.setCreateMain(true)

# disable OMP if not set by user since it is slower than pure MPI
if !haskey(ENV, "CGEN_NO_OMP")
    ParallelAccelerator.CGen.disableOMP()
end

function enableOMP()
    if !haskey(ENV, "CGEN_NO_OMP")
        ParallelAccelerator.CGen.enableOMP()
    end
end

OptFramework.setSaveOriginalFunction(false)

# smaller value means higher precedence
@enum Partitioning SEQ=1 ONE_D_VAR=2 TWO_D=3 ONE_D=4

# save distributed arrays of last compiled function for debugging
saved_array_partitionings = Dict{LHSVar,Partitioning}()
saved_linfo = nothing

function set_saved_array_partitionings(d, l)
    global saved_array_partitionings = d
    global saved_linfo = l
end

function get_saved_array_partitionings()
    return saved_array_partitionings
end

function get_saved_partitioning_for_symbol(name :: Symbol)
    lv = CompilerTools.LambdaHandling.toLHSVar(name, saved_linfo)
    return saved_array_partitionings[lv]
end

force_parallel = true

function setForceParallel(v::Bool)
    global force_parallel=v
end

# constant block size for 2D partitioning
BLOCK_SIZE = 128

function setBlockSize(v::Int)
    global BLOCK_SIZE=v
end

include("api.jl")
using HPAT.API
export data_source_HDF5, data_source_TXT
include("checkpoint.jl")
include("capture-api.jl")
include("domain-pass.jl")
include("datatable-pass.jl")
include("distributed-pass.jl")
include("cgen-hpat-pattern-match.jl")

# overwrite print so only rank 0 prints
#=function Base.print(xs...)
rank = MPI.Comm_rank(MPI.COMM_WORLD)
if rank==0 Base.print(STDOUT, xs) end
end

function Base.println(xs...)
rank = MPI.Comm_rank(MPI.COMM_WORLD)
if rank==0 Base.println(STDOUT, xs) end
end=#

# add HPAT pattern matching code generators to CGen
ParallelAccelerator.CGen.setExternalPatternMatchCall(CGenPatternMatch.pattern_match_call)
ParallelAccelerator.CGen.setExternalPatternMatchAssignment(CGenPatternMatch.from_assignment_match_dist)

HPAT_PATH = joinpath(dirname(@__FILE__), "..")

HPAT_INCLUDES = string("#include <ctime>\n#include <stdlib.h>\n#include \"", HPAT_PATH, "/deps/include/hpat.h\"\n")
ParallelAccelerator.CGen.addCgenUserOptions(ParallelAccelerator.CGen.CgenUserOptions(HPAT_INCLUDES,"",""))

function addHpatInclude(stmts, compileFlags, linkFlags)
    ParallelAccelerator.CGen.addCgenUserOptions(ParallelAccelerator.CGen.CgenUserOptions(stmts, compileFlags, linkFlags))
end

function ns_to_sec(x)
    x / 1000000000.0
end

function runDistributedPass(func :: GlobalRef, ast, signature :: Tuple)
    dir_start = time_ns()
    code = DistributedPass.from_root(string(func.name), ast)
    dir_time = time_ns() - dir_start
    @dprintln(2, "Distributed code = ", code)
    @dprintln(1, "accelerate: DistributedPass conversion time = ", ns_to_sec(dir_time))
    return code
end

function runDataTablePass(func :: GlobalRef, ast, signature :: Tuple)
    dir_start = time_ns()
    code = DataTablePass.from_root(string(func.name), ast)
    dir_time = time_ns() - dir_start
    @dprintln(2, "DataTable code = ", code)
    @dprintln(1, "accelerate: DataTablePass conversion time = ", ns_to_sec(dir_time))
    return code
end

function runDomainPass(func :: GlobalRef, ast, signature :: Tuple)
    dir_start = time_ns()
    code = DomainPass.from_root(string(func.name), ast)
    dir_time = time_ns() - dir_start
    @dprintln(2, "Domain code = ", code)
    @dprintln(1, "accelerate: DomainPass conversion time = ", ns_to_sec(dir_time))
    return code
end

function addCheckpointing(func :: GlobalRef, ast, signature :: Tuple)
    dir_start = time_ns()
    code = Checkpointing.from_root(string(func.name), ast, false)
    dir_time = time_ns() - dir_start
    @dprintln(2, "Added checkpointing = ", code)
    @dprintln(1, "accelerate: addCheckpointing conversion time = ", ns_to_sec(dir_time))
    return code
end

function addCheckpointingRestart(func :: GlobalRef, ast, signature :: Tuple)
    dir_start = time_ns()
    code = Checkpointing.from_root(string(func.name), ast, true)
    dir_time = time_ns() - dir_start
    @dprintln(2, "Added checkpointing and restart = ", code)
    @dprintln(1, "accelerate: addCheckpointingRestart conversion time = ", ns_to_sec(dir_time))
    return code
end

type MacroState
    # columns of each table
    tableCols::Dict{Symbol,Vector{Symbol}}
    # Types of each table column
    tableTypes::Dict{Symbol,Vector{Symbol}}
    # mapping table number to table name
    # domain-pass uses this to find tables of relational calls
    tableIds::Dict{Int,Symbol}
    # partitioning of arrays
    array_partitioning::Dict{Symbol,Partitioning}
    # unique id used in renaming
    unique_id::Int
    # symbol rename map for cases like t1=t1[cond] where output table needs to
    #  be rename for subsequent AST nodes.
    rename_map::Dict{Symbol,Symbol}

    function MacroState()
        new(Dict{Symbol,Vector{Symbol}}(), Dict{Symbol,Vector{Symbol}}(),
         Dict{Int,Symbol}(), Dict{Symbol,Partitioning}(), 0, Dict{Symbol,Symbol}())
    end
end

"""
    A macro pass that translates extensions such as DataSource()
"""
function captureHPAT(func, ast, sig)
    macro_state = MacroState()
    AstWalk(ast, CaptureAPI.process_node, macro_state)
    Base.pushmeta!(ast,:hpat_tables, macro_state.tableCols, macro_state.tableTypes, macro_state.tableIds)
    Base.pushmeta!(ast,:hpat_partitioning, macro_state.array_partitioning)
    return ast
end

"""
    Adds a duplicate of the function to be checkpointed with "`_`restart" appended to the name.
"""
function createCheckpointFunc(func, ast, sig)
    @dprintln(1, "createCheckpointFunc func = ", func)
    new_func = deepcopy(ast)
    new_func.args[1].args[1] = symbol(string(new_func.args[1].args[1],"_restart"))
    @dprintln(1, "createCheckpointFunc new_func = ", new_func.args[1].args[1])
    return CompilerTools.OptFramework.MoreWork(ast, CompilerTools.OptFramework.WorkItem[CompilerTools.OptFramework.WorkItem(hpat_checkpoint_internal, hpat_checkpoint_internal, Any[], new_func)])
end

const hpat_debug = [ OptPass(captureHPAT, PASS_MACRO) ]

# initialize set of compiler passes HPAT runs
const hpat =
    [ OptPass(captureHPAT, PASS_MACRO),
      OptPass(captureOperators, PASS_MACRO),
      OptPass(toCartesianArray, PASS_MACRO),
      OptPass(expandParMacro, PASS_MACRO),
      OptPass(runDomainPass, PASS_TYPED),
      OptPass(extractCallGraph, PASS_TYPED),
      OptPass(toDomainIR, PASS_TYPED),
      OptPass(runDataTablePass, PASS_TYPED),
      OptPass(toParallelIR, PASS_TYPED),
      OptPass(runDistributedPass, PASS_TYPED),
      OptPass(toFlatParfors, PASS_TYPED),
      OptPass(toCGen, PASS_TYPED) ]

const hpat_checkpoint =
    [ OptPass(createCheckpointFunc, PASS_MACRO),
      OptPass(captureHPAT, PASS_MACRO),
      OptPass(captureOperators, PASS_MACRO),
      OptPass(toCartesianArray, PASS_MACRO),
      OptPass(expandParMacro, PASS_MACRO),
      OptPass(runDomainPass, PASS_TYPED),
      OptPass(addCheckpointing, PASS_TYPED),
      OptPass(extractCallGraph, PASS_TYPED),
      OptPass(toDomainIR, PASS_TYPED),
      OptPass(runDataTablePass, PASS_TYPED),
      OptPass(toParallelIR, PASS_TYPED),
      OptPass(runDistributedPass, PASS_TYPED),
      OptPass(toFlatParfors, PASS_TYPED),
      OptPass(toCGen, PASS_TYPED) ]

const hpat_checkpoint_internal =
    [ OptPass(captureHPAT, PASS_MACRO),
      OptPass(captureOperators, PASS_MACRO),
      OptPass(toCartesianArray, PASS_MACRO),
      OptPass(expandParMacro, PASS_MACRO),
      OptPass(runDomainPass, PASS_TYPED),
      OptPass(addCheckpointingRestart, PASS_TYPED),
      OptPass(extractCallGraph, PASS_TYPED),
      OptPass(toDomainIR, PASS_TYPED),
      OptPass(runDataTablePass, PASS_TYPED),
      OptPass(toParallelIR, PASS_TYPED),
      OptPass(runDistributedPass, PASS_TYPED),
      OptPass(toFlatParfors, PASS_TYPED),
      OptPass(toCGen, PASS_TYPED) ]

append!(ParallelAccelerator.DomainIR.funcIgnoreList, map(x->GlobalRef(HPAT.API,x),[DomainPass.generatedFuncs; API.operators]))
append!(ParallelAccelerator.DomainIR.exprHeadIgnoreList, DomainPass.generatedExprHeads)
ParallelAccelerator.DomainIR.setExternalCallback(DomainPass.AstWalkCallback)
ParallelAccelerator.DomainIR.setExternalLiveCB(DomainPass.live_cb)
ParallelAccelerator.DomainIR.setExternalAliasCB(DomainPass.alias_cb)

HPAT_DEFAULT_DATAPATH = "input_data/"

if haskey(ENV, "HPAT_DEFAULT_DATA")
    global HPAT_DEFAULT_DATAPATH = ENV["HPAT_DEFAULT_DATA"]
elseif haskey(ENV, "SCRATCH")
    global HPAT_DEFAULT_DATAPATH = ENV["SCRATCH"]*"/input_data/"
else
    global HPAT_DEFAULT_DATAPATH = joinpath(dirname(@__FILE__), "..")*"/input_data/"
end

MPI.Comm_rank(MPI.COMM_WORLD)==0 && !isdir(HPAT_DEFAULT_DATAPATH) &&  mkdir(HPAT_DEFAULT_DATAPATH)

function getDefaultDataPath()
    return HPAT_DEFAULT_DATAPATH
end

function setDefaultDataPath(path::AbstractString)
    HPAT_DEFAULT_DATAPATH = path
end

function HPAT_finalize()
    MPI.Barrier(MPI.COMM_WORLD)
    MPI.Finalize()
end

atexit(HPAT_finalize)

function restart(func, args...)
    @dprintln(1, "HPAT restart func = ", func, " type = ", typeof(func), " args = ", args...)
    arg_type_tuple_expr = Expr(:tuple)
    arg_type_tuple_expr.args = map(x -> typeof(x), [args...])
    arg_type_tuple = eval(arg_type_tuple_expr)
    @dprintln(2, "arg_type_tuple = ", arg_type_tuple)
    gr = GlobalRef(Base.function_module(func,arg_type_tuple), symbol(func))
    @dprintln(2, "GlobalRef = ", gr)
    if haskey(CompilerTools.OptFramework.gOptFrameworkDict, gr)
        res = CompilerTools.OptFramework.gOptFrameworkDict[gr]
        @dprintln(2, "CompilerTools.OptFramework internally mapped ", gr, " to ", res)
        assert(typeof(res) == GlobalRef)
        new_gr = GlobalRef(res.mod, symbol(string(res.name,"_restart")))
        @dprintln(2, "Preparing to call the restart version of the function with name = ", new_gr)
        eval(new_gr)(args...)
    else
        throw(string("HPAT restart function could not find function to restart."))
    end
end

@noinline function __hpat_dist_broadcast(variable, size)
    convert(Int32, 1)
end

end # module


text:
# Space war game!
# A thing is anything that exist in the universe.
# Concrete type of Thing should always have the following fields:
#       1. position
#       2. size
abstract type Thing end

# Stores the coordinates of a game piece.
mutable struct Position # mutable because we expect the pieces to move around
    x::Int
    y::Int
end

# Stores the size of a game piece
struct Size
    width::Int
    height::Int
end

# Holds a single game piece
struct Widget
    name::String
    position::Position
    size::Size
end

# Functions that are applied for all Thing's
position(t::Thing) = t.position
size(t::Thing) = t.size
shape(t::Thing) = :unknown

# Type of weapons
@enum Weapon Laser Missile Gun

# Spaceship
struct Spaceship <: Thing
    position::Position
    size::Size
    weapon::Weapon
end
shape(s::Spaceship) = :saucer

struct Asteroid <: Thing
    position::Position
    size::Size
end

include("shape.jl")

# single-line functions to coordinate moving of the game pieces
move_up!(widget::Widget, to_move::Int) = widget.position.y -= to_move
move_down!(widget::Widget, to_move::Int) = widget.position.y += to_move
move_left!(widget::Widget, to_move::Int) = widget.position.x -= to_move
move_right!(widget::Widget, to_move::Int) = widget.position.x += to_move

function random_move()
    return rand([move_up!, move_down!, move_left!, move_right!])
end

function random_leap(w::Widget, move_func::Function, distance::Int)
    move_func(w, distance)
    return w
end

function clean_up_galaxy(asteroids, spaceships)
    ep = x -> println(x, " exploded!")
    foreach(ep, asteroids)
    foreach(ep, spaceships)
end

# Define pretty print functions
Base.show(io::IO, p::Position) = print(io, "(", p.x, ",", p.y, ")")
Base.show(io::IO, s::Size) = print(io, s.width, " x ", s.height)
Base.show(io::IO, w::Widget) = print(io, w.name, " at ", w.position, " size ", w.size)

function shoot(from::Widget, targets::Widget...)
    println("Type of targets: ", typeof(targets))
    for target in targets
        println(from.name, " --> ", target.name)
    end
end

# Function to make new asteroids
function make_asteroids(N::Int; pos_range=0:200, size_range=0:10:30)
    pos_rand() = rand(pos_range)
    sz_rand() = rand(size_range)
    return [Widget("Asteroid #$i", Position(pos_rand(), pos_rand()), Size(sz_rand(), sz_rand())) for i in 1:N]
end

# Function to make new spaceships
function make_spaceships(N::Int; pos_range=0:200, size_range=0:10:30)
    pos_rand() = rand(pos_range)
    sz_rand() = rand(size_range)
    return [Widget("Spaceship #$i", Position(pos_rand(), pos_rand()), Size(sz_rand(), sz_rand())) for i in 1:N]
end

# Special arrangement before attacks
function triangular_formation!(s1::Widget, s2::Widget, s3::Widget)
    x_offset = 30
    y_offset = 50
    s2.position.x = s1.position.x - x_offset
    s3.position.x =  s1.position.x + x_offset
    s2.position.y = s3.position.y = s1.position.y - y_offset
    # (s1, s2, s3)
end

# Random healthiness function for testing purposes
healthy(spaceship) = rand(Bool)

# make sure that the spaceship is healthy before any operation
function fire(f::Function, spaceship::Widget)
    if healthy(spaceship)
        f(spaceship)
    else
        println("Operation aborted as spaceship is not healthy")
    end
    return nothing
end

function process_file(func::Function, filename::AbstractString)
    ios = nothing
    try
        ios = open(filename)
        func(ios)
    finally
        close(ios)
    end
end

# explode an array of abjects function explode(things::AbstractVector{AnyObject})
function explode(things::AbstractVector{T}) where {T <: Thing}
    for t in things
        println("Exploding ", t)
    end
end

function explode(things::AbstractVector{T}) where {T}
    for t in things
        println("Exploding ", t)
    end
end

function tow(A::Spaceship, B::Thing)
    "tow 1"
end

# equivalent of parametric type
function tow(A::Spaceship, B::T) where {T <: Thing}
    "tow 2"
end

function group_anything(A::T, B::T) where {T <: Thing}
    println("Grouped ", A, " and ", B)
end


text:
using Oceananigans

# We'll want the make_vertical_slice_movie function from utils.jl to turn our
# rising thermal bubble output into a movie.
include("utils.jl")

# We'll set up a 2D model with an xz-slice so there's only 1 grid point in y.
Nx, Ny, Nz = 128, 128, 128     # Number of grid points in each dimension.
Lx, Ly, Lz = 2000, 2000, 2000  # Domain size (meters).
Nt, Δt = 5000, 10              # Number of time steps, time step size (seconds).

# Set up the model and use an artificially high viscosity ν and diffusivity κ.
model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), arch=GPU(), ν=4e-2, κ=4e-2)

# Get location of the cell centers in x, y, z and reshape them to easily
# broadcast over them when calculating hot_bubble_perturbation.
xC, yC, zC = model.grid.xC, model.grid.yC, model.grid.zC
xC, yC, zC = reshape(xC, (Nx, 1, 1)), reshape(yC, (Ny, 1, 1)), reshape(zC, (1, 1, Nz))

# Set a temperature perturbation with a Gaussian profile located at the center
# of the vertical slice. It roughly corresponds to a background temperature of
# T = 282.99 K and a bubble temperature of T = 283.01 K.
hot_bubble_perturbation = @. 0.01 * exp(-100 * ((xC - Lx/2)^2 + (yC - Ly/2)^2 + (zC + Lz/2)^2) / (Lx^2 + Ly^2 + Lz^2))
model.tracers.T.data .= 282.99 .+ 2 .* CuArray(hot_bubble_perturbation)

# Add a NetCDF output writer that saves NetCDF files to the current directory
# "." with a filename prefix of "thermal_bubble_3d_" every 10 iterations.
nc_writer = NetCDFOutputWriter(dir=".", prefix="thermal_bubble_3d_", frequency=25)
push!(model.output_writers, nc_writer)

time_step!(model, Nt, Δt)

make_vertical_slice_movie(model, nc_writer, "T", Nt, Δt, model.eos.T₀)


text:
"""
Overview: set of constants used over the project.
This helps minimizing literals throughout the work.
If you need to insert constants used througout the work,
then please add them here to help minimizing unmeaningful
numbers in the work.
Furthermore, add
"""

"""
Dimension of the treated problems (3D)
"""
const dim = 3

"""
Number of nodes of a vertex (1)
"""
const VE1n = 1

"""
Number of nodes of a linear line element (2)
"""
const BE2n = 2

"""
Number of nodes of a quadratic line element (3)
"""
const BE3n = 3

"""
Number of nodes of a linear triangular element (3)
"""
const TR3n = 3

"""
Number of nodes of a quadratic triangular element (6)
"""
const TR6n = 6

"""
Number of nodes of a linear tetrahedral element (4)
"""
const TE4n = 4

"""
Number of nodes of a quadratic tetrahedral element (10)
"""
const T10n = 10

"""
Number of nodes of a linear wedge element (6)
"""
const PE6n = 6

"""
Number of nodes of a quadratic 15-nodes wedge element (15)
"""
const P15n = 15

"""
Number of nodes of a quadratic 18-nodes wedge element (18)
"""
const P18n = 18

"""
Number of nodes of a linear square element (4)
"""
const QU4n = 4

"""
Number of nodes of a quadratic serendipity square element (8)
"""
const QU8n = 8

"""
Number of nodes of a quadratic Lagrange square element (9)
"""
const QU9n = 9

"""
Number of nodes of a linear hexahedral element (8)
"""
const HE8n = 8

"""
Number of nodes of a quadratic serendipity hexahedral element (20)
"""
const H20n = 20

"""
Number of nodes of a quadratic Lagrange hexahedral element (27)
"""
const H27n = 27

"""
Number of supported boundary element
"""
const Snse = 5

"""
Number of supported volume element
"""
const Vnse = 8

"""
Maximum number of elements connected one to the other
"""
const ncv_max = 30

"""
Maximum number nodes for an element
"""
const nne_max = 27

"""
resonances tolerance
"""
const ϵ_tol = 1e-1


text:
# TESTING THE BASIC METHODS FOR HdrImages #######################################################################################

img = HdrImage(3, 2)

@testset "HdrImages: Basic methods" begin
    @test img.width == 3
    @test img.height == 2
    @test valid_coordinates(img, 2, 1) == true
    @test valid_coordinates(img, 3, 2) == true #the last element of the matrix has to be (3,2) and not (2,1) as in C++
    @test valid_coordinates(img,-1, 4) == false
    @test valid_coordinates(img, 0, 4) == false
    @test pixel_offset(img,1,1) == 1
    @test pixel_offset(img,2,2) == 5
    @test get_pixel(img, 1, 2) == 4
    @test get_pixel(img, 3, 1) == 3

new_col = RGB(.012, 34, 999.9)
set_pixel(img, 1, 1, new_col)
set_pixel(img, 3, 2, new_col)

    @test img.pixels[1] == new_col
    @test img.pixels[get_pixel(img, 3, 2)] == new_col

end

# TESTING THE WRITING/READING METHODS FOR HdrImages ############################################################################

img2 = HdrImage(3, 2)

LE_REFERENCE_BYTES = [
    0x50, 0x46, 0x0a, 0x33, 0x20, 0x32, 0x0a, 0x2d, 0x31, 0x2e, 0x30, 0x0a,
    0x00, 0x00, 0xc8, 0x42, 0x00, 0x00, 0x48, 0x43, 0x00, 0x00, 0x96, 0x43,
    0x00, 0x00, 0xc8, 0x43, 0x00, 0x00, 0xfa, 0x43, 0x00, 0x00, 0x16, 0x44,
    0x00, 0x00, 0x2f, 0x44, 0x00, 0x00, 0x48, 0x44, 0x00, 0x00, 0x61, 0x44,
    0x00, 0x00, 0x20, 0x41, 0x00, 0x00, 0xa0, 0x41, 0x00, 0x00, 0xf0, 0x41,
    0x00, 0x00, 0x20, 0x42, 0x00, 0x00, 0x48, 0x42, 0x00, 0x00, 0x70, 0x42,
    0x00, 0x00, 0x8c, 0x42, 0x00, 0x00, 0xa0, 0x42, 0x00, 0x00, 0xb4, 0x42
]

BE_REFERENCE_BYTES = [
    0x50, 0x46, 0x0a, 0x33, 0x20, 0x32, 0x0a, 0x31, 0x2e, 0x30, 0x0a, 0x42,
    0xc8, 0x00, 0x00, 0x43, 0x48, 0x00, 0x00, 0x43, 0x96, 0x00, 0x00, 0x43,
    0xc8, 0x00, 0x00, 0x43, 0xfa, 0x00, 0x00, 0x44, 0x16, 0x00, 0x00, 0x44,
    0x2f, 0x00, 0x00, 0x44, 0x48, 0x00, 0x00, 0x44, 0x61, 0x00, 0x00, 0x41,
    0x20, 0x00, 0x00, 0x41, 0xa0, 0x00, 0x00, 0x41, 0xf0, 0x00, 0x00, 0x42,
    0x20, 0x00, 0x00, 0x42, 0x48, 0x00, 0x00, 0x42, 0x70, 0x00, 0x00, 0x42,
    0x8c, 0x00, 0x00, 0x42, 0xa0, 0x00, 0x00, 0x42, 0xb4, 0x00, 0x00
]
# Set pixels colors like the "reference_le.pfm" file
set_pixel(img2,1, 1, RGB(1.0e1, 2.0e1, 3.0e1))
set_pixel(img2,2, 1, RGB(4.0e1, 5.0e1, 6.0e1))
set_pixel(img2,3, 1, RGB(7.0e1, 8.0e1, 9.0e1))
set_pixel(img2,1, 2, RGB(1.0e2, 2.0e2, 3.0e2))
set_pixel(img2,2, 2, RGB(4.0e2, 5.0e2, 6.0e2))
set_pixel(img2,3, 2, RGB(7.0e2, 8.0e2, 9.0e2))

# Test the write function
buf = IOBuffer()
write(buf,img2)

@testset "HdrImages: Write" begin
    @test take!(buf) == LE_REFERENCE_BYTES
end

# Test the read function
line = IOBuffer(b"Hello\nWorld!")

@testset "HdrImages: Read" begin

    @test _read_line(line) == "Hello"
    @test _read_line(line) == "World!"
    @test _read_line(line) == ""

    @test _parse_img_size("3 2") == (3, 2)
    @test_throws InvalidPfmFileFormat _parse_img_size("-1 3")
    @test_throws InvalidPfmFileFormat _parse_img_size("1 2 3")

    @test _parse_endianness("1.0") == "BE"
    @test _parse_endianness("-1.0") == "LE"
    @test_throws InvalidPfmFileFormat _parse_endianness("abc")
    @test_throws InvalidPfmFileFormat _parse_endianness("2.0")

for reference_bytes in [ BE_REFERENCE_BYTES, LE_REFERENCE_BYTES ]
    img2_1 = read_pfm_image(IOBuffer(reference_bytes))

    @test img2.pixels[get_pixel(img2_1, 1, 1)] ≈ RGB(1.0e1, 2.0e1, 3.0e1)
    @test img2.pixels[get_pixel(img2_1, 2, 1)] ≈ RGB(4.0e1, 5.0e1, 6.0e1)
    @test img2.pixels[get_pixel(img2_1, 3, 1)] ≈ RGB(7.0e1, 8.0e1, 9.0e1)
    @test img2.pixels[get_pixel(img2_1, 1, 2)] ≈ RGB(1.0e2, 2.0e2, 3.0e2)
    @test img2.pixels[get_pixel(img2_1, 1, 1)] ≈ RGB(1.0e1, 2.0e1, 3.0e1)
    @test img2.pixels[get_pixel(img2_1, 2, 2)] ≈ RGB(4.0e2, 5.0e2, 6.0e2)
    @test img2.pixels[get_pixel(img2_1, 3, 2)] ≈ RGB(7.0e2, 8.0e2, 9.0e2)
end

buf = IOBuffer(transcode(UInt8, "PF\n3 2\n-1.0\nstop"))
    @test_throws InvalidPfmFileFormat read_pfm_image(buf)
end

# TESTING THE SAVING METHODS FOR LdrImages ############################################################################

img3 = HdrImage(2,1)

set_pixel(img3, 1, 1, RGB(5.0, 10.0, 15.0))  # Luminosity: 10.0
set_pixel(img3, 2, 1, RGB(500.0, 1000.0, 1500.0))  # Luminosity: 1000.0

@testset "LdrImages: Saving methods" begin
    @test 100 ≈ average_luminosity(img3, 0.0)
    normalize_image(img3, 1000.0, 100.0)
    @test isapprox(img3.pixels[get_pixel(img3, 1, 1)], RGB(0.5e2, 1.0e2, 1.5e2))
    @test isapprox(img3.pixels[get_pixel(img3, 2, 1)], RGB(0.5e4, 1.0e4, 1.5e4))

    clamp_image(img3)
    for cur_pixel in img3.pixels
        @test (cur_pixel.r >= 0) && (cur_pixel.r <= 1)
        @test (cur_pixel.g >= 0) && (cur_pixel.g <= 1)
        @test (cur_pixel.b >= 0) && (cur_pixel.b <= 1)
    end
end


text:
using UNet: minsize, outputsize

# Generate a fake batch of two images and return it with the expected ouput size.
function dummy_data(;inchannels,
    nclasses = 1,
    volume = false,
    padding = false,
    nlevels = 1)
    sz = minsize(padding, nlevels)
    is = volume ? (sz, sz, sz) : (sz, sz)
    os = UNet.outputsize(is, padding, nlevels)
    rand(Float32, (is..., inchannels, 2)), (os..., nclasses, 2)
end

@testset "unet" begin
    for ic ∈ [1 3],
        nc ∈ [1 5],
        vl ∈ [true false],
        bn ∈ [true false],
        pg ∈ [true false],
        up ∈ [:convt :nearest :bilinear]

        if vl && up == :bilinear
            continue
        end
        model = unet(inchannels = ic, nclasses = nc, volume = vl, basewidth = 4,
            batchnorm = bn, padding = pg, upsample = up, nlevels = 2)
        x, ys = dummy_data(inchannels = ic, nclasses = nc, volume = vl,
            padding = pg, nlevels = 2)
        @test (model(x) |> size) == ys
    end

    @test_throws ArgumentError unet(inchannels = 3, volume = true,
                                    upsample = :bilinear)
end