text:
module Checkpoint

using ..Ahorn, Maple

function getBgSprite(entity::Maple.ChapterCheckpoint)
    bg = get(entity, "bg", "")
    name = "objects/checkpoint/bg/$bg"

    if !isempty(bg)
        backgroundSprite = Ahorn.getSprite(name, "Gameplay")

        if backgroundSprite.width != 0 && backgroundSprite.height != 0
            return true, name
        end
    end

    return false, false
end

bgFallback = "objects/checkpoint/bg/1"
flashSprite = "objects/checkpoint/flash03"

function Ahorn.selection(entity::Maple.ChapterCheckpoint)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(bgFallback, x, y, jx=0.5, jy=1.0)
end

Ahorn.editingOptions(entity::Maple.ChapterCheckpoint) = Dict{String, Any}(
    "dreaming" => Dict{String, Any}(
        "Automatic" => nothing,
        "Dreaming" => true,
        "Awake" => false
    ),
    "inventory" => merge(
        Dict{String, Any}(
            "Automatic" => nothing
        ),
        Dict{String, Any}(
            inventory => inventory for inventory in Maple.inventories
        )
    ),
    "coreMode" => merge(
        Dict{String, Any}(
            "Automatic" => nothing
        ),
        Dict{String, Any}(
            mode => mode for mode in Maple.core_modes
        )
    ),
)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.ChapterCheckpoint, room::Maple.Room)
    exists, bgSprite = getBgSprite(entity)

    if exists
        Ahorn.drawSprite(ctx, bgSprite, 0, 0, jx=0.5, jy=1.0)
    end

    Ahorn.drawSprite(ctx, flashSprite, 0, 0, jx=0.5, jy=1.0)
end

end


text:
# ---
# title: 479. Largest Palindrome Product
# id: problem479
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories:
# link: <https://leetcode.com/problems/largest-palindrome-product/description/>
# hidden: true
# ---
#
# Find the largest palindrome made from the product of two n-digit numbers.
#
# Since the result could be very large, you should return the largest palindrome
# mod 1337.
#
#
#
# **Example:**
#
# Input: 2
#
# Output: 987
#
# Explanation: 99 x 91 = 9009, 9009 % 1337 = 987
#
#
#
# **Note:**
#
# The range of n is [1,8].
#
#
## @lc code=start
using LeetCode

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

ext: jl


text:
function rotate(key::Int, msg::String)
    return map(c -> rotate(key, c), msg)
end

function rotate(key::Int, c::Char)
    # https://en.wikipedia.org/wiki/Caesar_cipher
    # https://rosettacode.org/wiki/Caesar_cipher#Julia
    if c in 'a':'z'
        shft = 'a'
        # https://docs.julialang.org/en/v1/manual/strings/#man-characters
        shft + (c - shft + key) % 26
    elseif c in 'A':'Z'
        shft = 'A'
        shft + (c - shft + key) % 26
    else
        c
    end
end

# Source: https://exercism.org/tracks/julia/exercises/rotational-cipher/solutions/2c423955e21c45d28ad20d04af99cbbb
# https://docs.julialang.org/en/v1/manual/metaprogramming/#meta-non-standard-string-literals
macro R13_str(msg)
    return rotate(13, msg)
end

for i in 0:26
    # https://docs.julialang.org/en/v1/manual/metaprogramming/#Code-Generation
    # https://docs.julialang.org/en/v1/manual/metaprogramming/#Symbols
    # https://docs.julialang.org/en/v1/manual/metaprogramming/#man-expression-interpolation
    @eval macro $(Symbol('R', i, "_str"))(msg)
        return rotate($i, msg)
    end
end

max_stars_repo_name: joaopalmeiro/learning-julia

max_issues_repo_name: joaopalmeiro/learning-julia

max_forks_repo_name: joaopalmeiro/learning-julia


text:
import JSON: json

js_runtime() =
    readstring(open(joinpath(dirname(@__FILE__), "..", "runtime", "build.js")))

load_js_runtime() =
    display(MIME("text/html"), "")

pwid() = replace(string(gensym("pwid")), "#", "")

@compat show(io::IO, ::MIME"application/json", x::(@compat Union{Node, Patch})) =
    write(io, json(jsonfmt(x)))

@compat function show(io::IO, ::MIME"text/html", x::Node)
    id = pwid()
    write(io, """
""",
              """""",
              """
""")
end

function refdiff(a, b; label="")
    # Inspect a reference (vtree) diff of two nodes
    display(MIME("text/html"), string(""))
end


text:
file_dir = Pkg.dir("D4M")*"/examples/2Apps/1EntityAnalysis/Entity.jld";

using JLD

E = load(file_dir)["E"]

sum(logical(col2type(E,"/")),1);

En = sum(logical(E),1);

~,entity,count = find(En);

An = Assoc(count,entity,1);

using PyPlot
loglog(full(sum(Adj(An[:,StartsWith("LOCATION/,")]),2)) ,"o");


text:
export VoxelwiseConstant
export bounds

struct VoxelwiseConstant{K,AN<:AbsArrNumber{K}}
    val::AN  # scalar values in voxels
    lvxl::NTuple{K,VecFloat}  # locations of faces (not centers) of voxels
    ∆lvxl::NTuple{K,VecFloat}  # voxel edges

    function VoxelwiseConstant{K,AN}(val, lvxl) where {K,AN}
        size(val) .+ 1 == length.(lvxl) ||
            @error "length.(lvxl) = $(length.(lvxl)) should be greater by 1 in each dimension than size(val) = $(size(val))."

        for k = 1:K
            issorted(lvxl[k]) || @error "lvxl[$k] = $(lvxl[k]) should be sorted."
        end

        ∆lvxl = diff.(lvxl)
        @assert length.(∆lvxl)==size(val)

        return new(val, lvxl, ∆lvxl)
     end
end

VoxelwiseConstant(val::AN, lvxl::NTuple{K,AbsVecReal}) where {K,AN<:AbsArrNumber{K}} =
    VoxelwiseConstant{K,AN}(val, lvxl)

Base.size(vc::VoxelwiseConstant) = size(vc.val)
bounds(vc::VoxelwiseConstant) = (SVec(minimum.(vc.lvxl)), SVec(maximum.(vc.lvxl)))

(vc::VoxelwiseConstant{K})(x::AbsVecReal) where {K} = vc(SFloat{K}(x))

function (vc::VoxelwiseConstant{K,AN})(x::SReal{K}) where {K,AN<:AbsArrNumber{K}}
    lvxl = vc.lvxl
    for k = 1:K
        lvxl[k][1] ≤ x[k] ≤ lvxl[k][end] || @error "x[$k] = $(x[k]) should be between lvxl[$k][1] = $(lvxl[k][1]) and lvxl[k][end] = $(lvxl[k][end])."
    end

    # Find the x-containing voxel with an inclusive negative-end boundary and exclusive
    # positive-end boundary.
    ind_rng = MVec{K,UnitRange{Int64}}(undef)
    for k = 1:K
        xₖ = x[k]
        lvxlₖ = lvxl[k]

        indₑ = findlast(w -> w≤xₖ, @view(lvxlₖ[1:end-1]))

        at_bound = lvxlₖ[indₑ]==xₖ && indₑ≥2
        indₛ = indₑ - at_bound

        ind_rng[k] = indₛ:indₑ
    end
    CI = CartesianIndices(ind_rng.data)  # CartesianIndices{K}

    return mean(@view(vc.val[CI]))
end


text:
## Linear Solve

function linsolve(A::Operator,b;kwds...)
    if dimension(domain(A)) > 1
        pdesolve(A,b;kwds...)
    elseif isambiguous(domainspace(A))
        A=choosespaces(A,b)
        if isambiguous(domainspace(A))
            error("Cannot infer spaces")
        end
        linsolve(A,b;kwds...)
    else
        linsolve(qrfact(A),b;kwds...)
    end
end

function linsolve{OO<:Operator}(A::Array{OO},b;kwds...)
    if dimension(domain(A[end])) > 1
        pdesolve(A,b;kwds...)
    else
        linsolve(interlace(A),b;kwds...)
    end
end

for p in (1,2)
    @eval begin
        \{T<:Operator}(A::Array{T,$p},b::Union{Array,Number,Fun}) = linsolve(A,b)
        \{T<:Operator,F<:Fun}(A::Array{T,$p},b::Vector{F}) = linsolve(A,b)
        \{T<:Operator,F<:Number}(A::Array{T,$p},b::Vector{F}) = linsolve(A,b)
    end
end
\(A::Operator,b) = linsolve(A,b)


text:
# Pharmacokinetics
# Makoid C, Vuchetich J, Banakar V. 1996-1999. Basic Pharmacokinetics.

function firstobs(time::Vector{T}, obs::Vector, dosetime) where T <: Number
    @inbounds for i = 1:length(time)
        if time[i] >= dosetime && !isnanormissing(obs[i]) return i end
    end
    error("Observations not found")
end
function firstobs(time::Vector{T}, obs, vol, dosetime) where T <: Tuple
    @inbounds for i = 1:length(time)
        if time[i][1] >= dosetime && !isnanormissing(obs[i]) && !isnanormissing(vol[i]) return i end
    end
    error("Observations not found")
end

function ctaumin(time::AbstractVector, obs::AbstractVector, taulastp::Int)
    min = first(obs)
    if length(obs) == 1 return min end
    @inbounds for i = 2:taulastp
        if  obs[i] < min  min = obs[i] end
    end
    min
end

function ctmax(time::AbstractVector, obs::AbstractVector{T}, taulastp) where T
    cmax  = obs[1]
    tmaxn = 1
    if length(obs) == 1 return cmax, first(time), tmaxn end
    taulastp > length(obs) && error("taulastp > length(obs")
    @inbounds for i = 2:taulastp
        if obs[i] > cmax
            cmax  = obs[i]
            tmaxn = i
        end
    end
    return cmax, time[tmaxn], tmaxn
end
function ctmax(time::AbstractVector, obs::AbstractVector{T}) where T
    f = findfirst(!isnanormissing, obs)
    cmax  = obs[f]
    tmaxn = f
    if length(obs) - f == 0 return cmax, time[f], tmaxn end
    @inbounds for i = f+1:length(obs)
        if !isnanormissing(obs[i]) && obs[i] > cmax
            cmax  = obs[i]
            tmaxn = i
        end
    end
    return cmax, time[tmaxn], tmaxn
end
function ctmax(data::PKSubject)
    fobs = firstobs(data.time, data.obs, data.dosetime.time)
    if  data.dosetime.tau > 0
        taulastp = findlast(x -> x <= data.dosetime.time + data.dosetime.tau, data.time)
    else
        taulastp = length(data.obs)
    end
    cmax  = data.obs[fobs]
    tmaxn = fobs
    if length(data.obs) - fobs == 0 return cmax, data.time[fobs], tmaxn end
    @inbounds for i = fobs + 1:length(data.obs)
        if !isnanormissing(data.obs[i]) && data.obs[i] > cmax
            cmax  = data.obs[i]
            tmaxn = i
        end
    end
    return cmax, data.time[tmaxn], tmaxn
end

function logcpredict(t₁, t₂, tx, c₁, c₂)
    return exp(log(c₁) + (tx-t₁)/(t₂-t₁)*(log(c₂) - log(c₁)))
end

#Linear trapezoidal auc
function linauc(t₁, t₂, c₁, c₂)
    return (t₂-t₁)*(c₁+c₂)/2
end
#Linear trapezoidal aumc
function linaumc(t₁, t₂, c₁, c₂)
    return (t₂-t₁)*(t₁*c₁+t₂*c₂)/2
end
#Log trapezoidal auc
function logauc(t₁, t₂, c₁, c₂)
    return  (t₂-t₁)*(c₂-c₁)/log(c₂/c₁)
end
#Log trapezoidal aumc
function logaumc(t₁, t₂, c₁, c₂)
    return (t₂-t₁) * (t₂*c₂-t₁*c₁) / log(c₂/c₁) - (t₂-t₁)^2 * (c₂-c₁) / log(c₂/c₁)^2
end
#Intrapolation
#linear prediction bx from ax, a1 < ax < a2
function linpredict(a₁, a₂, ax, b₁, b₂)
    return (ax - a₁) / (a₂ - a₁)*(b₂ - b₁) + b₁
end

function slope(x, y)
    if length(x) != length(y) throw(ArgumentError("Unequal vector length!")) end
    n   = length(x)
    if n < 2 throw(ArgumentError("n < 2!")) end
    Σxy::Float64 = zero(Float64)
    Σx::Float64  = zero(Float64)
    Σy::Float64  = zero(Float64)
    Σx2::Float64 = zero(Float64)
    Σy2::Float64 = zero(Float64)
    @inbounds for i = 1:n
        xi = x[i]
        yi = y[i]
        Σxy += xi  * yi
        Σx  += xi
        Σy  += yi
        Σx2 += xi^2
        Σy2 += yi^2
    end
    a   = (n * Σxy - Σx * Σy)/(n * Σx2 - Σx^2)
    b   = (Σy * Σx2 - Σx * Σxy)/(n * Σx2 - Σx^2)
    r2  = (n * Σxy - Σx * Σy)^2/((n * Σx2 - Σx^2)*(n * Σy2 - Σy^2))

    n > 2 ? ar  = 1 - (1 - r2)*(n - 1)/(n - 2) : ar = NaN

    return a, b, r2, ar
end
function logslope(x, y)
    if length(x) != length(y) throw(ArgumentError("Unequal vector length!")) end
    n   = length(x)
    if n < 2 throw(ArgumentError("n < 2!")) end
    Σxy::Float64 = zero(Float64)
    Σx::Float64  = zero(Float64)
    Σy::Float64  = zero(Float64)
    Σx2::Float64 = zero(Float64)
    Σy2::Float64 = zero(Float64)
    @inbounds for i = 1:n
        xi = x[i]
        yi = log(y[i])
        Σxy += xi * yi
        Σx  += xi
        Σy  += yi
        Σx2 += xi^2
        Σy2 += yi^2
    end
    a   = (n * Σxy - Σx * Σy)/(n * Σx2 - Σx^2)
    b   = (Σy * Σx2 - Σx * Σxy)/(n * Σx2 - Σx^2)
    r2  = (n * Σxy - Σx * Σy)^2/((n * Σx2 - Σx^2)*(n * Σy2 - Σy^2))

    n > 2 ? ar  = 1 - (1 - r2)*(n - 1)/(n - 2) : ar = NaN

    return a, b, r2, ar
end #end slope

