text:
function _walkdown(f::Function, root::AbstractString;
        keepout::Function, onerr::Function
    )

    content = try; readdir(root)
        catch err; (onerr(root, err) === true) && return
    end

    # walk dir
    path::String = ""
    subi = firstindex(content)
    sub1 = lastindex(content)
    while subi <= sub1
        try
            for _ in subi:sub1
                name = content[subi]
                subi += 1

                path = joinpath(root, name)
                (f(path) === true) && return

                if isdir(path)
                    keepout(path) && continue
                    _walkdown(f, path; keepout, onerr)
                end
            end
        catch err;
            (onerr(path, err) === true) && return
        end
    end
end

function _walkdown_th(f::Function, root::AbstractString;
        keepout::Function, onerr::Function,
        endsig::Ref{Bool}, thfrec::Float64
    )

    endsig[] && return

    content = try; readdir(root)
        catch err; (onerr(root, err) === true) && (endsig[] = true; return)
    end

    # walk dir
    allth = (thfrec >= 1.0)
    path::String = ""
    subi = firstindex(content)
    sub1 = lastindex(content)
    @sync while subi <= sub1
        try
            for _ in subi:sub1
                name = content[subi]
                subi += 1

                path = joinpath(root, name)
                (f(path) === true) && (endsig[] = true)
                endsig[] && return

                if isdir(path)
                    keepout(path) && continue
                    if allth || (rand() < thfrec)
                        @spawn _walkdown_th(f, $path; keepout, onerr, endsig, thfrec)
                    else
                       _walkdown_th(f, path; keepout, onerr, endsig, thfrec)
                    end
                end
            end
        catch err
            (onerr(path, err) === true) && (endsig[] = true)
            endsig[] == true
        end
    end

end

"""
`walkdown_th(f, root; onerr::Function, keepout::Function, thfrec::Float64)`

walkdown the file tree applying `f` to all the founded paths.
The return value of `f` is consider a break flag, so if it returns `true`
the walk if over.
`keepout` is a filter that disallows walks a dir if returns `true`.
`thfrec` = [0,1] determines the frecuency of @spawn threads.
That is, `thfrec > 0` makes the method multithreaded.
[NOTE] This method do not waranty thread safetiness in any of it callbacks,
`f` or `keepout`. You must do it for yourself.

defaults:
    onerr = (path, err) -> (err isa InterruptException) ? rethrow(err) : nothing
    keepout = (dir) -> false
    thfrec = 0.0

"""
function walkdown(f::Function, root;
        keepout::Function = _default_keepout,
        onerr::Function = _default_onerr,
        thfrec::Float64 = 0.0
    )
    thfrec > 0.0 ?
        _walkdown_th(f, root; keepout, onerr, endsig = Ref(false), thfrec) :
        _walkdown(f, root; keepout, onerr)

    return nothing
end


text:
@reset
@testset "Moller-Plesset" begin
    @set {
          basis cc-pvdz
          scf_alg df
          jkfit cc-pvdz-jkfit
          rifit cc-pvdz-rifit
         }
    wf = @energy rmp2;
    @test isapprox(wf.CorrelationEnergy, -0.2062284914, rtol=tol)
    wf = @energy rmp3;
    @test isapprox(wf.CorrelationEnergy, -0.21283529233693912, rtol=tol)
end

ext: jl


text:
using Documenter, SearchMethods

makedocs(modules=[SearchMethods],
        doctest=true)

deploydocs(deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "https://github.com/EricAlcaide/SearchMethods.jl.git",
    julia  = "0.6.2",
    osname = "windows")


text:
module TestWater

using Base.Test
using SoilTracers.PhysChem: water_density, water_dissoc

@time @testset "Water properties" begin
    # test water density
    rho_4C = water_density(277.15)
    @test typeof(rho_4C) <: Real
    @test isapprox(rho_4C, 999.9, rtol=1e-4)

    # test water dissociation coefficient
    pK_w_25C = water_dissoc(298.15)
    @test typeof(pK_w_25C) <: Real
    @test isapprox(pK_w_25C, 14., rtol=1e-3)
end

end


text:
@testset "image graph construction" begin
  # these points are null, but we only need the size for the indices tests
  img = Array{SnPoint}(undef,2,3)
  K = constructImageGraph(img,"firstOrderDifference")
  @test K==[(1,2),(3,4),(5,6),(1,3),(2,4),(3,5),(4,6)]
  K2 = constructImageGraph(img,"secondOrderDifference")
  @test K2 == [(1,3,5),(2,4,6)]
  @test_throws ErrorException constructImageGraph(img,"K") # any other yields an error
end


text:
"""
    _parse_formula(f::String)::MetaboliteFormula

Parse a formula in format `C2H6O` into a [`MetaboliteFormula`](@ref), which is
basically a dictionary of atom counts in the molecule.
"""
function _parse_formula(f::String)::MetaboliteFormula
    res = Dict{String,Int}()
    pattern = @r_str "([A-Z][a-z]*)([1-9][0-9]*)?"

    for m in eachmatch(pattern, f)
        res[m.captures[1]] = isnothing(m.captures[2]) ? 1 : parse(Int, m.captures[2])
    end

    return res
end

"""
    _unparse_formula(f::MetaboliteFormula)::String

Format [`MetaboliteFormula`](@ref) to `String`.
"""
function _unparse_formula(f::MetaboliteFormula)::String
    return join(["$elem$n" for (elem, n) in f])
end


text:
module FakeTerminals

import REPL
using Test

export fake_terminal

function fake_terminal(f; options::REPL.Options=REPL.Options(confirm_exit=false))
    # Use pipes so we can easily do blocking reads
    # In the future if we want we can add a test that the right object
    # gets displayed by intercepting the display
    input = Pipe()
    output = Pipe()
    err = Pipe()
    Base.link_pipe!(input, reader_supports_async=true, writer_supports_async=true)
    Base.link_pipe!(output, reader_supports_async=true, writer_supports_async=true)
    Base.link_pipe!(err, reader_supports_async=true, writer_supports_async=true)

    term_env = get(ENV, "TERM", @static Sys.iswindows() ? "" : "dumb")
    term = REPL.Terminals.TTYTerminal(term_env, input.out, IOContext(output.in, :color=>true), err.in)
    f(term, input.in, output.out)
    t = @async begin
        close(input.in)
        close(output.in)
        close(err.in)
    end
    @test read(err.out, String) == ""
    wait(t)
end

end


text:
@testset "db" begin
    # drop database "test_db" if exists
    drop_db(DBNAME)

    # create db
    create_database(DBNAME)
    @test DBNAME in from_sql("SELECT datname FROM pg_database;")[!, 1]

    # enable uuid utils in "test_db"
    @test enable_uuid(DBNAME) == 0

    # table type
    struct TestTable end
    Base.string(::Type{TestTable}) = "test_table"

    # modify dbconfig
    dbconfig = copy(current_dbconfig())
    dbconfig[:dbname] = DBNAME

    # create table including uuid column into "test_db"
    column_names = [
        "ID",
        "FloatVal",
        "IntVal",
        "ByteVec",
    ]

    create_table(
        string(TestTable),
        """
        CREATE TABLE $(string(TestTable)) (
            $(column_names[1]) UUID DEFAULT uuid_generate_v4(),

            $(column_names[2]) FLOAT8 NOT NULL,
            $(column_names[3]) INT8 NOT NULL,
            $(column_names[4]) BYTEA COMPRESSION lz4 NOT NULL,

            PRIMARY KEY ($(column_names[1]))
        );
        """,
        dbconfig=dbconfig
    )

    column_names_from_sql = from_sql(
        """
        SELECT column_name
            FROM information_schema.columns
            WHERE table_schema = 'public' AND table_name = '$(string(TestTable))';
        """,
        dbconfig=dbconfig
    )[!, 1]

    @test all(lowercase(col) in column_names_from_sql for col in column_names)

    # test insert
    f = rand(Float64)
    i = rand(Int64)
    dim = 100
    b = rand(ComplexF64, dim, dim)

    df = DataFrame([
        :floatval => f,
        :intval => i,
        :ByteVec => hexbytes_str(b)
    ])

    to_sql(df, TestTable, dbconfig=dbconfig)

    df_from_sql = from_sql(TestTable, dbconfig=dbconfig)

    @test Float64(df_from_sql[1, :floatval]) == f
    @test Int64(df_from_sql[1, :intval]) == i
    @test all(reshape(hexbytes2array(ComplexF64, df_from_sql[1, :bytevec]), dim, dim) .== b)
end

@testset "insert" begin
    dbconfig = copy(current_dbconfig())
    dbconfig[:dbname] = DBNAME

    n = 10

    f = rand(Float64, n)
    i = rand(Int64, n)
    dim = 100
    b = [rand(ComplexF64, dim, dim) for _ in 1:n]

    df = DataFrame([
        :floatval => f,
        :intval => i,
        :ByteVec => hexbytes_str.(b)
    ])

    to_sql(df, TestTable, dbconfig=dbconfig)

    df_from_sql = from_sql(TestTable, dbconfig=dbconfig)

    for nᵢ in 1:n
        # index from db should +1 because we already insert 1 row in previous testset
        @test Float64(df_from_sql[nᵢ+1, :floatval]) == f[nᵢ]
        @test Int64(df_from_sql[nᵢ+1, :intval]) == i[nᵢ]
        @test all(reshape(hexbytes2array(ComplexF64, df_from_sql[nᵢ+1, :bytevec]), dim, dim) .== b[nᵢ])
    end
end


text:
call(f, xs...) = f(xs...)

# note for optimisers: set to zero
# p.Δ at the end of the weights update
function optimiser(ps, fs...)
  ps = [Param(p) for p in ps]
  fs = map(ps) do p
    os = map(f -> f(p), fs)
    () -> foreach(call, os)
  end
  () -> foreach(call, fs)
end

"""
    SGD(params, η = 0.1; decay = 0)

Classic gradient descent optimiser with learning rate `η`.
For each parameter `p` and its gradient `δp`, this runs `p -= η*δp`.

Supports inverse decaying learning rate if the `decay` argument is provided.
"""
SGD(ps, η = 0.1; decay = 0) =
  optimiser(ps, p -> invdecay(p, decay), p -> descent(p,η))

"""
    Momentum(params, η = 0.01; ρ = 0.9, decay = 0)

SGD with learning rate  `η`, momentum `ρ` and optional learning rate inverse decay.
"""
Momentum(ps, η = 0.01; ρ = 0.9, decay = 0) =
  optimiser(ps, p->invdecay(p,decay), p->momentum(p, ρ, η), p->descent(p,1))

"""
    Nesterov(params, η = 0.01; ρ = 0.9, decay = 0)

SGD with learning rate  `η`, Nesterov momentum `ρ` and optional learning rate inverse decay.
"""
Nesterov(ps, η = 0.01; ρ = 0.9, decay = 0) =
  optimiser(ps, p->invdecay(p,decay), p->nesterov(p, ρ, η), p->descent(p,1))

"""
    RMSProp(params, η = 0.001; ρ = 0.9, ϵ = 1e-8, decay = 0)

[RMSProp](http://www.cs.toronto.edu/~tijmen/csc321/slides/lecture_slides_lec6.pdf)
optimiser. Parameters other than learning rate don't need tuning. Often a good
choice for recurrent networks.
"""
RMSProp(ps, η = 0.001; ρ = 0.9, ϵ = 1e-8, decay = 0) =
  optimiser(ps, p->rmsprop(p; η=η, ρ=ρ, ϵ=ϵ), p->invdecay(p,decay), p->descent(p,1))

"""
    ADAM(params, η = 0.001; β1 = 0.9, β2 = 0.999, ϵ = 1e-08, decay = 0)

[ADAM](https://arxiv.org/abs/1412.6980v8) optimiser.
"""
ADAM(ps, η = 0.001; β1 = 0.9, β2 = 0.999, ϵ = 1e-08, decay = 0) =
  optimiser(ps, p->adam(p; η=η, β1=β1, β2=β2, ϵ=ϵ), p->invdecay(p,decay), p->descent(p,1))

"""
   ADAMW((params, η = 0.001; β1 = 0.9, β2 = 0.999, ϵ = 1e-08, decay = 0)

[ADAMW](https://arxiv.org/abs/1711.05101) fixing weight decay regularization in Adam.
"""
ADAMW(ps, η = 0.001; β1 = 0.9, β2 = 0.999, ϵ = 1e-08, decay = 0) =
  optimiser(ps, p->adam(p; η=η, β1=β1, β2=β2, ϵ=ϵ), p->descentweightdecay(p,1,decay))

"""
    AdaMax(params, η = 0.001; β1 = 0.9, β2 = 0.999, ϵ = 1e-08, decay = 0)

[AdaMax](https://arxiv.org/abs/1412.6980v9) optimiser. Variant of ADAM based on
the ∞-norm.
"""
AdaMax(ps, η = 0.002; β1 = 0.9, β2 = 0.999, ϵ = 1e-08, decay = 0) =
  optimiser(ps, p->adamax(p; η=η, β1=β1, β2=β2, ϵ=ϵ), p->invdecay(p,decay), p->descent(p,1))

"""
    ADAGrad(params, η = 0.01; ϵ = 1e-8, decay = 0)

[ADAGrad](http://www.jmlr.org/papers/volume12/duchi11a/duchi11a.pdf) optimiser.
Parameters don't need tuning.
"""
ADAGrad(ps, η = 0.01; ϵ = 1e-8, decay = 0) =
  optimiser(ps, p->adagrad(p; η=η, ϵ=ϵ), p->invdecay(p,decay), p->descent(p,1))

"""
    ADADelta(params; ρ = 0.9, ϵ = 1e-8, decay = 0)

[ADADelta](http://arxiv.org/abs/1212.5701) optimiser. Parameters don't need
tuning.
"""
ADADelta(ps; ρ = 0.9, ϵ = 1e-8, decay = 0) =
  optimiser(ps, p->adadelta(p; ρ=ρ, ϵ=ϵ), p->descent(p,1))

"""
    AMSGrad(params; η = 0.001, β1 = 0.9, β2 = 0.999, ϵ = 1e-08, decay = 0)

[AMSGrad](https://openreview.net/forum?id=ryQu7f-RZ) optimiser. Parameters don't need
tuning.
"""
AMSGrad(ps, η = 0.001; β1 = 0.9, β2 = 0.999, ϵ = 1e-08, decay = 0) =
  optimiser(ps, p -> amsgrad(p; η = η, β1 = β1, β2 = β2, ϵ = ϵ), p -> invdecay(p, decay), p -> descent(p, 1))

"""
    NADAM(params, η = 0.001; β1 = 0.9, β2 = 0.999, ϵ = 1e-08, decay = 0)

[NADAM](https://openreview.net/pdf?id=OM0jvwB8jIp57ZJjtNEZ) optimiser. Parameters other
than learning rate don't need tuning.
"""
NADAM(ps, η = 0.001; β1 = 0.9, β2 = 0.999, ϵ = 1e-08, decay = 0) =
  optimiser(ps, p->nadam(p; η=η, β1=β1, β2=β2, ϵ=ϵ), p->invdecay(p,decay), p->descent(p,1))


text:
using Test
using DataFrames
using CSVFiles
using CSV
using Query
using Statistics
#
# local imports; see ~/.julia/config/startup.jl
#
using GlobalDecls
using LAModelData: get_underlying_data, get_18s_level_from_doe, get_staying_put_rates,
    do_exit_rate_query, get_target_la, get_pay_class_from_region_code,
    UDATA
