text:
#using Ipopt: optimize!
using Base: Int64
using XLSX: length
#=using Pkg
Pkg.add("Tables")
Pkg.add("DataFrames")
Pkg.add("XLSX")
Pkg.add("ExcelReaders")
Pkg.add("JuMP")
Pkg.add("Ipopt")
=#
using XLSX, ExcelReaders, DataFrames, Tables, JuMP, Ipopt;
IOSource = XLSX.readdata("IO.xlsx", "io-table-5!A1:DV130");

#filepath cross system compatability code
if Sys.KERNEL === :NT || kern === :Windows
    pathmark = "\\"
else
    pathmark = "/"
end
#indexing vectors for initial data import groups
intermediaryTotalsCol = findall(x -> occursin("T4", x), string.(IOSource[3,:]));
intermediaryTotalsRow = findall(x -> occursin("T1", x), string.(IOSource[:,1]));
finalTotalsCol = findall(x -> occursin("T6", x), string.(IOSource[3,:]));
finalTotalsRow = findall(x -> occursin("Australian Production", x), string.(IOSource[:,2]));
finalDemandCol = findall(x -> occursin('Q', x), string.(IOSource[3,:]));
factorRow = findall(x -> occursin('P', x), string.(IOSource[:,1]));
IOSourceCol = vcat(intermediaryTotalsCol, finalDemandCol, finalTotalsCol);
IOSourceRow = vcat(intermediaryTotalsRow, factorRow, finalTotalsRow);
#initialising IO
IO = zeros(length(IOSourceRow), length(IOSourceCol));
#import numerical data into IO
IO[1:length(IOSourceRow), 1:length(IOSourceCol)] = IOSource[IOSourceRow, IOSourceCol];
#creating vectors of titles for IO
IOCodeRow = IOSource[IOSourceRow, 1];
IOCodeCol = IOSource[3, IOSourceCol];
IONameRow = IOSource[IOSourceRow, 2];
IONameCol = IOSource[2, IOSourceCol];

#code to sum public and private entities into one collumn
investment = findall(x -> occursin("Capital Formation", x), IONameCol);
IO[:, investment[1]]=sum(eachcol(IO[:, investment[1:2]]));
IO = IO[:,Not(investment[2])];
#alter title vectors accordingly (include Q in total investment collumn in IOcode)
IOCodeCol[investment[1]] = "Q3+Q4";
IOCodeCol = IOCodeCol[Not(investment[2])];
IONameCol[investment[1]] = "Private and Public Gross Fixed Capital Formation";
IONameCol = IONameCol[Not(investment[2])];
#creating a dictionary for the index of each collumn and row in IO by IOCode
IOColDict = Dict(IOCodeCol .=> [1:1:8;]);
IORowDict = Dict(IOCodeRow .=> [1:1:8;]);
IOCapForm = findall(x -> occursin("Capital Formation", x), IONameCol);
IOChangeInv = findall(x -> occursin("Changes in Inventories", x), IONameCol);

#importing relevant ASNA data for table 5
ASNAHouseCap = ExcelReaders.readxl("ASNAData"*pathmark*"5204039_Household_Capital_Account.xls", "Data1!A1:T71");
ASNANonFinCap = ExcelReaders.readxl("ASNAData"*pathmark*"5204018_NonFin_Corp_Capital_Account.xls", "Data1!A1:T71");
ASNAFinCap = ExcelReaders.readxl("ASNAData"*pathmark*"5204026_Fin_Corp_Capital_Account.xls", "Data1!A1:S71");
ASNAGovCap = ExcelReaders.readxl("ASNAData"*pathmark*"5204032_GenGov_Capital_Account.xls", "Data1!A1:AV71");
ASNAYearRow = findall(x -> occursin("2019", x), string.(ASNAHouseCap[:,1]));

#creating table 5a - allocation of investment expenditure (broken into subsections for dict referencing purposes)
#subsection a is fixed capital expenditure
table5aNameCol = ["Households", "Non-Financial Corporations", "Financial Corporations", "General Government", "Total"];
table5aNameRow = ["Domestic Commodities", "Imported Commodities, complementary", "Imported Commodities, competing",
"Taxes less subsidies on products", "Other taxes less subsidies on investment", "Total indirect taxes",
"Total fixed capital expenditure"];
table5aRowDict = Dict(table5aNameRow .=> [1:1:length(table5aNameRow);]);
table5aColDict = Dict(table5aNameCol .=> [1:1:length(table5aNameCol);]);
table5a = zeros(length(table5aNameRow), length(table5aNameCol));

#filling in totals collumn from corresponding IO data
table5a[table5aRowDict["Domestic Commodities"], table5aColDict["Total"]] = sum(IO[IORowDict["T1"],IOCapForm]);
table5a[table5aRowDict["Imported Commodities, complementary"], table5aColDict["Total"]] = sum(IO[IORowDict["P5"],IOCapForm]);
table5a[table5aRowDict["Imported Commodities, competing"], table5aColDict["Total"]] = sum(IO[IORowDict["P6"],IOCapForm]);
table5a[table5aRowDict["Taxes less subsidies on products"], table5aColDict["Total"]] = sum(IO[IORowDict["P3"],IOCapForm]);
table5a[table5aRowDict["Other taxes less subsidies on investment"], table5aColDict["Total"]] = sum(IO[IORowDict["P4"],IOCapForm]);
table5aTaxes = findall(x -> occursin("taxes", lowercase(x)), table5aNameRow);
table5aTaxes = table5aTaxes[Not(3)];
table5a[table5aRowDict["Total indirect taxes"], table5aColDict["Total"]] = sum(table5a[table5aTaxes,table5aColDict["Total"]]);
table5a[table5aRowDict["Total fixed capital expenditure"], table5aColDict["Total"]] = sum(table5a[Not(table5aRowDict["Total indirect taxes"]),table5aColDict["Total"]]);

#creating index variables for the measurements that we want
ASNAHouseCapTotCapForm = findall(x -> occursin("Gross fixed capital formation ;", x), string.(ASNAHouseCap[1,:]));
ASNANonFinCapTotCapForm = findall(x -> occursin("Gross fixed capital formation ;", x), string.(ASNANonFinCap[1,:]));
ASNAFinCapTotCapForm = findall(x -> occursin("Gross fixed capital formation ;", x), string.(ASNAFinCap[1,:]));
ASNAGenGovCapTotCapForm = findall(x -> occursin("General government ;  Gross fixed capital formation ;", x), string.(ASNAGovCap[1,:]));

#filling in totals row from ASNA Data
table5a[table5aRowDict["Total fixed capital expenditure"], table5aColDict["Households"]] = first(ASNAHouseCap[ASNAYearRow, ASNAHouseCapTotCapForm]);
table5a[table5aRowDict["Total fixed capital expenditure"], table5aColDict["Non-Financial Corporations"]] = first(ASNANonFinCap[ASNAYearRow, ASNANonFinCapTotCapForm]);
table5a[table5aRowDict["Total fixed capital expenditure"], table5aColDict["Financial Corporations"]] = first(ASNAFinCap[ASNAYearRow, ASNAFinCapTotCapForm]);
table5a[table5aRowDict["Total fixed capital expenditure"], table5aColDict["General Government"]] = first(ASNAGovCap[ASNAYearRow, ASNAGenGovCapTotCapForm]);

#filling in non-total values
for ring in [1:1:length(table5aColDict)-1;];
    table5a[table5aRowDict["Domestic Commodities"],ring] = table5a[table5aRowDict["Total fixed capital expenditure"],ring]*IO[IORowDict["T1"],IOCapForm[1]] / IO[IORowDict[missing],IOCapForm[1]];
    table5a[table5aRowDict["Imported Commodities, complementary"],ring] = table5a[table5aRowDict["Total fixed capital expenditure"],ring]*IO[IORowDict["P5"],IOCapForm[1]] / IO[IORowDict[missing],IOCapForm[1]];
    table5a[table5aRowDict["Imported Commodities, competing"],ring] = table5a[table5aRowDict["Total fixed capital expenditure"],ring]*IO[IORowDict["P6"],IOCapForm[1]] / IO[IORowDict[missing],IOCapForm[1]];
    table5a[table5aRowDict["Taxes less subsidies on products"],ring] = table5a[table5aRowDict["Total fixed capital expenditure"],ring]*IO[IORowDict["P3"],IOCapForm[1]] / IO[IORowDict[missing],IOCapForm[1]];
    table5a[table5aRowDict["Other taxes less subsidies on investment"],ring] = table5a[table5aRowDict["Total fixed capital expenditure"],ring]*IO[IORowDict["P4"],IOCapForm[1]] / IO[IORowDict[missing],IOCapForm[1]];
    table5a[table5aRowDict["Total indirect taxes"], ring] = sum(table5a[table5aTaxes, ring]);
end

#creating table 5b - allocation of investment expenditure (broken into subsections for dict referencing purposes)
#subsection b is fixed capital expenditure
table5bNameCol = ["Households", "Non-Financial Corporations", "Financial Corporations", "General Government", "Total"];
table5bNameRow = ["Domestic Commodities", "Imported Commodities, complementary", "Imported Commodities, competing",
"Taxes less subsidies on products", "Total change in inventories"];
table5bRowDict = Dict(table5bNameRow .=> [1:1:length(table5bNameRow);]);
table5bColDict = Dict(table5bNameCol .=> [1:1:length(table5bNameCol);]);
table5b = zeros(length(table5bNameRow), length(table5bNameCol));

#filling in totals collumn from corresponding IO data
table5b[table5bRowDict["Domestic Commodities"], table5bColDict["Total"]] = sum(IO[IORowDict["T1"],IOChangeInv]);
table5b[table5bRowDict["Imported Commodities, complementary"], table5bColDict["Total"]] = sum(IO[IORowDict["P5"],IOChangeInv]);
table5b[table5bRowDict["Imported Commodities, competing"], table5bColDict["Total"]] = sum(IO[IORowDict["P6"],IOChangeInv]);
table5b[table5bRowDict["Taxes less subsidies on products"], table5bColDict["Total"]] = sum(IO[IORowDict["P3"],IOChangeInv]);
table5b[table5bRowDict["Total change in inventories"], table5bColDict["Total"]] = sum(table5b[:,table5bColDict["Total"]]);

#creating index variables for the measurements that we want
ASNAHouseCapChangeInv = findall(x -> occursin("Changes in inventories ;", x), string.(ASNAHouseCap[1,:]));
ASNANonFinCapChangeInv = findall(x -> occursin("Changes in inventories ;", x), string.(ASNANonFinCap[1,:]));
ASNAFinCapChangeInv = findall(x -> occursin("Changes in inventories ;", x), string.(ASNAFinCap[1,:]));
ASNAGenGovCapChangeInv = findall(x -> occursin("General government ;  Changes in inventories ;", x), string.(ASNAGovCap[1,:]));

#filling in totals row from ASNA Data
table5b[table5bRowDict["Total change in inventories"], table5bColDict["Households"]] = first(ASNAHouseCap[ASNAYearRow, ASNAHouseCapChangeInv]);
table5b[table5bRowDict["Total change in inventories"], table5bColDict["Non-Financial Corporations"]] = first(ASNANonFinCap[ASNAYearRow, ASNANonFinCapChangeInv]);
table5b[table5bRowDict["Total change in inventories"], table5bColDict["Financial Corporations"]] = first(ASNAFinCap[ASNAYearRow, ASNAFinCapChangeInv]);
table5b[table5bRowDict["Total change in inventories"], table5bColDict["General Government"]] = first(ASNAGovCap[ASNAYearRow, ASNAGenGovCapChangeInv]);

#fill in non-total value with lagrangian optimisation
table5bscalingfact = abs(minimum(table5b)) * 1.5;
table5b = table5b .+ table5bscalingfact;
mod5b = Model(Ipopt.Optimizer);
@variable(mod5b, x[1:4, 1:4]);
@NLobjective(mod5b, Min, sum((x[i,j] - table5bscalingfact) ^ 2 for i in 1:4, j in 1:4));
for i in 1:4
    @constraint(mod5b, sum(x[:,i]) == table5b[table5bRowDict["Total change in inventories"],i]);
end
for i in 1:4
    @constraint(mod5b, sum(x[i,:]) == table5b[i,table5bColDict["Total"]]);
end
optimize!(mod5b);
table5b[1:4,1:4]=value.(x);

table5b[1:4,1:4] = table5b[1:4,1:4] .- table5bscalingfact/4;
table5b[:,5] = table5b[:,5] .- table5bscalingfact;
table5b[5,:] = table5b[5,:] .- table5bscalingfact;

#creating table 5c - allocation of investment expenditure (broken into subsections for dict referencing purposes)
#subsection c is totals
table5cNameCol = ["Households", "Non-Financial Corporations", "Financial Corporations", "General Government", "Total"];
table5cNameRow = ["Domestic Commodities", "Imported Commodities", "Taxes less subsidies on products", "Other taxes less subsidies on investment", "Total investment expenditure"];
table5cRowDict = Dict(table5cNameRow .=> [1:1:length(table5cNameRow);]);
table5cColDict = Dict(table5cNameCol .=> [1:1:length(table5cNameCol);]);
table5c = zeros(length(table5cNameRow), length(table5cNameCol));