#---------------------------------------------------------------------------
function aucpart(t₁, t₂, c₁::T, c₂::T, calcm, aftertmax) where T
    if calcm == :lint || c₁ <= zero(T) && c₂ <= zero(T)
        auc   =  linauc(t₁, t₂, c₁, c₂)
    elseif calcm == :logt && aftertmax && c₁ > zero(T) && c₂ > zero(T)
        auc   =  logauc(t₁, t₂, c₁, c₂)
    elseif calcm == :luld &&  c₁ > c₂ > zero(T)
        auc   =  logauc(t₁, t₂, c₁, c₂)
    elseif calcm == :luldt && aftertmax && c₁ > c₂ > zero(T)
        auc   =  logauc(t₁, t₂, c₁, c₂)
    #elseif calcm == :log && c₁ > zero(T) && c₂ > zero(T)
        #auc   =  logauc(t₁, t₂, c₁, c₂)
    else
        auc   =  linauc(t₁, t₂, c₁, c₂)
    end
    return auc
end
function aumcpart(t₁, t₂, c₁::T, c₂::T, calcm, aftertmax) where T
    if calcm == :lint || c₁ <= zero(T) && c₂ <= zero(T)
        aumc  = linaumc(t₁, t₂, c₁, c₂)
    elseif calcm == :logt && aftertmax && c₁ > zero(T) && c₂ > zero(T)
        aumc  = logaumc(t₁, t₂, c₁, c₂)
    elseif calcm == :luld &&  c₁ > c₂ > zero(T)
        aumc  = logaumc(t₁, t₂, c₁, c₂)
    elseif calcm == :luldt && aftertmax && c₁ > c₂ > zero(T)
        aumc  = logaumc(t₁, t₂, c₁, c₂)
    #elseif calcm == :log && c₁ > zero(T) && c₂ > zero(T)
        #aumc  = logaumc(t₁, t₂, c₁, c₂)
    else
        aumc  = linaumc(t₁, t₂, c₁, c₂)
    end
    return aumc
end
#---------------------------------------------------------------------------
function interpolate(t₁, t₂, tx, c₁::T, c₂::T, intpm, aftertmax) where T
    if intpm == :lint || c₁ <= zero(T) || c₂ <= zero(T)
        c = linpredict(t₁, t₂, tx, c₁, c₂)
    elseif intpm == :logt && aftertmax && c₁ > zero(T) && c₂ > zero(T)
        c = logcpredict(t₁, t₂, tx, c₁, c₂)
    elseif intpm == :luld && c₁ > c₂ > zero(T)
        c = logcpredict(t₁, t₂, tx, c₁, c₂)
    elseif intpm == :luldt && aftertmax && c₁ > c₂ > zero(T)
        c = logcpredict(t₁, t₂, tx, c₁, c₂)
    #elseif intpm == :log && c₁ > zero(T) && c₂ > zero(T)
        #c = logcpredict(t₁, t₂, tx, c₁, c₂)
    else
        c = linpredict(t₁, t₂, tx, c₁, c₂)
    end
    return c
end

################################################################################
# STEPS
################################################################################
# 1
function step_1_filterpksubj!(time, obs, dosetime)
    fobs     = firstobs(time, obs, dosetime)
    if fobs > 1
        inds = collect(1:fobs-1)
    else
        inds = Vector{Int}(undef, 0)
    end
    @inbounds for i = fobs:length(obs)
        if isnanormissing(obs[i])
            push!(inds, i)
        end
    end
    if length(inds) > 0
        deleteat!(time, inds)
        deleteat!(obs, inds)
    end
    time, obs
end
#=
function step_1_filterpksubj(time, obs, dosetime)
    fobs     = firstobs(time, obs, dosetime)
    ni = 0
    @inbounds for i = fobs:length(obs)
        isnanormissing(obs[i]) || begin ni += 1 end
    end
    inds = Vector{Int}(undef, ni)
    ni = 1
    @inbounds for i = fobs:length(obs)
        isnanormissing(obs[i]) || begin
        inds[ni] = i
        ni += 1
        end
    end
    time_cp = time[inds]
    obs_cp  = obs[inds]
    time_cp, obs_cp
end
=#
function dropkeldata!(keldata::KelData)
    if length(keldata.s) > 0 resize!(keldata.s, 0) end
    if length(keldata.e) > 0 resize!(keldata.e, 0) end
    if length(keldata.a) > 0 resize!(keldata.a, 0) end
    if length(keldata.b) > 0 resize!(keldata.b, 0) end
    if length(keldata.r) > 0 resize!(keldata.r, 0) end
    if length(keldata.ar) > 0 resize!(keldata.ar, 0) end
end
# 3
function step_3_elim!(result, data, adm, tmaxn, time_cp, obs_cp, time, keldata)
    dropkeldata!(keldata)
    obsnum = length(time_cp)
    excltime = time[data.kelrange.kelexcl]
    #keldata                = KelData()
    if data.kelauto
        if (adm != :iv && obsnum - tmaxn > 2) || (adm == :iv && obsnum - tmaxn > 1)
            if adm == :iv
                stimep = tmaxn
            else
                stimep = tmaxn + 1
            end
            timep = collect(stimep:obsnum)
            if length(data.kelrange.kelexcl) > 0
                exclinds = findall(x -> x in excltime, time_cp)
                filter!(x-> x ∉ exclinds, timep)
            end
            zcinds = findall(x -> x <= 0, obs_cp)
            filter!(x-> x ∉ zcinds, timep)
            if length(timep) > 2
                logconc    = log.(obs_cp)
                for i = length(timep)-2:-1:1
                    timepv = view(timep, i:length(timep))
                    sl = slope(view(time_cp, timepv), view(logconc, timepv))
                    if sl[1] < 0
                        push!(keldata, time_cp[timep[i]], time_cp[timep[end]], sl[1], sl[2], sl[3], sl[4])
                    end
                end
            end
        end
    else
        stimep = findfirst(x -> x >= time[data.kelrange.kelstart], time_cp)
        etimep = findlast(x -> x <= time[data.kelrange.kelend], time_cp)
        timep = collect(stimep:etimep)
        if length(data.kelrange.kelexcl) > 0
            @inbounds for i in data.kelrange.kelexcl
                filter!(x-> x ∉ findall(x -> x in excltime, time_cp), timep)
            end
        end
        zcinds = findall(x -> x <= 0, obs_cp)
        filter!(x-> x ∉ zcinds, timep)
        if length(timep) > 1
            sl = logslope(view(time_cp, timep), view(obs_cp, timep))
            push!(keldata, time_cp[stimep], time_cp[etimep], sl[1], sl[2], sl[3], sl[4])
        end
    end
    keldata, excltime
end
# 6
function step_6_areas(time_cp, obs_cp, calcm, tmaxn, tlastn)
    obsnum = length(time_cp)
    aucpartl  = Array{Float64, 1}(undef, obsnum - 1)
    aumcpartl = Array{Float64, 1}(undef, obsnum - 1)
    #Calculate all AUC/AUMC part based on data
    for i = 1:(obsnum - 1)
        aucpartl[i]  = aucpart(time_cp[i], time_cp[i + 1], obs_cp[i], obs_cp[i + 1], calcm, i >= tmaxn)
        aumcpartl[i] = aumcpart(time_cp[i], time_cp[i + 1], obs_cp[i], obs_cp[i + 1], calcm, i >= tmaxn)
    end

    #-----------------------------------------------------------------------
    #-----------------------------------------------------------------------
    auclast  = 0
    aumclast = 0
    @inbounds for i = 1:tlastn-1
        auclast  += aucpartl[i]
        aumclast += aumcpartl[i]
    end
    aucall  = auclast
    if tlastn < length(time_cp)
        @inbounds for i = tlastn:obsnum-1
            aucall  += aucpartl[i]
        end
    end
    aucpartl, aumcpartl, auclast, aumclast, aucall
end
"""
    nca(args...; kelauto = true,  elimrange = ElimRange(), dosetime = DoseTime(), kwargs...)

    nca(data, time, conc, sort; kelauto = true,  elimrange = ElimRange(), dosetime = DoseTime(), kwargs...)

    nca(data, time, conc; kelauto = true,  elimrange = ElimRange(), dosetime = DoseTime(), kwargs...)

    nca(time, conc; kelauto = true,  elimrange = ElimRange(), dosetime = DoseTime(), kwargs...)

Import data and perform NCA analysis.

Syntax simillar to [`pkimport`](@ref)

Applicable `kwargs` see  [`nca!`](@ref).
"""
function nca(args...; type::Symbol = :bps, th = NaN, bl = 0, kelauto = true,  elimrange = ElimRange(), dosetime = DoseTime(), kwargs...)
    if !(type in (:bps, :ur, :eff)) error("Unknown type") end
    if type == :bps
        pki    = pkimport(args...; kelauto = kelauto,  elimrange = elimrange, dosetime = dosetime)
    elseif type == :ur
        pki    = upkimport(args...; kelauto = kelauto,  elimrange = elimrange, dosetime = dosetime)
    elseif type == :eff
        pki    = pdimport(args...; th = th, bl = bl)
    end
    #kwargs = Dict{Symbol, Any}(kwargs)
    nca!(pki; kwargs...)
end
"""
    nca!(data::PKSubject{T,O}; adm = :ev, calcm = :lint, intpm = nothing, limitrule = nothing, verbose = 0, warn = true, io::IO = stdout, modify! = nothing) where T where O

* `adm` - administration:
    - `:ev` - extra vascular;
    - `:iv` - intravascular bolus;
* `calcm` - AUC/AUMC calculation method:
    - `:lint` - linear trapezoidal;
    - `:logt` - log-trapezoidal after Tmax;
    - `:luld` - linar up log down;
    - `:luldt` - linear up log down after Tmax;
* `intpm` - interpolation method:
    - `:lint` - linear trapezoidal;
    - `:logt` - log-trapezoidal after Tmax;
    - `:luld` - linar up log down;
    - `:luldt` - linear up log down after Tmax;
* `limitrule` use limitrule for data;
* `verbose` - print to `io`, 1: partial areas table, 2: 1, and results;
* `warn` - show warnings;
* `io` - output stream;
* `modify!` - function to modify output paramaters, call `modify!(data, result)` if difined.

Results:

* Cmax
* Tmax
* Cdose
* Tlag
* Clast
* AUClast
* AUMClast
* AUCall
* Rsq
* ARsq
* Kel
* HL
* LZint
* Clast_pred
* AUCinf
* AUCinf_pred
* AUMCinf
* AUMCinf_pred
* AUCpct
* MRTlast
* MRTinf
* MRTinf_pred
* Cllast
* Clinf
* Vzlast
* Vzinf
* Vssinf

Stable state (tau used):

* AUCtau
* AUMCtau
* Ctau
* Cavg
* Ctaumin
* Accind
* Fluc
* Fluctau
* Swing
* Swingtau
* MRTtauinf
* Cltau
* Vztau
"""
function nca!(data::PKSubject{T,O}; adm = :ev, calcm = :lint, intpm = nothing, limitrule::LimitRule = LimitRule(), verbose = 0, warn = true, io::IO = stdout, modify! = identity) where T where O

    result   = Dict{Symbol, Float64}()

    options =  Dict(:type => :bps, :adm => adm, :calcm => calcm, :intpm => intpm, :limitrule => limitrule, :verbose => verbose, :warn => warn, :modify! => modify!)

    if verbose > 0
        println(io, "  Non-compartmental Pharmacokinetic Analysis")
        if length(data.id) > 0
            print(io, "    Subject: ")
            for (k,v) in data.id
                print(io, "$(k) => $(v); ")
            end
            println(io, "")
        end
        println(io, "    Settings:")
        println(io, "    Method: $(calcm); Dose: $(data.dosetime.dose); Dose time: $(data.dosetime.time)")
        if data.dosetime.tau > 0
            println(io, "    Tau: $(data.dosetime.tau)")
        end
    end
    if isnothing(intpm) intpm = calcm end

    auctype  = promote_type(eltype(data.time), eltype(data.obs))

    if isapplicable(limitrule)
        time, obs = applylimitrule!(deepcopy(data.time), deepcopy(data.obs), limitrule)
    else
        time             = deepcopy(data.time)
        obs              = deepcopy(data.obs)
    end
################################################################################
    # STEP 1 FILTER ALL BEFORE DOSETIME AND ALL NAN OR MISSING VALUES
    time_cp, obs_cp = step_1_filterpksubj!(time, obs, data.dosetime.time)
    if length(obs_cp) < 2
        return NCAResult(data, options, result)
    end
################################################################################
    # STEP 2 - CMAX TMAX FOR TAU RANGE Clast Tlast
    result[:Obsnum] = obsnum = length(obs_cp)
    # If TAU set, calculates start and end timepoints for AUCtau
    if  data.dosetime.tau > zero(typeof(data.dosetime.tau))
        taulastp = findlast(x -> x <= data.dosetime.time + data.dosetime.tau, time_cp)
        result[:Ctaumin] = ctaumin(time_cp, obs_cp, taulastp)
    else
        taulastp = length(obs_cp)
    end
    result[:Cmax], result[:Tmax], tmaxn = ctmax(time_cp, obs_cp, taulastp)
    # C last and T last
    tlastn = 0
    @inbounds for i = length(obs_cp):-1:1
        if obs_cp[i] > zero(O)
            result[:Tlast]   = time_cp[i]
            result[:Clast]   = obs_cp[i]
            tlastn           = i
            break
        end
    end
################################################################################
    # STEP 3
    # Elimination
    keldata, excltime = step_3_elim!(result, data, adm, tmaxn, time_cp, obs_cp, data.time, data.keldata)
################################################################################
    # STEP 4
    if  data.dosetime.time > 0
        time_cp .-= data.dosetime.time
    end
################################################################################
    # STEP 5
    # Dose concentration
    # Dosetime is first point
    cdoseins = zero(Int)
    if  iszero(first(time_cp))
        result[:Cdose] = first(obs_cp)
    # Dosetime before first point
    else
        if adm == :iv
            if  first(obs_cp) > obs_cp[2] > zero(O)
                result[:Cdose] = logcpredict(first(time_cp), time_cp[2], 0, first(obs_cp), obs_cp[2])
            else
                result[:Cdose] = first(obs_cp)
            end
        else
            if  data.dosetime.tau > zero(typeof(data.dosetime.tau))
                result[:Cdose] = result[:Ctaumin]
            else
                result[:Cdose] = zero(O)
            end
        end
        cdoseins = 1
        pushfirst!(time_cp, zero(T))
        pushfirst!(obs_cp, result[:Cdose])
        # if time-zero point added some points shoud be shifted
        taulastp += 1
        tmaxn += 1
        tlastn += 1
    end
################################################################################
    # STEP 6
    # Areas
    aucpartl, aumcpartl, auclast, aumclast, aucall = step_6_areas(time_cp, obs_cp, calcm, tmaxn, tlastn)
    result[:AUClast]   = auclast
    result[:AUMClast]  = aumclast
    result[:AUCall]    = aucall