using StayingPutSim
using CareData
using ONSCodes
using StayingPutModelDriver

@testset "LADataTests" begin
    nt = get_underlying_data()
    @test nt.c19_21_2019[1,:cl_all_19to21] == 29930
    print(nt.c19_21_2019[!,:cl_stayput_ffc_20_pc] )
    @test nt.c19_21_2019[164,:cl_stayput_ffc_20_pc] == 48
    @test UDATA.c17_18_2017[1,Symbol("cl_all_17&18")]== 10_710
    @test UDATA.c19_21_2017[164,:cl_stayput_ffc_20] == 6
    # E09000028 is Southwark;
    # E12000007 is London
    @test get_18s_level_from_doe("E09000028") == 37
    @test do_exit_rate_query( 2019, "E09000028", 18 ) == (37,23)
    @test get_target_la( "E09000007", regional ) == "E12000007"
    spr = get_staying_put_rates( "E09000028", local_authority, false )
    println( "spr=$spr")
    @test spr[1] ≈ 23/37 # 18 yos in
    all = CareData.load_all(2020)
    @test all.ofdata[1,:council] == "England"
    @test all.ofdata[24,:sp_fee_2] ≈ 207.3314285714
end

@testset "Output Tests" begin
    outdir = "/home/graham_s/VirtualWorlds/projects/action_for_children/england/results/testcase/"
    num_systems = 3
    year = 2020
    StayingPutModelDriver.createmaintables( outdir, num_systems, :ccode, year )
end


text:
using GMAT
using Test

@testset "GMAT" begin
	@testset "Run Script" begin
		test_script = joinpath(@__DIR__, "test.script")
		output = run_script(test_script; output=tempname())
		@test length(output) == 2
	end
end


text:
# This file is a part of MGVI.jl, licensed under the MIT License (MIT).

__precompile__(true)

"""
    MGVI

An implementation of the Metric Gaussian Variational Inference algorithm.
"""
module MGVI

import ChainRulesCore
using Distributed
using LinearAlgebra
using Random
using SparseArrays
using Base.Iterators
using Distributions
using DistributionsAD
import ForwardDiff
using LinearMaps
using IterativeSolvers
using Optim
using PositiveFactorizations
import SparseArrays: blockdiag
using SparseArrays
using StaticArrays
using ValueShapes
import Zygote

using ChainRulesCore: AbstractTangent, Tangent, NoTangent, ZeroTangent, AbstractThunk, unthunk

include("custom_linear_maps.jl")
include("shapes.jl")
include("jacobian_maps.jl")
include("information.jl")
include("residual_samplers.jl")
include("mgvi_impl.jl")

end # module


text:
# # Chapter 1. Computer Arithmetic

# ## 1.1 Positional Systems
# $\beta \in \mathbb{Z}, \beta \ge 2$ : Integer base
#
# For $\sigma \in \{-1, 1\}$, $b_i \in \{ 0, 1, \ldots, \beta-1\}$,
#
# \begin{align*}
# (-1)^\sigma (b_n b_{n-1} \cdots b_0 . b_{-1} b_{-2} \cdots)_\beta
# &:= (-1)^\sigma \sum_{i=-\infty}^n b_i \beta^i  \\
# &= (-1)^\sigma (b_n \beta^n + b_{n-1} \beta^{n-1} + \cdots + b_0 + b_{-1} \beta^{-1} + b_{-2} \beta^{-2} + \cdots)
# \end{align*}
#
# Condition (a) $0 \le b_i \le \beta - 1$ for all $i$.
#
# Condition (b) $0 \le b_i \le \beta - 2$ for infinitely many $i$.

# **Prop**. For all nonzero $x \in \mathbb{R}$ has a unique representation above.
#
# **Proof**. (Uniqueness) Let $x = (-1)^\sigma \sum_{i=-\infty}^n b_i \beta^i = (-1)^{\sigma^\prime} \sum_{i=-\infty}^{n^\prime} b^\prime_i \beta^i$.
# We can assume $\sigma = \sigma^\prime = 1$ and $n \ge n^\prime$. Set $b_i = 0$ for $b_i^\prime = 0$ for $n < i \le n^\prime$.
# By (b) we have $|b_n - b^\prime_n|\beta^n =  |b_n \beta^n - b^\prime_n\beta^n| \le \sum_{i=-\infty}^{n-1}|b_i^\prime - b_i|\beta^i <
# \sum_{i=-\infty}^{n-1}(\beta-1)\beta^i = (\beta - 1)\beta^{n-1}\frac{1}{1-1/\beta} = \beta^n$.
# Since $b_n, b_n^\prime \in \mathbb{Z}$ and $|b_n - b_n^\prime| < 1 $, we get $b_n = b_n^\prime$.
# It is clear that $x - b_n$ also satisfies the conditions (a) and (b), we conclude $b_{n-1} = b_{n-1}^\prime$.
# By continuing this process, we see $b_i = b_i^\prime$ for all $i$.
#
# (Existance) Without loss of generality, we can assume $x > 0$.
# By dividing $x$ by some $\beta^i$, it suffices to show when $0 < x \le 1$.
# The existance is trivial if $x = 1$. We recursively construct for all $n > 0$,
# an expression $x_n^\prime = \Sigma_{i=1}^n b_{-i}\beta^{-i}$ such that $0 \le x - x_n^\prime < \beta^{-n}$.
# Let $b_{-1}$ be the integer $i \in [0, \beta - 1]$ such that $x \in [i \beta^{-1}, (i+1) \beta^{-1})$.
# It is clear that $0 \le x - x_1^\prime < \beta^{-1}$.
# Reaplacing $x$ by $(x - x_1^\prime)/\beta$ and follow the same procesure, we can find $b_{-2}$.
# Continuing this process, we obtain a sequence $b_{-i}$.
# If this sequance violates the condition (b), there exists $m$ such that $b_{-i} = \beta - 1$ for all $i > m$.
# By constuction, $x - x_{m-1}^\prime < \beta^{-m}$ but $x - x_{m-1}^\prime = \Sigma_{i=m}^{\infty} (\beta-1)\beta^{-i} = \beta^{-m}$, it leads a contradiction.

# ## 1.2 Floating Point Numbers
# \begin{align*}
# \mathbb{R} \supset &\mathbb{F}_\beta = \{ (-1)^\sigma m \times \beta^e \mid m = (b_0 . b_{-1} b_{-2} \cdots)_\beta, m \text{ satisfies (a) and (b)} \} \\
# \supset &\mathbb{F}_{\beta, p} = \{ x \in \mathbb{F}_\beta \mid m = (b_0 . b_{-1} b_{-2} \cdots b_{p-1})_\beta \} \\
# \supset &\mathbb{F}_{\beta, p}^{\check{e}, \hat{e}} = \{ x \in \mathbb{F}_{\beta, p} \mid \check{e} \le e \le \hat{e} \}
# \end{align*}
# $m$ : *mantiss*a of $x$, $e$ : *exponent* of $x$, $p$ : *precision*
#
# $\mathbb{F}_\beta$ : *floating point numbers*
#
# $\mathbb{F}_{\beta, p}^{\check{e}, \hat{e}} $ : *computer representable floating point numbers*
#

# **Def**. We say a floating point number $x = (-1)^\sigma (b_0 . b_{-1} b_{-2} \cdots)_\beta \times \beta^e$ is *normalized* if the leading digit $b_0$ is non-zero.

# Let's show the *largest positive normal number* $N_{\max}^n$ and
# the *smallest positive normal number* $N_{\min}^n$.
# The (union of) range $\{ x \in \mathbb{R} \mid N_{min}^n \le |x| \le N_{max}^n \}$ is said to be the *normalized range*.

#-
## Normalized range
const FloatTypes = (Float16, Float32, Float64)

for T in FloatTypes
    println("Normalized range of $(T) is ", floatmin(T), " <= |x| <= ", floatmax(T))
end

# ### 1.2.1 Subnormal numbers
# **Def**. $x \in \mathbb{F}_{\beta, p}^{\check{e}, \hat{e}}$ is said to be *subnormal* if $b_0 = 0$ and $e = \check{e}$.

# Let's find the *largest positive subnormal number* $N_{\max}^s$ and
# the *smallest positive subnormal nubmer* $N_{\min}^s$.

for T in FloatTypes
    println("$(T), N_min^s = ", nextfloat(zero(T)), ", N_max^s = ", prevfloat(floatmin(T)))
end

# ## 1.3. Rounding
# Let $\mathbb{R}^* = \mathbb{R} \cup \{ -\infty, +\infty \}, \mathbb{F}^* = \mathbb{F} \cup \{ -\infty, \infty \}$.
#
# **Def**. A *rounding* is a mapping $\bigcirc : \mathbb{R}^* \rightarrow \mathbb{F}^*$ satisfying the conditions (R1) and (R2):
#
# (R1) $x \in \mathbb{F}^* \Rightarrow \bigcirc (x) = x$.
#
# (R2, *Maximum quality*) $x, y \in \mathbb{R}^*, x \le y \Rightarrow \bigcirc (x) \le \bigcirc (y)$

# **Def**. Let $\bigcirc$ be a rounding.
#
# (R3) $\bigcirc$ is *odd* if $\bigcirc (-x) = - \bigcirc (x)$ for all $x \in \mathbb{R}^*$.
#
# (R4a) $\bigcirc$ is *round toward minus infinity* if $\bigcirc (x) \le x$ for all $x \in \mathbb{R}^*$.
#
# (R4a) $\bigcirc$ is *round toward plus infinity* if $\bigcirc (x) \ge x$ for all $x \in \mathbb{R}^*$.

# **Example**.
# *Rounded down:* $\triangledown (x) = \max \{y \in \mathbb{F}^* \mid y \le x\}$,
#
# *Rounded up:* $\triangle (x) = \min \{y \in \mathbb{F}^* \mid y \ge x\}$.

# If $x \in \mathbb{F}^*$, then $\triangledown (x) = x = \triangle (x)$, on the contrary
# if $x \notin \mathbb{F}^*$,  we have $\triangledown (x) = x = \triangle (x)$
#
# $\triangledown$ and $\triangle$ give the lower bound and upper bound of roundings: if $\bigcirc$ is a rounding,
# $\triangledown (x) \le \bigcirc (x) \le \triangle (x)$ for all $x$.

# Julia can create a float from a given `x` with a explicit rounding mode.
# See [`Core.Float32`](https://docs.julialang.org/en/v1/base/numbers/#Core.Float32-Tuple{Any}),
# [`Core.Float64`](https://docs.julialang.org/en/v1/base/numbers/#Core.Float64-Tuple{Any}).
# However, `0.1` is a floating-point literal and already converted into `Float64`,
# the expression `Float64(0.1, RoundDown)` is not equal to $\triangledown (0.1)$.

#-
using Printf

println(typeof(0.1))
@printf("%0.17f\n", Float64(0.1, RoundDown))

# Similarly, `1/10` is a floating-point literal.

#-
println(typeof(1/10))
@printf("%0.17f\n", Float64(1/10, RoundDown))

# Using a [rational number type](https://docs.julialang.org/en/v1/manual/complex-and-rational-numbers/#Rational-Numbers-1),
# we can calculate $\triangledown (0.1)$ properly.