#do totals calcuations to get all values in 5c

#=convert dataframe to dictionary
function increment!( d::Dict{S, T}, k::S, i::T) where {T<:Real, S<:Any}
    if haskey(d, k)
        d[k] += i
    else
        d[k] = i
    end
end
increment!(d::Dict{S, T}, k::S ) where {T<:Real, S<:Any} = increment!( d, k, one(T))

function df2dict( df::DataFrame, key_col::Symbol, val_col::Symbol=:null)
    keytype = typeof(df[1,key_col])
    if val_col == :null
        valtype = Int
    else
        valtype = typeof(df[1,val_col])
    end
    D = Dict{keytype, valtype}()
    for i=1:size(df,1)
        if !ismissing(df[i,key_col])
            if val_col == :null
                increment!( D, df[i,key_col] )
            elseif valtype <: Real
                increment!( D, df[i,key_col], df[i,val_col] )
            else
                if haskey(D, df[i,key_col])
                    @warn("non-unique entry: $(df[i,key_col])")
                else
                    D[df[i,key_col]] = df[i,val_col]
                end
            end
        end
    end
    return D
end
df[!, "IOCode"]=IOcode
insertcols!(df, 2, :name => vector)
D = df2dict(df, :IOCode, :x3)
=#


text:
# Mahalanobis distances

"""
    Mahalanobis(Q; skipchecks=false) <: Metric

Create a Mahalanobis distance (i.e., a bilinear form) with covariance matrix `Q`.
Upon construction, both symmetry/self-adjointness and positive semidefiniteness are checked,
where the latter check can be skipped by passing the keyword argument `skipchecks = true`.

# Example:
```julia
julia> A = collect(reshape(1:9, 3, 3)); Q = A'A;

julia> dist = Mahalanobis(Q)
Mahalanobis{Matrix{Int64}}([14 32 50; 32 77 122; 50 122 194])

julia> dist = Mahalanobis(A, skipchecks=true)
┌ Warning: matrix is not symmetric/Hermitian
└ @ Distances ...
Mahalanobis{Matrix{Int64}}([1 4 7; 2 5 8; 3 6 9])
"""
struct Mahalanobis{M<:AbstractMatrix} <: Metric
    qmat::M
    function Mahalanobis(Q::AbstractMatrix; skipchecks::Bool=false)
        # TODO: turn the warnings into errors in next breaking release
        ishermitian(Q) || @warn "matrix is not symmetric/Hermitian"
        if !skipchecks
            eigmin(Q) ≥ 0 || @warn "matrix is not positive semidefinite"
        end
        return new{typeof(Q)}(Q)
    end
end

"""
    SqMahalanobis(Q; skipchecks=false) <: Metric

Create a squared Mahalanobis distance (i.e., a bilinear form) with covariance matrix `Q`.
Upon construction, both symmetry/self-adjointness and positive semidefiniteness are checked,
where the latter check can be skipped by passing the keyword argument `skipchecks = true`.

# Example:
```julia
julia> A = collect(reshape(1:9, 3, 3)); Q = A'A;

julia> dist = SqMahalanobis(Q)
SqMahalanobis{Matrix{Int64}}([14 32 50; 32 77 122; 50 122 194])

julia> dist = SqMahalanobis(A, skipchecks=true)
┌ Warning: matrix is not symmetric/Hermitian
└ @ Distances ...
SqMahalanobis{Matrix{Int64}}([1 4 7; 2 5 8; 3 6 9])
"""
struct SqMahalanobis{M<:AbstractMatrix} <: SemiMetric
    qmat::M
    function SqMahalanobis(Q::AbstractMatrix; skipchecks::Bool=false)
        # TODO: turn the warnings into errors in next breaking release
        ishermitian(Q) || @warn "matrix is not symmetric/Hermitian"
        if !skipchecks
            eigmin(Q) ≥ 0 || @warn "matrix is not positive semidefinite"
        end
        return new{typeof(Q)}(Q)
    end
end

function result_type(d::Mahalanobis, ::Type{T1}, ::Type{T2}) where {T1,T2}
    z = zero(T1) - zero(T2)
    return typeof(sqrt(z * zero(eltype(d.qmat)) * z))
end

function result_type(d::SqMahalanobis, ::Type{T1}, ::Type{T2}) where {T1,T2}
    z = zero(T1) - zero(T2)
    return typeof(z * zero(eltype(d.qmat)) * z)
end

# SqMahalanobis

function (dist::SqMahalanobis)(a::AbstractVector, b::AbstractVector)
    if length(a) != length(b)
        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end

    Q = dist.qmat
    z = a - b
    return dot(z, Q * z)
end

sqmahalanobis(a::AbstractVector, b::AbstractVector, Q::AbstractMatrix) = SqMahalanobis(Q)(a, b)

function colwise!(r::AbstractArray, dist::SqMahalanobis, a::AbstractMatrix, b::AbstractMatrix)
    Q = dist.qmat
    get_colwise_dims(size(Q, 1), r, a, b)
    z = a - b
    dot_percol!(r, Q * z, z)
end

function colwise!(r::AbstractArray, dist::SqMahalanobis, a::AbstractVector, b::AbstractMatrix)
    Q = dist.qmat
    get_colwise_dims(size(Q, 1), r, a, b)
    z = a .- b
    dot_percol!(r, Q * z, z)
end