################################################################################
    # STEP 7
    # Other parameters
    #---------------------------------------------------------------------------
    result[:MRTlast]    = result[:AUMClast] / result[:AUClast]
    #---------------------------------------------------------------------------
    if data.dosetime.dose > 0
        result[:Cllast]           = data.dosetime.dose / result[:AUClast]
        result[:Dose]             = data.dosetime.dose
    end
    #-----------------------------------------------------------------------
    #-----------------------------------------------------------------------
    tlagn = findfirst(!iszero, obs_cp)
    if tlagn > 1 result[:Tlag] = time_cp[tlagn-1] else result[:Tlag] = zero(Float64) end

    if  length(data.keldata) > 0
        #data.keldata             = keldata
        result[:ARsq], rsqn      = findmax(keldata.ar)
        result[:Rsq]             = keldata.r[rsqn]
        result[:Kel]             = abs(keldata.a[rsqn])
        result[:LZ]              = keldata.a[rsqn]
        result[:LZint]           = keldata.b[rsqn]
        result[:Rsqn]            = rsqn
        result[:Clast_pred]      = exp(result[:LZint] + result[:LZ]*result[:Tlast])
        result[:HL]              = LOG2 / result[:Kel]
        result[:AUCinf]          = result[:AUClast] + result[:Clast] / result[:Kel]
        result[:AUCinf_pred]     = result[:AUClast] + result[:Clast_pred] / result[:Kel]
        result[:AUCpct]          = (result[:AUCinf] - result[:AUClast]) / result[:AUCinf] * 100
        result[:AUMCinf]         = result[:AUMClast] + result[:Tlast] * result[:Clast] / result[:Kel] + result[:Clast] / result[:Kel] ^ 2
        result[:MRTinf]          = result[:AUMCinf] / result[:AUCinf]
        if data.dosetime.dose > 0
            result[:Vzlast]          = data.dosetime.dose / result[:AUClast] / result[:Kel]
            result[:Vzinf]           = data.dosetime.dose / result[:AUCinf] / result[:Kel]
            result[:Clinf]           = data.dosetime.dose / result[:AUCinf]
            result[:Vssinf]          = result[:Clinf] * result[:MRTinf]
        end
    else
        result[:Kel] = NaN
    end
################################################################################
    # STEP 8
    # Steady-state parameters
    if data.dosetime.tau > 0
        eaucpartl  = eaumcpartl = 0.0
        if time_cp[taulastp] < data.dosetime.tau < time_cp[end]
            result[:Ctau] = interpolate(time_cp[taulastp], time_cp[taulastp + 1], data.dosetime.tau, obs_cp[taulastp], obs_cp[taulastp + 1], intpm, true)
            eaucpartl = aucpart(time_cp[taulastp], data.dosetime.tau, obs_cp[taulastp], result[:Ctau], calcm, true)
            eaumcpartl = aumcpart(time_cp[taulastp], data.dosetime.tau, obs_cp[taulastp], result[:Ctau], calcm, true)
                #remoove part after tau
        elseif data.dosetime.tau > time_cp[end] && result[:Kel] !== NaN
                #extrapolation
            result[:Ctau] = exp(result[:LZint] + result[:LZ] * (data.dosetime.tau + data.dosetime.time))
            eaucpartl = aucpart(time_cp[end], data.dosetime.tau, obs_cp[end], result[:Ctau], calcm, true)
            eaumcpartl = aumcpart(time_cp[end], data.dosetime.tau, obs_cp[end], result[:Ctau], calcm, true)
        else
            result[:Ctau] = obs_cp[taulastp]
        end

        auctau   = eaucpartl
        aumctau  = eaumcpartl
        @inbounds for i = 1:taulastp-1
            auctau  += aucpartl[i]
            aumctau += aumcpartl[i]
        end
        result[:AUCtau]   = auctau
        result[:AUMCtau]  = aumctau

        result[:Cavg]     = result[:AUCtau]/data.dosetime.tau
        if result[:Ctaumin] != 0
            result[:Swing]    = (result[:Cmax] - result[:Ctaumin])/result[:Ctaumin]
        end
        if result[:Ctau] != 0
            result[:Swingtau] = (result[:Cmax] - result[:Ctau])/result[:Ctau]
        end
        result[:Fluc]     = (result[:Cmax] - result[:Ctaumin])/result[:Cavg] * 100
        result[:Fluctau]  = (result[:Cmax] - result[:Ctau])/result[:Cavg] * 100
        #If Kel calculated
        result[:Cltau]           = data.dosetime.dose / result[:AUCtau]
        if result[:Kel] !== NaN
            result[:Accind]   = 1 / (1 - (exp(-result[:Kel] * data.dosetime.tau)))
            result[:MRTtauinf]       = (result[:AUMCtau] + data.dosetime.tau * (result[:AUCinf] - result[:AUCtau])) / result[:AUCtau]
            result[:Vztau]           = data.dosetime.dose / result[:AUCtau] / result[:Kel]
            result[:Vsstau]          = result[:Cltau] * result[:MRTtauinf]
        end

    end
################################################################################
    # Verbose output
    if verbose > 0
        aucpartlsum  = similar(aucpartl)
        aumcpartlsum = similar(aumcpartl)
        @inbounds for i = 1:length(aucpartl)
            aucpartlsum[i]  = sum(view(aucpartl, 1:i))
            aumcpartlsum[i] = sum(view(aumcpartl, 1:i))
        end
        if  data.dosetime.time > 0
            time_cp .+= data.dosetime.time
        end
        hnames = [:Time, :Concentrtion, :AUC, :AUC_cum, :AUMC, :AUMC_cum, :Info]
        mx = metida_table(collect(time_cp), collect(obs_cp), pushfirst!(aucpartl, 0.0),  pushfirst!(aucpartlsum, 0.0), pushfirst!(aumcpartl, 0.0),  pushfirst!(aumcpartlsum, 0.0), fill("", length(obs_cp));
        names = hnames)
        if cdoseins > 0
            mx[1, 7] = "D*"
        else
            mx[1, 7] = "D"
        end
        if !isnan(result[:Kel])
            @inbounds for i = 1:length(time_cp)
                if time_cp[i] >= keldata.s[rsqn] && time_cp[i] <= keldata.e[rsqn]
                    if length(data.kelrange.kelexcl) > 0
                        if time_cp[i] in excltime
                            mx[i, 7] = "Excl"
                        else
                            mx[i, 7] = "E"
                        end
                    else
                        mx[i, 7] = "E"
                    end
                end
            end
        end
        hnames = (["Time" "Concentrtion" "AUC"  "AUC" "AUMC" "AUMC" "Info"],
                  ["" "" "" "(cumulate)" "" "(cumulate)" ""])
        PrettyTables.pretty_table(io, mx; tf = PrettyTables.tf_compact, header = hnames)
        println(io, "")
        println(io, "    Cdose: $(result[:Cdose]), Dose time: $(data.dosetime.time)")
        println(io, "    Kel start: $(keldata.s[rsqn]); end: $(keldata.e[rsqn])")
        println(io, "")
        if data.dosetime.tau < time_cp[end] && data.dosetime.tau > 0
            println(io, "    Tau + dosetime is less then end time. Interpolation used.")
            println(io, "    Ctau: $(result[:Ctau])")
            println(io, "    AUC  final part: $(eaucpartl)")
            println(io, "    AUMC final part: $(eaumcpartl)")
            println(io, "")
        end
        if verbose > 1
            println(io, "    Results:")
            PrettyTables.pretty_table(io, result; tf = PrettyTables.tf_compact)
        end
    end
################################################################################

    ncares = NCAResult(data, options, result)
    modify!(ncares)

    #-----------------------------------------------------------------------
    return ncares
end

"""
    nca!(data::DataSet{T1}; adm = :ev, calcm = :lint, intpm = nothing, verbose = 0, warn = true, io::IO = stdout) where T1 <: PKSubject{T,O}  where T  where O

Non-compartmental (NCA) analysis of pharmacokinetic (PK) data.
"""
function nca!(data::DataSet{Subj}; adm = :ev, calcm = :lint, intpm = nothing, limitrule::LimitRule = LimitRule(), verbose = 0, warn = true, io::IO = stdout, modify! = identity) where Subj <: AbstractSubject
    result = Vector{NCAResult{Subj}}(undef, length(data))
    for i = 1:length(data)
        result[i] = nca!(data[i]; adm = adm, calcm = calcm, intpm = intpm, limitrule = limitrule, verbose = verbose, warn = warn, io = io, modify! = modify!)
    end
    DataSet(result)
end

function maxconc(subj::T) where T <: PKSubject
    maximum(subj.obs)
end
function minconc(subj::T) where T <: PKSubject
    minimum(subj.obs)
end

function exrate(time, conc, vol)
    length(time) == length(conc) == length(vol) || error("")
    er = Vector{Float64}(undef, length(time))
    @inbounds for i = 1:length(time)
        er[i] = conc[i]*vol[i]/(time[i][2] - time[i][1])
    end
    er
end

function step_1_filterupksubj(time, obs, vol, dosetime)
    fobs     = firstobs(time, obs, vol, dosetime)
    ni = 0
    @inbounds for i = fobs:length(obs)
        if !isnanormissing(obs[i]) && !isnanormissing(vol[i])
            ni += 1
        end
    end
    inds = Vector{Int}(undef, ni)
    ni = 1
    @inbounds for i = fobs:length(obs)
        if !isnanormissing(obs[i]) && !isnanormissing(vol[i])
            inds[ni] = i
            ni += 1
        end
    end
    time_cp = time[inds]
    obs_cp  = obs[inds]
    vol_cp  = vol[inds]
    time_cp, obs_cp, vol_cp
end

"""
    nca!(data::UPKSubject{T, O, VOL, V}; adm = :ev, calcm = :lint, intpm = nothing, limitrule::LimitRule = LimitRule(), verbose = 0, warn = true, io::IO = stdout, modify! = identity) where T where O where VOL where V

Non-compartmental (NCA) analysis of pharmacokinetic for urine data.

Results:

* AUCall
* AUClast
* Rlast
* Maxrate
* Tmax
* AR
* Vol
* Prec
* ARsq
* Rsq
* Kel
* LZ
* LZint
* Rsqn
* HL
* AUCinf
"""
function nca!(data::UPKSubject{T, O, VOL, V}; adm = :ev, calcm = :lint, intpm = nothing, limitrule::LimitRule = LimitRule(), verbose = 0, warn = true, io::IO = stdout, modify! = identity) where T where O where VOL where V

    result   = Dict{Symbol, Float64}()

    options =  Dict(:type => :urine, :adm => adm, :calcm => calcm, :intpm => intpm, :limitrule => limitrule, :verbose => verbose, :warn => warn, :modify! => modify!)

    if verbose > 0
        println(io, "  Non-compartmental Pharmacokinetic Analysis")
        println(io, "  Matrix: urine")
        if length(data.id) > 0
            print(io, "    Subject: ")
            for (k,v) in data.id
                print(io, "$(k) => $(v); ")
            end
            println(io, "")
        end
        println(io, "    Settings:")
        println(io, "    Method: $(calcm); Dose: $(data.dosetime.dose); Dose time: $(data.dosetime.time)")
    end
    if isnothing(intpm) intpm = calcm end

    time, obs, vol = step_1_filterupksubj(data.time, data.obs, data.vol, data.dosetime.time)

    mtime  = map(x-> (x[1]+x[2])/2, time)

    exr  = exrate(time, obs, vol)

    result[:AR]   = data.obs' * data.vol
    result[:Vol]  = sum(vol)

    if time[1][1] > data.dosetime.time
        pushfirst!(mtime, 0)
    else
        pushfirst!(mtime, time[1][1])
    end
    pushfirst!(obs, 0)
    pushfirst!(vol, 0)
    pushfirst!(exr, 0)

    result[:Maxrate], result[:Tmax], tmaxn = ctmax(mtime, exr)

    if data.dosetime.dose > 0
        result[:Prec]  = result[:AR]/data.dosetime.dose * 100
    end
    obsnum = length(exr)

    lastobs = length(exr)
    for i = length(exr):-1:1
        if exr[i] > 0
            break
        else
            lastobs = i
        end
    end

    # STEP 3
    # Elimination
    keldata, excltime = step_3_elim!(result, data, adm, tmaxn, mtime, exr, data.time, data.keldata)

    #result[:Kel]
    #result[:HL]
    aucpartl  = Array{Float64, 1}(undef, obsnum - 1)
    #Calculate all AUC part based on data
    for i = 1:(obsnum - 1)
        aucpartl[i] = aucpart(mtime[i], mtime[i + 1], exr[i], exr[i + 1], calcm, i >= tmaxn)
    end

    result[:AUCall]  = sum(aucpartl)
    result[:AUClast] = sum(aucpartl[1:lastobs-1])
    result[:Rlast]  = exr[end]

    if  length(data.keldata) > 0
        result[:ARsq], rsqn      = findmax(keldata.ar)
        result[:Rsq]             = keldata.r[rsqn]
        result[:Kel]             = abs(keldata.a[rsqn])
        result[:LZ]              = keldata.a[rsqn]
        result[:LZint]           = keldata.b[rsqn]
        result[:Rsqn]            = rsqn
        result[:HL]              = LOG2 / result[:Kel]

        result[:AUCinf]          = result[:AUClast] + result[:Rlast] / result[:Kel]
        result[:AUCpct]          = (result[:AUCinf] - result[:AUClast]) / result[:AUCinf] * 100
    end

    ncares = NCAResult(data, options, result)
    modify!(ncares)
    #-----------------------------------------------------------------------
    return ncares
end


text:
using Revise # reduce recompilation time
using Plots
using LinearAlgebra
using BenchmarkTools
using StaticArrays
using UnPack
using DelimitedFiles

push!(LOAD_PATH, "./src")
using CommonUtils
using Basis1D
using Basis2DTri
using UniformTriMesh

using SetupDG

push!(LOAD_PATH, "./examples/EntropyStableEuler.jl")
using EntropyStableEuler
using EntropyStableEuler.Fluxes2D

"Approximation parameters"
N = 2 # The order of approximation
K1D = 32
# N = 2 # The order of approximation
# K1D = 15
CFL = 0.05
T = 0.2 # endtime
inviscid_dissp = true
viscous_dissp = false
Re = 10000

# Becker viscous shocktube
const γ = 1.4
const M_0 = 3.0
const mu = 0.01
const lambda = 2/3*mu
const Pr = 3/4
const cp = γ/(γ-1)
const cv = 1/(γ-1)
const kappa = mu*cp/Pr

const v_inf = 0.2
const rho_0 = 1.0
const m_0 = 1.0
const v_0 = 1.0
const v_1 = (γ-1+2/M_0^2)/(γ+1)
const v_01 = sqrt(v_0*v_1)

const uL = v_0+v_inf
const uR = v_1+v_inf
const vL = 0.0
const vR = 0.0
const rhoL = m_0/v_0
const rhoR = m_0/v_1
const eL = 1/(2*γ)*((γ+1)/(γ-1)*v_01^2-v_0^2)
const eR = 1/(2*γ)*((γ+1)/(γ-1)*v_01^2-v_1^2)
const pL = (γ-1)*rhoL*eL
const pR = (γ-1)*rhoR*eR
const EL = pL/(γ-1)+0.5*rhoL*uL^2
const ER = pR/(γ-1)+0.5*rhoR*uR^2