#-
println(typeof(1//10))
@printf("%0.17f\n", Float64(1//10, RoundDown))

#-
## Rounding methods (in Float64)
rounding_modes = [RoundNearest,
## RoundNearestTiesAway,
## RoundNearestTiesUp,
                  RoundToZero,
                  RoundUp,
                  RoundDown]

for mode in rounding_modes
    print("$(mode)")
    for x in [1//10, -1//10, π]
        print(", ∘($(x)) = ")
        @printf("%0.17f", Float64(x, mode))
    end
    println()
end

# ### 1.3.1 Round to Zero
# **Def**. $\square_z : \mathbb{R}^* \rightarrow \mathbb{F}^*, \square_z (x) =
# \text{sign} (x) \max \{ y \in \mathbb{F}^* \mid y \le |x| \} =
# \text{sign} (x) \triangledown (|x|)$.

# ### 1.3.4 Rounding Errors

#-
for T in FloatTypes
    @printf("1//10, %s, Round down = %0.17f, Round up = %0.17f\n", T, T(1//10, RoundDown), T(1//10, RoundUp))
end

# ## 1.4. Floating Point Arithmetic
# For any one of the arithmetic operations $\star \in \{ +, -, \times, / \}$,
# let $⍟ \in \{ \oplus, \ominus, \otimes, \oslash \}$ denote the corresponding operator in $\mathbb{F}$.
#
# **Def**. The floating point arithmetis is said to have *maximum quality* if
# (R5) $x, y \in \mathbb{F}$ and $\star \in \{ +, -, \times, / \} \Rightarrow x ⍟ y = \bigcirc (x \star y)$.

# ### Julia Base
# Julia has a method [`Base.Rounding.setrounding`](https://docs.julialang.org/en/v1/base/numbers/#Base.Rounding.setrounding-Tuple{Type,Any}) but it only supports `BigFloat`.
# See [Deprecate setrounding #27166](https://github.com/JuliaLang/julia/pull/27166).

# ### FastRounding.jl
# [FastRounding.jl](https://github.com/JeffreySarnoff/FastRounding.jl) provides arithmetic operations
# with directed rounding modes.

#-
using FastRounding

a = 1/3
b = 3.0

println("a * b, rounded down = ", mul_round(a, b, RoundDown))

println("a * b, rounded down = ", ⊗₋(a, b)) # Equivalent exports

println("√2, rounded up = ", sqrt_round(2.0, RoundUp))

# **Example**. Define the tenary shift map
# $f : [0, 1] \rightarrow [0,1]$ by $f(x) = 3x \mod 1$.
# It has a cycle $\frac{1}{10} \rightarrow \frac{3}{10} \rightarrow \frac{9}{10} \rightarrow \frac{7}{10} \rightarrow \frac{1}{10}$ of four.

#-
## tenary shift map
tenary_shift(x) = mod(3x, 1)

let x = 0.1
    for i in 0:52
        @printf("x(%d) = %0.17f\n", i, x)
        x = tenary_shift(x)
    end
end

# **Example**. Calculate the values of the following three symbolically equivalent expressions at $x = 0.0031834$:
# 1. $r_1(x) = 1 - \frac{1}{\sqrt{1 + x}}$,
# 2. $r_2(x) = \frac{\sqrt{1 + x} - 1}{\sqrt{1 + x}}$,
# 3. $r_3(x) = \frac{x}{1 + x + \sqrt{1 + x}}$

#-
x = 0.0031834

r1(x) = 1 - 1 / sqrt(1 + x)
r2(x) = (sqrt(1 + x) - 1) / sqrt(1 + x)
r3(x) = x / (1 + x + sqrt(1 + x))

for T in FloatTypes
    @printf("%s, r1(x) = %0.17f, r2(x) = %0.17f, r3(x) = %0.17f\n", T, r1(T(x)), r2(T(x)), r3(T(x)))
end

# **Example**.  $\theta = 1.23456^\circ$,
# 1. $s_1(\theta) = \sin \left( \frac{\theta}{2} \right)^2$,
# 2. $s_2(\theta) = \frac{1 - \cos(\theta)}{2}$
#-
θ = 1.23456

s1(x) = sind(x / 2)^2
s2(x) = (1 - cosd(x)) / 2

for T in FloatTypes
    @printf("%s, s1(x) = %0.17f, s2(x) = %0.17f\n", T, s1(T(θ)), s2(T(θ)))
end

# **Def**. The machine epsilon $\epsilon_M$ is the smallest positive floating point number $x$ satisfying
# $\triangledown (1 + x) > 1$. i.e. $\min \{ x \in \mathbb{F} \mid \triangledown (1 + x) > 1 \}$.
#
# **Def**. The machine eta $\eta_M$ is the smallest positive floating point number $x$ satisfying
# $\triangledown (x) > 0$. It is equal to the smallest positive subnormal number $N_{\min}^s$.

#-
## Machine epsilon
for T in FloatTypes
    println("Machine epsilon of $(T) = ", eps(T))
end

## Machine epsilon
for T in FloatTypes
    println("Machine eta of $(T) = ", eps(zero(T)))
end

# ## 1.5. The IEEE standard

#-
using Base: exponent_bits, exponent_mask, exponent_one, significand_bits
for T in FloatTypes
    ## The followings are defined on >= 1.4
    @eval exponent_bias(::Type{$T}) = $(Int(exponent_one(T) >> significand_bits(T)))
    @eval exponent_max(::Type{$T}) = $(Int(exponent_mask(T) >> significand_bits(T)) - exponent_bias(T))

    println("$(T), exponent width in bits = ", exponent_bits(T), ", precision = ", precision(T),
            ", exponent bias = ", exponent_bias(T), ", maximal exponent = ", exponent_max(T))
end

#-
ieee_digits(::Type{T}) where {T<:Base.IEEEFloat} = (1, Base.exponent_bits(T), precision(T) - 1)

function print_float_format(x::T) where {T<:Base.IEEEFloat}
    bit_x = bitstring(x)
    digits = ieee_digits(T)
    n = sum(digits) + Base.length(digits) + 1
    print("-")
    i = 1
    for j in digits
        print(bit_x[i:i+j-1], "-")
        i += j
    end
    println()
end

#-
## Bitstring expression of normalized range
for T in FloatTypes
    n_min_n = floatmin(T)
    n_max_n = floatmax(T)
    println("The smallest positive normal number of $(T) = ", n_min_n)
    print_float_format(n_min_n)
    println()
    println("The largest normal number of $(T) = ", n_max_n)
    print_float_format(n_max_n)
    println()
end

#-
## Subnormal numbers
for T in FloatTypes
    n_min_s = nextfloat(zero(T))
    n_max_s = prevfloat(floatmin(T))
    println("The smallest positive subnormal number of $(T) = ", n_min_s)
    print_float_format(n_min_s)
    println()

    println("The largest subnormal number of $(T) = ", n_max_s)
    print_float_format(n_max_s)
    println()
end

#-
## NaN, -Inf, Inf
for T in FloatTypes
    ninf = T(-Inf)
    inf = T(Inf)
    println("NaN of $(T):")
    print_float_format(T(NaN))
    println()

    println("-Inf of $(T):")
    print_float_format(T(-Inf))
    println()

    println("Inf of $(T):")
    print_float_format(T(Inf))
    println()
end

# ## 1.6. Examples of Floating Point Computations
# **Example**. Print the first 40 factorials $n! = 1 \cdot 2 \cdots n$.

#-
factorial(n) = n == 0 ? 1 : n * factorial(n - 1)

for n in 0:40
    println(n, "! = ", factorial(n))
end

# **Example**. Evaluate $f(x, y) = 333.75y^6 + x^2(11x^2y^2 - y^6 - 121y^4 - 2) + 5.5y^8 + x/(2y)$ for $(x, y) = (77617, 33096)$

#-
rump(x::Real, y::Real) = oftype(x, 333.75) * y^6 + x^2 * (11x^2 * y^2 - y^6 - 121y^4 -2) + oftype(x, 5.5) * y^8 + x/(2y)

x, y = 77617, 33096
for T in FloatTypes
    println("rump(x, y) = ", rump(T(x), T(y)))
end

# **Example**. Plot
# \begin{align*}
# p(x) &= t^6 - 6t^5 + 15t^4 - 20t^3 + 15t^2 - 6t + 1 \ (\text{expanded}) \\
# &= (((((t - 6)t + 15)t - 20)t + 15)t - 6) + 1 \ (\text{Horner}) \\
# &= (t - 1)^6 \ (\text{factored})
# \end{align*}
#
# around $x = 1.0$.

#-
using Plots

xs = range(0.995, 1.005, length = 100)

## Expanded
p1(t) = t^6 - 6t^5 + 15t^4 - 20t^3 + 15t^2 - 6t + 1
plot(xs, p1.(xs), ls = :dash, label = "expanded")

## Horner
p2(t) = (((((t - 6) * t + 15) * t - 20) * t + 15) * t - 6) * t + 1
plot!(xs, p2.(xs), ls = :solid, label = "Horner")

## factored
p3(t) = (t - 1)^6
plot!(xs, p3.(xs), ls = :dashdot, label = "factored")

# ## 1.7 Computer Lab 1
# **Example**. Find the smallest $x \in (1, 2)$ such that $x \otimes \frac{1}{x} \neq 1$.
#-
for i in 1:1000000000
    x = nextfloat(1.0, i)
    y = 1/x
    if x * y != 1.0
        @printf("x = %0.17f\n", x)
        print_float_format(x)
        break
    end
end

# **Example**. Evaluate $f(x, y) = 9x^4 - y^4 +2y^2$ for $(x, y) = (40545, 70226)$.

#-
prob4(x, y) = 9x^4 - y^4 + 2y^2

x, y = 40545, 70226

for T in [Int64, Float16, Float32]
    println("f(x, y) = ", prob4(T(x), T(y)))
end

# Correct answer is (3x^2 - x^2)(3x^y + x^2) + 2y^2 = -1 * 9863382151 + 9863382152 = 1

max_stars_repo_name: matsueushi/study-validated-numerics

max_issues_repo_name: matsueushi/study-validated-numerics

max_forks_repo_name: matsueushi/study-validated-numerics


text:
include("Include.jl")

# Where are my files stored?
path_to_sim_files = "./best_ensemble"

# Load the ensemble for AhR -
time_AhR = [0.0 48.0]
(raw_data_AhR,scaled_data_AhR) = process_ensemble_data(path_to_sim_files,time_AhR,38);
writedlm("./plot_best/scaled_data_AhR.dat",scaled_data_AhR);

# Load the ensemble for CD38 -
time_CD38 = [0.0 24.0 48.0 72.0]
(raw_data_CD38,scaled_data_CD38) = process_ensemble_data(path_to_sim_files,time_CD38,41);
writedlm("./plot_best/scaled_data_CD38.dat",scaled_data_CD38);

# Load the ensemble for CD11b -
time_CD11b = [0.0 24.0 48.0 72.0]
(raw_data_CD11b,scaled_data_CD11b) = process_ensemble_data(path_to_sim_files,time_CD11b,39);
writedlm("./plot_best/scaled_data_CD11b.dat",scaled_data_CD11b);

# Load the ensemble for P47Phox -
time_P47Phox = [0.0 48.0]
(raw_data_P47Phox,scaled_data_P47Phox) = process_ensemble_data(path_to_sim_files,time_P47Phox,50);
writedlm("./plot_best/scaled_data_P47Phox.dat",scaled_data_P47Phox);

# Load the ensemble for EGR1 -
time_EGR1 = [0.0 48.0]
(raw_data_EGR1,scaled_data_EGR1) = process_ensemble_data(path_to_sim_files,time_EGR1,44);
writedlm("./plot_best/scaled_data_EGR1.dat",scaled_data_EGR1);

# Load the ensemble for IRF1 -
time_IRF1 = [0.0 48.0]
(raw_data_IRF1,scaled_data_IRF1) = process_ensemble_data(path_to_sim_files,time_IRF1,46);
writedlm("./plot_best/scaled_data_IRF1.dat",scaled_data_IRF1);

# Load the ensemble for GIF1 -
time_GIF1 = [0.0 48.0]
(raw_data_GIF1,scaled_data_GIF1) = process_ensemble_data(path_to_sim_files,time_GIF1,45);
writedlm("./plot_best/scaled_data_GIF1.dat",scaled_data_GIF1);

# Load the ensemble for PU1 -
time_PU1 = [0.0 48.0]
(raw_data_PU1,scaled_data_PU1) = process_ensemble_data(path_to_sim_files,time_PU1,52);
writedlm("./plot_best/scaled_data_PU1.dat",scaled_data_PU1);

# Load the ensemble for OCT4 -
time_OCT4 = [0.0 48.0]
(raw_data_OCT4,scaled_data_OCT4) = process_ensemble_data(path_to_sim_files,time_OCT4,48);
writedlm("./plot_best/scaled_data_OCT4.dat",scaled_data_OCT4);

# Load the ensemble for CEBPa -
time_CEBPa = [0.0 48.0]
(raw_data_PU1,scaled_data_CEBPa) = process_ensemble_data(path_to_sim_files,time_CEBPa,42);
writedlm("./plot_best/scaled_data_CEBPa.dat",scaled_data_CEBPa);


text:
module Post

export getTables, printTable, postMonteCarlo, getDivRatio, exportTable!
export collectElast!, reportElast!, fastElasticity, divRatio_ms15

using TypeDef, Helper, Tools, UserFileIO, Compute, Solve
using JLD, Printf, Statistics, LinearAlgebra

function postMonteCarlo(
  model  ::Int64, # ID for simulated model
  mcid   ::Int64, # ID for Monte Carlo estimate for simulation
  vInt   ::Bool,  # Vertical integration model (in Appendix)
  ms15   ::Bool,  # Market size * 1.5 (in Appendix)
  mc_mvno::Int64, # Margial cost of MVNO (different normalizations)
  saveTab::Bool=false, # Export CSV tables of div ratio & elasticities
  listBS ::Vector{Int64}=[1], # List of bootstrap batches (single by default)
  b1     ::Int64=1,
  b2     ::Int64=200)

  dat,dim,env,mat,par0,set,sim,indx,mask =
    readData(model,mcid,b1,b2,vInt,false,ms15,mc_mvno)
  eqb,pab = collectMC(listBS,model,dim,env,mat,par0)
  tab = getTables(eqb,vInt,dat,dim,env,mat,pab,set,sim,indx,mask)
  printTable(tab,vInt,dat,dim,env,eqb,mat,pab,set,sim,indx,mask,true)
  if saveTab==true
    exportTable!(tab,eqb,dat,dim,env,mat,pab,set,sim,indx,mask,model)
  end

  return dim,eqb,tab
end

function exportTable!(
  tab ::tabType,
  eqb ::Array{Array{equilType,1}},
  dat ::datType,
  dim ::dimType,
  env ::envType,
  mat ::matType,
  pab ::Array{parType,1},
  set ::setType,
  sim ::simType,
  indx::indxType,
  mask::maskType,
  model::Int64)

  lab = setupLab()

  ### Diversion ratio & elasticity
  @printf("\n Exporting Diversion ratio (Table A.3) to divRatio.csv and\n")
  @printf(" elasticity table (Table 3) to elasR.csv.\n")
  println("This will take about 10 minutes.")
  for b in eachindex(eqb)
    tic = time()
    initialize!(mat,dat,dim,pab[b],set)
    tab.divRatio[:,:,b] = getDivRatio(1,eqb[b],dim,lab,mat,pab[b],sim,indx,mask)
    collectElast!(tab,b,eqb[b][1],dim,lab,mat,pab[b],sim,indx,mask)
    @printf("Div. ratio for bootstrap %d complete in %.1f sec.\n",b,time()-tic)
  end

  fname = string(env.pathOut,"/tab.jld")
  jldopen(fname,"w") do file
    addrequire(file,TypeDef)
    write(file,"tab",tab)
  end

  csvname = string(env.pathOut,"/divRatio.csv")
  divRatio = mean(tab.divRatio,dims=3)
  fid = open(csvname,"w")
  @printf(fid,"%8s",lab.product[1])
  for i in lab.prodIdx[2:end]
    @printf(fid,", %8s",lab.product[i])
  end
  @printf(fid,"\n")
  for i=1:size(divRatio,1)
    @printf(fid,"%8.2f",divRatio[i,1])
    for j=2:size(divRatio,2)
      @printf(fid,", %8.2f",divRatio[i,j])
    end
    @printf(fid,"\n")
  end
  close(fid)

  fr = open(string(env.pathOut,"/elasR.csv"),"w")
  fw = open(string(env.pathOut,"/elasW.csv"),"w")
  elasR = mean(tab.elasR,dims=3)
  elasW = mean(tab.elasW,dims=3)
  @printf(fr,"%8s",lab.product[1])
  for i in lab.prodIdx[2:end]
    @printf(fr,", %8s",lab.product[i])
  end
  @printf(fr,"\n")
  for i=1:size(elasR,1)
    @printf(fr,"%8.3f",elasR[i,1])
    for j=2:size(elasR,2)
      @printf(fr,", %8.3f",elasR[i,j])
    end
    @printf(fr,"\n")
  end
  for i=1:dim.product
    @printf(fw,"%8.3f",elasW[i,1])
    for j=2:dim.mvno
      @printf(fw,", %8.3f",elasW[i,j])
    end
    @printf(fw,"\n")
  end
  close(fr)
  close(fw)

  # fx =
  # let sim=deepcopy(sim)
  #   b->getDivRatio(1,eqb[b],dim,mat,pab[b],sim,indx,mask)
  # end
  # pout = pmap(fx, 1:8)
  # return pout

end

function printTable(
  tab  ::tabType,
  vInt ::Bool,
  dat  ::datType,
  dim  ::dimType,
  env  ::envType,
  eqb  ::Array{Array{equilType,1}},
  mat  ::matType,
  pab  ::Array{parType,1},
  set  ::setType,
  sim  ::simType,
  indx ::indxType,
  mask ::maskType,
  elast::Bool=false)

  lab = setupLab()

  @printf("\n------------------ Tables for Model %2d ",env.model)
  @printf("------------------\n")

  if vInt==false
    ### Summary elasticity
    if elast==true
      reportElast!(tab,vInt,env)
    end

    ### Retail markup table
    @printf("\n1. Retail market (Table 6)\n")
    @printf("%10s  %10s  %10s  %8s  %8s\n","Network","Product","Tariff",
            "Price","Markup")

    for i in lab.prodIdx
      f = mat.firm[i,2]
      @printf("%10s  %10s  %10s  ",lab.network[f],lab.product[i],lab.prepost[i])
      @printf("%8.2f  %8.2f\n",tab.retail[i,1],tab.retail[i,2])
    end

    ### Wholesale markup table
    if vInt==false
      @printf("\n2. Wholesale market (Table A.5)\n")
      @printf("%8s  %10s  %28s  %28s  %28s\n","Upstream","Downstream",
              "Wholesale.price","Wholesale.markup","Wholesale.mcost")
      @printf("%8s  %10s  ","network","network")
      for k=1:3
        for j=1:3
          @printf("%8s  ",lab.prepost[j])
        end
      end
      @printf("\n")
      for i=1:3
        @printf("%8s  %10s  ",lab.network[i],"MVNO")
        whsale = tab.whsale[mat.mvno[:,i],:]
        mcMean = mean(tab.mcost[1][mat.mvno[:,i],:],dims=2)
        mcStd  = std(tab.mcost[1][mat.mvno[:,i],:],dims=2)
        for k=1:2
          for j=1:3
            @printf("%8.2f  ",whsale[j,k])
          end
        end
        for j=1:3
          @printf("%8.2f  ",mcMean[j])
        end
        @printf("\n%8s  %10s  ","","")
        for k=3:4
          for j=1:3
            @printf("%8.2f  ",whsale[j,k])
          end
        end
        for j=1:3
          @printf("%8.2f  ",mcStd[j])
        end
        @printf("\n")
      end
    end
  end # end of if vInt==false

  ### Fixed cost table
  @printf("\n3. Fixed Cost Bounds (Table 8/A.13)\n")
  @printf("%10s  %10s  %10s  %10s  %10s  %10s\n",
          "Mean","Lower","Upper","f̄̄","Diff1","Diff2")
  for i=1:3
    lb = mean(tab.fb[i,1,:])
    ub = mean(tab.fb[i,2,:])
    f̄̄  = mean(tab.fb[i,3,:])
    ls = std(tab.fb[i,1,:])
    us = std(tab.fb[i,2,:])
    f̄̄s = std(tab.fb[i,3,:])
    d1 = mean(tab.fb[i,2,:]-tab.fb[i,1,:])
    d2 = mean(tab.fb[i,3,:]-tab.fb[i,1,:])
    s1 = std(tab.fb[i,2,:]-tab.fb[i,1,:])
    s2 = std(tab.fb[i,3,:]-tab.fb[i,1,:])
    @printf("%10s  %10.0f  %10.0f  %10.0f  %10.0f  %10.0f\n",
            lab.network[i],lb,ub,f̄̄,d1,d2)
    @printf("%10s  %10.0f  %10.0f  %10.0f  %10.0f  %10.0f\n",
            "(se)",ls,us,f̄̄s,s1,s2)
  end

  if vInt==false
    @printf("\n4. Fixed Cost Bounds - last 2 quarters only (Table A.18)\n")
    printFixCostBounds(tab.fblast,1,lab)

    @printf("\n4. Fixed Cost Bounds - nonstationary model (Table A.19)\n")
    printFixCostBounds(tab.fbns,1,lab)

    @printf("\n4. Fixed Cost Bounds - nonstationary model (Table A.20)\n")
    printFixCostBounds(tab.fbns,2,lab)

    @printf("\n5. Background conditions (Table A.7)\n")
    @printf("%10s  %10s  %10s  %10s  %10s\n","Operator","Δ","Diff E","Diff N","f̄̄")
    t9mean = mean(tab.t9,dims=3)
    t9std  = std(tab.t9,dims=3)
    for i=1:3
      @printf("%10s  %10.2f",lab.network[i],t9mean[i,1])
      for j=2:4
        @printf("  %10.0f",t9mean[i,j])
      end
      @printf("\n%10s  %10.2f","(se)",t9std[i,1])
      for j=2:4
        @printf("  %10.0f",t9std[i,j])
      end
      @printf("\n")
    end
  end

  ### Joint incentives
  payoff = mean(tab.payoff,dims=3)
  ΔΠ0boot = tab.payoff[7:9,1,:] - tab.payoff[10:12,4,:]
  ΔΠ1boot = tab.payoff[1:3,1,:] - tab.payoff[ 4: 6,4,:]
  ΔΠ0 = mean(ΔΠ0boot,dims=2)
  ΔΠ1 = mean(ΔΠ1boot,dims=2)
  sd0 = std(ΔΠ0boot,dims=2)
  sd1 = std(ΔΠ1boot,dims=2)
  @printf("\n6. Joint incentives for fighting brand (Table A.9/A.11)\n")
  @printf("%10s  %10s  %10s\n","Network","No enter","Enter")
  for i=1:3
    @printf("%10s  %10.0f  %10.0f\n",lab.network[i],ΔΠ0[i],ΔΠ1[i])
    @printf("%10s  %10.0f  %10.0f\n","(se)",sd0[i],sd1[i])
  end
  # ΣΔΠ0 = mean(sum(ΔΠ0boot,dims=1))
  # ΣΔΠ1 = mean(sum(ΔΠ1boot,dims=1))
  # sdΣ0 = std(sum(ΔΠ0boot,dims=1))
  # sdΣ1 = std(sum(ΔΠ1boot,dims=1))
  ΣΔΠ0 = sum(ΔΠ0boot,dims=1)
  ΣΔΠ1 = sum(ΔΠ1boot,dims=1)
  @printf("%10s  %10.0f  %10.0f\n","TotalInc",mean(ΣΔΠ0),mean(ΣΔΠ1))
  @printf("%10s  %10.0f  %10.0f\n","(se)",std(ΣΔΠ0),std(ΣΔΠ1))
  ΔΠf = tab.Πchannel[4,1,:] - tab.Πchannel[4,14,:]
  @printf("%10s  %10.0f  %10.0f\n","Free",0,mean(ΔΠf)/1e6)
  @printf("%10s  %10.0f  %10.0f\n","(se)",0,std(ΔΠf)/1e6)

  ### Unilateral incentives
  @printf("\n7. Unilateral incentives for fighting brand (Table 7/A.12)\n")
  t8 = zeros(3,2,length(eqb));
  t8[1,1,:] = tab.payoff[10,3,:] - tab.payoff[10,4,:];
  t8[2,1,:] = tab.payoff[11,2,:] - tab.payoff[11,4,:];
  t8[3,1,:] = tab.payoff[ 9,4,:] - tab.payoff[12,4,:];
  t8[1,2,:] = tab.payoff[ 1,2,:] - tab.payoff[ 1,1,:];
  t8[2,2,:] = tab.payoff[ 2,3,:] - tab.payoff[ 2,1,:];
  t8[3,2,:] = tab.payoff[ 6,1,:] - tab.payoff[ 3,1,:];
  @printf("%10s  %10s  %10s  %10s\n","Network","No enter","Enter","Diff")
  for i=1:3
    ΔΠuni = t8[i,1,:]+t8[i,2,:]
    @printf("%10s  %10.0f  %10.0f  %10.0f\n",
            lab.network[i],mean(t8[i,1,:]),mean(t8[i,2,:]),mean(ΔΠuni))
    @printf("%10s  %10.0f  %10.0f  %10.0f\n","(se)",
            std(t8[i,1,:]),std(t8[i,2,:]),std(ΔΠuni))
  end
  corrOrange = cor(t8[1,2,:],-t8[1,1,:])
  @printf("Corr(Π(NE),Π(E)) for Orange: %.4f\n",corrOrange)

  ### Surplus
  @printf("\n8. Surplus impacts (Table 9/A.14)\n")
  labRow = ["Free's entry"," Variety"," Price","Fight brands","Total"]
  tsplus = mean(tab.tsplus,dims=3)
  tsplus_se = std(tab.tsplus,dims=3)
  @printf("%12s  %10s  %10s  %10s\n","Source","Consumer","Producer","Total")
  for i=1:5
    @printf("%12s",labRow[i])
    for j=1:3
      @printf("  %10.0f",tsplus[i,j]);
    end
    @printf("\n")
    @printf("%12s","(se)")
    for j=1:3
      @printf("  %10.0f",tsplus_se[i,j]);
    end
    @printf("\n")
  end
  @printf("\n8a. Producer surplus change by groups\n")
  idx_inc = [1,2,3,5,6,7]
  Δpsall  = tab.Πchannel[idx_inc,14,:] - tab.Πchannel[idx_inc,16,:]
  Δpsfree = tab.Πchannel[4,14,:] - tab.Πchannel[4,16,:]
  @printf("%10s  %10.2f\n","Inc+MVNO",sum(mean(Δpsall,dims=2))/1e6)
  @printf("%10s  %10.2f\n","Free",mean(Δpsfree)/1e6)

  ### Payoff matrix
  @printf("\n9. Payoff matrix (Table A.6/A.10)\n")
  payoff = mean(tab.payoff,dims=3)
  for i=1:12
    @printf("%10.0f",payoff[i,1])
    for j=2:4
      @printf("  %10.0f",payoff[i,j]);
    end
    @printf("\n")
  end

  if vInt==false
    ### Price effects
    @printf("\n10. Price effects in retail market (Table A.8)\n")
    Δprice = mean(tab.Δprice,dims=3)
    @printf("%10s  %10s  %10s  %10s  %10s\n","Operator","Product","Retail p",
            "Change(%)","Change")
    for i in lab.prodIdx[1:12]
      @printf("%10s  %10s",lab.product[i],lab.prepost[i])
      for j=1:3
        @printf("  %10.2f",Δprice[i,j])
      end
      @printf("\n")
    end

    ### Price effects
    @printf("\n11. Lerner index & Ebitda in retail market\n")
    lerner = mean(tab.lerner,dims=2)
    ebitda = mean(tab.ebitda,dims=2)
    @printf("%10s  %10s  %10s\n","Operator","Lerner","Ebitda")
    for i=1:4
      @printf("%10s  %10.4f  %10.4f\n",lab.network[i],lerner[i],ebitda[i])
    end
  end
end

function getTables(
  eqb ::Array{Array{equilType,1}},
  vInt::Bool,
  dat ::datType,
  dim ::dimType,
  env ::envType,
  mat ::matType,
  pab ::Array{parType,1},
  set ::setType,
  sim ::simType,
  indx::indxType,
  mask::maskType)

  nBS      = length(eqb)
  payoff   = zeros(4*3,4,nBS)
  Πperiod  = zeros(dim.network,dim.qtr,16,nBS)
  t9       = zeros(3,4,nBS)
  t9last   = zeros(3,4,2,nBS)
  t9ns     = zeros(3,4,2,nBS)
  fb       = zeros(3,3,nBS)
  fblast   = zeros(3,3,2,nBS)
  fbns     = zeros(3,3,2,nBS)
  bound    = zeros(3,2,nBS)
  divR     = zeros(dim.product,dim.product,nBS)
  elasR    = zeros(dim.product,dim.product,nBS)
  elasW    = zeros(dim.product,dim.mvno,nBS)
  Πchannel = zeros(dim.network,16,nBS)
  csplus   = zeros(16,nBS)
  tsplus   = zeros(5,3,nBS)
  Δprice   = zeros(dim.product,3,nBS)
  lerner   = zeros(4,nBS)
  ebitda   = zeros(4,nBS)

  pm = zeros(dim.product,2,nBS)
  wm = zeros(dim.product,2,nBS)
  mcost = [zeros(dim.product,nBS) for k=1:2]
  for b=1:nBS
    eq = eqb[b][1]
    for i=1:dim.product
      idx = findall(eq.price[i,:].>0)
      pm[i,1,b] = mean(eq.price[i,idx])
      if mask.mno[i]
        pm[i,2,b] = mean(eq.price[i,idx] - eq.mcost[2][i,idx])
      elseif mask.mvno[i] & (!vInt)
        pm[i,2,b] = mean(eq.price[i,idx] - eq.wprice[i,idx])
        wm[i,1,b] = mean(eq.wprice[i,idx])
        wm[i,2,b] = mean(eq.wprice[i,idx] - eq.mcost[1][i,idx])
      elseif mask.mvno[i] & vInt
        pm[i,2,b] = mean(eq.price[i,idx] - eq.mcost[2][i,idx])
      end
      mcost[1][i,b] = mean(eq.mcost[1][i,idx])
      mcost[2][i,b] = mean(eq.mcost[2][i,idx])
    end

    for t=1:dim.qtr
      ∑D  = zeros(4)
      ∑mD = zeros(4)
      ∑pD = zeros(4)
      entries = findall(eq.price[1:dim.mno,t].>0)
      for i in entries
        j = mat.firm[i,2]
        ∑D[j] += eq.demand[i,t]
        ∑mD[j] += eq.markup[i,t,2]*eq.demand[i,t]
        ∑pD[j] += eq.price[i,t]*eq.demand[i,t]
      end
      for i in entries
        j = mat.firm[i,2]
        weight = eq.demand[i,t]/∑D[j]
        lerner[j,b] += eq.markup[i,t,2]/eq.price[i,t] * weight / dim.qtr
      end
      for j=1:4
        if ∑pD[j]>0
          ebitda[j,b] += ∑mD[j]/∑pD[j] / dim.qtr
        end
      end
    end
  end
  retail = reshape([mean(pm,dims=3) std(pm,dims=3)],dim.product,4)
  whsale = reshape([mean(wm,dims=3) std(wm,dims=3)],dim.product,4)

  tab = tabType(retail,whsale,mcost,payoff,Πperiod,t9,t9last,t9ns,fb,fblast,
                fbns,bound,divR,elasR,elasW,Πchannel,csplus,tsplus,Δprice,
                lerner,ebitda)

  for b=1:nBS
    ### Equilibrium payoffs
    r      = view(tab.payoff,:,:,b)
    Πchanl = view(tab.Πchannel,:,:,b)
    Πper   = view(tab.Πperiod,:,:,:,b)
    loadPayoff!(r,Πchanl,Πper,eqb[b],vInt,dim,env,mat,sim,indx,mask)

    #### Table: Background condition
    t9 = view(tab.t9,:,:,b)
    CN = r[10:12,4]
    DN = [r[10,3],r[11,2],r[9,4]]
    NN = r[7:9,1]
    CE = r[4:6,4]
    DE = [r[4,3],r[5,2],r[3,4]]
    NE = r[1:3,1]

    @. t9[:,2] = DE - NE  # Diff E
    @. t9[:,3] = DN - NN  # Diff N
    @. t9[:,1] = (DE-CE)/(DE-NE) - (DN-CN)/(DN-NN) # Δ

    # f double bar
    num = zeros(3); den = zeros(3);
    @. num = CE*(DN-NN)-CN*(DE-NE)+DE*NN-DN*NE
    @. den = DE-DN-NE+NN
    @. t9[:,4] = num/den

    ### Table 9: bounds on fixed costs
    fb = view(tab.fb,:,:,b)
    fb[1,1] = r[7,1] - r[10,4]
    fb[2,1] = r[8,1] - r[11,4]
    fb[3,1] = r[9,1] - r[12,4]

    fb[1,2] = r[7,1] - r[ 7,2]
    fb[2,2] = r[8,1] - r[ 8,3]
    fb[3,2] = r[9,1] - r[12,1]

    fb[:,3] = t9[:,4]

    ### Table A.19/A.20: bounds on fixed costs (2-period model)
    nNS = [2,4]
    for t=1:2
      fbns   = view(tab.fbns,  :,:,t,b)
      fblast = view(tab.fblast,:,:,t,b)
      t9ns   = view(tab.t9ns,  :,:,t,b)
      t9last = view(tab.t9last,:,:,t,b)
      getFixCostBounds!(fbns,fblast,t9ns,t9last,Πper,nNS[t],dim)
    end

    ### Collect surplus table
    for i=1:16
      tab.csplus[i,b] = sum(eqb[b][i].csplus)
    end
    tsplus = view(tab.tsplus,:,:,b)
    surplusTable!(tsplus,Πchanl,vInt,eqb[b],dat,dim,env,mat,pab[b],set,sim,
                  indx,mask)

    ### Price difference between endog & exo pricing when Free exits
    Δp = view(tab.Δprice,:,:,b)
    peq = partialBertrand(eqb[b],dat,dim,env,mat,pab[b],set,sim,indx,mask,vInt)
    postEntry = 2:13
    p0    = eqb[b][3].price[:,postEntry]
    rate  = (p0-peq[1].price[:,postEntry])./p0 .* 100
    level = (p0-peq[1].price[:,postEntry])
    rate[rate[:].==0] .= NaN
    level[level[:].==0] .= NaN
    for j=1:dim.product
      Δp[j,1] = mean(skipmissing(p0[j,:]))
      Δp[j,2] = mean(skipmissing(rate[j,:]))
      Δp[j,3] = mean(skipmissing(level[j,:]))
    end

  end # end of b loop

  return tab
end

function surplusTable!(
  tsplus,
  Πchannel,
  vInt::Bool,
  eq  ::Array{equilType,1},
  dat ::datType,
  dim ::dimType,
  env ::envType,
  mat ::matType,
  par ::parType,
  set ::setType,
  sim ::simType,
  indx::indxType,
  mask::maskType)

  initialize!(mat,dat,dim,par,set)

  peq = deepcopy(eq[14])
  @. peq.pdline  = eq[16].pdline
  @. peq.price   *= peq.pdline
  @. peq.mcost[1] = eq[1].mcost[1] * peq.pdline
  @. peq.mcost[2] = eq[1].mcost[2] * peq.pdline
  @. peq.wprice  *= peq.pdline

  Πchanlpeq = zeros(dim.network)

  for t=1:dim.qtr
    setupIndex!(sim,indx,peq.pdline,dim,mat,mask,t)
    collectEquil!(peq,t,vInt,peq.mcost,dim,mat,par,sim,indx,mask)
    collectΠ!(Πchanlpeq,peq,vInt,t,dim,mat,indx,mask)
  end

  tsplus[2,1] = sum(eq[14].csplus-peq.csplus)
  tsplus[3,1] = sum(peq.csplus-eq[16].csplus)
  tsplus[2,2] = sum(Πchannel[:,14] - Πchanlpeq)
  tsplus[3,2] = sum(Πchanlpeq - Πchannel[:,16])
  tsplus[1,1] = sum(eq[14].csplus - eq[16].csplus)
  tsplus[1,2] = sum(Πchannel[:,14] - Πchannel[:,16])
  tsplus[4,1] = sum(eq[1].csplus - eq[14].csplus)
  tsplus[4,2] = sum(Πchannel[:,1] - Πchannel[:,14])
  tsplus[5,1] = sum(eq[1].csplus - eq[16].csplus)
  tsplus[5,2] = sum(Πchannel[:,1] - Πchannel[:,16])

  tsplus[:,3] = sum(tsplus[:,1:2],dims=2)
  @. tsplus /= 1e6

end

function loadPayoff!(
  payoff,
  ΠChannel,
  Πperiod,
  eq  ::Array{equilType,1},
  vInt::Bool,
  dim ::dimType,
  env ::envType,
  mat ::matType,
  sim ::simType,
  indx::indxType,
  mask::maskType)

  for x=1:16
    Πchanl = view(ΠChannel,:,x)
    for t=1:dim.qtr
    # for t=1:dim.qtr-2  ## Last 2 quarters
      Πper = view(Πperiod,:,t,x)
      setupIndex!(sim,indx,eq[x].pdline,dim,mat,mask,t)
      collectΠ!(Πchanl,eq[x],vInt,t,dim,mat,indx,mask)
      collectΠ!(Πper,  eq[x],vInt,t,dim,mat,indx,mask)
    end # end of t loop
  end # end of x loop
  payoff .= reshape(ΠChannel[1:3,:],4*3,4)/1e6
  # payoff .= reshape(mean(Πperiod[1:3,:,:],dims=2),4*3,4)/1e6 * 13

end

function collectElast!(
  tab ::tabType,
  b   ::Int64,
  eq  ::equilType,
  dim ::dimType,
  lab ::labType,
  mat ::matType,
  par ::parType,
  sim ::simType,
  indx::indxType,
  mask::maskType)

  eR = Array{Union{Missing,Float64}}(missing,dim.product,dim.product,dim.qtr)
  eW = Array{Union{Missing,Float64}}(missing,dim.product,dim.mvno,dim.qtr)

  for t=1:dim.qtr
    eRs = view(eR,:,:,t)
    eWs = view(eW,:,:,t)
    elasticity!(eRs,eWs,t,eq,dim,mat,par,sim,indx,mask)
  end

  elasR = zeros(dim.product,dim.product)
  for j=1:dim.product
    for i=1:dim.product
      elasR[i,j] = mean(skipmissing(eR[i,j,:]))
    end
  end

  elasW = zeros(dim.product,dim.mvno)
  for j=1:dim.mvno
    for i=1:dim.product
      elasW[i,j] = mean(skipmissing(eW[i,j,:]))
    end
  end

  elasR = elasR[:,lab.prodIdx]
  elasR = elasR[lab.prodIdx,:]

  tab.elasR[:,:,b] = elasR
  tab.elasW[:,:,b] = elasW

end

function elasticity!(
  elasR,
  elasW,
  t   ::Int64,
  eq  ::equilType,
  dim ::dimType,
  mat ::matType,
  par ::parType,
  sim ::simType,
  indx::indxType,
  mask::maskType)

  p      = eq.price[:,t]
  wprice = eq.wprice[:,t]
  mcost  = [eq.mcost[1][:,t], eq.mcost[2][:,t]]
  demand!(sim,p,dim,mat,par,indx)
  markupR,markupW,markupV = getMarkup(p,wprice,mcost)
  x = passThrough(markupR,markupW,markupV,dim,mat,sim,indx,mask)
  w = wprice[mask.mvno]

  prodIdx = findall(indx.anyActive)
  for j=1:dim.mvno
    for i in prodIdx #i=1:dim.product
      if sim.D[i]>0
        elasW[i,j] = 0
        for l in prodIdx #l=1:dim.product
          elasW[i,j] += sim.∂D∂p[i,l]*x[l,j]*w[j]/sim.D[i]
        end
      end
    end
  end

  for j in prodIdx #j=1:dim.product
    for i in prodIdx #i=1:dim.product
      if sim.D[i]>0
        elasR[i,j] = sim.∂D∂p[i,j]*p[j]/sim.D[i]
      else
        elasR[i,j] = 0
      end
    end
  end

end

function getDivRatio(
  free::Int64,
  eq  ::Array{equilType,1},
  dim ::dimType,
  lab ::labType,
  mat ::matType,
  par ::parType,
  sim ::simType,
  indx::indxType,
  mask::maskType)

  # dF = Array{Union{Missing,Float64}}(missing,dim.product,dim.product,dim.qtr)
  dF = zeros(dim.product,dim.product,dim.qtr)
  if free==1
    id = 1
  elseif free==0
    id = 3
  end
  for t=1:dim.qtr
    setupIndex!(sim,indx,eq[id].pdline,dim,mat,mask,t)
    p = eq[id].price[:,t]
    D,cs = demandOnly(sim,p,dim,mat,par,indx)
    prodIdx = findall(indx.anyActive)
    for j in prodIdx
      p1 = copy(p)
      p1[j] = p1[j] + 1.0
      D1,cs = demandOnly(sim,p1,dim,mat,par,indx)
      dDown = -(D1[j] - D[j])
      for i in prodIdx
        dD = D1[i] - D[i]
        dF[i,j,t] = dD/dDown
      end
      if sum(skipmissing(dF[:,j,t]))>0
        for i in eachindex(dF[:,j,t])
          @printf("%.4f\n",dF[i,j,t])
        end
      end
    end
    # println(sum(dF[:,:,t],dims=1))
  end

  divRatio = zeros(dim.product,dim.product)
  for j=1:dim.product
    for i=1:dim.product
      divRatio[i,j] = mean(skipmissing(dF[i,j,:]))*100
    end
    if divRatio[j,j] > -100.0
      dRatioOwn = abs(divRatio[j,j])
      for i=1:dim.product
        divRatio[i,j] /= (dRatioOwn/100)
      end
    end
  end
  divRatio = divRatio[:,lab.prodIdx]
  divRatio = divRatio[lab.prodIdx,:]
  return divRatio

end

function reportElast!(
  tab  ::tabType,
  vInt ::Bool,
  env  ::envType)

  m    = env.model
  mcid = env.mcid
  mcv  = env.mc_mvno
  opt  = false
  ms15 = false

  dat,dim,env,mat,par0,set,sim,indx,mask = readData(m,mcid,0,0,vInt,opt,ms15,mcv)
  eq = allocateEquil(dim)
  sim = allocateSim(dim)
  initialize!(mat,dat,dim,par0,set)
  marginalCost!(eq,dat,dim,env,mat,par0,set,sim,indx,mask,vInt)
  bertrandSolve!(eq,dat,dim,env,mat,par0,set,sim,indx,mask,vInt,1)

  lab = setupLab()
  collectElast!(tab,1,eq[1],dim,lab,mat,par0,sim,indx,mask)

  elasOwn = tr(tab.elasR[:,:,1])/dim.product
  elasReg = zeros(dim.region,dim.qtr)

  for t=1:dim.qtr
    p0 = copy(eq[1].price[:,t])
    setupIndex!(sim,indx,eq[1].pdline,dim,mat,mask,t)
    demandOnly(sim,p0,dim,mat,par0,indx)
    tms0 = sum(sim.ms,dims=1)

    p1 = copy(p0)
    @. p1[indx.anyActive] *= 1.01
    demandOnly(sim,p1,dim,mat,par0,indx)
    tms1 = sum(sim.ms,dims=1)
    for i=1:dim.region
      elasReg[i,t] = (tms1[i]-tms0[i])/tms0[i]*100
    end
  end
  elasMarket = mean(elasReg)

  @printf("\nElasticity table\n")
  @printf("%10s  %10s\n","Own","Market")
  @printf("%10.3f  %10.3f\n",elasOwn,elasMarket)

  return elasOwn, elasMarket
end

function fastElasticity(
  model::Int64,
  mcid ::Int64,
  vInt ::Bool,
  ver  ::String="824")

  dat,dim,env,mat,par0,set,sim,indx,mask =
    readData(model,mcid,0,0,vInt,false,false,0;ver)
  eq = allocateEquil(dim)
  sim = allocateSim(dim)
  initialize!(mat,dat,dim,par0,set)
  lab = setupLab()

  eR = Array{Union{Missing,Float64}}(missing,dim.product,dim.product,dim.qtr)

  for t=1:dim.qtr
    setupIndex!(sim,indx,eq[1].pdline,dim,mat,mask,t)
    demand!(sim,sim.p,dim,mat,par0,indx)
    elasRt = view(eR,:,:,t)

    prodIdx = findall(indx.anyActive)
    for j in prodIdx #j=1:dim.product
      for i in prodIdx #i=1:dim.product
        if sim.D[i]>0
          elasRt[i,j] = sim.∂D∂p[i,j]*sim.p[j]/sim.D[i]
        else
          elasRt[i,j] = 0
        end
      end
    end
  end

  elasR = zeros(dim.product,dim.product)
  for j=1:dim.product
    for i=1:dim.product
      elasR[i,j] = mean(skipmissing(eR[i,j,:]))
    end
  end

  elasR = elasR[:,lab.prodIdx]
  elasR = elasR[lab.prodIdx,:]
  elasOwn = tr(elasR)/dim.product

  @printf("\nOwn national demand elasticity (average):\n")
  @printf("%10.3f\n",elasOwn)

  return elasR, elasOwn

end

function getFixCostBounds!(
  fbns,
  fblast,
  t9ns,
  t9last,
  Πper,
  nNS::Int64, # length of nonstationary time periods
  dim::dimType)

  t0 = dim.qtr - nNS

  ### Bounds on fixed costs (2-period model)
  Π0 = reshape(mean(Πper[1:3,2:t0,:], dims=2)[:,1,:],3*4,4) * 12/1e6
  Π1 = reshape(mean(Πper[1:3,(t0+1):dim.qtr,:], dims=2)[:,1,:],3*4,4) * 12/1e6

  CN = hcat(Π0[10:12,4], Π1[10:12,4])
  DN = hcat([Π0[10,3],Π0[11,2],Π0[9,4]], [Π1[10,3],Π1[11,2],Π1[9,4]])
  NN = hcat(Π0[7:9,1], Π1[7:9,1])
  Π̂N = hcat([Π0[7,2],Π0[8,3],Π0[12,1]], [Π1[7,2],Π1[8,3],Π1[12,1]])

  CE = hcat(Π0[4:6,4], Π1[4:6,4])
  DE = hcat([Π0[4,3],Π0[5,2],Π0[3,4]], [Π1[4,3],Π1[5,2],Π1[3,4]])
  NE = hcat(Π0[1:3,1], Π1[1:3,1])

  fbns[:,1] = NN[:,2] - CN[:,2]
  fbns[:,2] = NN[:,2] - Π̂N[:,2]

  aE = CE[:,2] - CE[:,1]
  aN = CN[:,2] - CN[:,1]
  A  = DN[:,1] - NN[:,2] + aN
  B  = DE[:,1] - NE[:,2] + aE

  # δᴱ - δᴺ = Coef_f*f + Δ0 > 0
  Coef_f = zeros(3); Δ0 = zeros(3)
  @. Coef_f = 1/A - 1/B
  @. Δ0 = (DE[:,1]-CE[:,1])/B - (DN[:,1]-CN[:,1])/A
  @. fbns[:,3] = -Δ0./Coef_f

  ### Bounds on fixed costs (last period only)
  @. fblast = fbns

  @. t9last[:,2] = (DE - NE)[:,2]  # Diff E
  @. t9last[:,3] = (DN - NN)[:,2]  # Diff N
  @. t9last[:,1] = ((DE-CE)/(DE-NE) - (DN-CN)/(DN-NN))[:,2] # Δ
  num = zeros(3,2); den = zeros(3,2);
  @. num = CE*(DN-NN)-CN*(DE-NE)+DE*NN-DN*NE
  @. den = DE-DN-NE+NN
  @. t9last[:,4] = (num/den)[:,2]

  @. fblast[:,3] = t9last[:,4]

  #### Table: Background condition (2-period model)
  @. t9ns[:,1] = Δ0
  @. t9ns[:,2] = B
  @. t9ns[:,3] = A
  @. t9ns[:,4] = fbns[:,3]

end

end

max_stars_repo_name: yutec/FightBrand

max_issues_repo_name: yutec/FightBrand

max_forks_repo_name: yutec/FightBrand


text:
# texture and surface handling

# NOTE: the API for texture support is not final yet. some thoughts:
#
# - instead of CuTextureArray, use CuArray with an ArrayBuffer. This array could then
#   adapt to a CuTexture, or do the same for CuDeviceArray.

#
# Texture array
#

export CuTextureArray

"""
    CuTextureArray{T,N}(undef, dims)

`N`-dimensional dense texture array with elements of type `T`. These arrays are optimized
for texture fetching, and are only meant to be used as a source for [`CuTexture{T,N,P}`](@ref)
objects.

!!! warning
    Experimental API. Subject to change without deprecation.
"""
mutable struct CuTextureArray{T,N}
    buf::Mem.ArrayBuffer{T}
    dims::Dims{N}
    ctx::CuContext

    @doc """
        CuTextureArray{T,N}(undef, dims)

    Construct an uninitialized texture array of `N` dimensions specified in the `dims`
    tuple, with elements of type `T`. Use `Base.copyto!` to initialize this texture array,
    or use constructors that take a non-texture array to do so automatically.

    !!! warning
        Experimental API. Subject to change without deprecation.
    """
    function CuTextureArray{T,N}(::UndefInitializer, dims::Dims{N}) where {T,N}
        buf = Mem.alloc(Mem.Array{T}, dims)
        t = new{T,N}(buf, dims, context())
        finalizer(unsafe_destroy!, t)
        return t
    end
end

function unsafe_destroy!(t::CuTextureArray)
    @context! skip_destroyed=true t.ctx Mem.free(t.buf)
end

Base.unsafe_convert(T::Type{CUarray}, t::CuTextureArray) = Base.unsafe_convert(T, t.buf)

## array interface

Base.size(tm::CuTextureArray) = tm.dims
Base.length(tm::CuTextureArray) = prod(size(tm))

Base.eltype(tm::CuTextureArray{T,N}) where {T,N} = T

Base.sizeof(tm::CuTextureArray) = sizeof(eltype(tm)) * length(tm)

Base.pointer(t::CuTextureArray) = t.buf.ptr

## interop with other arrays

"""
    CuTextureArray(A::AbstractArray)

Allocate and initialize a texture buffer from host memory in `A`.

!!! warning
    Experimental API. Subject to change without deprecation.
"""
@inline function CuTextureArray{T,N}(xs::AbstractArray{<:Any,N}) where {T,N}
  A = CuTextureArray{T,N}(undef, size(xs))
  copyto!(A, convert(Array{T}, xs))
  return A
end

"""
    CuTextureArray(A::CuArray)

Allocate and initialize a texture buffer from device memory in `A`.

!!! warning
    Experimental API. Subject to change without deprecation.
"""
@inline function CuTextureArray{T,N}(xs::CuArray{<:Any,N}) where {T,N}
  A = CuTextureArray{T,N}(undef, size(xs))
  copyto!(A, convert(CuArray{T}, xs))
  return A
end

# idempotency
CuTextureArray{T,N}(xs::CuTextureArray{T,N}) where {T,N} = xs

CuTextureArray(A::AbstractArray{T,N}) where {T,N} = CuTextureArray{T,N}(A)

## memory operations

function Base.copyto!(dst::CuTextureArray{T,1}, src::Union{Array{T,1}, CuArray{T,1}}) where {T}
    size(dst) == size(src) || throw(DimensionMismatch("source and destination sizes must match"))
    Base.unsafe_copyto!(pointer(dst), pointer(src), length(dst))
    return dst
end

function Base.copyto!(dst::CuTextureArray{T,2}, src::Union{Array{T,2}, CuArray{T,2}}) where {T}
    size(dst) == size(src) || throw(DimensionMismatch("source and destination sizes must match"))
    Mem.unsafe_copy2d!(pointer(dst), Mem.Array,
                       pointer(src), isa(src, Array) ? Mem.Host : Mem.Device,
                       size(dst)...)
    return dst
end

function Base.copyto!(dst::CuTextureArray{T,3}, src::Union{Array{T,3}, CuArray{T,3}}) where {T}
    size(dst) == size(src) || throw(DimensionMismatch("source and destination sizes must match"))
    Mem.unsafe_copy3d!(pointer(dst), Mem.Array,
                       pointer(src), isa(src, Array) ? Mem.Host : Mem.Device,
                       size(dst)...)
    return dst
end

#
# Texture objects
#

export CuTexture

@enum_without_prefix CUaddress_mode CU_TR_

Base.convert(::Type{CUfilter_mode}, ::NearestNeighbour)      = CU_TR_FILTER_MODE_POINT
Base.convert(::Type{CUfilter_mode}, ::LinearInterpolation)   = CU_TR_FILTER_MODE_LINEAR
Base.convert(::Type{CUfilter_mode}, ::CubicInterpolation)    = CU_TR_FILTER_MODE_LINEAR

"""
    CuTexture{T,N,P}

`N`-dimensional texture object with elements of type `T`. These objects do not store data
themselves, but are bounds to another source of device memory. Texture objects can be passed
to CUDA kernels, where they will be accessible through the [`CuDeviceTexture`](@ref) type.

!!! warning
    Experimental API. Subject to change without deprecation.
"""
mutable struct CuTexture{T,N,P} <: AbstractArray{T,N}
    parent::P
    handle::CUtexObject

    interpolation::TextureInterpolationMode
    normalized_coordinates::Bool

    ctx::CuContext

    @doc """
        CuTexture{T,N,P}(parent::P; address_mode, filter_mode, normalized_coordinates)

    Construct a `N`-dimensional texture object with elements of type `T` as stored in
    `parent`.

    Several keyword arguments alter the behavior of texture objects:
    - `address_mode` (wrap, *clamp*, mirror): how out-of-bounds values are accessed. Can be
      specified as a value for all dimensions, or as a tuple of `N` entries.
    - `interpolation` (*nearest neighbour*, linear, bilinear): how non-integral indices are
      fetched. Nearest-neighbour fetches a single value, others interpolate between multiple.
    - `normalized_coordinates` (true, *false*): whether indices are expected to fall in the
      normalized `[0:1)` range.

    !!! warning Experimental API. Subject to change without deprecation.
    """
    function CuTexture{T,N,P}(parent::P;
                              address_mode::Union{CUaddress_mode,NTuple{N,CUaddress_mode}}=ADDRESS_MODE_CLAMP,
                              interpolation::TextureInterpolationMode=NearestNeighbour(),
                              normalized_coordinates::Bool=false) where {T,N,P}
        resDesc_ref = CUDA_RESOURCE_DESC(parent)

        flags = 0x0
        if normalized_coordinates
            flags |= CU_TRSF_NORMALIZED_COORDINATES
        end
        if eltype(T) <: Integer
            flags |= CU_TRSF_READ_AS_INTEGER
        end

        # we always need 3 address modes
        if !(address_mode isa Tuple)
            address_mode = ntuple(_->address_mode, N)
        end
        address_mode = tuple(address_mode..., ntuple(_->ADDRESS_MODE_CLAMP, 3 - N)...)

        if interpolation == CubicInterpolation()
            normalized_coordinates && throw(ArgumentError("Cubic interpolation is only supported with non-normalized coordinates"))
            N>2 && throw(ArgumentError("Cubic interpolation is only implemented for 1 and 2D textures"))
        end
        filter_mode = convert(CUfilter_mode, interpolation)

        texDesc_ref = Ref(CUDA_TEXTURE_DESC(
            address_mode, # addressMode::NTuple{3, CUaddress_mode}
            filter_mode, # filterMode::CUfilter_mode
            flags, # flags::UInt32
            1, # maxAnisotropy::UInt32
            filter_mode, # mipmapFilterMode::CUfilter_mode
            0, # mipmapLevelBias::Cfloat
            0, # minMipmapLevelClamp::Cfloat
            0, # maxMipmapLevelClamp::Cfloat
            ntuple(_->Cfloat(zero(eltype(T))), 4), # borderColor::NTuple{4, Cfloat}
            ntuple(_->Cint(0), 12)))

        texObject_ref = Ref{CUtexObject}(0)
        cuTexObjectCreate(texObject_ref, resDesc_ref, texDesc_ref, C_NULL)

        t = new{T,N,P}(parent, texObject_ref[], interpolation, normalized_coordinates, context())
        finalizer(unsafe_destroy!, t)
        return t
    end
end

function CUDA_RESOURCE_DESC(texarr::CuTextureArray)
    # FIXME: manual construction due to invalid padding (JuliaInterop/Clang.jl#238)
    # res = Ref{ANONYMOUS1_res}()
    # unsafe_store!(Ptr{CUarray}(pointer_from_objref(res)), texarr.handle)
    # resDesc_ref = Ref(CUDA_RESOURCE_DESC(
    #     CU_RESOURCE_TYPE_ARRAY, # resType::CUresourcetype
    #     res[], # res::ANONYMOUS1_res
    #     0 # flags::UInt32
    # ))
    resDesc_ref = Ref((CU_RESOURCE_TYPE_ARRAY, # resType::CUresourcetype
                       pointer(texarr), # 1 x UInt64
                       ntuple(_->Int64(0), 15), # 15 x UInt64
                       UInt32(0)))
    return resDesc_ref
end

function CUDA_RESOURCE_DESC(arr::CuArray{T,N}) where {T,N}
    # TODO: take care of allowed pitches
    1 <= N <= 2 || throw(ArgumentError("Only 1 or 2D CuArray objects can be wrapped in a texture"))

    format = convert(CUarray_format, eltype(T))
    channels = nchans(T)

    # FIXME: manual construction due to invalid padding (JuliaInterop/Clang.jl#238)
    resDesc_ref = Ref(((N == 1 ? CU_RESOURCE_TYPE_LINEAR : CU_RESOURCE_TYPE_PITCH2D), # resType::CUresourcetype
                       pointer(arr), # 1 x UInt64 (CUdeviceptr)
                       format, # 1/2 x UInt64 (CUarray_format)
                       UInt32(channels), # 1/2 x UInt64
                       (N == 2 ? size(arr, 1) : size(arr, 1) * sizeof(T)), # 1 x UInt64 nx
                       (N == 2 ? size(arr, 2) : 0), # 1 x UInt64 ny
                       (N == 2 ? size(arr, 1) * sizeof(T) : 0), # 1 x UInt64 pitch
                       ntuple(_->Int64(0), 11), # 11 x UInt64
                       UInt32(0)))
    return resDesc_ref
end

# allow passing mismatching refs where we're expecting a CUDA_RESOURCE_DESC
Base.unsafe_convert(::Type{Ptr{CUDA_RESOURCE_DESC}}, ref::Base.RefValue{T}) where {T} =
    convert(Ptr{CUDA_RESOURCE_DESC}, Base.unsafe_convert(Ptr{T}, ref))

function unsafe_destroy!(t::CuTexture)
    @context! skip_destroyed=true t.ctx cuTexObjectDestroy(t)
end

Base.convert(::Type{CUtexObject}, t::CuTexture) = t.handle

Base.parent(tm::CuTexture) = tm.parent

## array interface

Base.size(tm::CuTexture) = size(tm.parent)
Base.sizeof(tm::CuTexture) = Base.elsize(x) * length(x)

Base.show(io::IO, t::CuTexture{T,1}) where {T} =
    print(io, "$(length(t))-element $(nchans(T))-channel CuTexture(::$(typeof(parent(t)).name.name)) with eltype $T")
Base.show(io::IO, t::CuTexture{T}) where {T} =
    print(io, "$(join(size(t), '×')) $(nchans(T))-channel CuTexture(::$(typeof(parent(t)).name.name)) with eltype $T")

Base.show(io::IO, mime::MIME"text/plain", t::CuTexture) = show(io, t)

## interop with other arrays

CuTexture{T,N}(A::Union{CuTextureArray{T,N},CuArray{T,N}}; kwargs...) where {T,N} =
    CuTexture{T,N,typeof(A)}(A; kwargs...)

"""
    CuTexture(x::CuTextureArray{T,N})

Create a `N`-dimensional texture object withelements of type `T` that will be read from `x`.

!!! warning
    Experimental API. Subject to change without deprecation.
"""
CuTexture(x::CuTextureArray{T,N}; kwargs...) where {T,N} =
    CuTexture{T,N}(x; kwargs...)

"""
    CuTexture(x::CuArray{T,N})

Create a `N`-dimensional texture object that reads from a `CuArray`.

Note that it is necessary the their memory is well aligned and strided (good pitch).
Currently, that is not being enforced.

!!! warning
    Experimental API. Subject to change without deprecation.
"""
CuTexture(x::CuArray{T,N}; kwargs...) where {T,N} =
    CuTexture{T,N}(x; kwargs...)

memory_source(::Any) = error("Unknown texture source $(typeof(t))")
memory_source(::CuArray) = LinearMemory()
memory_source(::CuTextureArray) = ArrayMemory()

Adapt.adapt_storage(::Adaptor, t::CuTexture{T,N}) where {T,N} =
    CuDeviceTexture{T,N,typeof(memory_source(parent(t))),
                    t.normalized_coordinates, typeof(t.interpolation)}(size(t), t.handle)


text:
mutable struct OggEncoder
    streams::Dict{Clong,OggStreamState}
    pages::Dict{Clong,Vector{Vector{UInt8}}}

    function OggEncoder()
        enc = new(Dict{Clong,OggStreamState}(), Dict{Clong,Vector{Vector{UInt8}}}())

        # This seems to be causing problems.  :(
        #=
        finalizer(enc) do enc
            for serial in keys(enc.streams)
                ogg_stream_destroy(enc.streams[serial])
            end
        end
        =#

        return enc
    end
end

function show(io::IO, x::OggEncoder)
    num_streams = length(x.streams)
    if num_streams != 1
        write(io, "OggEncoder with $num_streams streams")
    else
        write(io, "OggEncoder with 1 stream")
    end
end

function ogg_stream_packetin(enc::OggEncoder, serial::Clong, data::Vector{UInt8}, packet_idx, last_packet::Bool, granulepos::Int64)
    if !haskey(enc.streams, serial)
        streamref = Ref{OggStreamState}(OggStreamState())
        status = ccall((:ogg_stream_init,libogg), Cint, (Ref{OggStreamState}, Cint), streamref, serial)
        if status != 0
            error("ogg_stream_init() failed: Unknown failure")
        end
        enc.streams[serial] = streamref[]

        # Also initialize enc.pages for this serial
        enc.pages[serial] = Vector{Vector{UInt8}}()
    end

    # Build ogg_packet structure
    packet = OggPacket( pointer(data), length(data), packet_idx == 0,
                        last_packet, granulepos, packet_idx )

    streamref = Ref{OggStreamState}(enc.streams[serial])
    packetref = Ref{OggPacket}(packet)
    status = ccall((:ogg_stream_packetin,libogg), Cint, (Ref{OggStreamState}, Ref{OggPacket}), streamref, packetref)
    enc.streams[serial] = streamref[]
    if status == -1
        error("ogg_stream_packetin() failed: Unknown failure")
    end
    return nothing
end

function ogg_stream_pageout(enc::OggEncoder, serial::Clong)
    if !haskey(enc.streams, serial)
        return nothing
    end
    streamref = Ref{OggStreamState}(enc.streams[serial])
    pageref = Ref{OggPage}(OggPage())
    status = ccall((:ogg_stream_pageout,libogg), Cint, (Ref{OggStreamState}, Ref{OggPage}), streamref, pageref)
    enc.streams[serial] = streamref[]
    if status == 0
        return nothing
    else
        return pageref[]
    end
end

function ogg_stream_flush(enc::OggEncoder, serial::Clong)
    if !haskey(enc.streams, serial)
        return nothing
    end
    streamref = Ref{OggStreamState}(enc.streams[serial])
    pageref = Ref{OggPage}(OggPage())
    status = ccall((:ogg_stream_flush,libogg), Cint, (Ref{OggStreamState}, Ref{OggPage}), streamref, pageref)
    enc.streams[serial] = streamref[]
    if status == 0
        return nothing
    else
        return pageref[]
    end
end

"""
    encode_all_packets(enc, packets, granulepos)

Feed all packets (with their corresponding granule positions) into encoder `enc`.
"""
function encode_all_packets(enc::OggEncoder, packets::Dict{Clong,Vector{Vector{UInt8}}}, granulepos::Dict{Clong,Vector{Int64}})
    pages = Vector{Vector{UInt8}}()

    # We're just going to chain streams together, not interleave them
    for serial in keys(packets)
        # Shove all packets into their own stream
        for packet_idx in 1:length(packets[serial])
            eos = packet_idx == length(packets[serial])
            ogg_stream_packetin(enc, serial, packets[serial][packet_idx], packet_idx - 1, eos, granulepos[serial][packet_idx])

            # A granulepos of zero signifies a header packet, which should be
            # flushed into its own page.  We know the header packets always
            # fit within a single page too, so we don't bother with the typical
            # while loop that would dump out excess data into its own page.
            if granulepos[serial][packet_idx] == 0
                page = ogg_stream_flush(enc, serial)
                push!(pages, read(page))
            end
        end

        # Flush and then pull all pages out.  Note that in a streaming situation,
        # normally you would call ogg_stream_pageout(enc, serial) periodically, as
        # you feed data into the stream with ogg_stream_packetin(), however in this
        # case all the data has already been written, so we can simply blithely call
        # ogg_stream_flush() and then read pages out until we run out of pages.
        # Note further that we could call ogg_stream_flush() from within the while
        # loop, however it is unnecessary as a single flush is enough to get the
        # pages moving, and this is better for code coverage purposes.
        page = ogg_stream_flush(enc, serial)
        while page != nothing
            push!(pages, read(page))
            page = ogg_stream_pageout(enc, serial)
        end
    end

    # Flatten pages
    return pages
end

function save(fio::IO, packets::Dict{Clong,Vector{Vector{UInt8}}}, granulepos::Dict{Clong,Vector{Int64}})
    enc = OggEncoder()
    pages = encode_all_packets(enc, packets, granulepos)
    for page in pages
        write(fio, page)
    end
end

function save(file_path::Union{File{format"OGG"},AbstractString}, packets::Dict{Clong,Vector{Vector{UInt8}}}, granulepos::Dict{Clong,Vector{Int64}})
    open(file_path, "w") do fio
        save(fio, packets, granulepos)
    end
end

# Convenience save() function for single-stream Ogg files, assigns an arbitrary stream ID
function save(file_path::Union{File{format"OGG"},AbstractString,IO}, packets::Vector{Vector{UInt8}}, granulepos::Vector{Int64})
    return save(file_path, Dict(Clong(1) => packets), Dict(Clong(1) => granulepos))
end


text:
gencsl(tbl; dw=1000, cf=:chrom, sf=:max_motif_start, ef=:max_motif_stop, strf=:max_motif_strand) = (tbl[cf], [div(a + b, 2) .+ (-dw:dw) for (a, b) in zip(tbl[sf], tbl[ef])], tbl[strf]);

@inline inc_meta_frag!(V, mpl, mpr, w) =  V[intersect(mpl:mpr, 1:size(V, 1))] .+= w

@inline inc_meta_cut!(V, mpl, mpr, w)  =  (inc_meta_index!(V, mpl, w) ; inc_meta_index!(V, mpr, w))
@inline inc_meta_mid!(V, mpl, mpr, w)  =  inc_meta_index!(V, div(mpl + mpr, 2), w)
@inline inc_meta_index!(V, i::Int, w)  =  (1 ≤ i ≤ size(V, 1)) && (V[i] += w)

@inline inc_heat_frag!(V, k, mpl, mpr, w) =  V[intersect(mpl:mpr, 1:size(V, 1)), k] .+= w
@inline inc_heat_cut!(V, k, mpl, mpr, w)  =  (inc_heat_index!(V, k, mpl, w) ; inc_heat_index!(V, k, mpr, w))
@inline inc_heat_left_cut!(V, k, mpl, mpr, w)  =  inc_heat_index!(V, k, mpl, w)
@inline inc_heat_mid!(V, k, mpl, mpr, w)  =  inc_heat_index!(V, k, div(mpl + mpr, 2), w)
@inline @inbounds inc_heat_index!(V, k, i::Int, w)  =  (1 ≤ i ≤ size(V, 1)) && (V[i, k] += w)

inc_meta_mid_width(bw) = (V, mpl, mpr, w) -> begin
    mid = div(mpl + mpr, 2)
    inc_meta_frag!(V, mid - bw, mid + bw, w)
end

inc_heat_mid_width(bw) = (V, k, mpl, mpr, w) -> begin
    mid = div(mpl + mpr, 2)
    bc = intersect((mid - bw):(mid + bw), 1:size(V, 1))
    V[bc, k] .+= w
end

inc_heat_atac_cut_width(bw) = (V, k, mpl, mpr, w, ao=4) -> begin
    inc_heat_frag!(V, k, mpl + ao - bw, mpl + ao + bw, w)
    inc_heat_frag!(V, k, mpr - ao - bw, mpr - ao + bw, w)
end

inc_meta_atac_cut_width(bw) = (V, mpl, mpr, w, ao=4) -> begin
    inc_meta_frag!(V, mpl + ao - bw, mpl + ao + bw, w)
    inc_meta_frag!(V, mpr - ao - bw, mpr - ao + bw, w)
end

@inline inc_heat_atac_cut!(V, k, mpl, mpr, w, ao=4)   =  (inc_heat_index!(V, k, mpl + ao, w) ; inc_heat_index!(V, k, mpr - ao, w))
@inline inc_heat_atac_frag!(V, k, mpl, mpr, w, ao=4)  =  inc_heat_frag!(V, k, mpl+ao, mpr-ao, w);
@inline inc_meta_atac_frag!(V, mpl, mpr, w, ao=4)     =  inc_meta_frag!(V, mpl+ao, mpr-ao, w);

function libsize(FM, dataentry, fn=identity, minfrag=0, maxfrag=1000000)

    σ = 0.0
    for i = 1:size(FM.FM, 2)
        fp = FM.FM[2, i] - FM.FM[1, i] + 1
        σ += ifelse(minfrag ≤ fp ≤ maxfrag, fn(FM.FM[dataentry, i]), 0.0)
    end
    σ
end

function libsize(FM, dataentry, fn=identity)
    σ = 0
    @inbounds for i = 1:size(FM.FM, 2)
        σ += fn(FM.FM[dataentry, i])
    end
    σ
end

function frag_heatmap_strand_norm(chroms, locations, strands, FM::Vector{FragMatrixSingle{T}}, inc_fun=inc_heat_mid!, fraglength=120, data_entry=4, fn=identity; norm_scale=1e+6, norm_type=:mpm, show_progress=true) where T

    fw = length(locations[1])
    H = zeros(Float64, fw, length(chroms))
    p = ProgressMeter.Progress(length(chroms)*length(FM))
    for F in FM
        TH = zeros(Float64, fw, length(chroms))
        frag_heatmap_strand!(TH, p, chroms, locations, strands, F, inc_fun, fraglength, data_entry, fn, show_progress=show_progress)
        if norm_type == :mpm
            H .+= norm_scale.*TH./F.totalecfrags[data_entry - 3]
        elseif norm_type == :mpm_f
            t = sum(fn, view(F.FM, data_entry, :))
            H .+= norm_scale.*TH./t
	    elseif norm_type == :frag_region
	        H .+= norm_scale.*TH/sum(TH)
        elseif norm_type == :total
            H += TH
        else
            H += TH
        end
    end
    if norm_type == :total
        rmul!(H, norm_scale/sum(F.totalecfrags[data_entry - 3] for F in FM))
    else
        rmul!(H, 1/length(FM))
    end
    H
end

function frag_heatmap_strand_norm(chroms, locations, strands, FM::Vector{FragMatrixPair{T}}, inc_fun=inc_heat_mid!, minfragsize=0, maxfragsize=500, data_entry=4, fn=identity; norm_scale=1e+6, norm_type=:mpm, show_progress=true) where T

    fw = length(locations[1])
    H = zeros(Float64, fw, length(chroms))
    p = ProgressMeter.Progress(length(chroms)*length(FM))
    for F in FM
        TH = zeros(Float64, fw, length(chroms))
        frag_heatmap_strand!(TH, p, chroms, locations, strands, F, inc_fun, minfragsize, maxfragsize, data_entry, fn, show_progress=show_progress)
        if norm_type == :mpm
            H .+= norm_scale.*TH./F.totalecfrags[data_entry - 3]
        elseif norm_type == :mpm_f
            t = sum(fn, view(F.FM, data_entry, :))
            H .+= norm_scale.*TH./t
        elseif norm_type == :frag
            σ = libsize(F, data_entry, fn, minfragsize, maxfragsize)
            H .+= norm_scale.*TH./σ
	elseif norm_type == :frag_region
	    H .+= norm_scale.*TH/sum(TH)
        elseif norm_type == :total
            H += TH
        else
            H += TH
        end
    end
    if norm_type == :total
        rmul!(H, 1e+6/sum(F.totalecfrags[data_entry - 3] for F in FM))
    else
        rmul!(H, 1/length(FM))
    end
    H
end

function frag_heatmap_strand(chroms, locations, strands, FM::Vector{FragMatrixPair{T}}, inc_fun=inc_heat_mid!, minfragsize=0, maxfragsize=500, data_entry=4, fn=identity; show_progress=true) where T

    fw = length(locations[1])
    H = zeros(Float64, fw, length(chroms))
    p = Progress(length(chroms)*length(FM))
    for F in FM
        frag_heatmap_strand!(H, p, chroms, locations, strands, F, inc_fun, minfragsize, maxfragsize, data_entry, fn, show_progress=show_progress)
    end
    H
end

function frag_heatmap_strand(chroms, locations, strands, FM::Vector{FragMatrixSingle{T}}, inc_fun=inc_heat_mid!, fraglength=120, data_entry=4, fn=identity; show_progress=true) where T

    fw = length(locations[1])
    H = zeros(Float64, fw, length(chroms))
    p = Progress(length(chroms)*length(FM))
    for F in FM
        frag_heatmap_strand!(H, p, chroms, locations, strands, F, inc_fun, fraglength, data_entry, fn, show_progress=show_progress)
    end
    H
end

function frag_heatmap_strand!(H, p, chroms, locations, strands, FM::FragMatrixPair{T}, inc_fun=inc_heat_mid!, minfragsize=0, maxfragsize=500, data_entry=4, fn=identity; show_progress=true) where T

    fw = length(locations[1])

    for (k, (c, l, s)) in enumerate(zip(chroms, locations, strands))
        show_progress && next!(p)
        V = get_frags(c, l, FM)
        isempty(V) && continue
        for i = 1:size(V, 2)
            fp = V[2, i] - V[1, i] + 1
            ((fp < minfragsize) || (fp > maxfragsize)) && continue
            if s == "+"
                mpl = V[1, i] - l.start + 1
                mpr = V[2, i] - l.start + 1
                strand = 1
            else
                mpl = l.stop - V[2, i] + 1
                mpr = l.stop - V[1, i] + 1
                strand = -1
            end

            ### flip missing?

            w = fn(V[data_entry, i])
            inc_fun(H, k, mpl, mpr,  w)
        end
    end
    H
end

function frag_heatmap_strand!(H, p, chroms, locations, strands, FM::FragMatrixSingle{T}, inc_fun=inc_heat_mid!, fraglength=120, data_entry=4, fn=identity; show_progress=true) where T

    fw = length(locations[1])

    for (k, (c, l, str)) in enumerate(zip(chroms, locations, strands))
        show_progress && next!(p)
        V = get_frags(c, (l.start - 2*fraglength):(l.stop + 2*fraglength), FM)
        isempty(V) && continue
        for i = 1:size(V, 2)
            strand, _ = GenomeFragments.get_strand_chrom_enc(V[3, i])

            if strand == STRAND_POS
                s = V[1, i] - l.start + 1
                e = s + fraglength - 1
            else
                e = V[2, i]  - l.start + 1
                s = e - fraglength + 1
            end
            if str == "-"
                s, e = fw - e + 1, fw - s + 1
            end
            w = fn(V[data_entry, i])
            inc_fun(H, k, s, e,  w)
        end
    end
    H
end

function average_heatmap(H::Vector{T}, δ) where {T}
    n = length(H)

    w = cld(n, δ)
    AH = zeros(Float64, w)
    te = zeros(Int, w)
    for i = 1:n
        wi = cld(i, δ)
        te[wi] += 1
        AH[wi] += H[i]
    end
    if te[end] < 0.75δ
        AH[end] += AH[end-1]
        te[end] += te[end-1]
    end
    AH./te
end

average_heatmap(H, δh, δv) = average_heatmap_vert(average_heatmap(H, δh), δv)
max_heatmap(H, δh, δv) = max_heatmap_vert(max_heatmap(H, δh), δv)

function average_heatmap_vert(H, δ=2)

    n, m = size(H)

    w = cld(n, δ)
    AH = zeros(Float64, w, m)
    te = zeros(Int, w)
    for i = 1:m
        for j = 1:n
            wj = cld(j, δ)
            (i == 1) && (te[wj] += 1)
            AH[wj, i] += H[j, i]
        end
    end

    if te[end] < 0.75δ
        for i = 1:m
            AH[end, i] += AH[end-1, i]
        end
        te[end] += te[end-1]
    end

    AH./te
end

function average_heatmap(H, δ=2)

    n, m = size(H)
    w = cld(m, δ)

    AH = zeros(Float64, n, w)
    te = zeros(Int, w)
    for i = 1:m
        wi = cld(i, δ)
        te[wi] += 1
        for j = 1:n
            AH[j, wi] += H[j, i]
        end
    end
    if te[end] < 0.75δ
        for j = 1:n
            AH[j, end] += AH[j, end-1]
        end
        te[end] += te[end-1]
    end
    AH./te'
end

function max_heatmap_vert(H, δ=2)

    n, m = size(H)

    w = cld(n, δ)
    AH = zeros(Float64, w, m)
    te = zeros(Int, w)
    for i = 1:m
        for j = 1:n
            wj = cld(j, δ)
            (i == 1) && (te[wj] += 1)
            AH[wj, i] = max(AH[wj, i], H[j, i])
        end
    end

    if te[end] < 0.75δ
        for i = 1:m

            AH[end, i] = max(AH[end, i], AH[end-1, i])
        end
        te[end] += te[end-1]
    end
    AH
end

function max_heatmap(H, δ=2)

    n, m = size(H)
    w = cld(m, δ)

    AH = zeros(Float64, n, w)
    te = zeros(Int, w)
    for i = 1:m
        wi = cld(i, δ)
        te[wi] += 1
        for j = 1:n
            AH[j, wi] = max(AH[j, wi], H[j, i])
        end
    end
    if te[end] < 0.75δ
        for j = 1:n
            AH[j, end] = max(AH[j, end], AH[j, end-1])
        end
        te[end] += te[end-1]
    end
    # AH./te'
    AH
end

function frag_vplot_strand_norm(chroms, locations, strands, FM::Vector{FragMatrixPair{T}}, inc_fun=inc_heat_mid!, minfragsize=0, maxfragsize=500, data_entry=4, fn=identity; norm_scale=1e+6, norm_type=:mpm, show_progress=true) where T

    fw = length(locations[1])
    fr = length(minfragsize:maxfragsize)
    H = zeros(Float64, fw, fr)
    p = ProgressMeter.Progress(length(chroms)*length(FM))
    for F in FM
        TH = zeros(Float64, fw, fr)
        frag_vplot_strand!(TH, p, chroms, locations, strands, F, inc_fun, minfragsize, maxfragsize, data_entry, fn, show_progress=show_progress)
        if norm_type == :mpm
            H .+= norm_scale.*TH./F.totalecfrags[data_entry - 3]
        elseif norm_type == :mpm_f
            t = sum(fn, view(F.FM, data_entry, :))
            H .+= norm_scale.*TH./t
        elseif norm_type == :frag
            σ = libsize(F, data_entry, fn, minfragsize, maxfragsize)
            H .+= norm_scale.*TH./σ
	elseif norm_type == :frag_region
	    H .+= norm_scale.*TH/sum(TH)
        elseif norm_type == :total
            H += TH
        else
            H += TH
        end
    end
    if norm_type == :total
        rmul!(H, 1e+6/sum(F.totalecfrags[data_entry - 3] for F in FM))
    else
        rmul!(H, 1/length(FM))
    end
    H
end

function frag_vplot_strand(chroms, locations, strands, FM::Vector{FragMatrixPair{T}}, inc_fun=inc_heat_mid!, minfragsize=0, maxfragsize=500, data_entry=4, fn=identity; show_progress=true) where T

    fw = length(locations[1])
    fr = length(minfragsize:maxfragsize)
    H = zeros(Float64, fw, fr)
    p = Progress(length(chroms)*length(FM))
    for F in FM
        frag_vplot_strand!(H, p, chroms, locations, strands, F, inc_fun, minfragsize, maxfragsize, data_entry, fn, show_progress=show_progress)
    end
    H
end

function frag_vplot_strand!(H, p, chroms, locations, strands, FM::FragMatrixPair{T}, inc_fun=inc_heat_mid!, minfragsize=0, maxfragsize=500, data_entry=4, fn=identity; show_progress=true) where T

    fw = length(locations[1])

    for (k, (c, l, s)) in enumerate(zip(chroms, locations, strands))
        show_progress && next!(p)
        V = get_frags(c, l, FM)
        isempty(V) && continue
        for i = 1:size(V, 2)
            fp = V[2, i] - V[1, i] + 1
            ((fp < minfragsize) || (fp > maxfragsize)) && continue
            if s == "+"
                mpl = V[1, i] - l.start + 1
                mpr = V[2, i] - l.start + 1
                strand = 1
            else
                mpl = l.stop - V[2, i] + 1
                mpr = l.stop - V[1, i] + 1
                strand = -1
            end

            ### flip missing?

            w = fn(V[data_entry, i])
            inc_fun(H, fp, mpl, mpr,  w)
        end
    end
    H
end

function average_heatmap(H::Vector{T}, δ) where {T}
    n = length(H)

    w = cld(n, δ)
    AH = zeros(Float64, w)
    te = zeros(Int, w)
    for i = 1:n
        wi = cld(i, δ)
        te[wi] += 1
        AH[wi] += H[i]
    end
    if te[end] < 0.75δ
        AH[end] += AH[end-1]
        te[end] += te[end-1]
    end
    AH./te
end

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

function frag_meta_strand_norm(chroms, locations, strands, FM::Vector{FragMatrixPair{T}}, inc_fun=inc_heat_mid!, minfragsize=0, maxfragsize=500, data_entry=4, fn=identity; norm_scale=1e+6, norm_type=:mpm, show_progress=true) where T

    fw = length(locations[1])
    H = zeros(Float64, fw)
    p = ProgressMeter.Progress(length(chroms)*length(FM))
    for F in FM
        TH = zeros(Float64, fw)
        frag_meta_strand!(TH, p, chroms, locations, strands, F, inc_fun, minfragsize, maxfragsize, data_entry, fn, show_progress=show_progress)
        if norm_type == :mpm
            H .+= norm_scale.*TH./F.totalecfrags[data_entry - 3]
        elseif norm_type == :mpm_f
            t = sum(fn, view(F.FM, data_entry, :))
            H .+= norm_scale.*TH./t
        elseif norm_type == :frag
            σ = libsize(F, data_entry, fn, minfragsize, maxfragsize)
            H .+= norm_scale.*TH./σ
	elseif norm_type == :frag_region
	    H .+= norm_scale.*TH/sum(TH)
        elseif norm_type == :total
            H += TH
        else
            H += TH
        end
    end
    if norm_type == :total
        rmul!(H, 1e+6/sum(F.totalecfrags[data_entry - 3] for F in FM))
    else
        rmul!(H, 1/length(FM))
    end

    rmul!(H, 1/length(chroms))
    H
end

function frag_meta_strand(chroms, locations, strands, FM::Vector{FragMatrixPair{T}}, inc_fun=inc_heat_mid!, minfragsize=0, maxfragsize=500, data_entry=4, fn=identity; show_progress=true) where T

    fw = length(locations[1])
    H = zeros(Float64, fw)
    p = Progress(length(chroms)*length(FM))
    for F in FM
        frag_meta_strand!(H, p, chroms, locations, strands, F, inc_fun, minfragsize, maxfragsize, data_entry, fn, show_progress=show_progress)
    end
    rmul!(H, 1/length(chroms))
    H
end

function frag_meta_strand!(H, p, chroms, locations, strands, FM::FragMatrixPair{T}, inc_fun=inc_heat_mid!, minfragsize=0, maxfragsize=500, data_entry=4, fn=identity; show_progress=true) where T

    fw = length(locations[1])

    for (k, (c, l, s)) in enumerate(zip(chroms, locations, strands))
        show_progress && next!(p)
        V = get_frags(c, l, FM)
        isempty(V) && continue
        for i = 1:size(V, 2)
            fp = V[2, i] - V[1, i] + 1
            ((fp < minfragsize) || (fp > maxfragsize)) && continue
            if s == "+"
                mpl = V[1, i] - l.start + 1
                mpr = V[2, i] - l.start + 1
                strand = 1
            else
                mpl = l.stop - V[2, i] + 1
                mpr = l.stop - V[1, i] + 1
                strand = -1
            end

            w = fn(V[data_entry, i])
            inc_fun(H, mpl, mpr,  w)
        end
    end
    H
end

function fragregion(chrom, location, strand, FM::FragMatrixSingle{T}, inc_fun=inc_meta_frag!, fraglength=120, dataentry=4, fn=identity; pos=true, neg=true) where {T}
    fw = length(location)
    P = zeros(fw)
    V = get_frags(chrom, location, FM)

    for i = 1:size(V, 2)
        rstrand, _ = GenomeFragments.get_strand_chrom_enc(V[3, i])
        !pos && (rstrand == STRAND_POS) && continue
        !neg && (rstrand == STRAND_NEG) && continue
        if rstrand == STRAND_POS
            s = V[1, i] - location.start + 1
            e = s + fraglength - 1
        else
            e = V[2, i]  - location.start + 1
            s = e - fraglength + 1
        end
        if strand == "-"
            s, e = fw - e + 1, fw - s + 1
        end
        w = fn(V[dataentry, i])
        inc_fun(P, s, e,  w)
    end
    P
end

@inline function fragregion(chrom, location, FM::FragMatrixPair{T}, inc_fun=inc_meta_frag!, minfrag=0, maxfrag=1000, dataentry=4, fn=identity) where {T}
    fw = length(location)
    P = zeros(fw)
    V = get_frags(chrom, location, FM)

    for col in eachcol(V)
        fp = col[2] - col[1] + 1
        s = col[1] - location.start + 1
        e = col[2] - location.start + 1
        w = fn(col[dataentry])
        (minfrag ≤ fp ≤ maxfrag) && inc_fun(P, s, e, w)
    end

    P
end

@inline function fragregion(chrom, location, FMS::Vector{FragMatrixPair{T}}, inc_fun=inc_meta_frag!, minfrag=0, maxfrag=1000, dataentry=4, fn=identity) where {T}
    fw = length(location)
    P = zeros(fw)

    for FM in FMS
        V = get_frags(chrom, location, FM)

        for col in eachcol(V)
            fp = col[2] - col[1] + 1
            s = col[1] - location.start + 1
            e = col[2] - location.start + 1
            w = fn(col[dataentry])
            (minfrag ≤ fp ≤ maxfrag) && inc_fun(P, s, e, w)
        end

        P
    end
    P
end

function count_cut_sites(chroms, locations, strands, FM, minfragsize=0, maxfragsize=1000, data_entry=4, fn=identity) where T

    cutsites = zeros(Int, length(chroms))
    for (k, (c, l, s)) in enumerate(zip(chroms, locations, strands))

        V = get_frags(c, l, FM)
        isempty(V) && continue
        for i = 1:size(V, 2)
            fp = V[2, i] - V[1, i] + 1
            ((fp < minfragsize) || (fp > maxfragsize)) && continue

            w = fn(V[data_entry, i])
            cutsites[k] += ifelse(V[1, i] ∈ l, w, 0)
            cutsites[k] += ifelse(V[2, i] ∈ l, w, 0)
        end
    end
    cutsites

end

function readsintersecting(chroms, locations, FM, dataentry=4, fn=identity)
    total = zeros(Int, length(chroms))
    for (k, (c, l)) in enumerate(zip(chroms, locations))

        v = get_frags(c, l, FM)
        for i = 1:size(v, 2)
            total[k] += fn(v[dataentry, i])
        end
    end
    total
end

function peakheights(chroms, locations, FM::Vector{FragMatrixPair{T}} ; minfragsize=0, maxfragsize=1000, dataentry=4, fn=one, scale=1e+6) where {T}

    PH = zeros(length(chroms))
    σ = 0
    σ = sum(libsize(F, dataentry, fn) for F in FM)

    for (k, (c, l)) in enumerate(zip(chroms, locations))
        PH[k] = maximum(mapreduce(F -> fragregion(c, l, F, inc_meta_frag!, minfragsize, maxfragsize, dataentry, fn), +, FM))
    end

    NH = scale.*PH./σ
    NH
end

function peakheightsloc(chroms, locations, FM ; minfragsize=0, maxfragsize=1000, dataentry=4, fn=one, scale=1e+6)

    PH = zeros(length(chroms))
    PS = zeros(Int, length(chroms))

    for (k, (c, l)) in enumerate(zip(chroms, locations))
        PH[k], PS[k] = findmax(fragregion(c, l, FM, inc_meta_frag!, minfragsize, maxfragsize, dataentry, fn))
    end

    PH, PS
end

function frag_heat_cut_pair_strand_norm(chroms, locations, strands, FM::Vector{FragMatrixPair{T}}, minfragsize=0, maxfragsize=500, data_entry=4, fn=identity; norm_scale=1e+6, norm_type=:mpm, show_progress=true) where T

    fw = length(locations[1])
    # fr = length(minfragsize:maxfragsize)
    H = zeros(Float64, fw, fw)
    p = ProgressMeter.Progress(length(chroms)*length(FM))
    for F in FM
        TH = zeros(Float64, fw, fw)
        frag_heat_cut_pair_strand!(TH, p, chroms, locations, strands, F, minfragsize, maxfragsize, data_entry, fn, show_progress=show_progress)
        if norm_type == :mpm
            H .+= norm_scale.*TH./F.totalecfrags[data_entry - 3]
        elseif norm_type == :mpm_f
            t = sum(fn, view(F.FM, data_entry, :))
            H .+= norm_scale.*TH./t
        elseif norm_type == :frag
            σ = libsize(F, data_entry, fn, minfragsize, maxfragsize)
            H .+= norm_scale.*TH./σ
	elseif norm_type == :frag_region
	    H .+= norm_scale.*TH/sum(TH)
        elseif norm_type == :total
            H += TH
        else
            H += TH
        end
    end
    if norm_type == :total
        rmul!(H, 1e+6/sum(F.totalecfrags[data_entry - 3] for F in FM))
    else
        rmul!(H, 1/length(FM))
    end
    H
end

function frag_heat_cut_pair_strand!(H, p, chroms, locations, strands, FM::FragMatrixPair{T}, minfragsize=0, maxfragsize=500, data_entry=4, fn=identity; show_progress=true) where T

    fw = length(locations[1])

    for (k, (c, l, s)) in enumerate(zip(chroms, locations, strands))
        show_progress && next!(p)
        V = get_frags(c, l, FM)
        isempty(V) && continue
        for i = 1:size(V, 2)
            fp = V[2, i] - V[1, i] + 1
            ((fp < minfragsize) || (fp > maxfragsize)) && continue
            if s == "+"
                mpl = V[1, i] - l.start + 1
                mpr = V[2, i] - l.start + 1
                strand = 1
            else
                mpl = l.stop - V[2, i] + 1
                mpr = l.stop - V[1, i] + 1
                strand = -1
            end

            w = fn(V[data_entry, i])
            if (1 <= mpl <= fw) && (1 <= mpr <= fw)
                H[mpl, mpr] += w
            end
        end
    end
    H
end

function frag_heat_cut_pair_strand_norm(chroms, locations, strands, FM::Vector{FragMatrixPair{T}}, minfragsize=0, maxfragsize=500, data_entry=4, fn=identity; norm_scale=1e+6, norm_type=:mpm, show_progress=true) where T

    fw = length(locations[1])
    # fr = length(minfragsize:maxfragsize)
    H = zeros(Float64, fw, fw)
    p = ProgressMeter.Progress(length(chroms)*length(FM))
    for F in FM
        TH = zeros(Float64, fw, fw)
        frag_heat_cut_pair_strand!(TH, p, chroms, locations, strands, F, minfragsize, maxfragsize, data_entry, fn, show_progress=show_progress)
        if norm_type == :mpm
            H .+= norm_scale.*TH./F.totalecfrags[data_entry - 3]
        elseif norm_type == :mpm_f
            t = sum(fn, view(F.FM, data_entry, :))
            H .+= norm_scale.*TH./t
        elseif norm_type == :frag
            σ = libsize(F, data_entry, fn, minfragsize, maxfragsize)
            H .+= norm_scale.*TH./σ
	elseif norm_type == :frag_region
	    H .+= norm_scale.*TH/sum(TH)
        elseif norm_type == :total
            H += TH
        else
            H += TH
        end
    end
    if norm_type == :total
        rmul!(H, 1e+6/sum(F.totalecfrags[data_entry - 3] for F in FM))
    else
        rmul!(H, 1/length(FM))
    end
    H
end

function frag_tensor_cut_pair_strand_norm(chroms, locations, strands, FM::Vector{FragMatrixPair{T}}, minfragsize=0, maxfragsize=500, data_entry=4, fn=identity; norm_scale=1e+6, norm_type=:mpm, show_progress=true) where T

    fw = length(locations[1])
    fr = length(chroms)
    H = [spzeros(Float64, fw, fw) for i = 1:fr]

    p = ProgressMeter.Progress(length(chroms)*length(FM))
    for F in FM
        TH = [spzeros(Float64, fw, fw) for i = 1:fr]
        frag_tensor_cut_pair_strand!(TH, p, chroms, locations, strands, F, minfragsize, maxfragsize, data_entry, fn, show_progress=show_progress)
        if norm_type == :mpm
            H .+= norm_scale.*TH./F.totalecfrags[data_entry - 3]
        elseif norm_type == :mpm_f
            t = sum(fn, view(F.FM, data_entry, :))
            H .+= norm_scale.*TH./t
        elseif norm_type == :frag
            σ = libsize(F, data_entry, fn, minfragsize, maxfragsize)
            H .+= norm_scale.*TH./σ
	    elseif norm_type == :frag_region
	        H .+= norm_scale.*TH/sum(TH)
        elseif norm_type == :total
            H .+= TH
        else
            H .+= TH
        end
    end
    if norm_type == :total
        #rmul!(H,)
        f = 1e+6/sum(F.totalecfrags[data_entry - 3] for F in FM)
        H .*= f
    else
        H ./= length(FM)
    end
    H
end

function frag_tensor_cut_pair_strand!(H, p, chroms, locations, strands, FM::FragMatrixPair{T}, minfragsize=0, maxfragsize=500, data_entry=4, fn=identity; show_progress=true) where T

    fw = length(locations[1])

    for (k, (c, l, s)) in enumerate(zip(chroms, locations, strands))
        show_progress && next!(p)
        V = get_frags(c, l, FM)
        isempty(V) && continue
        for i = 1:size(V, 2)
            fp = V[2, i] - V[1, i] + 1
            ((fp < minfragsize) || (fp > maxfragsize)) && continue
            if s == "+"
                mpl = V[1, i] - l.start + 1
                mpr = V[2, i] - l.start + 1
                strand = 1
            else
                mpl = l.stop - V[2, i] + 1
                mpr = l.stop - V[1, i] + 1
                strand = -1
            end

            w = fn(V[data_entry, i])
            if (1 <= mpl <= fw) && (1 <= mpr <= fw)
                H[k][mpl, mpr] += w
            end
        end
    end
    H
end