function _pairwise!(r::AbstractMatrix, dist::SqMahalanobis, a::AbstractMatrix, b::AbstractMatrix)
    Q = dist.qmat
    m, na, nb = get_pairwise_dims(size(Q, 1), r, a, b)

    Qa = Q * a
    Qb = Q * b
    sa2 = dot_percol(a, Qa)
    sb2 = dot_percol(b, Qb)
    mul!(r, a', Qb)

    for j = 1:nb
        @simd for i = 1:na
            @inbounds r[i, j] = max(sa2[i] + sb2[j] - 2 * r[i, j], 0)
        end
    end
    r
end

function _pairwise!(r::AbstractMatrix, dist::SqMahalanobis, a::AbstractMatrix)
    Q = dist.qmat
    m, n = get_pairwise_dims(size(Q, 1), r, a)

    Qa = Q * a
    sa2 = dot_percol(a, Qa)
    mul!(r, a', Qa)

    for j = 1:n
        for i = 1:(j - 1)
            @inbounds r[i, j] = r[j, i]
        end
        r[j, j] = 0
        for i = (j + 1):n
            @inbounds r[i, j] = max(sa2[i] + sa2[j] - 2 * r[i, j], 0)
        end
    end
    r
end

# Mahalanobis

function (dist::Mahalanobis)(a::AbstractVector, b::AbstractVector)
    sqrt(SqMahalanobis(dist.qmat, skipchecks = true)(a, b))
end

mahalanobis(a::AbstractVector, b::AbstractVector, Q::AbstractMatrix) = Mahalanobis(Q)(a, b)

function colwise!(r::AbstractArray, dist::Mahalanobis, a::AbstractMatrix, b::AbstractMatrix)
    sqrt!(colwise!(r, SqMahalanobis(dist.qmat, skipchecks = true), a, b))
end

function colwise!(r::AbstractArray, dist::Mahalanobis, a::AbstractVector, b::AbstractMatrix)
    sqrt!(colwise!(r, SqMahalanobis(dist.qmat, skipchecks = true), a, b))
end

function _pairwise!(r::AbstractMatrix, dist::Mahalanobis, a::AbstractMatrix, b::AbstractMatrix)
    sqrt!(_pairwise!(r, SqMahalanobis(dist.qmat, skipchecks = true), a, b))
end

function _pairwise!(r::AbstractMatrix, dist::Mahalanobis, a::AbstractMatrix)
    sqrt!(_pairwise!(r, SqMahalanobis(dist.qmat, skipchecks = true), a))
end

ext: jl


text:
cd(@__DIR__)
lines = readlines("../monitor.txt")
lines = split.(lines, ":")

filterlines(name) = [eval(Meta.parse(l[2])) for l in filter(x->x[1]==name, lines)]

rws, wws = filterlines("readweight"), filterlines("writeweight")

rws = hcat(rws...)
wws = hcat(wws...)


text:
"""
    defines_strides(::Type{T}) -> Bool

Is strides(::T) defined? It is assumed that types returning `true` also return a valid
pointer on `pointer(::T)`.
"""
defines_strides(x) = defines_strides(typeof(x))
_defines_strides(::Type{T}, ::Type{T}) where {T} = false
_defines_strides(::Type{P}, ::Type{T}) where {P,T} = defines_strides(P)
defines_strides(::Type{T}) where {T} = _defines_strides(parent_type(T), T)
defines_strides(::Type{<:StridedArray}) = true
function defines_strides(::Type{<:SubArray{T,N,P,I}}) where {T,N,P,I}
    stride_preserving_index(I) === True()
end
defines_strides(::Type{<:BitArray}) = true

#=
    stride_preserving_index(::Type{T}) -> StaticBool

Returns `True` if strides between each element can still be derived when indexing with an
instance of type `T`.
=#
stride_preserving_index(::Type{T}) where {T<:AbstractRange} = True()
stride_preserving_index(::Type{T}) where {T<:Int} = True()
stride_preserving_index(::Type{T}) where {T} = False()
function stride_preserving_index(::Type{T}) where {N,T<:Tuple{Vararg{Any,N}}}
    if all(eachop(_stride_preserving_index, ntuple(static, StaticInt(N)), T))
        return True()
    else
        return False()
    end
end
function _stride_preserving_index(::Type{T}, i::StaticInt) where {T}
    return stride_preserving_index(field_type(T, i))
end

"""
    known_offsets(::Type{T}) -> Tuple
    known_offsets(::Type{T}, dim) -> Union{Int,Nothing}

Returns a tuple of offset values known at compile time. If the offset of a given axis is
not known at compile time `nothing` is returned its position.
"""
known_offsets(x, dim) = known_offsets(typeof(x), dim)
known_offsets(::Type{T}, dim) where {T} = known_offsets(T, to_dims(T, dim))
function known_offsets(::Type{T}, dim::CanonicalInt) where {T}
    if ndims(T) < dim
        return 1
    else
        return known_offsets(T)[dim]
    end
end

known_offsets(x) = known_offsets(typeof(x))
function known_offsets(::Type{T}) where {T}
    return eachop(_known_offsets, ntuple(static, StaticInt(ndims(T))), axes_types(T))
end
_known_offsets(::Type{T}, dim::StaticInt) where {T} = known_first(field_type(T, dim))

known_offsets(::Type{<:StrideIndex{N,R,C,S,O}}) where {N,R,C,S,O} = known(O)

"""
    offsets(A) -> Tuple
    offsets(A, dim) -> Union{Int,StaticInt}

Returns offsets of indices with respect to 0. If values are known at compile time,
it should return them as `Static` numbers.
For example, if `A isa Base.Matrix`, `offsets(A) === (StaticInt(1), StaticInt(1))`.
"""
@inline offsets(x, i) = static_first(indices(x, i))
offsets(::Tuple) = (One(),)
offsets(x::StrideIndex) = getfield(x, :offsets)
offsets(x) = eachop(_offsets, ntuple(static, StaticInt(ndims(x))), x)
function _offsets(x::X, dim::StaticInt{D}) where {X,D}
    start = known_first(axes_types(X, dim))
    if start === nothing
        return first(axes(x, dim))
    else
        return static(start)
    end
end
# we can't generate an axis for `StrideIndex` so this is performed manually here
@inline offsets(x::StrideIndex, dim::Int) = getfield(offsets(x), dim)
@inline offsets(x::StrideIndex, ::StaticInt{dim}) where {dim} = getfield(offsets(x), dim)

"""
    known_offset1(::Type{T}) -> Union{Int,Nothing}

Returns the linear offset of array `x` if known at compile time.
"""
@inline known_offset1(x) = known_offset1(typeof(x))
@inline function known_offset1(::Type{T}) where {T}
    if ndims(T) === 0
        return 1
    else
        return known_offsets(T, 1)
    end
end

"""
    offset1(x) -> Union{Int,StaticInt}

Returns the offset of the linear indices for `x`.
"""
@inline function offset1(x::X) where {X}
    o1 = known_offset1(X)
    if o1 === nothing
        if ndims(X) === 0
            return 1
        else
            return offsets(x, 1)
        end
    else
        return static(o1)
    end
end

"""
    contiguous_axis(::Type{T}) -> StaticInt{N}

Returns the axis of an array of type `T` containing contiguous data.
If no axis is contiguous, it returns a `StaticInt{-1}`.
If unknown, it returns `nothing`.
"""
contiguous_axis(x) = contiguous_axis(typeof(x))
contiguous_axis(::Type{<:StrideIndex{N,R,C}}) where {N,R,C} = static(C)
contiguous_axis(::Type{<:StrideIndex{N,R,nothing}}) where {N,R} = nothing
function contiguous_axis(::Type{T}) where {T}
    is_forwarding_wrapper(T) ? contiguous_axis(parent_type(T)) : nothing
end
contiguous_axis(::Type{<:DenseArray}) = One()
contiguous_axis(::Type{<:BitArray}) = One()
contiguous_axis(::Type{<:AbstractRange}) = One()
contiguous_axis(::Type{<:Tuple}) = One()
function contiguous_axis(::Type{T}) where {T<:VecAdjTrans}
    c = contiguous_axis(parent_type(T))
    if c === nothing
        return nothing
    elseif c === One()
        return StaticInt{2}()
    else
        return -One()
    end
end
function contiguous_axis(::Type{T}) where {T<:MatAdjTrans}
    c = contiguous_axis(parent_type(T))
    if c === nothing
        return nothing
    elseif isone(-c)
        return c
    else
        return StaticInt(3) - c
    end
end
function contiguous_axis(::Type{T}) where {T<:PermutedDimsArray}
    c = contiguous_axis(parent_type(T))
    if c === nothing
        return nothing
    elseif isone(-c)
        return c
    else
        return from_parent_dims(T, c)
    end
end
function contiguous_axis(::Type{<:Base.ReshapedArray{T, N, A, Tuple{}}}) where {T, N, A}
    c = contiguous_axis(A)
    if c !== nothing && isone(-c)
        return StaticInt(-1)
    elseif dynamic(is_column_major(A) & is_dense(A))
        return StaticInt(1)
    else
        return nothing
    end
end
# we're actually looking at the original vector indices before transpose/adjoint
function contiguous_axis(::Type{<:Base.ReshapedArray{T, 1, A, Tuple{}}}) where {T, A <: VecAdjTrans}
    contiguous_axis(parent_type(A))
end
function contiguous_axis(::Type{T}) where {T<:SubArray}
    _contiguous_axis(T, contiguous_axis(parent_type(T)))
end

_contiguous_axis(::Type{A}, ::Nothing) where {T,N,P,I,A<:SubArray{T,N,P,I}} = nothing
_contiguous_axis(::Type{A}, c::StaticInt{-1}) where {T,N,P,I,A<:SubArray{T,N,P,I}} = c
function _contiguous_axis(::Type{A}, c::StaticInt{C}) where {T,N,P,I,A<:SubArray{T,N,P,I},C}
    if field_type(I, c) <: AbstractUnitRange
        return from_parent_dims(A)[C]
    elseif field_type(I, c) <: AbstractArray
        return -One()
    elseif field_type(I, c) <: CanonicalInt
        return -One()
    else
        return nothing
    end
end

# contiguous_if_one(::StaticInt{1}) = StaticInt{1}()
# contiguous_if_one(::Any) = StaticInt{-1}()
function contiguous_axis(::Type{R}) where {T,N,S,A<:Array{S},R<:ReinterpretArray{T,N,S,A}}
    if isbitstype(S)
        return One()
    else
        return nothing
    end
end

"""
    contiguous_axis_indicator(::Type{T}) -> Tuple{Vararg{StaticBool}}

Returns a tuple boolean `Val`s indicating whether that axis is contiguous.
"""
function contiguous_axis_indicator(::Type{A}) where {D,A<:AbstractArray{<:Any,D}}
    return contiguous_axis_indicator(contiguous_axis(A), Val(D))
end
contiguous_axis_indicator(::A) where {A<:AbstractArray} = contiguous_axis_indicator(A)
contiguous_axis_indicator(::Nothing, ::Val) = nothing
function contiguous_axis_indicator(c::StaticInt{N}, dim::Val{D}) where {N,D}
    map(eq(c), ntuple(static, dim))
end

function rank_to_sortperm(R::Tuple{Vararg{StaticInt,N}}) where {N}
    sp = ntuple(zero, Val{N}())
    r = ntuple(n -> sum(R[n] .≥ R), Val{N}())
    @inbounds for n = 1:N
        sp = Base.setindex(sp, n, r[n])
    end
    return sp
end

stride_rank(::Type{<:StrideIndex{N,R}}) where {N,R} = static(R)
stride_rank(x) = stride_rank(typeof(x))
function stride_rank(::Type{T}) where {T}
    is_forwarding_wrapper(T) ? stride_rank(parent_type(T)) : nothing
end
stride_rank(::Type{<:DenseArray{T,N}}) where {T,N} = ntuple(static, StaticInt(N))
stride_rank(::Type{BitArray{N}}) where {N} = ntuple(static, StaticInt(N))
stride_rank(::Type{<:AbstractRange}) = (One(),)
stride_rank(::Type{<:Tuple}) = (One(),)

stride_rank(::Type{T}) where {T<:VecAdjTrans} = (StaticInt(2), StaticInt(1))
stride_rank(::Type{T}) where {T<:MatAdjTrans} = _stride_rank(T, stride_rank(parent_type(T)))
_stride_rank(::Type{T}, ::Nothing) where {T<:MatAdjTrans} = nothing
function _stride_rank(::Type{T}, rank) where {T<:MatAdjTrans}
    return (getfield(rank, 2), getfield(rank, 1))
end

function stride_rank(::Type{T},) where {T<:PermutedDimsArray}
    return _stride_rank(T, stride_rank(parent_type(T)))
end
_stride_rank(::Type{T}, ::Nothing) where {T<:PermutedDimsArray} = nothing
_stride_rank(::Type{T}, r) where {T<:PermutedDimsArray} = permute(r, to_parent_dims(T))

stride_rank(::Type{T}) where {T<:SubArray} = _stride_rank(T, stride_rank(parent_type(T)))
_stride_rank(::Any, ::Any) = nothing
_stride_rank(::Type{T}, r::Tuple) where {T<:SubArray} = permute(r, to_parent_dims(T))

stride_rank(x, i) = stride_rank(x)[i]
function stride_rank(::Type{R}) where {T,N,S,A<:Array{S},R<:Base.ReinterpretArray{T,N,S,A}}
    return ntuple(static, StaticInt(N))
end
@inline function stride_rank(::Type{A}) where {NB,NA,B<:AbstractArray{<:Any,NB},A<:Base.ReinterpretArray{<:Any,NA,<:Any,B,true}}
    NA == NB ? stride_rank(B) : _stride_rank_reinterpret(stride_rank(B), gt(StaticInt{NB}(), StaticInt{NA}()))
end
@inline function stride_rank(::Type{A}) where {N,B<:AbstractArray{<:Any,N},A<:Base.ReinterpretArray{<:Any,N,<:Any,B,false}}
    stride_rank(B)
end

@inline _stride_rank_reinterpret(sr, ::False) = (One(), map(Base.Fix2(+, One()), sr)...)
@inline _stride_rank_reinterpret(sr::Tuple{One,Vararg}, ::True) = map(Base.Fix2(-, One()), tail(sr))
function contiguous_axis(::Type{R}) where {T,N,S,B<:AbstractArray{S,N},R<:ReinterpretArray{T,N,S,B,false}}
    contiguous_axis(B)
end
# if the leading dim's `stride_rank` is not one, then that means the individual elements are split across an axis, which ArrayInterface
# doesn't currently have a means of representing.
@inline function contiguous_axis(::Type{A}) where {NB,NA,B<:AbstractArray{<:Any,NB},A<:Base.ReinterpretArray{<:Any,NA,<:Any,B,true}}
    _reinterpret_contiguous_axis(stride_rank(B), dense_dims(B), contiguous_axis(B), gt(StaticInt{NB}(), StaticInt{NA}()))
end
@inline _reinterpret_contiguous_axis(::Any, ::Any, ::Any, ::False) = One()
@inline _reinterpret_contiguous_axis(::Any, ::Any, ::Any, ::True) = Zero()
@generated function _reinterpret_contiguous_axis(t::Tuple{One,Vararg{StaticInt,N}}, d::Tuple{True,Vararg{StaticBool,N}}, ::One, ::True) where {N}
    for n in 1:N
        if t.parameters[n+1].parameters[1] === 2
            if d.parameters[n+1] === True
                return :(StaticInt{$n}())
            else
                return :(Zero())
            end
        end
    end
    :(Zero())
end

function stride_rank(::Type{Base.ReshapedArray{T, N, P, Tuple{Vararg{Base.SignedMultiplicativeInverse{Int},M}}}}) where {T,N,P,M}
    _reshaped_striderank(is_column_major(P), Val{N}(), Val{M}())
end
function stride_rank(::Type{Base.ReshapedArray{T, 1, A, Tuple{}}}) where {T, A}
    IfElse.ifelse(is_column_major(A) & is_dense(A), (static(1),), nothing)
end
function stride_rank(::Type{Base.ReshapedArray{T, 1, LinearAlgebra.Adjoint{T, A}, Tuple{}}}) where {T, A <: AbstractVector{T}}
    IfElse.ifelse(is_dense(A), (static(1),), nothing)
end
function stride_rank(::Type{Base.ReshapedArray{T, 1, LinearAlgebra.Transpose{T, A}, Tuple{}}}) where {T, A <: AbstractVector{T}}
    IfElse.ifelse(is_dense(A), (static(1),), nothing)
end

_reshaped_striderank(::True, ::Val{N}, ::Val{0}) where {N} = ntuple(static, StaticInt(N))
_reshaped_striderank(_, __, ___) = nothing

"""
    contiguous_batch_size(::Type{T}) -> StaticInt{N}

Returns the Base.size of contiguous batches if `!isone(stride_rank(T, contiguous_axis(T)))`.
If `isone(stride_rank(T, contiguous_axis(T)))`, then it will return `StaticInt{0}()`.
If `contiguous_axis(T) == -1`, it will return `StaticInt{-1}()`.
If unknown, it will return `nothing`.
"""
contiguous_batch_size(x) = contiguous_batch_size(typeof(x))
contiguous_batch_size(::Type{T}) where {T} = _contiguous_batch_size(contiguous_axis(T), stride_rank(T))
_contiguous_batch_size(_, __) = nothing
function _contiguous_batch_size(::StaticInt{D}, ::R) where {D,R<:Tuple}
    if R.parameters[D].parameters[1] === 1
        return Zero()
    else
        return nothing
    end
end
_contiguous_batch_size(::StaticInt{-1}, ::R) where {R<:Tuple} = -One()

contiguous_batch_size(::Type{Array{T,N}}) where {T,N} = Zero()
contiguous_batch_size(::Type{BitArray{N}}) where {N} = Zero()
contiguous_batch_size(::Type{<:AbstractRange}) = Zero()
contiguous_batch_size(::Type{<:Tuple}) = Zero()
function contiguous_batch_size(::Type{T}) where {T<:Union{Transpose,Adjoint}}
    return contiguous_batch_size(parent_type(T))
end
function contiguous_batch_size(::Type{T}) where {T<:PermutedDimsArray}
    return contiguous_batch_size(parent_type(T))
end
function contiguous_batch_size(::Type{S}) where {N,NP,T,A<:AbstractArray{T,NP},I,S<:SubArray{T,N,A,I}}
    return _contiguous_batch_size(S, contiguous_batch_size(A), contiguous_axis(A))
end
_contiguous_batch_size(::Any, ::Any, ::Any) = nothing
function _contiguous_batch_size(::Type{<:SubArray{T,N,A,I}}, b::StaticInt{B}, c::StaticInt{C}) where {T,N,A,I,B,C}
    if I.parameters[C] <: AbstractUnitRange
        return b
    else
        return -One()
    end
end
contiguous_batch_size(::Type{<:Base.ReinterpretArray{T,N,S,A}}) where {T,N,S,A} = Zero()

"""
    is_column_major(A) -> StaticBool

Returns `True()` if elements of `A` are stored in column major order. Otherwise returns `False()`.
"""
is_column_major(A) = is_column_major(stride_rank(A), contiguous_batch_size(A))
is_column_major(sr::Nothing, cbs) = False()
is_column_major(sr::R, cbs) where {R} = _is_column_major(sr, cbs)
is_column_major(::AbstractRange) = False()

# cbs > 0
_is_column_major(sr::R, cbs::StaticInt) where {R} = False()
# cbs <= 0
_is_column_major(sr::R, cbs::Union{StaticInt{0},StaticInt{-1}}) where {R} = is_increasing(sr)

"""
    dense_dims(::Type{<:AbstractArray{N}}) -> Tuple{Vararg{StaticBool,N}}

Returns a tuple of indicators for whether each axis is dense.
An axis `i` of array `A` is dense if `stride(A, i) * Base.size(A, i) == stride(A, j)`
where `stride_rank(A)[i] + 1 == stride_rank(A)[j]`.
"""
dense_dims(x) = dense_dims(typeof(x))
function dense_dims(::Type{T}) where {T}
    is_forwarding_wrapper(T) ? dense_dims(parent_type(T)) : nothing
end
_all_dense(::Val{N}) where {N} = ntuple(_ -> True(), Val{N}())

dense_dims(::Type{<:DenseArray{T,N}}) where {T,N} = _all_dense(Val{N}())
dense_dims(::Type{BitArray{N}}) where {N} = _all_dense(Val{N}())
dense_dims(::Type{<:AbstractRange}) = (True(),)
dense_dims(::Type{<:Tuple}) = (True(),)
function dense_dims(::Type{T}) where {T<:VecAdjTrans}
    dense = dense_dims(parent_type(T))
    if dense === nothing
        return nothing
    else
        return (True(), first(dense))
    end
end
function dense_dims(::Type{T}) where {T<:MatAdjTrans}
    dense = dense_dims(parent_type(T))
    if dense === nothing
        return nothing
    else
        return (last(dense), first(dense))
    end
end
function dense_dims(::Type{T}) where {T<:PermutedDimsArray}
    dense = dense_dims(parent_type(T))
    if dense === nothing
        return nothing
    else
        return permute(dense, to_parent_dims(T))
    end
end
function dense_dims(::Type{S}) where {N,NP,T,A<:AbstractArray{T,NP},I,S<:SubArray{T,N,A,I}}
    return _dense_dims(S, dense_dims(A), Val(stride_rank(A)))
end
@inline function dense_dims(::Type{A}) where {NB, NA, B <: AbstractArray{<:Any,NB},A<: Base.ReinterpretArray{<:Any, NA, <:Any, B, true}}
    ddb = dense_dims(B)
    IfElse.ifelse(Static.le(StaticInt(NB), StaticInt(NA)), (True(), ddb...), Base.tail(ddb))
end
@inline function dense_dims(::Type{A}) where {NB, NA, B <: AbstractArray{<:Any,NB},A<: Base.ReinterpretArray{<:Any, NA, <:Any, B, false}}
    dense_dims(B)
end

_dense_dims(::Type{S}, ::Nothing, ::Val{R}) where {R,N,NP,T,A<:AbstractArray{T,NP},I,S<:SubArray{T,N,A,I}} = nothing
@generated function _dense_dims(
    ::Type{S},
    ::D,
    ::Val{R},
) where {D,R,N,NP,T,A<:AbstractArray{T,NP},I,S<:SubArray{T,N,A,I}}
    still_dense = true
    sp = rank_to_sortperm(R)
    densev = Vector{Bool}(undef, NP)
    for np = 1:NP
        spₙ = sp[np]
        if still_dense
            still_dense = D.parameters[spₙ] <: True
        end
        densev[spₙ] = still_dense
        # a dim not being complete makes later dims not dense
        still_dense &= (I.parameters[spₙ] <: Base.Slice)::Bool
    end
    dense_tup = Expr(:tuple)
    for np = 1:NP
        Iₙₚ = I.parameters[np]
        if Iₙₚ <: AbstractUnitRange
            if densev[np]
                push!(dense_tup.args, :(True()))
            else
                push!(dense_tup.args, :(False()))
            end
        elseif Iₙₚ <: AbstractVector
            push!(dense_tup.args, :(False()))
        end
    end
    # If n != N, then an axis was indexed by something other than an integer or `AbstractUnitRange`, so we return `nothing`.
    if length(dense_tup.args) === N
        return dense_tup
    else
        return nothing
    end
end

function dense_dims(T::Type{<:Base.ReshapedArray})
    d = dense_dims(parent_type(T))
    if d === nothing
        return nothing
    elseif all(d)
        return n_of_x(StaticInt(ndims(T)), True())
    else
        return n_of_x(StaticInt(ndims(T)), False())
    end
end

is_dense(A) = is_dense(typeof(A))
is_dense(::Type{A}) where {A} = _is_dense(dense_dims(A))
_is_dense(::Tuple{False,Vararg}) = False()
_is_dense(t::Tuple{True,Vararg}) = _is_dense(Base.tail(t))
_is_dense(t::Tuple{True}) = True()
_is_dense(t::Tuple{}) = True()
_is_dense(::Nothing) = False()

"""
    known_strides(::Type{T}) -> Tuple
    known_strides(::Type{T}, dim) -> Union{Int,Nothing}

Returns the strides of array `A` known at compile time. Any strides that are not known at
compile time are represented by `nothing`.
"""
known_strides(x, dim) = known_strides(typeof(x), dim)
known_strides(::Type{T}, dim) where {T} = known_strides(T, to_dims(T, dim))
function known_strides(::Type{T}, dim::CanonicalInt) where {T}
    # see https://github.com/JuliaLang/julia/blob/6468dcb04ea2947f43a11f556da9a5588de512a0/base/reinterpretarray.jl#L148
    if ndims(T) < dim
        return known_length(T)
    else
        return known_strides(T)[dim]
    end
end
known_strides(::Type{<:StrideIndex{N,R,C,S,O}}) where {N,R,C,S,O} = known(S)

known_strides(x) = known_strides(typeof(x))
known_strides(::Type{T}) where {T<:Vector} = (1,)
function known_strides(::Type{T}) where {T<:MatAdjTrans}
    permute(known_strides(parent_type(T)), to_parent_dims(T))
end
@inline function known_strides(::Type{T}) where {T<:VecAdjTrans}
    strd = first(known_strides(parent_type(T)))
    return (strd, strd)
end
@inline function known_strides(::Type{T}) where {T<:PermutedDimsArray}
    permute(known_strides(parent_type(T)), to_parent_dims(T))
end
@inline function known_strides(::Type{T}) where {T<:SubArray}
    permute(known_strides(parent_type(T)), to_parent_dims(T))
end
function known_strides(::Type{T}) where {T}
    if ndims(T) === 1
        return (1,)
    else
        return size_to_strides(known_size(T), 1)
    end
end

"""
    strides(A) -> Tuple{Vararg{Union{Int,StaticInt}}}
    strides(A, dim) -> Union{Int,StaticInt}

Returns the strides of array `A`. If any strides are known at compile time,
these should be returned as `Static` numbers. For example:
```julia
julia> A = rand(3,4);

julia> ArrayInterface.strides(A)
(static(1), 3)
```

Additionally, the behavior differs from `Base.strides` for adjoint vectors:

```julia
julia> x = rand(5);

julia> ArrayInterface.strides(x')
(static(1), static(1))
```

This is to support the pattern of using just the first stride for linear indexing, `x[i]`,
while still producing correct behavior when using valid cartesian indices, such as `x[1,i]`.
```
"""
strides(A::StrideIndex) = getfield(A, :strides)
@inline strides(A::Vector{<:Any}) = (StaticInt(1),)
@inline strides(A::Array{<:Any,N}) where {N} = (StaticInt(1), Base.tail(Base.strides(A))...)
@inline function strides(x::X) where {X}
    if is_forwarding_wrapper(X)
        return strides(parent(x))
    elseif defines_strides(X)
        return size_to_strides(size(x), One())
    else
        return Base.strides(x)
    end
end

_is_column_dense(::A) where {A<:AbstractArray} =
    defines_strides(A) &&
    (ndims(A) == 0 || Bool(is_dense(A)) && Bool(is_column_major(A)))

# Fixes the example of https://github.com/JuliaArrays/ArrayInterfaceCore.jl/issues/160
function strides(A::ReshapedArray)
    if _is_column_dense(parent(A))
        return size_to_strides(size(A), One())
    else
        pst = strides(parent(A))
        psz = size(parent(A))
        # Try dimension merging in order (starting from dim1).
        # `sz1` and `st1` are the `size`/`stride` of dim1 after dimension merging.
        # `n` indicates the last merged dimension.
        # note: `st1` should be static if possible
        sz1, st1, n = merge_adjacent_dim(psz, pst)
        n == ndims(A.parent) && return size_to_strides(size(A), st1)
        return _reshaped_strides(size(A), One(), sz1, st1, n, Dims(psz), Dims(pst))
    end
end

@inline function _reshaped_strides(::Dims{0}, reshaped, msz::Int, _, ::Int, ::Dims, ::Dims)
    reshaped == msz && return ()
    throw(ArgumentError("Input is not strided."))
end
function _reshaped_strides(asz::Dims, reshaped, msz::Int, mst, n::Int, apsz::Dims, apst::Dims)
    st = reshaped * mst
    reshaped = reshaped * asz[1]
    if length(asz) > 1 && reshaped == msz && asz[2] != 1
        msz, mst′, n = merge_adjacent_dim(apsz, apst, n + 1)
        reshaped = 1
    else
        mst′ = Int(mst)
    end
    sts = _reshaped_strides(tail(asz), reshaped, msz, mst′, n, apsz, apst)
    return (st, sts...)
end

merge_adjacent_dim(::Tuple{}, ::Tuple{}) = 1, One(), 0
merge_adjacent_dim(szs::Tuple{Any}, sts::Tuple{Any}) = Int(szs[1]), sts[1], 1
function merge_adjacent_dim(szs::Tuple, sts::Tuple)
    if szs[1] isa One # Just ignore dimension with size 1
        sz, st, n = merge_adjacent_dim(tail(szs), tail(sts))
        return sz, st, n + 1
    elseif szs[2] isa One # Just ignore dimension with size 1
        sz, st, n = merge_adjacent_dim((szs[1], tail(tail(szs))...), (sts[1], tail(tail(sts))...))
        return sz, st, n + 1
    elseif (szs[1], szs[2], sts[1], sts[2]) isa NTuple{4,StaticInt} # the check could be done during compiling.
        if sts[2] == sts[1] * szs[1]
            szs′ = (szs[1] * szs[2], tail(tail(szs))...)
            sts′ = (sts[1], tail(tail(sts))...)
            sz, st, n = merge_adjacent_dim(szs′, sts′)
            return sz, st, n + 1
        else
            return Int(szs[1]), sts[1], 1
        end
    else # the check can't be done during compiling.
        sz, st, n = merge_adjacent_dim(Dims(szs), Dims(sts), 1)
        if (szs[1], sts[1]) isa NTuple{2,StaticInt} && szs[1] != 1
            # But the 1st stride might still be static.
            return sz, sts[1], n
        else
            return sz, st, n
        end
    end
end

function merge_adjacent_dim(psz::Dims{N}, pst::Dims{N}, n::Int) where {N}
    sz, st = psz[n], pst[n]
    while n < N
        szₙ, stₙ = psz[n+1], pst[n+1]
        if sz == 1
            sz, st = szₙ, stₙ
        elseif stₙ == st * sz
            sz *= szₙ
        elseif szₙ != 1
            break
        end
        n += 1
    end
    return sz, st, n
end

# `strides` for `Base.ReinterpretArray`
function strides(A::Base.ReinterpretArray{T,<:Any,S,<:AbstractArray{S},IsReshaped}) where {T,S,IsReshaped}
    _is_column_dense(parent(A)) && return size_to_strides(size(A), One())
    stp = strides(parent(A))
    ET, ES = static(sizeof(T)), static(sizeof(S))
    ET === ES && return stp
    IsReshaped && ET < ES && return (One(), _reinterp_strides(stp, ET, ES)...)
    first(stp) == 1 || throw(ArgumentError("Parent must be contiguous in the 1st dimension!"))
    if IsReshaped
        # The wrapper tell us `A`'s parent has static size in dim1.
        # We can make the next stride static if the following dim is still dense.
        sr = stride_rank(parent(A))
        dd = dense_dims(parent(A))
        stp′ = _new_static(stp, sr, dd, ET ÷ ES)
        return _reinterp_strides(tail(stp′), ET, ES)
    else
        return (One(), _reinterp_strides(tail(stp), ET, ES)...)
    end
end
_new_static(P,_,_,_) = P # This should never be called, just in case.
@generated function _new_static(p::P, ::SR, ::DD, ::StaticInt{S}) where {S,N,P<:NTuple{N,Union{Int,StaticInt}},SR<:NTuple{N,StaticInt},DD<:NTuple{N,StaticBool}}
    sr = fieldtypes(SR)
    j = findfirst(T -> T() == sr[1]()+1, sr)
    if !isnothing(j) && !(fieldtype(P, j) <: StaticInt) && fieldtype(DD, j) === True
        return :(tuple($((i == j ? :(static($S)) : :(p[$i]) for i in 1:N)...)))
    else
        return :(p)
    end
end
@inline function _reinterp_strides(stp::Tuple, els::StaticInt, elp::StaticInt)
    if elp % els == 0
        N = elp ÷ els
        return map(Base.Fix2(*, N), stp)
    else
        return map(stp) do i
            d, r = divrem(elp * i, els)
            iszero(r) || throw(ArgumentError("Parent's strides could not be exactly divided!"))
            d
        end
    end
end

strides(::AbstractRange) = (One(),)
function strides(x::VecAdjTrans)
    st = first(strides(parent(x)))
    return (st, st)
end
@inline strides(B::MatAdjTrans) = permute(strides(parent(B)), to_parent_dims(B))
@inline strides(B::PermutedDimsArray) = permute(strides(parent(B)), to_parent_dims(B))

getmul(x::Tuple, y::Tuple, ::StaticInt{i}) where {i} = getfield(x, i) * getfield(y, i)
function strides(A::SubArray)
    return eachop(getmul, to_parent_dims(A), map(maybe_static_step, A.indices), strides(parent(A)))
end

maybe_static_step(x::AbstractRange) = static_step(x)
maybe_static_step(_) = nothing

@generated function size_to_strides(sz::S, init) where {N,S<:Tuple{Vararg{Any,N}}}
    out = Expr(:block, Expr(:meta, :inline))
    t = Expr(:tuple, :init)
    prev = :init
    i = 1
    while i <= (N - 1)
        if S.parameters[i] <: Nothing || (i > 1 &&  t.args[i - 1] === :nothing)
            push!(t.args, :nothing)
        else
            next = Symbol(:val_, i)
            push!(out.args, :($next = $prev * getfield(sz, $i)))
            push!(t.args, next)
            prev = next
        end
        i += 1
    end
    push!(out.args, t)
    return out
end

strides(a, dim) = strides(a, to_dims(a, dim))
function strides(a::A, dim::CanonicalInt) where {A}
    if is_forwarding_wrapper(A)
        return strides(parent(a), to_parent_dims(A, dim))
    else
        return Base.stride(a, Int(dim))
    end
end

@inline stride(A::AbstractArray, ::StaticInt{N}) where {N} = strides(A)[N]
@inline stride(A::AbstractArray, ::Val{N}) where {N} = strides(A)[N]
stride(A, i) = Base.stride(A, i) # for type stability


text:
using Combinatorics

include("library.jl")

f = open("input-01-a")
lines = parse.(Int, readlines(f))

println("\nSolutions for Day 1")
@assert (@show sum(fuelForModule.(lines))) == 3380880
@assert (@show sum(rocketEqFuel.(lines))) == 5068454

println("\nSolutions for Day 2")

inputDay2 = parse.(Int, split(read("input-02", String), ","))
vm2 = ElfVM(inputDay2)
@assert 3654868 == @show runWithArguments(copy(vm2), 12, 2)
@assert (70, 14) == @show findArguments(copy(vm2), 19690720)

println("\nSolutions for Day 3")

# Puzzle
wires = createWires.(parseDirectionString.(readlines(open("input-03"))))
intersections_puzzle = intersections(wires[1], wires[2])
@assert 352 == @show minimum(manhattan.(intersections_puzzle))
@assert 43848 == @show minimum(intersection_offset.(intersections_puzzle))

println("\nSolutions for Day 4")

@assert 1019 == @show length(filter(check_password, 248345:746315))
@assert 660 == @show length(filter(check_password_stronger, 248345:746315))

println("\nSolutions for Day 5")
inputDay5 = parse.(Int, split(read("input-05", String), ","))
vm5 = ElfVM(copy(inputDay5))
push!(vm5.stdin, 1)
run!(vm5)
@assert [0, 0, 0, 0, 0, 0, 0, 0, 0, 13818007] == @show vm5.stdout

inputDay5 = parse.(Int, split(read("input-05", String), ","))
vm5_b = ElfVM(copy(inputDay5))
push!(vm5_b.stdin, 5)
run!(vm5_b)
@assert 3176266 == @show vm5_b.stdout[1]

println("\nSolutions for Day 6")
orbitMap = OrbitMap(readlines(open("input-06")))
@assert 271151 == @show sum(values(orbitMap.depth))
@assert 388 == @show distance(orbitMap, orbitMap.centers["YOU"], orbitMap.centers["SAN"])

println("\nSolutions for Day 7")
# Amplifier Controller Software
inputDay7 = parse.(Int, split(read("input-07", String), ","))

function day7_vm1_runner(config::Vector)::Int
    # Set up the VM, link the threads.
    vm = MultiThreadElfVM(
        inputDay7,
        [1, 2, 3, 4, 5],
        Dict( 1 => 2, 2 => 3, 3 => 4, 4 => 5 ))

    # Provide the threads with input.
    map(enumerate(config)) do (pid, val) push_stdin!(vm, pid, val) end
    push_stdin!(vm, 1, 0)

    # Run and extract result.
    run!(vm)
    vm.threads[5].stdout[1]
end

@assert 116680 == @show maximum(day7_vm1_runner.(permutations([0, 1, 2, 3, 4])))

function day7_vm2_runner(config::Vector)::Int
    # Set up the VM, link the threads.
    day7_vm2 = MultiThreadElfVM(
        inputDay7,
        [1, 2, 3, 4, 5],
        Dict( 1 => 2, 2 => 3, 3 => 4, 4 => 5, 5 => 1 ))

    # Provide the threads with input.
    map(enumerate(config)) do (pid, val) push_stdin!(day7_vm2, pid, val) end
    push_stdin!(day7_vm2, 1, 0)

    # Run and extract result.
    run!(day7_vm2)
    day7_vm2.threads[1].stdin[1]
end

@assert 89603079 == @show maximum(day7_vm2_runner.(permutations([5, 6, 7, 9, 8])))

println("\nSolutions for Day 8")

# Read data, remove terminating newline, convert to Vector{Int}
input_8 = parse.(Int, collect(read("input-08", String)[1:end-1]))
# Convert to a Matrix{Int, 3}. Note that the julia REPL outputs x and y flipped.
image_depth = length(input_8) ÷ 25 ÷ 6
image_data = reshape(input_8, (25, 6, image_depth ))

# [..] the Elves would like you to find the layer that contains the fewest 0 digits.
layer_digit_count(data, layer, digit) = length(filter(x -> x == digit, data[:,:,layer]))
zeros = layer_digit_count.(Ref(image_data), 1:image_depth, 0)
layer_with_least_zeros = findall(zeros .== minimum(zeros))[1]

# On that layer, what is the number of 1 digits multiplied by the number of 2 digits?
@assert 1820 == @show layer_digit_count(image_data, layer_with_least_zeros, 1) * layer_digit_count(image_data, layer_with_least_zeros, 2)

# What message is produced after decoding your image?
pixel_color(data, x, y)::Int = filter(c -> c != 2, data[x, y, :])[1]
image_decoded = pixel_color.(Ref(image_data), transpose(1:25), 1:6)
# Solution read manually from image output: ZUKCJ

println("\nSolutions for Day 9")

cpy_code = [109,1,204,-1,1001,100,1,100,1008,100,16,101,1006,101,0,99]
cpy_vm = ElfVM(copy(cpy_code))
run!(cpy_vm)
@assert cpy_code == cpy_vm.stdout

# Load the boost programm into a vm
boost_code = parse.(Int, split(read("input-09", String), ","))
boost_vm = ElfVM(copy(boost_code))
push!(boost_vm.stdin, 1)
run!(boost_vm)
@assert 3638931938 == @show boost_vm.stdout[1]

boost_vm = ElfVM(copy(boost_code))
push!(boost_vm.stdin, 2)
run!(boost_vm)
@assert 86025 == @show boost_vm.stdout[1]

println("\nSolutions for Day 10")

file_10 = filter(c -> c == '#' || c == '.', read("input-10", String))
vec_10 = map(c -> c .== '#', collect(file_10))
# The asteroid field is square
side_length = Int(sqrt(length(vec_10)))
asteroid_field = reshape(vec_10, (side_length, side_length))

function free_view(field, ax, ay, bx, by)::Bool
    steps = gcd(bx-ax, by-ay)
    # We don't count the asteroid itself, only other asteroids.
    # This also prevents us from dividing by zero.
    if steps == 0 return false end
    dx = (bx-ax) ÷ steps
    dy = (by-ay) ÷ steps
    for i in 1:(steps - 1)
        # Check if the possible location has an asteroid.
        if field[ax + i * dx, ay + i * dy]
            return false
        end
    end
    true
end

function count_visible_asteroids(field, cx, cy)::Int
    counter = 0
    for x in 1:side_length, y in 1:side_length
        if field[x, y] && free_view(field, cx, cy, x, y)
            counter += 1
        end
    end
    counter
end

function get_visible_asteroids(field, cx, cy)::Vector
    asteroid_list = []
    for x in 1:side_length, y in 1:side_length
        if field[x, y] && free_view(field, cx, cy, x, y)
            push!(asteroid_list, (x, y))
        end
    end
    asteroid_list
end

# Test how many asteroids you can see from this position
function day_10()
    max_count = 0
    best = (0, 0)
    for x in 1:side_length, y in 1:side_length
        if asteroid_field[x, y]
            new_count = count_visible_asteroids(asteroid_field, x, y)
            if new_count > max_count
                max_count = new_count
                best = (x, y)
            end
        end
    end
    (max_count, best)
end

choosen_asteroid = day_10()

@assert 267 == @show choosen_asteroid[1]

# Part 2 of day 10, now we need to determine the order in which the asteroids
# are vaporized.

# I am able to see 267 asteroids, this means I only have to order those.
# I won't even get to a second rotation.

visible_asteroids = get_visible_asteroids(asteroid_field, choosen_asteroid[2]...)

sweep_angle((x, y)) = atan(-((x)+0.000001), y)
sweep_angle_from_laser((x, y)) = sweep_angle((x - choosen_asteroid[2][1], y - choosen_asteroid[2][2]))

sorted_asteroids = sort(visible_asteroids, by=sweep_angle_from_laser)
@assert (14, 10) == @show sorted_asteroids[200]


text:
################################################################################

"Backtrace on hidden Markov model object"
function backTrace(self::HMM)
  state = 0
  min = self.tbM[1][end]

  for ix in eachindex(self.tbM)
    if self.tbM[ix][end] <= min
      state = ix
      min = self.tbM[ix][end]
    end

    tb = fill(-1, size(self.tbM[1], 1))

    for ix in reverse(eachindex(self.tbM[1]))
      vlocal = min
      global sig = -1

      for jx in eachindex(self.tbM)
        dob = self.tbM[jx][ix]
        if dob <= vlocal
          vlocal = dob
          sig = jx
        end
      end

      tb[ix] = state
      state = sig

    end

    tb = tb[1:end - 1]

    return tb
  end
end

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

max_stars_repo_name: grabherr/HiddenMarkovModelReaders

max_issues_repo_name: grabherr/HiddenMarkovModelReaders

max_forks_repo_name: grabherr/HiddenMarkovModelReaders


text:
#Gram-Schmidt Orthogonalization

function normalized_orthogonal_projection(b,Z::AbstractArray)
    annihilator =  I - Z * inv(Z'*Z) *Z'
    projection = annihilator*b
    return projection / norm(projection)
end

function gram_schmidt(x::AbstractArray)

    U = similar(x,Float64)

    for col in 1:size(U,2)
        b = X[:,col]       # vector we're going to project
        Z = X[:,1:col - 1] # first i-1 columns of X
        U[:,col] = normalized_orthogonal_projection(b, Z)
   end
   return U
end

y = [1, 3, -3]
X = [1 0; 0 -6; 2 2];

Py1 = X * inv(X'X) * X' * y #ordinary projection

U = gram_schmidt(X) #construct orthogonal set

Py2 = U*U'*y #use orthonormal projection result

#perform the same projection using QR decomposition
Q, R = qr(X)
Q = Matrix(Q)
R = Matrix(R)

Py3 = Q*Q'*y # Q is formed of columns from orthonormal basis of x
beta = inv(R)*Q'*y

max_stars_repo_name: shanemcmiken/quantecon-notebooks-julia

max_issues_repo_name: shanemcmiken/quantecon-notebooks-julia

max_forks_repo_name: shanemcmiken/quantecon-notebooks-julia


text:
function moduloexponentiation(a, b, m)
    ret = 1
    a = a%m
    while b > 0
        if b&1 == 1
            ret *= a
            ret %= m
        end
        b = b >> 1
        a = (a*a) % m
    end
    ret
end

function solve()
    n, k = [parse(Int, x) for x in split(readline())]
    m = 1000000007
    ret = 0
    numsequence = zeros(Int, k)
    for c=k:-1:1
        ns = moduloexponentiation(k÷c, n, m)
        for i=2:(k÷c)
            ns -= numsequence[c*i]
            while ns < 0
                ns += m
            end
        end
        ns %= m
        numsequence[c] = ns
        ret += ns*c
        ret %= m
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
# File extensions. GDAL is the catch-all for everything else
const EXT = Dict(
    GRDfile => (".grd", ".gri"),
    NCDfile => (".nc",),
    SMAPfile => (".h5",)
)
const REV_EXT = Dict(
    ".grd" => GRDfile,
    ".gri" => GRDfile,
    ".nc" => NCDfile,
    ".h5" => SMAPfile
)

# Get the source backend for a file extension, falling back to GDALfile
_sourcetype(filename::AbstractString) = get(REV_EXT, splitext(filename)[2], GDALfile)
_sourcetype(filenames::NamedTuple) = _sourcetype(first(filenames))

# Internal read method
function _open(f, filename::AbstractString; kw...)
    _open(f, _sourcetype(filename), filename; kw...)
end


text:
# coding: utf-8

# In[1]:

#fractal is one of the interesting topics in geometry
#it is usually described by a recursive function
#voila,here we are!
using Plots

# In[2]:

#initialize
x1=0
x2=3
y=0
bar_height=5
between_interval=10
n=6;

# In[3]:

#create rectangle shape
rectangle(w,h,x,y)=Shape(x.+[0,w,w,0],y.-[0,0,h,h])

# In[4]:

#cantor set is one of the simplest fractal shape
#at each iteration,we divide each bar into three equal parts
#we remove the mid part from each bar and keep the rest
#this effectively creates a binary tree
#check the link below for more details on math
# https://www.math.stonybrook.edu/~scott/Book331/Cantor_sets.html
function cantor_set(x1,x2,y,n,
               bar_height=5,between_interval=10)

    #base case
    if n==0
        return

    else

        #viz the first 1/3 and the last 1/3
        plot!(rectangle((x2-x1)/3,bar_height,x1,y-between_interval))
        plot!(rectangle((x2-x1)/3,bar_height,x2-(x2-x1)/3,y-between_interval))

        #recursion
        cantor_set(x1,x1+(x2-x1)/3,
                   y-between_interval,
                   n-1)
        cantor_set(x2-(x2-x1)/3,x2,
                   y-between_interval,
                   n-1)

    end

end

# In[5]:

#viz
#as n increases
#the bar gets too slim to be visible
gr(size=(250,250))
fig=plot(legend=false,grid=false,axis=false,ticks=false)
plot!(rectangle((x2-x1),bar_height,x1,y))
cantor_set(x1,x2,y,n)
fig

max_stars_repo_path: cantor set.jl

max_issues_repo_path: cantor set.jl

max_forks_repo_path: cantor set.jl


text:
export MISSolverType

"""
# Description
    MISSolverType(;
        splitting_type = SlowFastSplitting(),
        fast_model = AtmosAcousticGravityLinearModel,
        mis_method = MIS2,
        fast_method = LSRK54CarpenterKennedy,
        nsubsteps = 50,
    )

This solver type constructs an ODE solver using a generalization
of the split-explicit Runge-Kutta method. Known as the Multirate
Infinitesimal Step (MIS) method, this solver solves ODEs with
the partitioned form:

```math
    \\dot{Q} = f_{fast}(Q, t) + f_{slow}(Q, t)
```

where the right-hand-side functions `f_fast` and `f_slow` denote
fast and slow dynamics respectively, depending on the state `Q`.

# Arguments
- `splitting_type` (DiscreteSplittingType): The type of discrete
    splitting to apply to the right-hand side.
    Default: `SlowFastSplitting()`
- `fast_model` (Type): The model describing fast dynamics.
    Default: `AtmosAcousticGravityLinearModel`
- `mis_method` (Function): Function defining the particular MIS
    method to be used.
    Default: `MIS2`
- `fast_method` (Function): Function defining the fast solver.
    Default: `LSRK54CarpenterKennedy`
- `nsubsteps` (Tuple): Tuple denoting the total number of times
    to substep the fast process.
    Default: `(50,)`
- `discrete_splitting` (Boolean): Boolean denoting whether a PDE level or
    discretized level splitting should be used. If `true` then the PDE is
    discretized in such a way that `f_fast + f_slow` is equivalent to
    discretizing the original PDE directly.
    Default: `false`

### References
 - [KnothWensch2014](@cite)
"""
struct MISSolverType{DS} <: AbstractSolverType
    # The type of discrete splitting to apply to the right-hand side
    splitting_type::DS
    # The model describing fast dynamics
    fast_model::Type
    # Main MIS function
    mis_method::Function
    # Fast RK solver
    fast_method::Function
    # Substepping parameter for the fast processes
    nsubsteps::Tuple
    # Whether to use a PDE level or discrete splitting
    discrete_splitting::Bool

    function MISSolverType(;
        splitting_type = SlowFastSplitting(),
        fast_model = AtmosAcousticGravityLinearModel,
        mis_method = MIS2,
        fast_method = LSRK54CarpenterKennedy,
        nsubsteps = (50,),
        discrete_splitting = false,
    )

        DS = typeof(splitting_type)

        return new{DS}(
            splitting_type,
            fast_model,
            mis_method,
            fast_method,
            nsubsteps,
            discrete_splitting,
        )
    end
end

"""
    getdtmodel(ode_solver::MISSolverType, bl)

A function which returns a model representing the dynamics
with the most restrictive time-stepping requirements.
"""
function getdtmodel(ode_solver::MISSolverType, bl)
    # Most restrictive dynamics are part of the fast model
    return ode_solver.fast_model(bl)
end

"""
# Description
    function solversetup(
        ode_solver::MISSolverType{SlowFastSplitting},
        dg,
        Q,
        dt,
        t0,
        diffusion_direction,
    )

Creates an ODE solver for the partition slow-fast ODE
using an MIS method with explicit time-integration.
The splitting of the fast (acoustic and gravity waves)
dynamics is done in _all_ spatial directions.
"""
function solversetup(
    ode_solver::MISSolverType{SlowFastSplitting},
    dg,
    Q,
    dt,
    t0,
    diffusion_direction,
)

    # Extract fast model and define a DG model
    # for the fast processes (acoustic/gravity waves
    # in all spatial directions)
    fast_model = ode_solver.fast_model(dg.balance_law)
    fast_dg = DGModel(
        fast_model,
        dg.grid,
        dg.numerical_flux_first_order,
        dg.numerical_flux_second_order,
        dg.numerical_flux_gradient,
        state_auxiliary = dg.state_auxiliary,
        state_gradient_flux = dg.state_gradient_flux,
        states_higher_order = dg.states_higher_order,
        direction = EveryDirection(),
    )

    # Using the RemainderModel, we subtract away the
    # fast processes and define a DG model for the
    # slower processes (advection and diffusion)
    if ode_solver.discrete_splitting
        numerical_flux_first_order =
            (dg.numerical_flux_first_order, (dg.numerical_flux_first_order,))
    else
        numerical_flux_first_order = dg.numerical_flux_first_order
    end
    slow_dg = remainder_DGModel(
        dg,
        (fast_dg,);
        numerical_flux_first_order = numerical_flux_first_order,
    )

    solver = ode_solver.mis_method(
        slow_dg,
        fast_dg,
        ode_solver.fast_method,
        ode_solver.nsubsteps[1],
        Q;
        dt = dt,
        t0 = t0,
    )

    return solver
end

function solversetup(
    ode_solver::MISSolverType{HEVISplitting},
    dg,
    Q,
    dt,
    t0,
    diffusion_direction,
)

    # Extract fast model and define a DG model
    # for the fast processes (acoustic/gravity waves
    # in all spatial directions)
    fast_model = ode_solver.fast_model(dg.balance_law)
    fast_method = ode_solver.fast_method
    fast_dg = DGModel(
        fast_model,
        dg.grid,
        dg.numerical_flux_first_order,
        dg.numerical_flux_second_order,
        dg.numerical_flux_gradient,
        state_auxiliary = dg.state_auxiliary,
        state_gradient_flux = dg.state_gradient_flux,
        states_higher_order = dg.states_higher_order,
        direction = EveryDirection(),
    )

    # Using the RemainderModel, we subtract away the
    # fast processes and define a DG model for the
    # slower processes (advection and diffusion)
    if ode_solver.discrete_splitting
        numerical_flux_first_order =
            (dg.numerical_flux_first_order, (dg.numerical_flux_first_order,))
    else
        numerical_flux_first_order = dg.numerical_flux_first_order
    end
    slow_dg = remainder_DGModel(
        dg,
        (fast_dg,);
        numerical_flux_first_order = numerical_flux_first_order,
    )

    # Defining horizontal and vertical fast tendencies
    fast_dg_h = DGModel(
        fast_model,
        dg.grid,
        dg.numerical_flux_first_order,
        dg.numerical_flux_second_order,
        dg.numerical_flux_gradient,
        state_auxiliary = dg.state_auxiliary,
        state_gradient_flux = dg.state_gradient_flux,
        states_higher_order = dg.states_higher_order,
        direction = HorizontalDirection(),
    )
    fast_dg_v = DGModel(
        fast_model,
        dg.grid,
        dg.numerical_flux_first_order,
        dg.numerical_flux_second_order,
        dg.numerical_flux_gradient,
        state_auxiliary = dg.state_auxiliary,
        state_gradient_flux = dg.state_gradient_flux,
        states_higher_order = dg.states_higher_order,
        direction = VerticalDirection(),
    )

    fast_dg_tendencies = (fast_dg_h, fast_dg_v)

    if length(ode_solver.nsubsteps) == 1
        nsubsteps = getnsubsteps(
            ode_solver.mis_method,
            ode_solver.nsubsteps[1],
            real(eltype(Q)),
        )
        fast_method =
            (dg, Q) -> ode_solver.fast_method(
                dg,
                Q,
                dt / ode_solver.nsubsteps[1],
                nsubsteps,
            )
    elseif length(ode_solver.nsubsteps) == 2
        fast_method =
            (dg, Q) -> ode_solver.fast_method(dg, Q, ode_solver.nsubsteps[2])
    end

    solver = ode_solver.mis_method(
        slow_dg,
        fast_dg_tendencies,
        fast_method,
        ode_solver.nsubsteps[1],
        Q;
        dt = dt,
        t0 = t0,
    )
    return solver
end

max_forks_repo_name: jm-c/CLIMA


text:
module IEEE_754_2008

import Base: min, max, minmax, precision, ldexp, frexp

include("extensions.jl")
include("ulpufp.jl")

# include bypasses precompilation
# modulename = :QNaN; @eval begin import $modulename; using $modulename end
modulename = :QNaN; @eval importall $modulename
modulename = :SignificantBits; @eval importall $modulename

export negate, copySign,
    isSignMinus, isNormal, isFinite, is Zero, isSubnormal, isInfinite, isNaN,
    scaleB, logB,
    compareQuietEqual, compareQuietNotEqual,
    totalOrder, totalOrderMag,
    minNum, maxNum, minNumMag, maxNumMag,
    # extensions
    min, max, minmax, maxmin, fast_minmax, fast_maxmin,
    minmaxMag, maxmin, maxminMag,
    # ulpufp
    epsilon, eta, ulp, ufp, sulp, sufp,
    # QNaN
    qnan,
    # SignificantBits
    significand_bits, significant_bits

end # module


text:
using ReadCoverage

path_bam = "examples/data/RamDA_72h_A09.uniq.q40.bam"
path_bed12 = "examples/data/gencode.vM15.primary_assembly.annotation.protein_coding.head.bed"
out_prefix = "examples/out/RamDA_72h_A09.uniq.q40"

############
chrom = "chr19"
leftpos = 3205000
rightpos = 3207000
cov = readcoverage_bam(path_bam, chrom, leftpos, rightpos)
cov = readcoverage_bam(path_bam, chrom, leftpos, rightpos, output_prefix=out_prefix)
plot_read_coverage(cov)

################################
# Absolute genebody coverage
abcov = absolute_genebodycoverage(path_bam, path_bed12)
abcov = absolute_genebodycoverage(path_bam, path_bed12, output_prefix=out_prefix, bin_size=100)
plot_absolute_coverage(abcov)

################################
# Absolute genebody coverage
relcov = relative_genebodycoverage(path_bam, path_bed12);

max_depth = 0
relcov = relative_genebodycoverage(path_bam, path_bed12, output_prefix=out_prefix, max_depth=max_depth);
plot_relative_coverage(relcov)


text:
push!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

using ClimaCore.Geometry, LinearAlgebra, UnPack
import ClimaCore: Fields, Domains, Topologies, Meshes, Spaces
import ClimaCore: slab
import ClimaCore.Operators
import ClimaCore.Geometry
using LinearAlgebra, IntervalSets
using OrdinaryDiffEq: ODEProblem, solve, SSPRK33

using ClimaCore.RecursiveApply
using ClimaCore.RecursiveApply: rdiv, rmap

using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())

const parameters = (
    ϵ = 0.1,  # perturbation size for initial condition
    l = 0.5, # Gaussian width
    k = 0.5, # Sinusoidal wavenumber
    ρ₀ = 1.0, # reference density
    c = 2,
    g = 10,
)

numflux_name = get(ARGS, 1, "rusanov")
boundary_name = get(ARGS, 2, "")

domain = Domains.RectangleDomain(
    Geometry.XPoint(-2π)..Geometry.XPoint(2π),
    Geometry.YPoint(-2π)..Geometry.YPoint(2π),
    x1periodic = true,
    x2periodic = boundary_name != "noslip",
    x2boundary = boundary_name != "noslip" ? nothing : (:south, :north),
)

n1, n2 = 16, 16
Nq = 4
Nqh = 7
mesh = Meshes.EquispacedRectangleMesh(domain, n1, n2)
grid_topology = Topologies.GridTopology(mesh)
quad = Spaces.Quadratures.GLL{Nq}()
space = Spaces.SpectralElementSpace2D(grid_topology, quad)

Iquad = Spaces.Quadratures.GLL{Nqh}()
Ispace = Spaces.SpectralElementSpace2D(grid_topology, Iquad)

function init_state(coord, p)
    @unpack x, y = coord
    # set initial state
    ρ = p.ρ₀

    # set initial velocity
    U₁ = cosh(y)^(-2)

    # Ψ′ = exp(-(x2 + p.l / 10)^2 / 2p.l^2) * cos(p.k * x) * cos(p.k * y)
    # Vortical velocity fields (u₁′, u₂′) = (-∂²Ψ′, ∂¹Ψ′)
    gaussian = exp(-(y + p.l / 10)^2 / 2p.l^2)
    u₁′ = gaussian * (y + p.l / 10) / p.l^2 * cos(p.k * x) * cos(p.k * x)
    u₁′ += p.k * gaussian * cos(p.k * x) * sin(p.k * y)
    u₂′ = -p.k * gaussian * sin(p.k * x) * cos(p.k * y)

    u = Cartesian12Vector(U₁ + p.ϵ * u₁′, p.ϵ * u₂′)
    # set initial tracer
    θ = sin(p.k * y)

    return (ρ = ρ, ρu = ρ * u, ρθ = ρ * θ)
end

y0 = init_state.(Fields.coordinate_field(space), Ref(parameters))

function flux(state, p)
    @unpack ρ, ρu, ρθ = state
    u = ρu / ρ
    return (ρ = ρu, ρu = ((ρu ⊗ u) + (p.g * ρ^2 / 2) * I), ρθ = ρθ * u)
end

function energy(state, p)
    @unpack ρ, ρu = state
    u = ρu / ρ
    return ρ * (u.u1^2 + u.u2^2) / 2 + p.g * ρ^2 / 2
end

function total_energy(y, parameters)
    sum(state -> energy(state, parameters), y)
end

# numerical fluxes
wavespeed(y, parameters) = sqrt(parameters.g)

roe_average(ρ⁻, ρ⁺, var⁻, var⁺) =
    (sqrt(ρ⁻) * var⁻ + sqrt(ρ⁺) * var⁺) / (sqrt(ρ⁻) + sqrt(ρ⁺))

function roeflux(n, (y⁻, parameters⁻), (y⁺, parameters⁺))
    Favg = rdiv(flux(y⁻, parameters⁻) ⊞ flux(y⁺, parameters⁺), 2)

    λ = sqrt(parameters⁻.g)

    ρ⁻, ρu⁻, ρθ⁻ = y⁻.ρ, y⁻.ρu, y⁻.ρθ
    ρ⁺, ρu⁺, ρθ⁺ = y⁺.ρ, y⁺.ρu, y⁺.ρθ

    u⁻ = ρu⁻ / ρ⁻
    θ⁻ = ρθ⁻ / ρ⁻
    uₙ⁻ = u⁻' * n

    u⁺ = ρu⁺ / ρ⁺
    θ⁺ = ρθ⁺ / ρ⁺
    uₙ⁺ = u⁺' * n

    # in general thermodynamics, (pressure, soundspeed)
    p⁻ = (λ * ρ⁻)^2 * 0.5
    c⁻ = λ * sqrt(ρ⁻)

    p⁺ = (λ * ρ⁺)^2 * 0.5
    c⁺ = λ * sqrt(ρ⁺)

    # construct roe averges
    ρ = sqrt(ρ⁻ * ρ⁺)
    u = roe_average(ρ⁻, ρ⁺, u⁻, u⁺)
    θ = roe_average(ρ⁻, ρ⁺, θ⁻, θ⁺)
    c = roe_average(ρ⁻, ρ⁺, c⁻, c⁺)

    # construct normal velocity
    uₙ = u' * n

    # differences
    Δρ = ρ⁺ - ρ⁻
    Δp = p⁺ - p⁻
    Δu = u⁺ - u⁻
    Δρθ = ρθ⁺ - ρθ⁻
    Δuₙ = Δu' * n

    # constructed values
    c⁻² = 1 / c^2
    w1 = abs(uₙ - c) * (Δp - ρ * c * Δuₙ) * 0.5 * c⁻²
    w2 = abs(uₙ + c) * (Δp + ρ * c * Δuₙ) * 0.5 * c⁻²
    w3 = abs(uₙ) * (Δρ - Δp * c⁻²)
    w4 = abs(uₙ) * ρ
    w5 = abs(uₙ) * (Δρθ - θ * Δp * c⁻²)

    # fluxes!!!

    fluxᵀn_ρ = (w1 + w2 + w3) * 0.5
    fluxᵀn_ρu =
        (w1 * (u - c * n) + w2 * (u + c * n) + w3 * u + w4 * (Δu - Δuₙ * n)) *
        0.5
    fluxᵀn_ρθ = ((w1 + w2) * θ + w5) * 0.5

    Δf = (ρ = -fluxᵀn_ρ, ρu = -fluxᵀn_ρu, ρθ = -fluxᵀn_ρθ)
    rmap(f -> f' * n, Favg) ⊞ Δf
end

numflux = if numflux_name == "central"
    Operators.CentralNumericalFlux(flux)
elseif numflux_name == "rusanov"
    Operators.RusanovNumericalFlux(flux, wavespeed)
elseif numflux_name == "roe"
    roeflux
end

function rhs!(dydt, y, (parameters, numflux), t)

    # ϕ' K' W J K dydt =  -ϕ' K' I' [DH' WH JH flux.(I K y)]
    #  =>   K dydt = - K inv(K' WJ K) K' I' [DH' WH JH flux.(I K y)]

    # where:
    #  ϕ = test function
    #  K = DSS scatter (i.e. duplicates points at element boundaries)
    #  K y = stored input vector (with duplicated values)
    #  I = interpolation to higher-order space
    #  D = derivative operator
    #  H = suffix for higher-order space operations
    #  W = Quadrature weights
    #  J = Jacobian determinant of the transformation `ξ` to `x`
    #
    Nh = Topologies.nlocalelems(y)

    F = flux.(y, Ref(parameters))
    dydt .= Operators.slab_weak_divergence(F)

    Operators.add_numerical_flux_internal!(numflux, dydt, y, parameters)

    Operators.add_numerical_flux_boundary!(
        dydt,
        y,
        parameters,
    ) do normal, (y⁻, parameters)
        y⁺ = (ρ = y⁻.ρ, ρu = y⁻.ρu - dot(y⁻.ρu, normal) * normal, ρθ = y⁻.ρθ)
        numflux(normal, (y⁻, parameters), (y⁺, parameters))
    end

    # 6. Solve for final result
    dydt_data = Fields.field_values(dydt)
    dydt_data .= rdiv.(dydt_data, space.local_geometry.WJ)

    M = Spaces.Quadratures.cutoff_filter_matrix(
        Float64,
        space.quadrature_style,
        3,
    )
    Operators.tensor_product!(dydt_data, M)

    return dydt
end

dydt = Fields.Field(similar(Fields.field_values(y0)), space)
rhs!(dydt, y0, (parameters, numflux), 0.0);

# Solve the ODE operator
prob = ODEProblem(rhs!, y0, (0.0, 200.0), (parameters, numflux))
sol = solve(
    prob,
    SSPRK33(),
    dt = 0.02,
    saveat = 1.0,
    progress = true,
    progress_message = (dt, u, p, t) -> t,
)

ENV["GKSwstype"] = "nul"
import Plots
Plots.GRBackend()

dirname = "dg_$(numflux_name)"
if boundary_name != ""
    dirname = "$(dirname)_$(boundary_name)"
end
path = joinpath(@__DIR__, "output", dirname)
mkpath(path)

anim = Plots.@animate for u in sol.u
    Plots.plot(u.ρθ, clim = (-1, 1))
end
Plots.mp4(anim, joinpath(path, "tracer.mp4"), fps = 10)

Es = [total_energy(u, parameters) for u in sol.u]
Plots.png(Plots.plot(Es), joinpath(path, "energy.png"))

function linkfig(figpath, alt = "")
    # buildkite-agent upload figpath
    # link figure in logs if we are running on CI
    if get(ENV, "BUILDKITE", "") == "true"
        artifact_url = "artifact://$figpath"
        print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n")
    end
end

linkfig("output/$(dirname)/energy.png", "Total Energy")


text:
using DataStructures

function solve()
    s = readline()
    l = length(s)
    if l ≤ 2
        n = parse(Int, s)
        r = parse(Int, reverse(s))
        n % 8 == 0 || r % 8 == 0
    else
        r = false
        c = counter(s)
        for i=112:8:992
            rr = true
            j = counter(string(i))
            for (d, v) in j
                if c[d] < v
                    rr = false
                end
            end
            if rr
                r = true
                break
            end
        end
        r
    end
end

println(solve() ? "Yes" : "No")


text:
### A Pluto.jl notebook ###
# v0.17.1

using Markdown
using InteractiveUtils

# ╔═╡ d3b9397f-5673-47ca-8a01-62fd985f121f
begin

	# load packages -
	using BSON
	using PlutoUI
	using DataFrames
	using GLPK
	using PrettyTables

	# setup paths -
	_PATH_TO_ROOT = pwd()
	_PATH_TO_SRC = joinpath(_PATH_TO_ROOT,"src")
	_PATH_TO_MODEL = joinpath(_PATH_TO_ROOT, "model")

	# include my code lib -
	include(joinpath(_PATH_TO_SRC, "Include.jl"))

	# show -
	nothing
end

# ╔═╡ 7bf3ac90-637d-4992-8835-6798d46f71f5
md"""
### Example: Using Flux Balance Analysis to Compute Fractional Conversion

We want to produce a product _P_ (desired product) by converting feedstock A₁ using a cell-free biochemical process operating in a well-mixed continuous microfluidic chip with a single input and a single output, and a liquid reaction volume of V = 100 μL. The reaction network is shown in Fig 1.

__Assumptions__
* Microfluidic chip is well-mixed and operates at steady-state
* Constant T, P on the chip
* Liquid phase is ideal

__Compute__
* Use flux balance analysis to compute the optimal open extent of reaction $\dot{\epsilon}_{i}$, where the objective function is to maximize _P_
* Compute the state table for the mol flow rates into (you decide) and from the chip (requires FBA)
* Compute the fractional conversion of the feedstock A₁

"""

# ╔═╡ 7aa00ab3-b105-415a-a7cf-111afc1829d7
md"""
__Fig. 1__ Schematic of the reaction network operating on the steady-state chip."
"""

# ╔═╡ f2316c41-0c45-4dd2-8749-3d67256b08b4
PlutoUI.LocalResource("./figs/Fig-FBA-ToyNetwork.png")

# ╔═╡ 6a0c8767-6de3-44b2-9d84-04e504e43c1b
md"""
##### Compute the optimal open reaction extents using Flux Balance Analysis (FBA)

The FBA problem is (traditinally) a Linear Programming (LP) problem in which a linear objective function $\mathcal{O}$:

$$\mathcal{O} = \sum_{j=1}^{\mathcal{R}}c_{j}\dot{\epsilon}_{j}$$

is minimized (or maximized) subject to a collection of linear constraints, and bound constraints on the permissible values of the extents $\dot{\epsilon}_{i}$. In this case, we know that the steady-state mol flow rate for species $i$ leaving the chip given (where $s=1$ denotes the inlet and $s=2$ denotes the outlet):

$$\dot{n}_{i,2} = \dot{n}_{i,1} + \sum_{j=1}^{\mathcal{R}}\sigma_{ij}\dot{\epsilon}_{j}$$

Thus, because $\dot{n}_{i,2}\geq{0}$, the FBA problem is subject to the mol constraints:

$$\dot{n}_{i,1} + \sum_{j=1}^{\mathcal{R}}\sigma_{ij}\dot{\epsilon}_{j}\geq{0}\qquad\forall{i}$$

In other words, when searching for the optimal set of $\dot{\epsilon}_{j}$ we have to select values that give physically realistic answers (we can't have a negative mol flow rate). Next, the $\dot{\epsilon}_{j}$ terms are bounded from above and below:

$$\mathcal{L}_{j}\leq\dot{\epsilon}_{j}\leq\mathcal{U}_{j}\qquad{j=1,2\dots,\mathcal{R}}$$

where the $\mathcal{L}_{j}$ and $\mathcal{U}_{j}$ denote the lower and upper bounds that $\dot{\epsilon}_{j}$ can take. Remember that the open extents are just reaction rates times the volume. Thus, the lower and upper bounds describe the permissible range that we expect the rate _could_ obtain.

"""

# ╔═╡ 78eb6166-ac8a-4c28-8f8f-ae00b5c78b3f
begin

	# setup the volume -
	V = 100.0*(1.0/1.0e6)

	# setup flow rate in
	n_dot_in = [
		10.0 	; # 1 A₁
		3.0 	; # 2 A₂
		0.0 	; # 3 B
		0.0 	; # 4 P
		1.0 	; # 5 C
		0.0 	; # 6 x
		0.0 	; # 7 y
	]

	# Setup STM -
	S = [

		# r₁ r₂ r₃
		-1.0 0.0 0.0 ; # 1 A₁
		0.0 0.0 -1.0 ; # 2 A₂
		1.0 -1.0 0.0 ; # 3 B
		0.0 1.0 0.0  ; # 4 P
		0.0 0.0 1.0  ; # 5 C
		-1.0 0.0 1.0 ; # 6 x
		1.0 0.0 -1.0 ; # 7 y
	];
	(number_of_states, number_of_reactions) = size(S)

	# set the flux bounds array -
	flux_bounds_array = [

		# ℒ 𝒰
		0.0 10.0 	; # 1 r₁
		0.0 10.0  	; # 2 r₂
		0.0 20.0 	; # 3 r₃
	];

	# set the species bounds array -
	species_bounds_array = [

		# ℒ lower     𝒰 upper
		-n_dot_in[1] 1000.0 				; # 1 A₁
		-n_dot_in[2] 1000.0 				; # 2 A₂
		-n_dot_in[3] 1000.0 				; # 3 B
		-n_dot_in[4] 1000.0 				; # 4 P
		-n_dot_in[5] 1000.0 				; # 5 C
		-n_dot_in[6] 1000.0 				; # 6 x
		-n_dot_in[7] 1000.0 				; # 7 y
	];

	# max P formation -
	obj_vector = zeros(3)
	obj_vector[2] = -1.0

	# show -
	nothing
end

# ╔═╡ f3efe4fa-c78d-42ab-bf9e-ecd7011874d9
begin

	# compute the optimal flux -
	result = calculate_optimal_flux_distribution(S, flux_bounds_array, species_bounds_array, obj_vector)

	# initialize some storage -
	flux_table = Array{Any,2}(undef,number_of_reactions,4)

	# reactions -
	reaction_strings = [
		"A₁ + x => B + y" 	; # 1 r₁
		"B => P" 			; # 2 r₂
		"A₂ + y => C + x" 	; # 3 r₃
	]

	# populate the state table -
	for reaction_index = 1:number_of_reactions
		flux_table[reaction_index,1] = reaction_strings[reaction_index]
		flux_table[reaction_index,2] = result.calculated_flux_array[reaction_index]
		flux_table[reaction_index,3] = flux_bounds_array[reaction_index,1]
		flux_table[reaction_index,4] = flux_bounds_array[reaction_index,2]
	end

	with_terminal() do

		# header row -
		flux_table_header_row = (["Reaction","ϵᵢ_dot", "ϵ₁_dot LB", "ϵ₁_dot UB"],["","mol/time", "mol/time", "mol/time"]);

		# write the table -
		pretty_table(flux_table; header=flux_table_header_row)
	end

end

# ╔═╡ 45740c44-a846-4b7d-acc1-9105aa67e786
begin

	# compute output -
	ϵ_vector = result.calculated_flux_array
	n_dot_out = n_dot_in .+ S*ϵ_vector

	# make a pretty table -
	state_table = Array{Any,2}(undef, number_of_states,3)
	species_array = ["A₁", "A₂", "B", "P", "C", "x","y"]
	for state_index = 1:number_of_states

		# populate -
		state_table[state_index,1] = species_array[state_index]
		state_table[state_index,2] = n_dot_in[state_index]
		state_table[state_index,3] = n_dot_out[state_index]
	end

	with_terminal() do

		# header row -
		state_table_header_row = (["Species","nᵢ_dot_in","nᵢ_dot_out"],["","mol/time", "mol/time"]);

		# write the table -
		pretty_table(state_table; header=state_table_header_row)
	end
end

# ╔═╡ 05e4f243-4a1b-4962-a36a-7a931b8e7936
md"""

##### Compute the fractional conversion of A₁ and A₂
The fractional conversion of a _reactant_ species for an _open_ system (with a single input and a single output) is given by:

$$f_{i} = \frac{\dot{n}_{i,j-1} - \dot{n}_{i,j}}{\dot{n}_{i,j-1}}$$

where $\dot{n}_{i,j-1}$ denotes the mol flow rate of species $i$ in stream $j-1$ (inlet), and $\dot{n}_{i,j}$ denotes the mol flow rate of species $i$ in stream $j$ (outlet).

"""

# ╔═╡ 85db6161-58f6-4fea-ad00-7eb5aae85fd7
with_terminal() do

	# compute fᵢ
	f₁ = (n_dot_in[1] - n_dot_out[1])/(n_dot_in[1])
	f₂ = (n_dot_in[2] - n_dot_out[2])/(n_dot_in[2])

	# display -
	println("Fraction conversion A₁ = $(f₁) and A₂ = $(f₂)")

end

# ╔═╡ 6c747a6e-4bf7-11ec-2b43-af47c5975080
html"""
<style>
main {
    max-width: 1200px;
    width: 85%;
    margin: auto;
    font-family: "Roboto, monospace";
}

a {
    color: blue;
    text-decoration: none;
}

.H1 {
    padding: 0px 30px;
}
</style>"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"

[compat]
BSON = "~0.3.4"
DataFrames = "~1.2.2"
GLPK = "~0.15.1"
PlutoUI = "~0.7.19"
PrettyTables = "~1.2.3"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "0bc60e3006ad95b4bb7497698dd7c6d649b9bc06"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.1"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[BSON]]
git-tree-sha1 = "ebcd6e22d69f21249b7b8668351ebf42d6dc87a1"
uuid = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
version = "0.3.4"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[BenchmarkTools]]
deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
git-tree-sha1 = "61adeb0823084487000600ef8b1c00cc2474cd47"
uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
version = "1.2.0"

[[BinaryProvider]]
deps = ["Libdl", "Logging", "SHA"]
git-tree-sha1 = "ecdec412a9abc8db54c0efc5548c64dfce072058"
uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
version = "0.5.10"

[[Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.8+0"

[[CEnum]]
git-tree-sha1 = "215a9aa4a1f23fbd05b92769fdd62559488d70e9"
uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
version = "0.4.1"

[[CodecBzip2]]
deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"]
git-tree-sha1 = "2e62a725210ce3c3c2e1a3080190e7ca491f18d7"
uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd"
version = "0.7.2"

[[CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "dce3e3fea680869eaa0b774b2e8343e9ff442313"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.40.0"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "d785f42445b63fc86caa08bb9a9351008be9b765"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.2.2"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[GLPK]]
deps = ["BinaryProvider", "CEnum", "GLPK_jll", "Libdl", "MathOptInterface"]
git-tree-sha1 = "8c9b0ce6d476800c7a7e8cc0db5a60e241cf107a"
uuid = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
version = "0.15.1"

[[GLPK_jll]]
deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "01de09b070d4b8e3e1250c6542e16ed5cad45321"
uuid = "e8aa6df9-e6ca-548a-97ff-1f85fc5b8b98"
version = "5.0.0+0"

[[GMP_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d"

[[Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MathOptInterface]]
deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "JSON", "LinearAlgebra", "MutableArithmetics", "OrderedCollections", "Printf", "SparseArrays", "Test", "Unicode"]
git-tree-sha1 = "afa62f733d78f63c2292730994dcb076835cf1d2"
uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
version = "0.10.5"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[MutableArithmetics]]
deps = ["LinearAlgebra", "SparseArrays", "Test"]
git-tree-sha1 = "7bb6853d9afec54019c1397c6eb610b9b9a19525"
uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
version = "0.3.1"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "ae4bbcadb2906ccc085cf52ac286dc1377dceccc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.2"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "e071adf21e165ea0d904b595544a8e514c8bb42c"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.19"

[[PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "db3a23166af8aebf4db5ef87ac5b00d36eb771e2"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.4.0"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "d940010be611ee9d67064fe559edbb305f8cc0eb"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.2.3"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Profile]]
deps = ["Printf"]
uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "fed34d0e71b91734bf0a7e10eb1bb05296ddbcd0"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─7bf3ac90-637d-4992-8835-6798d46f71f5
# ╟─7aa00ab3-b105-415a-a7cf-111afc1829d7
# ╟─f2316c41-0c45-4dd2-8749-3d67256b08b4
# ╟─6a0c8767-6de3-44b2-9d84-04e504e43c1b
# ╠═78eb6166-ac8a-4c28-8f8f-ae00b5c78b3f
# ╟─f3efe4fa-c78d-42ab-bf9e-ecd7011874d9
# ╟─45740c44-a846-4b7d-acc1-9105aa67e786
# ╟─05e4f243-4a1b-4962-a36a-7a931b8e7936
# ╟─85db6161-58f6-4fea-ad00-7eb5aae85fd7
# ╟─d3b9397f-5673-47ca-8a01-62fd985f121f
# ╟─6c747a6e-4bf7-11ec-2b43-af47c5975080
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using HypothesisTests
import HypothesisTests: HypothesisTest, pvalue

"""
Goodness of Fit Test.
"""
abstract type KernelGoodnessOfFitTest <: HypothesisTest end

# implement this for whatever gof-test
"""
Performs the specified goodness-of-fit test.

# Arguments
- `t::KSDTest|FSSDTest`: the parameters for the test.
- `q`: model which needs to have an implementation for `gradlogpdf`.
- `x::AbstractArray`: samples to test against.
"""
function perform(::KernelGoodnessOfFitTest, q, xs) end

# HypothesisTests.jl interface
function pvalue(::KernelGoodnessOfFitTest) end

abstract type KernelGoodnessOfFitResult end


text:
using BinDeps

arch_name = "msys2-x86_64-latest.tar.xz"
depsdir = Pkg.dir("MSys2", "deps")
downloadsdir = joinpath(depsdir, "download")
mkpath(downloadsdir)
download(joinpath("http://repo.msys2.org/distrib", arch_name),
          joinpath(downloadsdir, arch_name))
(name, ext, ext2) = BinDeps.splittarpath(arch_name)
run(unpack_cmd(joinpath(downloadsdir, arch_name), depsdir, ext, ext2))
using MSys2
MSys2.update()
MSys2.install("pacman")
MSys2.upgrade()
MSys2.upgrade()


text:
#!/usr/bin/env julia

module PE12

trig(n::Int)=n*(n+1)÷2

function divisorlength(n::Int)
    local s=0
    local ed=floor(√n)
    for i in 1:ed
        if n%i==0
            if i≠ed
                s+=1
            end
            s+=1
        end
    end
    s
end

end

export pe12

function pe12()
    local n=1
    while PE12.divisorlength(PE12.trig(n))≤500
        n+=1
    end
    PE12.trig(n)
end

if !haskey(ENV,"PROJECT_EULER_WITHOUT_RESULT")
    pe12()|>print
end

max_stars_repo_path: 12. Highly divisible triangular number/12.jl

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_name: miRoox/ProjectEuler


text:
macro search_and_replace(ex, search, replace)

    ex = Expr(:quote, ex)
    search = Expr(:quote, search)
    replace = Expr(:quote, replace)
    esc(quote
        subst(x) = MLStyle.@match x begin
            $search => $replace
            _ => x
        end
        MacroTools.prewalk(subst, $ex)
    end)
end

# I would need to use DocumentFormat.jl style code to do this
function search_and_replace(ex::String, search::String, replace::String, eval_fun = ex -> @eval Main $ex)

end


text:
module Watchtower
using ..Ahorn, Maple

@mapdef Entity "CommunalHelper/NoOverlayLookout" NoOverlayLookout(
    x::Integer,
    y::Integer,
    summit::Bool=false,
    onlyY::Bool=false,
)

const placements = Ahorn.PlacementDict(
    "Watchtower (No Overlay) (Communal Helper)" => Ahorn.EntityPlacement(
        NoOverlayLookout,
    ),
)

Ahorn.nodeLimits(entity::NoOverlayLookout) = 0, -1

sprite = "objects/lookout/lookout05.png"

function Ahorn.selection(entity::NoOverlayLookout)
    nodes = get(entity.data, "nodes", ())
    x, y = Ahorn.position(entity)

    res = Ahorn.Rectangle[Ahorn.getSpriteRectangle(sprite, x, y, jx=0.5, jy=1.0)]

    for node in nodes
        nx, ny = Int.(node)

        push!(res, Ahorn.getSpriteRectangle(sprite, nx, ny, jx=0.5, jy=1.0))
    end

    return res
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::NoOverlayLookout)
    px, py = Ahorn.position(entity)

    for node in get(entity.data, "nodes", [])
        nx, ny = Int.(node)

        Ahorn.drawArrow(ctx, px, py - 8, nx, ny - 8, Ahorn.colors.selection_selected_fc, headLength=6)
        Ahorn.drawSprite(ctx, sprite, nx, ny, jx=0.5, jy=1.0)

        px, py = nx, ny
    end
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::NoOverlayLookout, room::Maple.Room)
    x, y = Ahorn.position(entity)

    Ahorn.drawSprite(ctx, sprite, x, y, jx=0.5, jy=1.0)
end

end


text:
export AbstractHLAData

abstract type AbstractHLAData end