"Time integration"
rk4a,rk4b,rk4c = rk45_coeffs()
CN = (N+1)*(N+2)/2  # estimated trace constant for CFL
h  = 2/K1D
dt = CFL * h / CN
Nsteps = convert(Int,ceil(T/dt))
dt = T/Nsteps
dt0 = dt

"Mesh related variables"
VX, VY, EToV = uniform_tri_mesh(Int(K1D/2*3),K1D)
@. VX = (VX/4*3)+1/4
@. VY = (VY+1)/2

# initialize ref element and mesh
rd = init_reference_tri(N)
md = init_mesh((VX,VY),EToV,rd)

# Make domain periodic
@unpack Nfaces,Vf = rd
@unpack xf,yf,K,mapM,mapP,mapB = md
LX,LY = (x->maximum(x)-minimum(x)).((VX,VY)) # find lengths of domain
mapPB = build_periodic_boundary_maps(xf,yf,LX,LY,Nfaces*K,mapM,mapP,mapB)
mapP[mapB] = mapPB
@pack! md = mapP

# construct hybridized SBP operators
@unpack M,Dr,Ds,Vq,Pq,Vf,wf,nrJ,nsJ = rd
Qr = Pq'*M*Dr*Pq
Qs = Pq'*M*Ds*Pq
Ef = Vf*Pq
Br = diagm(wf.*nrJ)
Bs = diagm(wf.*nsJ)
Qrh = .5*[Qr-Qr' Ef'*Br;
         -Br*Ef  Br]
Qsh = .5*[Qs-Qs' Ef'*Bs;
         -Bs*Ef  Bs]

Vh = [Vq;Vf]
Ph = M\transpose(Vh)
VhP = Vh*Pq

# make sparse skew symmetric versions of the operators"
# precompute union of sparse ids for Qr, Qs
Qrhskew = .5*(Qrh-transpose(Qrh))
Qshskew = .5*(Qsh-transpose(Qsh))

# interpolate geofacs to both vol/surf nodes
@unpack rxJ,sxJ,ryJ,syJ = md
rxJ,sxJ,ryJ,syJ = (x->Vh*x).((rxJ,sxJ,ryJ,syJ)) # interp to hybridized points
@pack! md = rxJ,sxJ,ryJ,syJ

# pack SBP operators into tuple
@unpack LIFT = rd
ops = (Qrhskew,Qshskew,VhP,Ph,LIFT,Vq)

function euler_fluxes_2D(UL,UR)
    rhoL = first(UL); betaL = last(UL)
    rhoR = first(UR); betaR = last(UR)
    logL = (log.(rhoL),log.(betaL))
    logR = (log.(rhoR),log.(betaR))
    return euler_fluxes(UL,UR,logL,logR)
end

function euler_fluxes_2D(UL,UR,logL,logR)
    return euler_fluxes(UL...,UR...,logL...,logR...)
end

# 2d version
function euler_fluxes_2D(rhoL,uL,vL,betaL,rhoR,uR,vR,betaR,
                      rhologL,betalogL,rhologR,betalogR)

    rholog = logmean.(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean.(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = (@. .5*(rhoL+rhoR))
    uavg   = (@. .5*(uL+uR))
    vavg   = (@. .5*(vL+vR))

    unorm = (@. uL*uR + vL*vR)
    pa    = (@. rhoavg/(betaL+betaR))
    f4aux = (@. rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm)

    FxS1 = (@. rholog*uavg)
    FxS2 = (@. FxS1*uavg + pa)
    FxS3 = (@. FxS1*vavg)
    FxS4 = (@. f4aux*uavg)

    FyS1 = (@. rholog*vavg)
    FyS2 = (@. FxS3)
    FyS3 = (@. FyS1*vavg + pa)
    FyS4 = (@. f4aux*vavg)
    return (FxS1,FxS2,FxS3,FxS4),(FyS1,FyS2,FyS3,FyS4)
end

function init_BC_funs(md::MeshData)
    @unpack xf,yf,mapP,mapB,nxJ,nyJ,sJ = md
    xb,yb = (x->x[mapB]).((xf,yf))

    leftwall     = mapB[findall(@. abs(xb+0.5) < 1e-12)]
    rightwall    = mapB[findall(@. abs(xb-1.0) < 1e-12)]
    topwall      = mapB[findall(@. abs(yb-1.0) < 1e-12)]
    bottomwall   = mapB[findall(@. abs(yb) < 1e-12)]
    vwall = [leftwall;rightwall]

    function impose_BCs_inviscid!(QP,Qf,md::MeshData)
        QP[1][leftwall] .= rhoL
        QP[2][leftwall] .= uL
        QP[3][leftwall] .= vL
        QP[4][leftwall] .= rhoL/(2*pL)

        QP[1][rightwall] .= Qf[1][rightwall]#rhoR
        QP[2][rightwall] .= Qf[2][rightwall]#rhoR*uR
        QP[3][rightwall] .= Qf[3][rightwall]#rhoR*vR
        QP[4][rightwall] .= Qf[4][rightwall]#rhoR/(2*pR)
    end

    function impose_BCs_lam!(lamP,lam,md)
        lam[leftwall]   .= 0.0
        lam[rightwall]  .= 0.0
        lamP[leftwall]  .= 0.0
        lamP[rightwall] .= 0.0
    end

    function impose_BCs_entropyvars!(VUP,VUf,md::MeshData)
        VL = v_ufun(rhoL,rhoL*uL,rhoL*vL,EL)
        VR = v_ufun(rhoR,rhoR*uR,rhoR*vR,ER)

        VUP[1][leftwall] .= VL[1]
        VUP[2][leftwall] .= VL[2]
        VUP[3][leftwall] .= VL[3]
        VUP[4][leftwall] .= VL[4]

        VUP[1][rightwall] .= VUf[1][rightwall]#VR[1]
        VUP[2][rightwall] .= VUf[2][rightwall]#VR[2]
        VUP[3][rightwall] .= VUf[3][rightwall]#VR[3]
        VUP[4][rightwall] .= VUf[4][rightwall]#VR[4]
    end

    function impose_BCs_stress!(σxP,σyP,σxf,σyf,VUf,md::MeshData)
        # Adiabatic no-slip BC
        @. σxP[1][vwall] = σxf[1][vwall]
        @. σyP[1][vwall] = σyf[1][vwall]
        @. σxP[2][vwall] = σxf[2][vwall]
        @. σyP[2][vwall] = σyf[2][vwall]
        @. σxP[3][vwall] = σxf[3][vwall]
        @. σyP[3][vwall] = σyf[3][vwall]
        @. σxP[4][vwall] = σxf[4][vwall]
        @. σyP[4][vwall] = σyf[4][vwall]
    end
    return impose_BCs_inviscid!,impose_BCs_lam!,impose_BCs_entropyvars!,impose_BCs_stress!
end
impose_BCs_inviscid!,impose_BCs_lam!,impose_BCs_entropyvars!,impose_BCs_stress! = init_BC_funs(md)

function update_flux!(rhsQ,QM,QP,Uf,UP,LFc,nxJ,nyJ,mapP,Nfq,K,Nd,inviscid_dissp)
    for k = 1:K
        for i = 1:Nfq
            fx,fy = euler_fluxes_2D((x->x[i,k]).(QP),(x->x[i,k]).(QM))
            tmp_nxJ = nxJ[i,k]
            tmp_nyJ = nyJ[i,k]
            tmp_LF = LFc[i,k]
            for d = 1:Nd
                if inviscid_dissp
                    rhsQ[d][i,k] = fx[d]*tmp_nxJ+fy[d]*tmp_nyJ-tmp_LF*(UP[d][i,k]-Uf[d][i,k])
                else
                    rhsQ[d][i,k] = fx[d]*tmp_nxJ+fy[d]*tmp_nyJ
                end
            end
        end
    end
end

function flux_differencing!(QF,Qh,Qrh,Qsh,rxJ,ryJ,sxJ,syJ,Nh,Nq,Nd,K)
    for k = 1:K
        rx = rxJ[1,k]
        ry = ryJ[1,k]
        sx = sxJ[1,k]
        sy = syJ[1,k]
        for j = 1:Nh # col idx
            for i = j:Nh # row idx
                if i <= Nq || j <= Nq # Skip lower right block
                    Fx,Fy = euler_fluxes_2D((x->x[i,k]).(Qh),(x->x[j,k]).(Qh))
                    Qr = Qrh[i,j]
                    Qs = Qsh[i,j]
                    for d = 1:Nd
                        update_val = 2*((rx*Qr+sx*Qs)*Fx[d]+(ry*Qr+sy*Qs)*Fy[d])
                        QF[d][i,k] += update_val
                        QF[d][j,k] -= update_val
                    end
                end
            end
        end
    end
end

function rhs_inviscid!(Q,md::MeshData,ops,flux_fun,compute_rhstest,inviscid_dissp,VU,Qh,QF,QM,QP,Uf,UP,rhsQ,tmp,tmp2,lam,LFc)
    @unpack rxJ,sxJ,ryJ,syJ,sJ,J,wJq = md
    @unpack nxJ,nyJ,sJ,mapP,mapB,K = md
    Qrh,Qsh,VhP,Ph,Lf,Vq = ops
    Nh,Nq = size(VhP)
    Nfq = size(Lf,2)
    Nd = length(Q)

    fill!.(QF,0.0)

    # entropy var projection
    # VU = v_ufun((x->Vq*x).(Q)...)
    Q = (x->Vq*x).(Q)
    vector_norm(U) = sum((x->x.^2).(U))
    @. VU[3] = Q[2]^2+Q[3]^2 # rhoUnorm
    @. VU[4] = Q[4]-.5*VU[3]/Q[1] # rhoe
    @. VU[1] = log(0.4*VU[4]/(Q[1]^1.4)) # sU #TODO: hardcoded gamma
    @. VU[1] = (-Q[4]+VU[4]*(2.4-VU[1]))/VU[4]
    @. VU[2] = Q[2]/VU[4]
    @. VU[3] = Q[3]/VU[4]
    @. VU[4] = -Q[1]/VU[4]

    # VU = (x->VhP*x).(VU)
    Uh = (x->VhP*x).(VU)

    # Uh = u_vfun(VU...)
    @. tmp = Uh[2]^2+Uh[3]^2 #vUnorm
    @. tmp2 = (0.4/((-Uh[4])^1.4))^(1/0.4)*exp(-(1.4 - Uh[1] + tmp/(2*Uh[4]))/0.4) # rhoeV
    @. Uh[1] = tmp2.*(-Uh[4])
    @. Uh[2] = tmp2.*Uh[2]
    @. Uh[3] = tmp2.*Uh[3]
    @. Uh[4] = tmp2.*(1-tmp/(2*Uh[4]))

    # # convert to rho,u,v,beta vars
    # (rho,rhou,rhov,E) = Uh
    # beta = betafun(rho,rhou,rhov,E)
    # Qh = (rho, rhou./rho, rhov./rho, beta) # redefine Q = (rho,u,v,β)
    @. tmp = Uh[1]/(2*0.4*(Uh[4]-.5*(Uh[2]^2+Uh[3]^2)/Uh[1])) #beta
    @. Qh[1] = Uh[1]
    @. Qh[2] = Uh[2]./Uh[1]
    @. Qh[3] = Uh[3]./Uh[1]
    @. Qh[4] = tmp

    # # compute face values
    # QM = (x->x[Nq+1:Nh,:]).(Qh)
    # QP = (x->x[mapP]).(QM)
    # impose_BCs_inviscid!(QP,QM,md)
    for i = 1:Nd
        QM[i] .= Qh[i][Nq+1:Nh,:]
        QP[i] .= QM[i][mapP]
    end
    impose_BCs_inviscid!(QP,QM,md)

    # simple lax friedrichs dissipation
    for i = 1:Nd
        Uf[i] .= Uh[i][Nq+1:Nh,:]
    end
    (rhoM,rhouM,rhovM,EM) = Uf
    rhoUM_n = @. (rhouM*nxJ + rhovM*nyJ)/sJ
    @. lam  = abs(sqrt(abs(rhoUM_n/rhoM))+sqrt(γ*(γ-1)*(EM-.5*rhoUM_n^2/rhoM)/rhoM))
    lamP = lam[mapP]
    impose_BCs_lam!(lamP,lam,md)
    @. LFc = .25*max(lam,lamP)*sJ

    for i = 1:Nd
        UP[i] .= Uf[i][mapP]
    end
    update_flux!(rhsQ,QM,QP,Uf,UP,LFc,nxJ,nyJ,mapP,Nfq,K,Nd,inviscid_dissp)
    rhsQ = (x->Lf*x).(rhsQ)

    flux_differencing!(QF,Qh,Qrh,Qsh,rxJ,ryJ,sxJ,syJ,Nh,Nq,Nd,K)
    rhsQ = (x->Ph*x).(QF).+rhsQ
    rhsQ = (x -> -x./J).(rhsQ)

    rhstest = 0.0
    if compute_rhstest
        for fld in eachindex(rhsQ)
            VUq = VU[fld][1:Nq,:]
            rhstest += sum(wJq.*VUq.*(Vq*rhsQ[fld]))
        end
    end
    return rhsQ,rhstest
end

function dg_grad!(Qx,Qy,Q,Qf,QP,md::MeshData,rd::RefElemData)
    @unpack rxJ,sxJ,ryJ,syJ,nxJ,nyJ,mapP,J = md
    @unpack Dr,Ds,LIFT,Vf = rd
    Np = size(Dr,1)
    Nd = size(Qx,1)

    rxj,sxj,ryj,syj = (x->@view x[1:Np,:]).((rxJ,sxJ,ryJ,syJ))
    Qr = (x->Dr*x).(Q)
    Qs = (x->Ds*x).(Q)

    volx(ur,us) = @. rxj*ur + sxj*us
    voly(ur,us) = @. ryj*ur + syj*us
    surfx(uP,uf) = LIFT*(@. .5*(uP-uf)*nxJ)
    surfy(uP,uf) = LIFT*(@. .5*(uP-uf)*nyJ)
    rhsx = volx.(Qr,Qs) .+ surfx.(QP,Qf)
    rhsy = voly.(Qr,Qs) .+ surfy.(QP,Qf)
    for i = 1:Nd
        Qx[i] .= rhsx[i]./J
        Qy[i] .= rhsy[i]./J
    end
    #return (x->x./J).(rhsx),(x->x./J).(rhsy)
end

function dg_div!(rhs,Qx,Qxf,QxP,Qy,Qyf,QyP,md::MeshData,rd::RefElemData)

    @unpack rxJ,sxJ,ryJ,syJ,nxJ,nyJ,mapP,J = md
    @unpack Dr,Ds,LIFT,Vf = rd
    Nd = length(Qx)

    Np = size(Dr,1)
    rxj,sxj,ryj,syj = (x->@view x[1:Np,:]).((rxJ,sxJ,ryJ,syJ))
    Qxr = (x->Dr*x).(Qx)
    Qxs = (x->Ds*x).(Qx)
    Qyr = (x->Dr*x).(Qy)
    Qys = (x->Ds*x).(Qy)

    vol(uxr,uxs,uyr,uys) = @. rxj*uxr + sxj*uxs + ryj*uyr + syj*uys
    surf(uxP,uxf,uyP,uyf) = LIFT*(@. .5*((uxP-uxf)*nxJ + (uyP-uyf)*nyJ))
    #rhs = vol.(Qxr,Qxs,Qyr,Qys) .+ surf.(QxP,Qxf,QyP,Qyf)
    for i = 1:Nd
        rhs[i] .= (vol(Qxr[i],Qxs[i],Qyr[i],Qys[i]) .+ surf(QxP[i],Qxf[i],QyP[i],Qyf[i]))./J
    end
    #rhs .= (x->x./J).(rhs)
    #return (x->x./J).(rhs)
end

function init_visc_fxn(λ,μ,Pr)
    let λ=-λ,μ=μ,Pr=Pr
        function viscous_matrices!(Kxx,Kxy,Kyy,v)
            v1,v2,v3,v4 = v
            inv_v4_cubed = @. 1/(v4^3)
            λ2μ = (λ+2.0*μ)
            Kxx[2,2] = inv_v4_cubed*-λ2μ*v4^2
            Kxx[2,4] = inv_v4_cubed*λ2μ*v2*v4
            Kxx[3,3] = inv_v4_cubed*-μ*v4^2
            Kxx[3,4] = inv_v4_cubed*μ*v3*v4
            Kxx[4,2] = inv_v4_cubed*λ2μ*v2*v4
            Kxx[4,3] = inv_v4_cubed*μ*v3*v4
            Kxx[4,4] = inv_v4_cubed*-(λ2μ*v2^2 + μ*v3^2 - γ*μ*v4/Pr)

            Kxy[2,3] = inv_v4_cubed*-λ*v4^2
            Kxy[2,4] = inv_v4_cubed*λ*v3*v4
            Kxy[3,2] = inv_v4_cubed*-μ*v4^2
            Kxy[3,4] = inv_v4_cubed*μ*v2*v4
            Kxy[4,2] = inv_v4_cubed*μ*v3*v4
            Kxy[4,3] = inv_v4_cubed*λ*v2*v4
            Kxy[4,4] = inv_v4_cubed*(λ+μ)*(-v2*v3)

            Kyy[2,2] = inv_v4_cubed*-μ*v4^2
            Kyy[2,4] = inv_v4_cubed*μ*v2*v4
            Kyy[3,3] = inv_v4_cubed*-λ2μ*v4^2
            Kyy[3,4] = inv_v4_cubed*λ2μ*v3*v4
            Kyy[4,2] = inv_v4_cubed*μ*v2*v4
            Kyy[4,3] = inv_v4_cubed*λ2μ*v3*v4
            Kyy[4,4] = inv_v4_cubed*-(λ2μ*v3^2 + μ*v2^2 - γ*μ*v4/Pr)
        end
        return viscous_matrices!
    end
end
viscous_matrices! = init_visc_fxn(lambda,mu,Pr)

function rhs_viscous!(Q,md::MeshData,rd::RefElemData,Re,viscous_dissp,rhs,VU,VUx,VUy,sigma_x,sigma_y,penalization,Kxx,Kyy,Kxy)
    @unpack Pq,Vq,Vf,LIFT = rd
    @unpack K,mapP,mapB,J,wJq = md
    Nd = length(Q)
    Nq = size(Vq,1)
    Np = size(Pq,1)
    Nfq = size(Vf,1)

    fill!.(rhs,0.0)
    fill!.(sigma_x,0.0)
    fill!.(sigma_y,0.0)

    # entropy var projection
    # VU = v_ufun((x->Vq*x).(Q)...)
    Q = (x->Vq*x).(Q)
    vector_norm(U) = sum((x->x.^2).(U))
    @. VU[3] = Q[2]^2+Q[3]^2 # rhoUnorm
    @. VU[4] = Q[4]-.5*VU[3]/Q[1] # rhoe
    @. VU[1] = log(0.4*VU[4]/(Q[1]^1.4)) # sU #TODO: hardcoded gamma
    @. VU[1] = (-Q[4]+VU[4]*(2.4-VU[1]))/VU[4]
    @. VU[2] = Q[2]/VU[4]
    @. VU[3] = Q[3]/VU[4]
    @. VU[4] = -Q[1]/VU[4]
    VU = (x->Pq*x).(VU)

    # compute and interpolate to quadrature
    VUf = (x->Vf*x).(VU)
    VUP = (x->x[mapP]).(VUf)
    impose_BCs_entropyvars!(VUP,VUf,md)

    dg_grad!(VUx,VUy,VU,VUf,VUP,md,rd)
    VUx = (x->Vq*x).(VUx)
    VUy = (x->Vq*x).(VUy)
    VUq = (x->Vq*x).(VU)

    # initialize sigma_x,sigma_y = viscous rhs
    vqi = zeros(Float64,Nd)
    for e = 1:K
        # mult by matrices and perform local projections
        for i = 1:Nq
            @. vqi = getindex(VUq,i,e)
            viscous_matrices!(Kxx,Kxy,Kyy,vqi)
            for col = 2:Nd
                vxi_col = VUx[col][i,e]
                vyi_col = VUy[col][i,e]
                for row = 2:Nd
                    sigma_x[row][i,e] += Kxx[row,col]*vxi_col + Kxy[row,col]*vyi_col
                    sigma_y[row][i,e] += Kxy[col,row]*vxi_col + Kyy[row,col]*vyi_col
                end
            end
        end
    end

    rhstest = 0.0
    for fld in 1:4
        rhstest += sum(wJq.*VUx[fld].*sigma_x[fld])
        rhstest += sum(wJq.*VUy[fld].*sigma_y[fld])
    end

    sigma_x = (x->Pq*x).(sigma_x)
    sigma_y = (x->Pq*x).(sigma_y)

    sxf = (x->Vf*x).(sigma_x)
    syf = (x->Vf*x).(sigma_y)
    sxP = (x->x[mapP]).(sxf)
    syP = (x->x[mapP]).(syf)
    impose_BCs_stress!(sxP,syP,sxf,syf,VUf,md)

    if viscous_dissp
        # # add penalty
        # #tau = .5
        # tau = @. -1/Re/VUf[4]
        # dV = ((xP,x)->xP-x).(VUP,VUf)
        # avgV = ((xP,x)->1/2*(xP+x)).(VUP,VUf)
        # @. penalization[2] = tau*dV[2]
        # @. penalization[3] = tau*dV[3]
        # @. penalization[4] = tau*dV[4]

        # tau = tau[mapB]
        # @. penalization[2][mapB] = tau*dV[2][mapB]
        # @. penalization[3][mapB] = tau*dV[3][mapB]
        # if BCTYPE == 1
        #     @. penalization[4][mapB] = -tau*(avgV[2][mapB]*dV[2][mapB]
        #                                     +avgV[3][mapB]*dV[3][mapB])/VUf[4][mapB]
        # else
        #     @. penalization[4][mapB] = -tau*(avgV[2][mapB]*dV[2][mapB]
        #                                     +avgV[3][mapB]*dV[3][mapB]
        #                                     +  dV[4][mapB]*dV[4][mapB]/2)/VUf[4][mapB]
        # end

        # penalization = (x->LIFT*x).(penalization)
    end

    dg_div!(rhs,sigma_x,sxf,sxP,sigma_y,syf,syP,md,rd)

    if viscous_dissp
        return rhs .+ penalization, rhstest
    else
        return rhs, rhstest
    end
end

#####################
### Time Stepping ###
#####################
bcopy!(x,y) = x .= y

# define initial conditions at nodes
@unpack x,y = md
Ma = .3
rho = ones(size(x))
u = zeros(size(x))
v = zeros(size(x))
p = (1/(Ma^2*γ))*ones(size(x))

"""Initial condition"""
function bisection_solve_velocity(x,max_iter,tol)
    v_L = v_1
    v_R = v_0
    num_iter = 0

    L_k = kappa/m_0/cv
    f(v) = -x+2*L_k/(γ+1)*(v_0/(v_0-v_1)*log((v_0-v)/(v_0-v_01))-v_1/(v_0-v_1)*log((v-v_1)/(v_01-v_1)))

    v_new = (v_L+v_R)/2
    while num_iter < max_iter
        v_new = (v_L+v_R)/2

        if abs(f(v_new)) < tol
            return v_new
        elseif sign(f(v_L)) == sign(f(v_new))
            v_L = v_new
        else
            v_R = v_new
        end
        num_iter += 1
    end

    return v_new
end

const max_iter = 100
const tol = 1e-14

function exact_sol_viscous_shocktube(x,t)
    u   = bisection_solve_velocity(x-v_inf*t,max_iter,tol)
    rho = m_0/u
    e   = 1/(2*γ)*((γ+1)/(γ-1)*v_01^2-u^2)
    return rho, rho*(v_inf+u), 0.0, rho*(e+1/2*(v_inf+u)^2)
end

Q = exact_sol_viscous_shocktube.(x,0.0)
Q = ([x[1] for x in Q], [x[2] for x in Q], [x[3] for x in Q], [x[4] for x in Q])

#Q = primitive_to_conservative(rho,u,v,p)
Q = collect(Q) # make Q,resQ arrays of arrays for mutability
resQ = [zeros(size(x)) for i in eachindex(Q)]

@unpack Pq,Vq,Vf,LIFT = rd
@unpack K,mapP,mapB,J = md
Qrh,Qsh,VhP,Ph,Lf,Vq = ops
Nh,Nq = size(VhP)
Nfq = size(Lf,2)
Nd = length(Q)
Nq = size(Vq,1)
Np = size(Pq,1)
Nfq = size(Vf,1)
rhs = (zeros(Float64,Np,K),zeros(Float64,Np,K),zeros(Float64,Np,K),zeros(Float64,Np,K))
VU = (zeros(Float64,Nq,K),zeros(Float64,Nq,K),zeros(Float64,Nq,K),zeros(Float64,Nq,K))
VUx = (zeros(Float64,Np,K),zeros(Float64,Np,K),zeros(Float64,Np,K),zeros(Float64,Np,K))
VUy = (zeros(Float64,Np,K),zeros(Float64,Np,K),zeros(Float64,Np,K),zeros(Float64,Np,K))
sigma_x = (zeros(Float64,Nq,K),zeros(Float64,Nq,K),zeros(Float64,Nq,K),zeros(Float64,Nq,K))
sigma_y = (zeros(Float64,Nq,K),zeros(Float64,Nq,K),zeros(Float64,Nq,K),zeros(Float64,Nq,K))
penalization = (zeros(Float64,Nfq,K),zeros(Float64,Nfq,K),zeros(Float64,Nfq,K),zeros(Float64,Nfq,K))
Kxx = MMatrix{4,4,Float64}(zeros(Nd,Nd))
Kyy = MMatrix{4,4,Float64}(zeros(Nd,Nd))
Kxy = MMatrix{4,4,Float64}(zeros(Nd,Nd))
Qh = (zeros(Float64,Nh,K),zeros(Float64,Nh,K),zeros(Float64,Nh,K),zeros(Float64,Nh,K))
QF = (zeros(Float64,Nh,K),zeros(Float64,Nh,K),zeros(Float64,Nh,K),zeros(Float64,Nh,K))
QM = (zeros(Float64,Nfq,K),zeros(Float64,Nfq,K),zeros(Float64,Nfq,K),zeros(Float64,Nfq,K))
QP = (zeros(Float64,Nfq,K),zeros(Float64,Nfq,K),zeros(Float64,Nfq,K),zeros(Float64,Nfq,K))
Uf = (zeros(Float64,Nfq,K),zeros(Float64,Nfq,K),zeros(Float64,Nfq,K),zeros(Float64,Nfq,K))
UP = (zeros(Float64,Nfq,K),zeros(Float64,Nfq,K),zeros(Float64,Nfq,K),zeros(Float64,Nfq,K))
rhsQ_tmp = (zeros(Float64,Nfq,K),zeros(Float64,Nfq,K),zeros(Float64,Nfq,K),zeros(Float64,Nfq,K))
tmp = zeros(Float64,Nh,K)
tmp2 = zeros(Float64,Nh,K)
lam = zeros(Float64,Nfq,K)
LFc = zeros(Float64,Nfq,K)

function dopri45_coeffs()
    rk4a = [0.0             0.0             0.0             0.0             0.0             0.0         0.0
            0.2             0.0             0.0             0.0             0.0             0.0         0.0
            3.0/40.0        9.0/40.0        0.0             0.0             0.0             0.0         0.0
            44.0/45.0      -56.0/15.0       32.0/9.0        0.0             0.0             0.0         0.0
            19372.0/6561.0 -25360.0/2187.0  64448.0/6561.0  -212.0/729.0    0.0             0.0         0.0
            9017.0/3168.0  -355.0/33.0      46732.0/5247.0  49.0/176.0      -5103.0/18656.0 0.0         0.0
            35.0/384.0      0.0             500.0/1113.0    125.0/192.0     -2187.0/6784.0  11.0/84.0   0.0 ]

    rk4c = vec([0.0 0.2 0.3 0.8 8.0/9.0 1.0 1.0 ])

    # coefficients to evolve error estimator = b1-b2
    rk4E = vec([71.0/57600.0  0.0 -71.0/16695.0 71.0/1920.0 -17253.0/339200.0 22.0/525.0 -1.0/40.0 ])

    return rk4a,rk4E,rk4c
end

function rhsRK!(Q,rd,md,Re,ops,euler_fluxes,inviscid_dissp,viscous_dissp,rhs,VU,Qh,QF,QM,QP,Uf,UP,rhsQ,tmp,tmp2,lam,LFc,VUx,VUy,sigma_x,sigma_y,penalization,Kxx,Kyy,Kxy)
    rhsQ,_ = rhs_inviscid!(Q,md,ops,euler_fluxes,false,inviscid_dissp,VU,Qh,QF,QM,QP,Uf,UP,rhsQ,tmp,tmp2,lam,LFc)
    visc_rhsQ,visc_test = rhs_viscous!(Q,md,rd,Re,viscous_dissp,rhs,VU,VUx,VUy,sigma_x,sigma_y,penalization,Kxx,Kyy,Kxy)
    bcopy!.(rhsQ, @. rhsQ + visc_rhsQ)

    let Pq=rd.Pq, Vq=rd.Vq, wJq=md.wJq
        rhstest = 0.0
        rhstest_visc = 0.0
        VU = v_ufun((x->Vq*x).(Q)...)
        VUq = (x->Vq*Pq*x).(VU)
        for field in eachindex(rhsQ)
            rhstest += sum(wJq.*VUq[field].*(Vq*rhsQ[field]))
            rhstest_visc += sum(wJq.*VUq[field].*(Vq*visc_rhsQ[field]))
        end
        rhstest_visc = rhstest_visc + visc_test
        return rhsQ,rhstest,rhstest_visc
    end
end

rka,rkE,rkc = dopri45_coeffs()

# DOPRI storage
Qtmp = similar.(Q)
rhsQrk = ntuple(x->zero.(Q),length(rkE))

errEst = 0.0
prevErrEst = 0.0

t = 0.0
i = 0
interval = 5

dthist = Float64[dt]
thist = Float64[0.0]
errhist = Float64[0.0]
vischist = Float64[0.0]
rhstesthist = Float64[0.0]
wsJ = diagm(rd.wf)*md.sJ

prevQ = [zeros(size(x)) for i in eachindex(Q)]

rhsQ,_,_ = rhsRK!(Q,rd,md,Re,ops,euler_fluxes_2D,inviscid_dissp,viscous_dissp,rhs,VU,Qh,QF,QM,QP,Uf,UP,rhsQ_tmp,tmp,tmp2,lam,LFc,VUx,VUy,sigma_x,sigma_y,penalization,Kxx,Kyy,Kxy)
bcopy!.(rhsQrk[1],rhsQ) # initialize DOPRI rhs (FSAL property)
@time begin
while t < T
    # DOPRI step and
    rhstest = 0.0
    rhstest_visc = 0.0
    for INTRK = 2:7
        k = zero.(Qtmp)
        for s = 1:INTRK-1
            bcopy!.(k, @. k + rka[INTRK,s]*rhsQrk[s])
        end
        bcopy!.(Qtmp, @. Q + dt*k)
        rhsQ,rhstest,rhstest_visc = rhsRK!(Qtmp,rd,md,Re,ops,euler_fluxes_2D,inviscid_dissp,viscous_dissp,rhs,VU,Qh,QF,QM,QP,Uf,UP,rhsQ_tmp,tmp,tmp2,lam,LFc,VUx,VUy,sigma_x,sigma_y,penalization,Kxx,Kyy,Kxy)
        bcopy!.(rhsQrk[INTRK],rhsQ)
    end
    errEstVec = zero.(Qtmp)
    for s = 1:7
        bcopy!.(errEstVec, @. errEstVec + rkE[s]*rhsQrk[s])
    end

    errTol = 1e-5
    errEst = 0.0
    for field = 1:length(Qtmp)
        errEstScale = @. abs(errEstVec[field]) / (errTol*(1+abs(Q[field])))
        errEst += sum(errEstScale.^2) # hairer seminorm
    end
    errEst = sqrt(errEst/(length(Q[1])*4))
    if errEst < 1.0 # if err small, accept step and update
            bcopy!.(Q, Qtmp)
            global t += dt
            bcopy!.(rhsQrk[1], rhsQrk[7]) # use FSAL property
    end
    order = 5
    dtnew = .8*dt*(.9/errEst)^(.4/(order+1)) # P controller
    if i > 0 # use PI controller if prevErrEst available
            dtnew *= (prevErrEst/max(1e-14,errEst))^(.3/(order+1))
    end
    global dt = max(min(10*dt0,dtnew),1e-9) # max/min dt
    global prevErrEst = errEst

    push!(dthist,dt)
    push!(thist,t)
    push!(vischist,rhstest_visc)
    push!(rhstesthist,rhstest)

    global i = i + 1  # number of total steps attempted
    if i%interval==0
        global prevQ
        preverr = sum(norm.(Q .- prevQ))
        bcopy!.(prevQ, Q)
        println("i = $i, t = $t, dt = $dtnew, errEst = $errEst, \nrhstest = $rhstest, rhstest_visc = $rhstest_visc, preverr = $preverr")
    end
end

end

################
### Plotting ###
################

#plotting nodes
@unpack Vp = rd
gr(aspect_ratio=:equal,legend=false,
   markerstrokewidth=0,markersize=2)
#vv = Vp*(@. (Q[2]/Q[1])^2 + (Q[3]/Q[1])^2)

exact_U = @. exact_sol_viscous_shocktube.(x,T)
exact_rho = [x[1] for x in exact_U]
exact_rhou = [x[2] for x in exact_U]
exact_rhov = [x[3] for x in exact_U]
exact_E = [x[4] for x in exact_U]

rho = Q[1]
rhou = Q[2]
rhov = Q[3]
E = Q[4]

exact_rho_p = Vp*exact_rho

# # ignore right boundary TODO: how to enforce BC correctly?
# domain = findall(x .<= 0.5)
# exact_rho,exact_rhou,exact_rhov,exact_E,rho,rhou,rhov,E = (x->x[domain]).((exact_rho,exact_rhou,exact_rhov,exact_E,rho,rhou,rhov,E))

Linferr = maximum(abs.(exact_rho-rho))/maximum(abs.(rho)) +
          maximum(abs.(exact_rhou-rhou))/maximum(abs.(rhou)) +
          #maximum(abs.(exact_rhov-rhov))/maximum(abs.(rhov).+1e-14) +
          maximum(abs.(exact_E-E))/maximum(abs.(E))

J = J[1] # TODO: assume uniform mesh
L1err = sum(J.*abs.(exact_rho-rho))/sum(J.*abs.(rho)) +
        sum(J.*abs.(exact_rhou-rhou))/sum(J.*abs.(rhou)) +
        #sum(J.*abs.(exact_rhov-rhov))/sum(J.*abs.(rhov).+1e-14) +
        sum(J.*abs.(exact_E-E))/sum(J.*abs.(E))
println("N = $N, K = $K")
println("L1 error is $L1err")
println("Linf error is $Linferr")

xp = Vp*x
yp = Vp*y
vv = Vp*Q[1]
scatter(xp,yp,vv,zcolor=vv,camera=(0,90),colorbar=:right)
#scatter(xp,yp,exact_rho_p,zcolor=exact_rho_p,camera=(0,90),colorbar=:right)

max_stars_repo_name: jlchan/ESDG-CNS

max_issues_repo_name: jlchan/ESDG-CNS

max_forks_repo_name: jlchan/ESDG-CNS


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

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

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

module MapReduceTest
using ParallelAccelerator

#ParallelAccelerator.DomainIR.set_debug_level(3)
#ParallelAccelerator.ParallelIR.set_debug_level(3)
#ParallelAccelerator.CGen.set_debug_level(3)
#ParallelAccelerator.set_debug_level(3)

@acc function map_1(A)
    a = map(x -> x + 1, A)
    b = map(x -> x - 1, A)
    c = a .* b
    return c
end

@acc function map_2(A, B)
    map((x, y) -> x * y, A, B)
end

@acc function map_3(A)
    map!(x -> x - 1, A)
    return A
end

@acc function map_4(A, B)
    map!(x -> x - 1, A, B)
    sum(A.+B)
end

@acc function reduce_1(A)
    reduce((x,y) -> x + y, 0, A)
end

@acc function reduce_2(A)
    reduce((x,y) -> x + y, 0, map(x -> x * x, A))
end

@acc function sum_A(col)
    A = ones(Int, col, col, col)
    sum(A)
end

@acc function sum_A_1(col)
    A = ones(Int, col, col, col)
    B = sum(A, 1)
    C = sum(B, 2)
end

@acc function sum_A_cond_1(col)
    A = ones(Int, col)
    A[1] = -1
    sum(A[A .> 0])
end

@acc function sum_A_cond_2(col)
    A = ones(Int, col, col)
    A[1,1] = -1
    B = [ i for i = 1:col ]
    C = B .> 2
    cartesianarray(Tuple{Int}, (col,)) do j
        sum(A[C, j:j])
    end
end

@acc function sum_A_range_1(col)
    A = ones(Int, col, col)
    A[1,2] = 0
    sum(A[:, 2])
end

function test1()
    A = Int[x for x = 1:10]
    map_1(A) == @noacc map_1(A)
end

function test2()
    A = Int[x for x = 1:10]
    map_2(A, A) == @noacc map_2(A, A)
end

function test3()
    A = Int[x for x = 1:10]
    map_3(A)
    A[1] == 0 && A[10] == 9
end

function test4()
    A = Int[0 for x = 1:10]
    B = Int[x for x = 1:10]
    map_4(A, B) == @noacc map_4(A, B)
end

function test5()
    A = Int[x for x = 1:10]
    reduce_1(A) == sum(A)
end

function test6()
    A = Int[x for x = 1:10]
    reduce_2(A) == sum(abs2, A)
end

function test7()
    return sum_A(5) == @noacc sum_A(5)
end

function test8()
    return sum_A_1(5) == @noacc sum_A_1(5)
end

function test9()
    return sum_A_cond_1(5) == @noacc sum_A_cond_1(5)
end

function test10()
    return sum_A_cond_2(5) == @noacc sum_A_cond_2(5)
end

function test11()
    return sum_A_range_1(5) == @noacc sum_A_range_1(5)
end

end

using Base.Test
using MapReduceTest
println("Testing map and reduce...")

@test MapReduceTest.test1()
@test MapReduceTest.test2()
@test MapReduceTest.test3()
@test MapReduceTest.test4()
@test MapReduceTest.test5()
@test MapReduceTest.test6()
@test MapReduceTest.test7()
@test MapReduceTest.test8()
@test MapReduceTest.test9()
@test MapReduceTest.test10()
@test MapReduceTest.test11()

println("Done testing map and reduce.")


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

"""
    hull(::Type{T}, a::ValidInterval)

    Returns the hull of the interval `a` with boundary type `T`.

    Note that hull() is the fundamental conversion operation in IEEE 1788-2015.
"""
hull(::Type{T}, a::ValidInterval{T}) where T = a

hull(::Type{T}, x::Interval) where T <: AbstractFloat =
    unsafe_interval( T(x.lo, RoundDown), T(x.hi, RoundUp) )
hull(::Type{T}, ::EmptyInterval) where T = emptyinterval(T)

promote_rule(::Type{Interval{T}}, ::Type{Interval{S}}) where {T, S} =
    Interval{promote_type(T, S)}

promote_rule(::Type{Interval{T}}, ::Type{S}) where {T, S} =
    Interval{promote_type(T, S)}

# convert numbers to intervals
convert(::Type{Interval{T}}, x::T) where T = interval(x)
convert(::Type{Interval}, x::T) where T = convert(Interval{float(T)}, x)

# convert between interval subtypes
convert(::Type{Interval{T}}, x::Interval) where T = hull(T, x)
convert(::Type{Interval}, x::Interval) = x

real(a::ValidInterval) = a
float(x::ValidInterval) = x
big(x::ValidInterval) = hull(BigFloat, x)

"""
    unsafe_atomic(::Type{<:AbstractFloat}, x)

Construct the tightest interval of a given type that contains the value `x`.

# Examples

Construct an `Interval{Float64}` containing a given `BigFloat`:

```julia
julia> x = big"0.1"
1.000000000000000000000000000000000000000000000000000000000000000000000000000002e-01

julia> i = IntervalArithmetic.unsafe_atomic(Float64, x)
[0.0999999, 0.100001]

julia> i isa Interval{Float64}
true

julia> i.lo <= x <= i.hi
true

julia> i.hi === nextfloat(i.lo)
true
```

Construct an `Interval{Float32}` containing a the real number 0.1 in two ways:

```julia
julia> i2 = IntervalArithmetic.unsafe_atomic(Float32, 1//10)
[0.0999999, 0.100001]

julia> i1 === i2
true

julia> i.lo <= 1//10 <= i.hi
true
```
"""
function unsafe_atomic end

unsafe_atomic(::Type{T}, x::T) where T <: AbstractFloat =
    isfinite(x) ? unsafe_interval(x, x) : signed_infinity(x)

unsafe_atomic(::Type{T}, x::S) where {T<:AbstractFloat, S<:Real} =
    isfinite(x) ?
        unsafe_interval(T(x, RoundDown), T(x, RoundUp)) :
        signed_infinity(T(x))


text:
mutable struct ModelBlock

    # TempField means the variable will not be in the
    # snapshot. It can still be output in the hist file
    # but only for debugging usage.

    ev   :: Env
    fi   :: Field
    tmpfi :: TempField
    co   :: Union{Core, Nothing}
    dt   :: Union{DataTable, Nothing}

    function ModelBlock(
        ev :: Env;
        init_core :: Bool = false,
    )

        fi = Field(ev)
        tmpfi = TempField(ev)

        mb = new(
            ev,
            fi,
            tmpfi,
            nothing,
            nothing,
        )

        co = (init_core) ? Core(ev, tmpfi) : nothing
        dt = DataTable(
            Nz = ev.Nz,
            Nx = ev.Nx,
            Ny = ev.Ny,
            mask_sT=ev.topo.sfcmask_sT,
            mask_T=ev.topo.mask_T,
        )

        mb.co = co
        mb.dt = dt

        for (k, varinfo) in getDynamicVariableList(mb; varsets=["ALL",])
            writeLog("Register variable: {:s}", string(k))
            varref, grid, mask = varinfo
            regVariable!(dt, k, grid, mask, varref)
        end

        for (k, varinfo) in getCompleteVariableList(mb, :static)
            writeLog("Register variable: {:s}", string(k))
            varref, grid, mask = varinfo
            regVariable!(dt, k, grid, mask, varref)
        end

        return mb
    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
module ResortLantern

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Resort Lantern" => Ahorn.EntityPlacement(
        Maple.ResortLantern
    )
)

holderSprite = "objects/resortLantern/holder.png"
lanternSprite = "objects/resortLantern/lantern00.png"

function isLeftSide(entity::Maple.ResortLantern, room::Maple.Room)
    x, y = Ahorn.position(entity) .+ (0, 2)
    tx, ty = floor(Int, x / 8), floor(Int, y / 8)

    return get(room.fgTiles.data, (ty, tx + 2), '0') == '0'
end

function Ahorn.selection(entity::Maple.ResortLantern, room::Maple.Room)
    x, y = Ahorn.position(entity)

    leftSide = isLeftSide(entity, room)

    return Ahorn.coverRectangles([
        Ahorn.getSpriteRectangle(holderSprite, x, y, sx=leftSide ? 1 : -1),
        Ahorn.getSpriteRectangle(lanternSprite, x, y)
    ])
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.ResortLantern, room::Maple.Room)
    leftSide = isLeftSide(entity, room)

    Ahorn.drawSprite(ctx, holderSprite, 0, 0, sx=leftSide ? 1 : -1)
    Ahorn.drawSprite(ctx, lanternSprite, 0, 0)
end

end


text:
mutable struct Cons{T}
    car::T
    cdr::Any
end

struct Nil{T} end

const List{T} = Union{Cons{T},Nil{T}}

listof(::Type{T}) where {T} = Nil{T}()
listof(::Type{T}, x::T, xs::T...) where {T} = Cons{T}(x, listof(T, xs...))
Base.empty(::List{T}) where {T} = Nil{T}()

car(x::Cons) = x.car
cdr(x::Cons{T}) where {T} = x.cdr::Union{Cons{T},Nil{T}}

Base.eltype(::Type{List{T}}) where {T} = T
Base.IteratorSize(::Type{<:List}) = Base.SizeUnknown()
Base.iterate(::Nil) = nothing
Base.iterate(lst::Cons{T}, x = lst) where {T} =
    if x isa Cons{T}
        car(x), cdr(x)
    else
        nothing
    end

mutable struct LockedLinkedList{T}
    list::List{T}
    lock::ReentrantLock
    isempty::Threads.Atomic{Bool}
end

LockedLinkedList{T}() where {T} =
    LockedLinkedList{T}(Nil{T}(), ReentrantLock(), Threads.Atomic{Bool}(true))

setlist(l::LockedLinkedList{T}, list::List{T}) where {T} = @set l.list = list

function locked(f, l::LockedLinkedList)
    lock(l.lock) do
        f(l.list)
    end
end

function Base.empty!(l::LockedLinkedList)
    lock(l.lock) do
        l.list = empty(l.list)
        l.isempty[] = true
    end
end

Base.isempty(l::LockedLinkedList) = l.isempty[]

function setcdr!(lst::Cons{T}, cdr::List{T}) where {T}
    lst.cdr = cdr
    return lst
end

function setcdr!(l::LockedLinkedList{T}, cdr::List{T}) where {T}
    lock(l.lock) do
        if l.list isa Nil
            l.list = cdr
            if !(cdr isa Nil)
                l.isempty[] = false
            end
        else
            setcdr!(l.list, cdr)
            l.isempty[] = false
        end
    end
    return
end

function trypopfirst!(l::LockedLinkedList)
    # isempty(l) && return nothing  # TODO: is it OK?
    lock(l.lock) do
        list = l.list
        if list isa Nil
            return nothing
        else
            l.list = newlist = cdr(list)
            if newlist isa Nil
                l.isempty[] = true
            end
            return Some(car(list))
        end
    end
end


text:
function verified_prediction_lags(lags)
    ηs = sort(lags[lags .!= 0])

    if length(ηs) % 2 != 0
        throw(ArgumentError("Need an even number of lags (as many negative as positive prediction lags)"))
    end

    ηs_neg = ηs[ηs .< 0]
    ηs_pos = ηs[ηs .> 0]

    if length(ηs_neg) != length(ηs_pos)
        throw(ArgumentError("There must be as many negative as positive prediction lags. Got $ηs_neg and $ηs_pos"))
    end

    ηs_neg_sort = sort(abs.(ηs_neg))
    ηs_pos_sort = sort(ηs_pos)
    if !(all(ηs_neg_sort .== ηs_pos_sort))
        throw(ArgumentError("Negative amd positive prediction lags must be symmetric. Got $ηs_neg and $ηs_pos"))
    end

    return ηs
end

function verified_prediction_lags(lag::Int)
    # If only one η is provided (either positive or negative), just make a prediction
    # for that lag.
    ηs = [-abs(lag), abs(lag)]
end

"""
    PredictiveAsymmetryTest(predictive_test::CausalityTest)

The parameters for a predictive asymmetry causality test [1].

## Mandatory keywords

- **`predictive_test`**: An instance of a predictive causality test that explicitly
    uses prediction lags (e.g. [`VisitationFrequencyTest`](@ref) or
    [`TransferOperatorGridTest`](@ref)).

## About the prediction lags

The prediction lags in the predictive causality test must consist of `n` negative
integers and `n` positive integers that are symmetric around zero.

In other words, negative lags  must exactly match the positive lags but with opposite
sign. The zero lag can be included, but will be ignored, so it is possible to give
ranges too.

## Examples

```julia
test_visitfreq = VisitationFrequencyTest(ηs = [-5, -4, -2, -1, 0, 1, 2, 4, 5])
test_transferoperator = TransferOperatorGridTest(ηs = -3:3)

# Note that `predictive_test` is a *mandatory* keyword.
PredictiveAsymmetryTest(predictive_test = test_visitfreq)
PredictiveAsymmetryTest(predictive_test = test_transferoperator)
```

## References

1. Diego, David, Kristian Agasøster Haaga, Jo Brendryen, and Bjarte Hannisdal.
    A simple test for causal asymmetry in complex systems. In prep.
"""
Base.@kwdef struct PredictiveAsymmetryTest{TEST, N} <: CausalityTest where TEST
    predictive_test::TEST

    function PredictiveAsymmetryTest(test::T) where {T <: TransferEntropyCausalityTest}
        # Check that prediction lags are okay
        verified_prediction_lags(test.ηs)
        N = length(test.ηs[test.ηs .> 0])
        return new{T, N}(test)
    end
end

function causalbalance(lags, tes)
    posinds = findall(lags .> 0)
    neginds = findall(lags .< 0)
    sum(tes[posinds]) - sum(tes[neginds])
end

function update_ηs(test::VisitationFrequencyTest)

    VisitationFrequencyTest(
        k = test.k, l = test.l, m = test.m, n = test.n, b = test.b,
        τ = test.τ,
        binning_summary_statistic = test.binning_summary_statistic,
        binning = test.binning,
        ηs = verified_prediction_lags(test.ηs))
end

function update_ηs(test::TransferOperatorGridTest)

    TransferOperatorGridTest(
        k = test.k, l = test.l, m = test.m, n = test.n, b = test.b,
        τ = test.τ,
        binning_summary_statistic = test.binning_summary_statistic,
        binning = test.binning,
        ηs = verified_prediction_lags(test.ηs))
end

return_predictive_asymmetry(ηs, As) = As
return_predictive_asymmetry(η::Int, As) = As[1]

function predictive_asymmetry(source, target,
        p::PredictiveAsymmetryTest{T, N}) where {T <: TransferEntropyCausalityTest, N}

    # Update the test parameters so that we have symmetric prediction lags
    test = update_ηs(p.predictive_test)

    # Predictions from source to target.
    preds = causality(source, target, test)

    # The number of predictive asymmetries is half the number of prediction lags,
    # which is encoded in the type parameter `N`
    ηs = test.ηs
    As = zeros(Float64, N)

    for (i, η) in enumerate(ηs[ηs .> 0])
        lag_idxs = findfirst(ηs .== -η):findfirst(ηs .== η)
        As[i] = causalbalance(ηs[lag_idxs], preds[lag_idxs])
    end

    return return_predictive_asymmetry(p.predictive_test.ηs, As)
end

function causality(source::AbstractVector{T}, target::AbstractVector{T}, p::PredictiveAsymmetryTest{CT}) where {T<:Real, CT}
    predictive_asymmetry(source, target, p)
end

# there is no need to define custom causality method for a uncertain data
# as long as the predictive tests that it uses supports them

export
    PredictiveAsymmetryTest,
    predictive_asymmetry,
    update_ηs,
    verified_prediction_lags


text:
######
# This file is part of the BivariateCopulas.jl package.
#
#   Module definition
#
#           University of Liverpool, Institute for Risk and Uncertainty
#
#                                           Authors: Ander Gray & Liam T. Henry
#                                           Email:  ander.gray@liverpool.ac.uk
######

module BivariateCopulas

using Reexport
using Requires

@reexport using Distributions             #   https://github.com/JuliaStats/Distributions.jl
using LinearAlgebra                       #   Main Library
using Interpolations                      #   https://github.com/JuliaMath/Interpolations.jl

import Distributions: cdf, quantile, rand

import Base: rand

global n = 200                       # Number of descretizations
global pn = 50                       # Number of descretization for plotting
global bOt = 0.0001                    # smallest quantile to use if left tail is unbounded. Should we use these to define the copula ranges?
global tOp = 0.9999                    # largest quantile to use if right tail is unbounded
#global bOt = 0;
#global tOp = 1;

abstract type AbstractCopula <: Real end
abstract type AbstractJoint <: Real end
abstract type AbstractMarginal <: ContinuousUnivariateDistribution end

include("archimedean/archimedean.jl")
include("copulas.jl")
include("NormalDistribution.jl")

function __init__()
    @require PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" begin
        @require PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" begin
            using PyCall
            using PyPlot
            import PyPlot: plot
            include("plotting.jl")
            return using3D()
        end
    end
end

export copula,
    density,
    cdf,
    sample,
    rand,
    CholeskyGaussian,
    conditional,
    calcCdf,
    calcDensity,
    M,
    W,
    Pi,
    Frank,
    Clayton,
    Gaussian,
    norm_cdf,
    bivariate_cdf,
    mvNormCdf,
    Joint,
    invCdf,
    marginal,
    quantile,
    ρCopula,
    SpearmanCopula,
    KendalCopula,
    τCopula,
    samplePlot,
    plotDensity,
    plotDen,
    plotCdf,
    plot,
    plotContourCdf,
    plotContourDen,
    scatter,
    rotate90,
    rotate180,
    rotate270,
    ArchimedeanCopula,
    Clayton,
    Independence,
    rosenblatt,
    inverse_rosenblatt,
    τ

end #BivariateCopulas


text:
module JuliaPythonSnippets

export @pyexample_str, @pyshow_str, @named

using PyCall
using Markdown

macro pyexample_str(str)
    code = """
    ```python
    $str
    ```
    """
    quote
        PyCall.@py_str $str
        Markdown.parse($code)
    end
end

macro pyshow_str(str)
    quote
        pyshow_impl(PyCall.@py_str $str 'o')
    end
end

function pyshow_impl(obj::PyObject)
    str = pybuiltin("repr")(obj)
    return Markdown.parse("""
    ```
    $str
    ```
    """)
end

const _scopes = Dict{Symbol,Module}()

macro named(name::Symbol, block)
    scope = get!(() -> Module(name), _scopes, name)
    Base.eval(scope, :(using JuliaPythonSnippets))
    Base.eval(scope, :(using PyCall))
    QuoteNode(Base.eval(scope, block))
end

end # module

max_stars_repo_name: tkf/julia-python-snippets

max_issues_repo_name: tkf/julia-python-snippets

max_forks_repo_name: tkf/julia-python-snippets


text:
module globals

# Free space impedance
imp0 = 377.0

# Courant number
cdtds = 1.0

end


text:
# import some useful things

import Base: +, -, *, <, ==
import Base: copy, show
# provides polynomial ring and polynomials interface
using AbstractAlgebra
# provides dot product operation
import LinearAlgebra: dot

#-----------------------------------------------------------------------------
# Signature2 stuff

struct Signature2{Poly}
    index::Int
    monom::Poly
end

function <ₛ(s1::Signature2, s2::Signature2)
    s1.index < s2.index || (s1.index == s2.index && s1.monom < s2.monom)
end

==(a::Signature2, b::Signature2) = a.index == b.index && a.monom == b.monom

<(a::Signature2, b::Signature2) = a <ₛ b
Base.isless(a::Signature2, b::Signature2) = a < b

#-----------------------------------------------------------------------------

# our module element will store one thing: a tuple of polynomials
mutable struct LabeledPolynomial{Poly}
    ev::Poly
    sgn::Signature2{Poly}
end

function normalize(f::LabeledPolynomial)
    lead = leading_coefficient(f.ev)
    LabeledPolynomial(map_coefficients(c -> c // lead, f.ev), f.sgn)
end

function normalize(f)
    lead = leading_coefficient(f)
    map_coefficients(c -> c // lead, f)
end

function normalize!(f::LabeledPolynomial)
    lead = leading_coefficient(f.ev)
    f.ev = map_coefficients(c -> c // lead, f.ev)
    f
end

# some magic to prin2t ModuleElement nicely
Base.show(io::IO, ::MIME"text/plain", m::LabeledPolynomial) = println(io, "($(m.ev), $(m.sgn))")
Base.show(io::IO, m::LabeledPolynomial) = print(io, "($(m.ev), $(m.sgn))")

Base.show(io::IO, m::Signature2) = print(io, "<$(m.index), $(m.monom)>")
Base.show(io::IO, ::MIME"text/plain", m::Signature2) = println(io, "<$(m.index), $(m.monom)>")

#-----------------------------------------------------------------------------
# Arithmetic operations

*(t::Poly, s::Signature2{Poly}) where {Poly} = Signature2(s.index, s.monom*t)
*(s::Signature2{Poly}, t::Poly) where {Poly} = t * s


text:
## INSPECTING LEARNING NETWORKS

"""
    tree(N)

Return a named-tuple respresentation of the ancestor tree `N`
(including training edges)

"""
function tree(W::Node)
    mach = W.machine
    if mach === nothing
        value2 = nothing
        endkeys = []
        endvalues = []
    else
        value2 = mach.model
        endkeys = (Symbol("train_arg", i) for i in eachindex(mach.args))
        endvalues = (tree(arg) for arg in mach.args)
    end
    keys = tuple(:operation,  :model,
                 (Symbol("arg", i) for i in eachindex(W.args))...,
                 endkeys...)
    values = tuple(W.operation, value2,
                   (tree(arg) for arg in W.args)...,
                   endvalues...)
    return NamedTuple{keys}(values)
end
tree(s::Source) = (source = s,)

# """
#    args(tree; train=false)

# Return a vector of the top level args of the tree associated with a node.
# If `train=true`, return the `train_args`.
# """
# function args(tree; train=false)
#     keys_ = filter(keys(tree) |> collect) do key
#         match(Regex("^$("train_"^train)arg[0-9]*"), string(key)) !== nothing
#     end
#     return [getproperty(tree, key) for key in keys_]
# end

"""
    models(N::AbstractNode)

A vector of all models referenced by a node `N`, each model appearing
exactly once.

"""
function models(W::AbstractNode)
    models_ = filter(flat_values(tree(W)) |> collect) do model
        model isa Model
    end
    return unique(models_)
end

"""
    sources(N::AbstractNode)

A vector of all sources referenced by calls `N()` and `fit!(N)`. These
are the sources of the ancestor graph of `N` when including training
edges.

Not to be confused with `origins(N)`, in which training edges are
excluded.

See also: [`origins`](@ref), [`source`](@ref).
"""
function sources(W::AbstractNode; kind=:any)
    if kind == :any
        sources_ = filter(flat_values(tree(W)) |> collect) do value
            value isa Source
        end
    else
        sources_ = filter(flat_values(tree(W)) |> collect) do value
            value isa Source && value.kind == kind
        end
    end
    return unique(sources_)
end

"""
    machines(N::AbstractNode [, model::Model])

List all machines in the ancestor graph of node `N`, optionally
restricting to those machines whose corresponding model matches the
specifed `model`.

Here two models *match* if they have the same, possibly nested
hyperparameters, or, more precisely, if
`MLJModelInterface.is_same_except(m1, m2)` is `true`.

See also `MLJModelInterface.is_same_except`.

"""
function machines(W::Node, model=nothing)
    if W.machine === nothing
        machs = vcat((machines(arg) for arg in W.args)...) |> unique
    else
        machs = vcat(Machine[W.machine, ],
                    (machines(arg) for arg in W.args)...,
                    (machines(arg) for arg in W.machine.args)...) |> unique
    end
    model == nothing && return machs
    return filter(machs) do mach
        mach.model == model
    end
end

args(::Source) = []
args(N::Node) = N.args
train_args(::Source) = []
train_args(N::Node{<:Machine}) = N.machine.args
train_args(N::Node{Nothing}) = []

"""
    children(N::AbstractNode, y::AbstractNode)

List all (immediate) children of node `N` in the ancestor graph of `y`
(training edges included).

"""
children(N::AbstractNode, y::AbstractNode) = filter(nodes(y)) do W
    N in args(W) || N in train_args(W)
end |> unique

"""
    lower_bound(type_itr)

Return the minimum type in the collection `type_itr` if one exists
(mininum in the sense of `<:`). If `type_itr` is empty, return `Any`,
and in all other cases return the universal lower bound `Union{}`.

"""
function lower_bound(Ts)
    isempty(Ts) && return Any
    sorted = sort(collect(Ts), lt=<:)
    candidate = first(sorted)
    all(T -> candidate <: T, sorted[2:end]) && return candidate
    return Union{}
end

function _lower_bound(Ts)
    Unknown in Ts && return Unknown
    return lower_bound(Ts)
end

MLJModelInterface.input_scitype(N::Node) = Unknown
MLJModelInterface.input_scitype(N::Node{<:Machine}) =
    input_scitype(N.machine.model)

## MANIPULATING LEARNING NETWORKS

"""
    reset!(N::Node)

Place the learning network terminating at node `N` into a state in
which `fit!(N)` will retrain from scratch all machines in its
dependency tape. Does not actually train any machine or alter
fit-results. (The method simply resets `m.state` to zero, for every
machine `m` in the network.)

"""
function reset!(W::Node)
    for mach in machines(W)
        mach.state = 0 # to do: replace with dagger object
    end
end


text:
# Example how to use the residual block
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using Test, LinearAlgebra, InvertibleNetworks

# Input
nx = 64
ny = 64
nz = 64
n_in = 10
n_hidden = 20
batchsize = 2

k1 = 3; p1 = 1; s1 = 1  # standard res-net block
k2 = 3; p2 = 1; s2 = 1

# Input image
X = glorot_uniform(nx, ny, nz, n_in, batchsize)

# Residual blocks
RB = ResidualBlock3D(n_in, n_hidden; k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2)
RB0 = ResidualBlock3D(n_in, n_hidden; k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2)

# Observed data
Y = RB.forward(X)

# Predicted data
Y0 = RB0.forward(X)
@test isnothing(RB0.W1.grad) # after forward pass, gradients are not set

# Residual
ΔY = Y0 - Y

# Backward pass: need to pass data residual and original input X (as layer is not invertible)
ΔX = RB0.backward(ΔY, X)   # returns derivative w.r.t input
@test ~isnothing(RB0.W1.grad)    # after inverse pass, gradients are set
@test ~isnothing(RB0.b1.grad)


text:
using DelimitedFiles
using MultivariateStats
import Statistics: mean, median

writecsv(filename::AbstractString, data) = writedlm(filename, data, ',')
readcsv(filename::AbstractString) = readdlm(filename, ',')
readcsv(filename::AbstractString, ::Type{T}) where {T} = readdlm(filename, ',', T)

input = ARGS[1]
output1 = ARGS[2]
output2 = ARGS[3]
dim = ARGS[4]

# input = "/data2/koki/ICCIPCA/Data/Cortical_SMART/Data.csv"
# output1 = "/data2/koki/ICCIPCA/Data/Cortical_SMART/MultivariateStats/Eigen_vectors.csv"
# output2 = "/data2/koki/ICCIPCA/Data/Cortical_SMART/MultivariateStats/Eigen_values.csv"
# dim = "10"

dim = parse(Int64, dim)
println(input)
println(output1)
println(output2)
println(dim)

# Data loading
out = readcsv(input)
# Library size
libsize = sum(out, dims=1)
# CPMED
out = median(libsize) .* (out ./ libsize)
# Log transformation
out = log10.(out .+ 1)
# PCA
mval = vec(mean(out, dims=2))
res = MultivariateStats.fit(PCA, out, mean=mval, maxoutdim=dim)
# Eigenvectors
writecsv(output1, MultivariateStats.transform(res, out)')
# Eigenvalues
writecsv(output2, principalvars(res))

max_stars_repo_name: rikenbit/onlinePCA-experiments

max_issues_repo_name: rikenbit/onlinePCA-experiments

max_forks_repo_name: rikenbit/onlinePCA-experiments


text:
using Pumas, Plots, Random
Random.seed!(0);

repeated_dose_regimen = DosageRegimen(100, time=0, ii=4, addl=2)

choose_covariates() = (Wt = rand(55:80),)

pop = Population(map(i -> Subject(id = i,
                                  events = repeated_dose_regimen,
                                  observations =(dv=Float64[],),
                                  covariates = choose_covariates()),
                                  1:24))

mymodel = @model begin
  @param   begin
    cl ∈ RealDomain(lower = 0.0, init = 1.0)
    tv ∈ RealDomain(lower = 0.0, init = 10.0)
    ka ∈ RealDomain(lower = 0.0, init = 1.0)
    q  ∈ RealDomain(lower = 0.0, init = 0.5)
    Ω  ∈ PDiagDomain(init = [0.9,0.07, 0.05])
    σ_prop ∈ RealDomain(lower = 0,init = 0.03)
  end

  @random begin
    η ~ MvNormal(Ω)
  end

  @covariates Wt

  @pre begin
    CL = cl * (Wt/70)^0.75 * exp(η[1])
    Vc = tv * (Wt/70) * exp(η[2])
    Ka = ka * exp(η[3])
    Vp = 30.0
    Q  = q
  end

  @dynamics Depots1Central1Periph1

  @derived begin
      cp := @. 1000*(Central / Vc) # We use := because we don't want simobs to store the variable
      dv ~ @. Normal(cp, abs(cp)*σ_prop)
    end
end

param = init_param(mymodel)
obs = simobs(mymodel, pop, param, obstimes=1:1:72)
plot(obs)

result = fit(mymodel, Subject.(obs), param, Pumas.FOCEI())

alternative_param = (
    cl = 0.5,
    tv = 9.0,
    ka = 1.3,
    q  = 0.3,
    Ω  = Diagonal([0.18,0.04, 0.03]),
    σ_prop = 0.04)

fit(mymodel, Subject.(obs), alternative_param, Pumas.FOCEI())

infer(result)

pop_big = Population(map(i -> Subject(id = i,
                                  events = repeated_dose_regimen,
                                  observations =(dv=Float64[],),
                                  covariates = choose_covariates()),
                                  1:100))
obs_big = simobs(mymodel, pop_big, param, obstimes=1:1:72)
result_big = fit(mymodel, Subject.(obs_big), param, Pumas.FOCEI())
infer(result_big)

mymodel_misspec = @model begin
  @param   begin
    cl ∈ RealDomain(lower = 0.0, init = 1.0)
    tv ∈ RealDomain(lower = 0.0, init = 20.0)
    ka ∈ RealDomain(lower = 0.0, init = 1.0)
    Ω  ∈ PDiagDomain(init = [0.12, 0.05, 0.08])
    σ_prop ∈ RealDomain(lower = 0, init = 0.03)
  end

  @random begin
    η ~ MvNormal(Ω)
  end

  @pre begin
    CL = cl * (Wt/70)^0.75 * exp(η[1])
    Vc = tv * (Wt/70) * exp(η[2])
    Ka = ka * exp(η[3])
  end
  @covariates Wt

  @dynamics Depots1Central1

  @derived begin
      cp = @. 1000*(Central / Vc)
      dv ~ @. Normal(cp, abs(cp)*σ_prop)
    end
end

alternative_param_no_q = (
    cl = 0.5,
    tv = 9.0,
    ka = 1.3,
    Ω  = Diagonal([0.18,0.04, 0.03]),
    σ_prop = 0.04)

result_misspec = fit(mymodel_misspec, Subject.(obs), alternative_param_no_q, Pumas.FOCEI())

wres = wresiduals(result)
wres_misspec = wresiduals(result_misspec)
p1 = plot([w.wres.dv for w in wres], title="Correctly specified", legend=false)
p2 = plot([w.wres.dv for w in wres_misspec], title = "Misspecified", legend=false)
plot(p1, p2)


text:
using Test
using ParallelStencil
import ParallelStencil: @reset_parallel_stencil, @is_initialized, SUPPORTED_PACKAGES, PKG_CUDA, PKG_THREADS, INDICES
import ParallelStencil: @require, longnameof, @prettyexpand, prettystring
import ParallelStencil: checkargs_parallel, validate_body, parallel
using ParallelStencil.Exceptions
using ParallelStencil.FiniteDifferences3D
ix, iy, iz = INDICES[1], INDICES[2], INDICES[3]
TEST_PACKAGES = SUPPORTED_PACKAGES
@static if PKG_CUDA in TEST_PACKAGES
    import CUDA
    if !CUDA.functional() TEST_PACKAGES = filter!(x->x≠PKG_CUDA, TEST_PACKAGES) end
end

@static for package in TEST_PACKAGES  eval(:(
    @testset "$(basename(@__FILE__)) (package: $(nameof($package)))" begin
        @testset "1. parallel macros" begin
            @init_parallel_stencil($package, Float64, 3)
            @require @is_initialized()
            @testset "@parallel kernelcall" begin # NOTE: calls must go to ParallelStencil.ParallelKernel.parallel and must therefore give the same result as in ParallelKernel (tests copied from there).
                @static if $package == PKG_CUDA
                    call = prettystring(parallel(:(f(A))))
                    @test occursin("CUDA.@cuda blocks = ParallelStencil.ParallelKernel.compute_nblocks(length.(ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A))), ParallelStencil.ParallelKernel.compute_nthreads(length.(ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A))))) threads = ParallelStencil.ParallelKernel.compute_nthreads(length.(ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A)))) f(A, ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A)))", call)
                    @test occursin("CUDA.synchronize()", call)
                    call = prettystring(parallel(:ranges, :(f(A))))
                    @test occursin("CUDA.@cuda blocks = ParallelStencil.ParallelKernel.compute_nblocks(length.(ParallelStencil.ParallelKernel.promote_ranges(ranges)), ParallelStencil.ParallelKernel.compute_nthreads(length.(ParallelStencil.ParallelKernel.promote_ranges(ranges)))) threads = ParallelStencil.ParallelKernel.compute_nthreads(length.(ParallelStencil.ParallelKernel.promote_ranges(ranges))) f(A, ParallelStencil.ParallelKernel.promote_ranges(ranges))", call)
                    call = prettystring(parallel(:nblocks, :nthreads, :(f(A))))
                    @test occursin("CUDA.@cuda blocks = nblocks threads = nthreads f(A, ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.compute_ranges(nblocks .* nthreads)))", call)
                    call = prettystring(parallel(:ranges, :nblocks, :nthreads, :(f(A))))
                    @test occursin("CUDA.@cuda blocks = nblocks threads = nthreads f(A, ParallelStencil.ParallelKernel.promote_ranges(ranges))", call)
                    call = prettystring(parallel(:nblocks, :nthreads, :(stream=mystream), :(f(A))))
                    @test occursin("CUDA.@cuda blocks = nblocks threads = nthreads stream = mystream f(A, ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.compute_ranges(nblocks .* nthreads)))", call)
                elseif $package == PKG_THREADS
                    @test prettystring(parallel(:(f(A)))) == "f(A, ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A)))"
                    @test prettystring(parallel(:ranges, :(f(A)))) == "f(A, ParallelStencil.ParallelKernel.promote_ranges(ranges))"
                    @test prettystring(parallel(:nblocks, :nthreads, :(f(A)))) == "f(A, ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.compute_ranges(nblocks .* nthreads)))"
                    @test prettystring(parallel(:ranges, :nblocks, :nthreads, :(f(A)))) == "f(A, ParallelStencil.ParallelKernel.promote_ranges(ranges))"
                    @test prettystring(parallel(:(stream=mystream), :(f(A)))) == "f(A, ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A)))"
                end;
            end;
            @testset "@parallel kernel" begin
                @require @is_initialized()
                A  = @zeros(4, 5, 6)
                @parallel function write_indices!(A)
                    @all(A) = $ix + ($iy-1)*size(A,1) + ($iz-1)*size(A,1)*size(A,2); # NOTE: $ix, $iy, $iz come from ParallelStencil.INDICES.
                    return
                end
                @parallel write_indices!(A);
                @test all(Array(A) .== [ix + (iy-1)*size(A,1) + (iz-1)*size(A,1)*size(A,2) for ix=1:size(A,1), iy=1:size(A,2), iz=1:size(A,3)])
            end;
            @testset "apply masks" begin
                @require @is_initialized()
                expansion = string(@prettyexpand 1 @parallel sum!(A, B) = (@all(A) = @all(A) + @all(B); return))
                @test occursin("if @within(\"@all\", A)", expansion)
                @test string(@prettyexpand @within("@all", A)) == string(:($ix <= size(A, 1) && ($iy <= size(A, 2) && $iz <= size(A, 3))))
            end;
            @reset_parallel_stencil()
        end;
        @testset "2. Exceptions" begin
            @init_parallel_stencil($package, Float64, 3)
            @require @is_initialized
            @testset "arguments @parallel" begin
                @test_throws ArgumentError checkargs_parallel();                                                  # Error: isempty(args)
                @test_throws ArgumentError checkargs_parallel(:(f()), :(something));                              # Error: last arg is not function or a kernel call.
                @test_throws ArgumentError checkargs_parallel(:(f()=99), :(something));                           # Error: last arg is not function or a kernel call.
                @test_throws ArgumentError checkargs_parallel(:(something), :(f()=99));                           # Error: if last arg is a function, then there cannot be any other argument.
                @test_throws ArgumentError checkargs_parallel(:ranges, :nblocks, :nthreads, :something, :(f()));  # Error: if last arg is a kernel call, then we have an error if length(posargs) > 3
                @test_throws ArgumentError validate_body(:(a = b + 1))
                @test_throws ArgumentError validate_body(:(a = b + 1; @all(A) = @all(B) + 1))
                @test_throws ArgumentError validate_body(:(A = @all(B) + 1; @all(A) = @all(B) + 1))
            end;
            @reset_parallel_stencil()
        end;
    end;
)) end == nothing || true;


text:
using Test
using CompScienceMeshes

fn = joinpath(@__DIR__, "assets", "world.msh")
border = CompScienceMeshes.read_gmsh_mesh(fn, physical="WorldsEnd", dimension=1)
coast  = CompScienceMeshes.read_gmsh_mesh(fn, physical="Coast", dimension=1)
sea    = CompScienceMeshes.read_gmsh_mesh(fn, physical="Sea", dimension=2)

border_vertices = skeleton(border, 0)
coast_vertices  = skeleton(coast, 0)
sea_vertices    = skeleton(sea, 0)

interior_vertices = submesh(v -> !(v in border_vertices || v in coast_vertices), sea_vertices)

@test length(border_vertices) + length(coast_vertices) + length(interior_vertices) == length(sea_vertices)