text:
#lags a variable in the time series
#assumes the order field is intended to be in ascending order
function lagwithin!(df::DataFrame,
  targets::Vector{<:DField},
  groups::Vector{<:DField},
  period::DField;
  lags::Int=1,
  laggednames::Vector{<:DField} = (lags ≠ 1 ?
    (s::DField->Symbol(:L, lags, s)).(targets) : (s::DField->Symbol(:L, s)).(targets)))::Nothing

  T::Type = eltype(df[!,period]) #the type of the period column

  issorted(df,[groups; period]) || error("df must be sorted in lagwithin")

  Ntargets::Int = length(targets)   #pre-allocate the space
  for t ∈ 1:Ntargets
    df[!,laggednames[t]] =
      Vector{Union{eltype(df[!,targets[t]]), Missing}}(undef, size(df, 1))
  end

  for subdf ∈ groupby(df, groups)
    Nsub::Int = size(subdf, 1)

    #only run these routines if we need to
    if Nsub > lags
      for i::Int ∈ (lags+1):Nsub #iterate for all values
        for t::Int ∈ 1:Ntargets
          cur::T = subdf[i, period]

          #first see if we can lag the easy way
          if subdf[i-lags, period] == cur - lags
            subdf[i,laggednames[t]] = subdf[i-lags, targets[t]]
          elseif lags ≠ 1 #now check the hard way (finding the lagged year), pointless if lags==1
            location::NInt = findfirst(isequal(cur-lags), subdf[1:(i-1),period])

            #record the lagged value if it is available
            if !isnothing(location)
              subdf[i,laggednames[t]] = subdf[location, targets[t]]
            end
          end
        end #targets for loop
      end #periods for loop
    end
  end

  return nothing
end

#helper method to handle the case of a single target and/or group
lagwithin!(df::DataFrame,
  targets::Union{DField, Vector{<:DField}},
  groups::Union{DField, Vector{<:DField}},
  period::DField;
  lags::Int=1,
  laggedname::NDField = nothing,
  laggednames::Vector{<:DField} = (lags ≠ 1 ?
    (s::DField->Symbol(:L, lags, s)).([targets;]) : (s::DField->Symbol(:L, s)).([targets;]))
    )::Nothing = (lagwithin!(df, [targets;], [groups;], period,
      lags=lags, laggednames=[something(laggedname, laggednames);]))

#creates a differenced column
function differencewithin!(df::DataFrame,
  targets::Vector{<:DField},
  groups::Vector{<:DField},
  period::DField;
  differencednames::Vector{<:DField} = (s::DField->Symbol(:D, s)).(targets),
  createlag::Bool=true,
  deletelag::Bool=true,
  laggednames::Vector{<:DField} = (deletelag ?
    (s::DField->Symbol(:L, s, :_temp)).(targets) : (s->Symbol(:L, s)).(targets))
  )::Nothing

  createlag && lagwithin!(df, targets, groups, period,
    laggednames=laggednames)
  for t ∈ 1:length(targets)
    df[!, differencednames[t]] = df[!, targets[t]] .- df[!, laggednames[t]]
    deletelag && select!(df, Not(laggednames[t]))
  end

  return nothing
end

#helper method to handle the case of a single target and/or group
differencewithin!(df::DataFrame,
  targets::Union{DField, Vector{<:DField}},
  groups::Union{DField, Vector{<:DField}},
  period::DField;
  differencedname::NDField = nothing,
  differencednames::Vector{<:DField} = (s::DField->Symbol(:D, s)).([targets;]),
  createlag::Bool=true,
  deletelag::Bool=true,
  laggedname::NDField = nothing,
  laggednames::Vector{<:DField} = (deletelag ?
    (s::DField->Symbol(:L, s, :_temp)).([targets;]) : (s::DField->Symbol(:L, s)).([targets;]))
  )::Nothing = differencewithin!(df, [targets;], [groups;], period,
    createlag=createlag, deletelag=deletelag,
    differencednames=[something(differencedname, differencednames);],
    laggednames=[something(laggedname, laggednames);])

##################Updated versions

###if these work, consider moving to finometrics
#lag within a group
function lagwithin2sorted(vals::AbstractVector{T},
  group::AbstractVector{G}, laggedgroup::Vector{<:Union{G, Missing}}, ::Type{T}, ::Type{G};
  ) where {T<:Any, G<:Any}

  #idea here is to lag everyting, then account for the boundaries between groups
  #lag everything
  lagged::Vector{Union{T, Missing}} = [missing; vals[1:(end-1)]]

  #fixes lags across groups
  missingindicies::Vector{Bool} = ((lg,g)->
    ismissing(lg) || ismissing(g) || lg≠g).(laggedgroup,group)
  lagged[missingindicies] .= missing

  return lagged
end

#no dataframe version - just lags the vector
#not called in the dataframe versions
function lagwithin2sorted(vals::TV,group::TG) where {TV, TG}
  laggedgroup = Vector{Union{eltype(group), Missing}}([missing; group[1:(end-1)]])
  return lagwithin2sorted(vals, group, laggedgroup, eltype(TV), eltype(TG))
end

function eliminatestaledates(lagged::AbstractVector{T}, date::AbstractVector{D},
  ::Type{T}, ::Type{D};
  laggeddate::AbstractVector{Union{D,Missing}}=error("laggeddate is required"),
  maxnotstale::Any=error("maxnotstale is required"),
  ) where {T<:Any, D<:Any}

  local missingindicies::Vector{Bool}

  if all((laggeddate .≤ date) .| (laggeddate .=== missing))
    missingindicies = ((ld, d)->
      ismissing(ld) || ismissing(d) || d - maxnotstale > ld).(laggeddate, date)
  #  println(DataFrame(laggeddate=laggeddate, date=date)[1:50,:])
  # with reversed dates, we are effectually computing a lead instead of a lag
  elseif all((laggeddate .≥ date) .| (laggeddate .=== missing))
    missingindicies = ((nd, d)->
      ismissing(nd) || ismissing(d) || nd - maxnotstale > d).(laggeddate, date)
    #@assert sum(missingindicies) ≠ length(lagged)
  else
    @assert false #something is wrong (df not sorted properly?)
  end

  lagged[missingindicies] .= missing

  #@assert !(all(missing .=== lagged))

  return lagged
end

function lagwithin2sorted(
  vals::AbstractVector,
  group::AbstractVector,
  date::AbstractVector{D};
  laggedgroup::Vector{Union{G, Missing}}=error("lagged group is required with lagwithin2sorted"),
  laggeddate::Vector{Union{D,Missing}}=error(
    "lagged group is required with lagwithin2sorted if maxnotstale is provided"),
  maxnotstale::Any = error(
    "maxnotstale is required if date vec is provided in lagwithin2sorted")) where {G<:Any, D<:Any}

  local lagged = lagwithin2sorted(vals, group, laggedgroup, eltype(vals), eltype(group))
  lagged = eliminatestaledates(lagged, date, eltype(lagged), eltype(date); laggeddate, maxnotstale)
  return lagged
end

#lags multiple columns within a dataframe
function lagwithin2sorted!(df::DataFrame, vals::Vector{<:DField}, group::DField;
  date::NDField = nothing,
  maxnotstale::Any = nothing,
  laggedvals::Vector{<:DField} = (s->Symbol(:L, s)).(vals))

  local laggedgroup::Vector{Union{eltype(df[!,group]), Missing}}
  local laggeddate::Vector{Union{eltype(df[!,date]), Missing}}

  #lag the groups, and dates if needed, once
  laggedgroup = Vector{Union{eltype(group), Missing}}([missing; df[!, group][1:(end-1)]])
  #allocate the taskes for multithreading
  tasks::Vector{Task} = Vector{Task}(undef, length(vals))

  #deploy the lags
  if maxnotstale === nothing
    tasks = @sync((s::DField -> Threads.@spawn lagwithin2sorted(
      df[!, s], df[!, group], laggedgroup, eltype(df[!, s]), eltype(df[!, group]),
      )).(vals))
  else
    laggeddate = [missing; df[!,date][1:(end-1)]]
    laggeddate[laggedgroup .!== df[!,group]] .= missing
    tasks = @sync((s::DField -> Threads.@spawn lagwithin2sorted(
      df[!, s], df[!, group], df[!,date]; laggedgroup, laggeddate, maxnotstale)).(vals))
  end

  for i ∈ 1:length(vals)
    df[!, laggedvals[i]] = fetch(tasks[i])
  end

  return nothing
end

function lagwithin2!(df::DataFrame, vals::Vector{<:DField}, group::DField;
  date::DField = error("date field is now required"),
  maxnotstale::Any = nothing,
  laggedvals::Vector{<:DField} = (s->Symbol(:L, s)).(vals))

  issorted(df, [group, date]) || error("df must be sorted by $group, $date")

  return lagwithin2sorted!(df, vals, group; date, maxnotstale, laggedvals)
end

function differencewithin2!(df::DataFrame,
  vals::Vector{<:DField},
  group::DField;
  date::DField = error("date field is now required"),
  differencedvals::Vector{<:DField} = (s::DField->Symbol(:D, s)).(vals),
  createlag::Bool=true,
  deletelag::Bool=true,
  maxnotstale::Any = nothing,
  laggedvals::Vector{<:DField} = (deletelag ?
    (s::DField->Symbol(:L, s, :_temp)).(vals) : (s->Symbol(:L, s)).(vals))
  )::Nothing

  (date ∈ vals) && error("date field cannot be differenced. workaround- create a copy")

  createlag && lagwithin2!(df, vals, group; date, laggedvals, maxnotstale)
  #println(df[1:50,:])
  for t ∈ 1:length(vals)
    df[!, differencedvals[t]] = df[!, vals[t]] .- df[!, laggedvals[t]]
  end
  deletelag && select!(df, Not(laggedvals))

  return nothing
end

function leadwithin2!(df::DataFrame, vals::Vector{<:DField}, group::DField;
  date::DField = error("date field is now required"),
  maxnotstale::Any = nothing,
  leadvals::Vector{<:DField} = (s->Symbol(:N, s)).(vals))

  issorted(df, [group, date]) || error("df must be sorted by $group, $date")

  #we reverse the target fields without copying
  inputfields = [date; vals; group;] |> unique!
  Nrows::Int = size(df,1)
  revdf = (f->f=>view(df[!,f], Nrows:-1:1)).(inputfields) |> (flds)->DataFrame(flds, copycols=false)
  #revdf = (f->f=>df[Nrows:-1:1,f]).(inputfields) |> DataFrame! (similar but copies columns)

  #now "lag" the reversed df (which is effectivelly a lead)
  lagwithin2sorted!(revdf, vals, group, date=date, maxnotstale=maxnotstale, laggedvals=leadvals)
  revrevdf = revdf[Nrows:-1:1,Not(inputfields)]
  for f ∈ propertynames(revrevdf)
    df[!, f] = revrevdf[!, f]
  end

  return nothing
end

max_stars_repo_name: clintonTE/Finometrics

max_issues_repo_name: clintonTE/Finometrics

max_forks_repo_name: clintonTE/Finometrics


text:
@testset "OpenCL.Macros" begin
    @testset "OpenCL.Macros version platform" begin
        for platform in cl.platforms()

            version = cl.opencl_version(platform)

            v11 = cl.@min_v11? platform true : false
            v12 = cl.@min_v12? platform true : false
            v20 = cl.@min_v20? platform true : false

            if version == v"1.0"
                @test v11 == false
                @test v12 == false
                @test v20 == false
            elseif version == v"1.1"
                @test v11 == true
                @test v12 == false
                @test v20 == false
            elseif version == v"1.2"
                @test v11 == true
                @test v12 == true
                @test v20 == false
            elseif version == v"2.0"
                @test v11 == true
                @test v12 == true
                @test v20 == true
            end
        end
    end

    @testset "OpenCL.Macros version device" begin
        for platform in cl.platforms()
            for device in cl.devices(platform)
                version = cl.opencl_version(device)

                v11 = cl.@min_v11? device true : false
                v12 = cl.@min_v12? device true : false
                v20 = cl.@min_v20? device true : false

                if version == v"1.0"
                    @test v11 == false
                    @test v12 == false
                    @test v20 == false
                elseif version == v"1.1"
                    @test v11 == true
                    @test v12 == false
                    @test v20 == false
                elseif version == v"1.2"
                    @test v11 == true
                    @test v12 == true
                    @test v20 == false
                elseif version == v"2.0"
                    @test v11 == true
                    @test v12 == true
                    @test v20 == true
                end
            end
        end
    end
end

ext: jl


text:
function do_gen_tests()
    println("testing url: hdfs://user:password@hdfshost:9000/root/folder/file.csv")
    up = urlparse("hdfs://user:password@hdfshost:9000/root/folder/file.csv")
    @test up.scheme == "hdfs"
    @test up.netloc == "user:password@hdfshost:9000"
    @test up.url == "/root/folder/file.csv"
    @test up.params == ""
    @test up.query == ""
    @test up.fragment == ""
    @test string(hostname(up)) == "hdfshost"
    @test string(port(up)) == "9000"
    @test string(username(up)) == "user"
    @test string(password(up)) == "password"
    @test urlunparse(up) == "hdfs://user:password@hdfshost:9000/root/folder/file.csv"

    println("testing url: https://user:password@httphost:9000/path1/path2;paramstring?q=a&p=r#frag")
    up = urlparse("https://user:password@httphost:9000/path1/path2;paramstring?q=a&p=r#frag")
    @test up.scheme == "https"
    @test up.netloc == "user:password@httphost:9000"
    @test up.url == "/path1/path2"
    @test up.params == "paramstring"
    @test up.query == "q=a&p=r"
    @test up.fragment == "frag"
    @test string(hostname(up)) == "httphost"
    @test string(port(up)) == "9000"
    @test string(username(up)) == "user"
    @test string(password(up)) == "password"
    @test urlunparse(up) == "https://user:password@httphost:9000/path1/path2;paramstring?q=a&p=r#frag"

    println("testing url: ftp://ftp.is.co.za/rfc/rfc1808.txt")
    up = urlparse("ftp://ftp.is.co.za/rfc/rfc1808.txt")
    @test up.scheme == "ftp"
    @test up.netloc == "ftp.is.co.za"
    @test up.url == "/rfc/rfc1808.txt"
    @test up.params == ""
    @test up.query == ""
    @test up.fragment == ""
    @test string(hostname(up)) == "ftp.is.co.za"
    @test string(port(up)) == "nothing"
    @test string(username(up)) == "nothing"
    @test string(password(up)) == "nothing"
    @test urlunparse(up) == "ftp://ftp.is.co.za/rfc/rfc1808.txt"

    println("testing url: http://www.ietf.org/rfc/rfc2396.txt")
    up = urlparse("http://www.ietf.org/rfc/rfc2396.txt")
    @test up.scheme == "http"
    @test up.netloc == "www.ietf.org"
    @test up.url == "/rfc/rfc2396.txt"
    @test up.params == ""
    @test up.query == ""
    @test up.fragment == ""
    @test string(hostname(up)) == "www.ietf.org"
    @test string(port(up)) == "nothing"
    @test string(username(up)) == "nothing"
    @test string(password(up)) == "nothing"
    @test urlunparse(up) == "http://www.ietf.org/rfc/rfc2396.txt"

    println("testing url: ldap://[2001:db8::7]/c=GB?objectClass?one")
    up = urlparse("ldap://[2001:db8::7]/c=GB?objectClass?one")
    @test up.scheme == "ldap"
    @test up.netloc == "[2001:db8::7]"
    @test up.url == "/c=GB?objectClass?one"
    @test up.params == ""
    @test up.query == ""
    @test up.fragment == ""
    @test string(hostname(up)) == "2001:db8::7"
    @test string(port(up)) == "nothing"
    @test string(username(up)) == "nothing"
    @test string(password(up)) == "nothing"
    @test urlunparse(up) == "ldap://[2001:db8::7]/c=GB?objectClass?one"

    println("testing url: mailto:John.Doe@example.com")
    up = urlparse("mailto:John.Doe@example.com")
    @test up.scheme == "mailto"
    @test up.netloc == ""
    @test up.url == "John.Doe@example.com"
    @test up.params == ""
    @test up.query == ""
    @test up.fragment == ""
    @test string(hostname(up)) == "nothing"
    @test string(port(up)) == "nothing"
    @test string(username(up)) == "nothing"
    @test string(password(up)) == "nothing"
    @test urlunparse(up) == "mailto:John.Doe@example.com"

    println("testing url: news:comp.infosystems.www.servers.unix")
    up = urlparse("news:comp.infosystems.www.servers.unix")
    @test up.scheme == "news"
    @test up.netloc == ""
    @test up.url == "comp.infosystems.www.servers.unix"
    @test up.params == ""
    @test up.query == ""
    @test up.fragment == ""
    @test string(hostname(up)) == "nothing"
    @test string(port(up)) == "nothing"
    @test string(username(up)) == "nothing"
    @test string(password(up)) == "nothing"
    @test urlunparse(up) == "news:comp.infosystems.www.servers.unix"

    println("testing url: tel:+1-816-555-1212")
    up = urlparse("tel:+1-816-555-1212")
    @test up.scheme == "tel"
    @test up.netloc == ""
    @test up.url == "+1-816-555-1212"
    @test up.params == ""
    @test up.query == ""
    @test up.fragment == ""
    @test string(hostname(up)) == "nothing"
    @test string(port(up)) == "nothing"
    @test string(username(up)) == "nothing"
    @test string(password(up)) == "nothing"
    @test urlunparse(up) == "tel:+1-816-555-1212"

    println("testing url: telnet://192.0.2.16:80/")
    up = urlparse("telnet://192.0.2.16:80/")
    @test up.scheme == "telnet"
    @test up.netloc == "192.0.2.16:80"
    @test up.url == "/"
    @test up.params == ""
    @test up.query == ""
    @test up.fragment == ""
    @test string(hostname(up)) == "192.0.2.16"
    @test string(port(up)) == "80"
    @test string(username(up)) == "nothing"
    @test string(password(up)) == "nothing"
    @test urlunparse(up) == "telnet://192.0.2.16:80/"

    println("testing url: urn:oasis:names:specification:docbook:dtd:xml:4.1.2")
    up = urlparse("urn:oasis:names:specification:docbook:dtd:xml:4.1.2")
    @test up.scheme == "urn"
    @test up.netloc == ""
    @test up.url == "oasis:names:specification:docbook:dtd:xml:4.1.2"
    @test up.params == ""
    @test up.query == ""
    @test up.fragment == ""
    @test string(hostname(up)) == "nothing"
    @test string(port(up)) == "nothing"
    @test string(username(up)) == "nothing"
    @test string(password(up)) == "nothing"
    @test urlunparse(up) == "urn:oasis:names:specification:docbook:dtd:xml:4.1.2"

end


text:
module Condition

using ..Space, ..Tagging, ..Traits
export |ᶜ, cnd, conditions, cond!, condf, BinaryPointwise, Conditional, ConditionException
export ==ₚ

# # Conditioning
# Conditioning a variable restricts the output to be consistent with some proposition.

"`x` given `y` is true"
struct Conditional{X, Y}
  x::X
  y::Y
end

@inline x |ᶜ y = cnd(x, y)
@inline cnd(x, y) = Conditional(x, y)

"Conditions variable was conditioned on are not satisfied"
struct ConditionException <: Exception end

@inline condf(traits, ω, x, y) = Bool(y(ω)) ? x(ω) : throw(ConditionException())

@inline condf(ω::Ω, x, y) where Ω = condf(traits(Ω), ω, x, y)

@inline tagcondition(ω, condition) = tag(ω, (condition = condition,))

# If error are violated then throw error
@inline (c::Conditional)(ω) = condf(ω, c.x, c.y)

"Conditions on `xy`"
conditions(xy::Conditional) = xy.y

# Implement x(\omega) when x is conditioned
# Imolement logpdf when `x` is conditioned

"""
`cond!(ω::Ω, bool)`

Condition intermediate values from within the functional definition of a `RandVar`

```
function x_(ω)
  x = 0.0
  xs = Float64[]
  while bernoulli(ω, 0.8, Bool)
    x += uniform(ω, -5.0, 5.0)
    cond!(ω, x <=ₛ 1.0)
    cond!(ω, x >=ₛ -1.0)
    push!(xs, x)
  end
  xs
end

x = ciid(x_)
samples = rand(x, 100; alg = SSMH)
```
"""
cond!(ω::Ω, bool) where Ω = cond!(traits(Ω), ω, bool)
@inline cond!(traits, ω, bool) = nothing

"`selfcond(x, pred)` conditioned on `predₚ(x)`"
selfcond(x, pred) = cnd(x, ω -> pred(x(ω)))

end


text:
include("../common/PWGrid_v02.jl")
include("../common/ortho_gram_schmidt.jl")
include("../common/wrappers_fft.jl")

include("EnergiesT.jl")
include("PotentialsT.jl")
include("op_K.jl")
include("op_V_loc.jl")
include("op_H.jl")
include("calc_rho.jl")
include("calc_grad_v2.jl")
include("calc_Energies.jl")
include("KS_solve_Emin_sd.jl")
include("KS_solve_Emin_cg.jl")
include("Poisson_solve.jl")
include("LDA_VWN.jl")
include("Kprec.jl")

include("../common/gen_dr_center.jl")
include("../common/calc_strfact_v2.jl")
include("../common/calc_ewald_v2.jl")

include("diag_lobpcg.jl")
include("KS_solve_SCF.jl")

function test_main( Ns )

    const LatVecs = 16.0*diagm( ones(3) )

    pw = PWGrid( Ns, LatVecs )

    const Ω  = pw.Ω
    const r  = pw.r
    const G  = pw.gvec.G
    const G2 = pw.gvec.G2
    const Npoints = prod(Ns)
    const Ngwx = pw.gvecw.Ngwx

    @printf("Ns   = (%d,%d,%d)\n", Ns[1], Ns[2], Ns[3])
    @printf("Ngwx = %d\n", Ngwx)

    const actual = Npoints/Ngwx
    const theor = 1/(4*pi*0.25^3/3)
    @printf("Compression: actual, theor: %f , %f\n", actual, theor)

    # Lithium atom
    Zion = 3.0
    const Nstates = 2
    Focc = [2.0, 1.0]

    Xpos = reshape( [0.0, 0.0, 0.0], (3,1) )

    Sf = calc_strfact( Xpos, 1, [1], pw.gvec.G )
    E_nn = calc_ewald( pw, Sf, Xpos, 1, [1], [Zion] )

    # Set up potential, using HGH pseudopotential for H
    # contains only local pseudopotential
    const rloc = 0.4
    const c1 = -14.034868
    const c2 = 9.553476
    const c3 = -1.766488
    const c4 = 0.084370
    #
    Vg = zeros(ComplexF64,Npoints)
    pre1 = -4*pi*Zion/Ω
    pre2 = sqrt(8*pi^3)*rloc^3/Ω
    #
    for ig=2:Npoints
        Gr = sqrt(G2[ig])*rloc
        expGr2 = exp(-0.5*Gr^2)
        Vg[ig] = pre1/G2[ig]*expGr2 + pre2*expGr2 * (c1 + c2*(3-Gr^2) +
                 c3*(15 - 10*Gr^2 + Gr^4) + c4*(105 - 105*Gr^2 + 21*Gr^4 - Gr^6) )
    end
    # limiting value
    Vg[1] = 2*pi*Zion*rloc^2 + (2*pi)^1.5 * rloc^3 * (c1 + 3.0*c2 + 15*c3 + 105*c4)
    V_ionic = real( G_to_R(Ns, Vg .* Sf) ) * Npoints

    # Need to sum up over Nspecies for more than one species type
    # We simply need reshape because we only have one species type here.
    V_ionic = reshape( V_ionic, (Npoints) )

    psi, Energies, Potentials = KS_solve_Emin_cg( pw, V_ionic, Focc, Nstates, NiterMax=1000, E_NN = E_nn )

    Y = ortho_gram_schmidt(psi)
    mu = Y' * op_H( pw, Potentials, Y )
    evals, evecs = eig(mu)
    psi = Y*evecs

    #@printf("Solution by self-consistent field method\n")
    #Energies, Potentials, psi, evals = KS_solve_SCF( pw, V_ionic, Focc, Nstates, β=0.3 )

    for st = 1:Nstates
        @printf("State # %d, Energy = %f\n", st, real(evals[st]))
    end

    @printf("E_nn    = %18.10f\n", E_nn)
    @printf("E total = %18.10f\n", E_nn + Energies.Total)

end

@time test_main( [64,64,64] )
#@time test_main( [100,100,100] )


text:
@testset "Threaded loading" begin
    if Threads.nthreads() > 1
        pkgid = Base.PkgId(FileIO.idCSVFiles.second, "CSVFiles")
        @assert !Base.root_module_exists(pkgid) "threaded loaidng test requires CSVFiles not loaded"
        # When threaded, loading a new backend should be locked
        # https://github.com/JuliaIO/FileIO.jl/issues/336
        files = map(["data.csv", "file.csv"]) do filename
            joinpath(@__DIR__, "files", filename)
        end
        @test_nowarn Threads.@threads for file in files
            load(file)
        end
    end
end


text:
DIFFEQ_DEP_WARN = """
Direct propagation of keyword arguments to DifferentialEquations.jl is deprecated.
From now on pass any DiffEq-related keywords as a `NamedTuple` using the
explicit keyword `diffeq` instead.
"""

@deprecate numericallyapunov(args...; kwargs...) lyapunov_from_data(args...; kwargs...)
@deprecate grassberger(args...; kwargs...) grassberger_dim(args...; kwargs...)


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

export read_step, set_topology!, set_cell!, path

__ptr(trajectory::Trajectory) = __ptr(trajectory.__handle)
__const_ptr(trajectory::Trajectory) = __const_ptr(trajectory.__handle)

"""
The ``Trajectory`` function opens a trajectory file, using the file at the given
``path``. The opening ``mode`` can be ``'r'`` for read, ``'w'`` for write or
``'a'`` for append, and defaults to ``'r'``. The optional ``format`` parameter
give a specific file format to use when opening the file.
"""
function Trajectory(path::AbstractString, mode::Char='r', format::AbstractString="")
    ptr = @__check_ptr(lib.chfl_trajectory_with_format(
        pointer(path), Int8(mode), pointer(format),
    ))
    return Trajectory(CxxPointer(ptr, is_const=false))
end

"""
Apply the function `f` to the result of `Trajectory(args...)` and close the
resulting trajectory upon completion, similar to ``open(f, args...)``.
"""
function Trajectory(f::Function, args...)
    tr = Trajectory(args...)
    try
        f(tr)
    finally
        close(tr)
    end
end

"""
Read the next step of the ``trajectory``, and return the corresponding ``Frame``.
"""
function Base.read(trajectory::Trajectory)
    @assert isopen(trajectory)
    frame = Frame()
    __check(lib.chfl_trajectory_read(
        __ptr(trajectory), __ptr(frame)
    ))
    return frame
end

"""
Read the given ``step`` of the ``trajectory``, and return the corresponding
``Frame``.
"""
function read_step(trajectory::Trajectory, step::Integer)
    @assert isopen(trajectory)
    frame = Frame()
    __check(lib.chfl_trajectory_read_step(
        __ptr(trajectory), UInt64(step), __ptr(frame)
    ))
    return frame
end

"""
Write the given ``frame`` to the ``trajectory``.
"""
function Base.write(trajectory::Trajectory, frame::Frame)
    @assert isopen(trajectory)
    __check(lib.chfl_trajectory_write(
        __ptr(trajectory), __const_ptr(frame)
    ))
    return nothing
end

"""
Set the ``Topology`` associated with a ``trajectory``. This topology will be
used when reading and writing the file, replacing any topology in the file.
"""
function set_topology!(trajectory::Trajectory, topology::Topology)
    @assert isopen(trajectory)
    __check(lib.chfl_trajectory_set_topology(
        __ptr(trajectory), __const_ptr(topology)
    ))
    return nothing
end

"""
Set the ``Topology`` associated with a ``trajectory`` by reading the first frame
of the file at ``path``; and extracting the topology of this frame. The optional
``format`` parameter can be used to specify the file format.
"""
function set_topology!(trajectory::Trajectory, path::AbstractString, format::AbstractString = "")
    @assert isopen(trajectory)
    __check(lib.chfl_trajectory_topology_file(
        __ptr(trajectory), pointer(path), pointer(format)
    ))
    return nothing
end

"""
Set the ``cell`` associated with a ``trajectory``. This cell will be used when
reading and writing the file, replacing any unit cell in the file.
"""
function set_cell!(trajectory::Trajectory, cell::UnitCell)
    @assert isopen(trajectory)
    __check(lib.chfl_trajectory_set_cell(__ptr(trajectory), __const_ptr(cell)))
    return nothing
end

"""
Get the number of steps (the number of frames) in a ``trajectory``.
"""
function Base.size(trajectory::Trajectory)
    @assert isopen(trajectory)
    result = Ref{UInt64}(0)
    __check(lib.chfl_trajectory_nsteps(
        __const_ptr(trajectory), result
    ))
    return result[]
end

"""
Get the number of steps (the number of frames) in a ``trajectory``.
"""
function Base.length(trajectory::Trajectory)
    size(trajectory)
end

"""
Get the path used to open a ``trajectory``.
"""
function path(trajectory::Trajectory)
    @assert isopen(trajectory)
    result = Ref{Ptr{UInt8}}(0)
    Base.cconvert(Ptr{Ptr{UInt8}}, result)
    __check(
        lib.chfl_trajectory_path(__const_ptr(trajectory), Base.unsafe_convert(Ptr{Ptr{UInt8}}, result))
    )
    return unsafe_string(result[])
end

"""
Close a ``trajectory``. This function flushes any buffer content to the hard
drive, and frees the associated memory. Necessary when running on the REPL to
finish  writing.
"""
function Base.close(trajectory::Trajectory)
    # Manually free and set the pointer to NULL
    lib.chfl_trajectory_close(trajectory.__handle.__ptr)
    trajectory.__handle.__ptr = 0
    return nothing
end

"""
Check if the ``trajectory`` is open.
"""
function Base.isopen(trajectory::Trajectory)
    return Int(trajectory.__handle.__ptr) != 0
end

# Iteration support
Base.iterate(trajectory::Trajectory, state=0) = state >= length(trajectory) ? nothing : (read_step(trajectory, state), state+1)
Base.eltype(trajectory::Trajectory) = Frame


text:
download("https://gr-framework.org/downloads/gr-fonts.tgz", "fonts.tgz")
run(`tar xzf fonts.tgz`)
rm("fonts.tgz")

ENV["GKS_FONT_DIRS"] = joinpath(pwd(), "fonts", "urw-base35")

const fonts = ("Times Roman", "Times Italic", "Times Bold", "Times Bold Italic", "Helvetica", "Helvetica Oblique", "Helvetica Bold", "Helvetica Bold Oblique", "Courier", "Courier Oblique", "Courier Bold", "Courier Bold Oblique", "Bookman Light", "Bookman Light Italic", "Bookman Demi", "Bookman Demi Italic", "New Century Schoolbook Roman", "New Century Schoolbook Italic", "New Century Schoolbook Bold", "New Century Schoolbook Bold Italic", "Avantgarde Book", "Avantgarde Book Oblique", "Avantgarde Demi", "Avantgarde Demi Oblique", "Palatino Roman", "Palatino Italic", "Palatino Bold", "Palatino Bold Italic", "Zapf Chancery Medium Italica", "Zapf Dingbats")

using GR

selntran(0)
setcharheight(0.018)
settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_HALF)

y = 0.95
for fontname in fonts
  global y
  font = loadfont(fontname)
  settextfontprec(font, GR.TEXT_PRECISION_OUTLINE)
  text(0.5, y, "The quick brown fox jumps over the lazy dog")
  tbx, tby = inqtext(0.5, y, "The quick brown fox jumps over the lazy dog")
  fillrect(tbx[1], tbx[2], tby[1], tby[3])
  y -= 0.03
end

updatews()


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

compute the capacity of a classical-quantum channel
adapted from https://github.com/hfawzi/cvxquad/blob/master/examples/cq_channel_capacity.m
and listing 1 in "Efficient optimization of the quantum relative entropy" by H. Fawzi and O. Fawzi
=#

struct ClassicalQuantum{T <: Real} <: ExampleInstanceJuMP{T}
    n::Int
    complex::Bool
    use_EF::Bool
end

function build(inst::ClassicalQuantum{T}) where {T <: Float64}
    n = inst.n
    @assert !(inst.complex && inst.use_EF)
    rt2 = sqrt(T(2))
    R = (inst.complex ? Complex{T} : T)
    function hermtr1()
        ρ = randn(R, n, n)
        ρ = ρ * ρ'
        ρ ./= tr(ρ)
        return Hermitian(ρ)
    end
    ρs = [hermtr1() for _ in 1:n]
    Hs = [dot(ρ, log(ρ)) for ρ in ρs]

    model = JuMP.Model()
    JuMP.@variable(model, prob[1:n] >= 0)
    JuMP.@constraint(model, sum(prob) == 1)
    JuMP.@variable(model, epi)
    JuMP.@objective(model, Max, -epi + dot(prob, Hs))

    entr = zeros(JuMP.AffExpr, Cones.svec_length(R, n))
    ρ_vec = zeros(T, length(entr))
    for (ρ, p) in zip(ρs, prob)
        Cones.smat_to_svec!(ρ_vec, ρ, rt2)
        entr += p * ρ_vec
    end

    aff = vcat(epi, 1, entr)
    if inst.use_EF
        add_spectral(MatNegEntropyEigOrd(), n, aff, model)
    else
        JuMP.@constraint(model, aff in Hypatia.EpiPerSepSpectralCone{Float64}(
            Cones.NegEntropySSF(), Cones.MatrixCSqr{T, R}, n))
    end

    # save for use in tests
    model.ext[:epi] = epi
    model.ext[:entr] = entr

    return model
end

function test_extra(inst::ClassicalQuantum{T}, model::JuMP.Model) where T
    stat = JuMP.termination_status(model)
    @test stat == MOI.OPTIMAL
    (stat == MOI.OPTIMAL) || return

    # check constraint
    tol = eps(T)^0.2
    epi_opt = JuMP.value(model.ext[:epi])
    entr_opt = JuMP.value.(model.ext[:entr])
    Entr_opt = zeros(inst.complex ? Complex{T} : T, inst.n, inst.n)
    Cones.svec_to_smat!(Entr_opt, entr_opt, sqrt(T(2)))
    λ = eigvals(Hermitian(Entr_opt, :U))
    @test minimum(λ) >= -tol
    qe_result = get_val(pos_only(λ), MatNegEntropy())
    @test epi_opt ≈ qe_result atol=tol rtol=tol
    return
end


text:
# Import Base functions (copied from JUDI)
import Base.*, Base./, Base.+, Base.-
import Base.copy!, Base.copy
import Base.sum, Base.ndims, Base.reshape, Base.fill!, Base.axes, Base.dotview
import Base.eltype, Base.length, Base.size, Base.iterate, Base.show, Base.display, Base.showarg
import Base.maximum, Base.minimum, Base.push!
import Base.Broadcast.broadcasted, Base.BroadcastStyle, Base.Broadcast.DefaultArrayStyle, Base.Broadcast, Base.broadcast!
import Base.getindex, Base.setindex!, Base.firstindex, Base.lastindex
import Base.similar, Base.isapprox, Base.isequal
import Base.materialize!, Base.materialize
import Base.promote_shape, Base.diff, Base.cumsum, Base.cumsum!

# h5geo includes
include("H5Seis.jl")
include("H5SeisCon.jl")
include("read_H5SeisCon.jl")
include("H5Geometry.jl")
include("H5judiVector.jl")
include("auxiliaryFunctions.jl")


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module ModiaMath.Frames

This module contains functions for **frames** that is coordinate systems in 3D.
The orientation of a frame is described either with a 3x3 **rotation matrix**
or with a **quaternion vector** and its origin is described with a **Vector3D**:

- `const `[`ModiaMath.RotationMatrix`](@ref) = SMatrix{3,3,Float64,9}`:
  Type of a Rotation matrix to rotate from a frame 1 into a frame 2.

- `const `[`ModiaMath.Quaternion`](@ref) = SVector{4,Float64}`:
  Type of a Quaternion vector to rotate from a frame 1 into a frame 2.

- `const `[`ModiaMath.Vector3D`](@ref) = SVector{3,Float64}`:
  Type of a vector in 3D space (e.g. position vector of the origin of a frame).

The following constants are defined

- `const `[`ModiaMath.NullRotation`](@ref):
  RotationMatrix with no rotation from a frame 1 into a frame 2.

- `const `[`ModiaMath.NullQuaternion`](@ref):
  Quaternion vector with no rotation from a frame 1 into a frame 2.

- `const `[`ModiaMath.ZeroVector3D`](@ref):
  Vector3D with only zero elements.

If an angle is given as an argument to one of the functions below, it might be a
number (interpreted as having unit `rad`) or a number with a unit
(for example: `using Unitful; angle = 90u"°"`).

# Constructors for a RotationMatrix R

The following functions return a [`ModiaMath.RotationMatrix`](@ref)` R`
to rotate a frame 1 into a frame 2.

| Function                                           | Description                                |
|:---------------------------------------------------|:-------------------------------------------|
| [`ModiaMath.rot1`](@ref)(angle)                    | Rotate around `angle` along x-axis         |
| [`ModiaMath.rot2`](@ref)(angle)                    | Rotate around `angle` along y-axis         |
| [`ModiaMath.rot3`](@ref)(angle)                    | Rotate around `angle` along z-axis         |
| [`ModiaMath.rot123`](@ref)(angle1, angle2, angle3) | Rotate around angles along x,y,z-axes      |
| [`ModiaMath.rot_e`](@ref)(e, angle)                | Rotate around `angle` along unit vector `e`|
| [`ModiaMath.rot_nxy`](@ref)(nx, ny)                | `nx`/`ny` are in x/y-direction of frame 2  |
| [`ModiaMath.from_q`](@ref)(q)                      | Return `R` from Quaternion `q`             |

# Constructors for a Quaternion q

The following functions return a [`ModiaMath.Quaternion`](@ref)` q`
to rotate a frame 1 into a frame 2.
Since `q` and `-q` define the same rotation the constructor functions have
a keyword argument `q_guess::Quaternion = NullQuaternion`.
From the two possible solutions `q` the one is returned that is closer to `q_guess`.

| Function                                            | Description                                |
|:----------------------------------------------------|:-------------------------------------------|
| [`ModiaMath.qrot1`](@ref)(angle)                    | Rotate around `angle` along x-axis         |
| [`ModiaMath.qrot2`](@ref)(angle)                    | Rotate around `angle` along y-axis         |
| [`ModiaMath.qrot3`](@ref)(angle)                    | Rotate around `angle` along z-axis         |
| [`ModiaMath.qrot123`](@ref)(angle1, angle2, angle3) | Rotate around angles along x,y,z-axes      |
| [`ModiaMath.qrot_e`](@ref)(e, angle)                | Rotate around `angle` along unit vector `e`|
| [`ModiaMath.qrot_nxy`](@ref)(nx, ny)                | `nx`/`ny` are in x/y-direction of frame 2  |
| [`ModiaMath.from_R`](@ref)(R)                       | Return `q` from RotationMatrix `R`         |

# Operations on Frames

The following functions provide operations on frames. The orientation of a frame is
defined with argument `Rq` meaning it can be either a
[`ModiaMath.RotationMatrix`](@ref) ` R` or a
[`ModiaMath.Quaternion`](@ref) ` q` (to rotate a frame 1 into a frame 2).

| Function                                         | Description                                   |
|:-------------------------------------------------|:----------------------------------------------|
| [`ModiaMath.resolve1`](@ref)(Rq, v2)             | Transform vector `v` from frame 2 to frame 1  |
| [`ModiaMath.resolve2`](@ref)(Rq, v1)             | Transform vector `v` from frame 1 to frame 2  |
| [`ModiaMath.absoluteRotation`](@ref)(Rq01, Rq12) | Return rotation 0->2 from rot. 0->1 and 1->2  |
| [`ModiaMath.relativeRotation`](@ref)(Rq01, Rq02) | Return rotation 1->2 from rot. 0->1 and 0->2  |
| [`ModiaMath.inverseRotation`](@ref)(Rq01)        | Return rotation 1->0 from rot, 0->1           |
| [`ModiaMath.planarRotationAngle`](@ref)(e,v1,v2) | Return angle of planar rotation along `e`     |
| [`ModiaMath.eAxis`](@ref)(axis)                  | Return unit vector `e` in direction of `axis` |
| [`ModiaMath.skew`](@ref)(v)                      | Return skew-symmetric matrix of vector v      |

# Interpolation of Frames

Given a set of frames by a vector `r` of position vectors to their origins and
and an optional vector `q` of Quaternions of their absolute orientations, then
the following functions interpolate linearly in these frames:

| Function                                          | Description                                    |
|:--------------------------------------------------|:-----------------------------------------------|
| [`ModiaMath.Path`](@ref)(r,q)                     | Return path defined by a vector of frames      |
| [`ModiaMath.t_pathEnd`](@ref)(path)               | Return path parameter `t_end` of last frame    |
| [`ModiaMath.interpolate`](@ref)(path,t)           | Return `(rt,qt)` of Path at path parameter `t` |
| [`ModiaMath.interpolate_r`](@ref)(path,t)         | Return `rt` of Path at path parameter `t`      |

# Examples

```julia
using ModiaMath
using Unitful

# R1,R2,R3 are the same RotationMatrices
R1 = ModiaMath.rot1(pi/2)
R2 = ModiaMath.rot1(90u"°")
R3 = ModiaMath.rot_e([1,0,0], pi/2)
```

# Main developer

[Martin Otter](https://rmc.dlr.de/sr/de/staff/martin.otter/),
[DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)

*The functions of this module are mostly a mapping of some of the functions
of the Modelica Standard Library from Modelica
([Modelica.Mechanics.MultiBody.Frames](https://doc.modelica.org/help/Modelica_Mechanics_MultiBody_Frames.html#Modelica.Mechanics.MultiBody.Frames))
to Julia (taking advantage of Julia features
such as multiple dispatch and unit package Unitful)*.
"""
module Frames

export RotationMatrix, NullRotation  , assertRotationMatrix
export Quaternion    , NullQuaternion, assertQuaternion
export Vector3D      , ZeroVector3D

export  rot1,  rot2,  rot3,  rot123,  rot_e,  rot_nxy, from_q
export qrot1, qrot2, qrot3, qrot123, qrot_e, qrot_nxy, from_R

export resolve1, resolve2, absoluteRotation, relativeRotation, inverseRotation
export planarRotationAngle, eAxis

export Path, t_pathEnd, interpolate, interpolate_r
export skew

using StaticArrays
import ModiaMath

@static if VERSION >= v"0.7.0-DEV.2005"
    using LinearAlgebra
end

include("vector3D.jl")
include("rotationMatrix.jl")
include("quaternion.jl")
include("interpolation.jl")

end


text:
###############################################################################
#
# Update leaf allocation
#
###############################################################################
"""
    leaf_allocation!(
                node::SPACSimple{FT},
                laba::FT
    ) where {FT<:AbstractFloat}
    leaf_allocation!(
                node::SPACSimple{FT},
                photo_set::AbstractPhotoModelParaSet{FT},
                vmax::FT
    ) where {FT<:AbstractFloat}
    leaf_allocation!(
                node::SPACSimple{FT},
                photo_set::AbstractPhotoModelParaSet{FT},
                laba::FT,
                vmax::FT
    ) where {FT<:AbstractFloat}

Update leaf area and maximal carboxylation rate, given
- `node` [`SPACSimple`] type struct
- `photo_set` [`AbstractPhotoModelParaSet`] type struct
- `laba` Given leaf area per basal area
- `vmax` Given Vcmax25
"""
function leaf_allocation!(
            node::SPACSimple{FT},
            laba::FT
) where {FT<:AbstractFloat}
    node.laba         = laba;
    node.hs.leaf.area = node.hs.stem.area * laba;
    node.lai          = laba / node.gaba;

    return nothing
end

function leaf_allocation!(
            node::SPACSimple{FT},
            photo_set::AbstractPhotoModelParaSet{FT},
            vmax::FT
) where {FT<:AbstractFloat}
    node.ps.Vcmax25 = vmax;
    node.ps.Jmax25  = vmax * node.vtoj;
    node.ps.Rd25    = vmax * photo_set.VR;

    return nothing
end

function leaf_allocation!(
            node::SPACSimple{FT},
            photo_set::AbstractPhotoModelParaSet{FT},
            laba::FT,
            vmax::FT
) where {FT<:AbstractFloat}
    leaf_allocation!(node, laba);
    leaf_allocation!(node, photo_set, vmax);

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
#!/usr/bin/env julia

using EAGOParametricInterval

# write your own tests here
println("Testing Checks and Utilities...")
t = @elapsed include("ParamChk_Tests.jl")
println("done (took $t seconds).")

println("Testing Sparse Preconditioner...")
t = @elapsed include("SparseCntr_Tests.jl")
println("done (took $t seconds).")

println("Testing Parametric Contractors...")
t = @elapsed include("ParametricContractor_Tests.jl")
println("done (took $t seconds).")


text:
##################################################################
# READ TROUBLESOME FILES WITH TOO MUCH WHITESPACE IN PLACES
##################################################################
read_and_unspace() = strip(replace(readline(), r"\s+", " "))
read_ints() = [parse(Int, x) for x in split(read_and_unspace(), " ")]

##################################################################
# PRIMES AND DIVISORS
##################################################################

# prime sieve with prime factors dictionary
function primesAndFactors(n::Int)
    sieve = [true for i in 1:n]
    sieve[1] = false
    prime_factors = [Dict{Int,Int}() for i in 1:n]
    for i in 2:n
        if sieve[i]
            prime_factors[i][i] = 1
            for j in (2i):i:n
                sieve[j] = false
                mul = 1
                k = j/i
                while k%i == 0
                    mul += 1
                    k/=i
                end
                prime_factors[j][i] = mul
            end
        end
    end
    primes = [i for i in 1:n if sieve[i]]
    return (sieve, primes, prime_factors)
end

# get the divisors of a number given its prime factorisation
function divisors_from_factors(factors::Dict{Int,Int})
    if length(factors)==0
        return Array{Int,1}()
    end
    ndivisors = prod(1+collect(values(factors)))

    bases = collect(keys(factors))
    exponents = [factors[x] for x in bases]
    divisors = zeros(Int, ndivisors)
    offset = 1
    N = length(bases)

    function enumerate_divs(cur_idx, cur_val)
        if cur_idx > N
            divisors[offset] = cur_val
            offset+=1
        else
            mul = cur_val
            for i in 0:exponents[cur_idx]
                enumerate_divs(cur_idx+1, mul)
                mul *= bases[cur_idx]
            end
        end
    end
    enumerate_divs(1, 1)
    return divisors
end

##################################################################
# Least unused natural number structure
##################################################################
type least_unused
    start_pos::Int
    cur_pos::Int
    unused::Base.Collections.PriorityQueue{Int, Int, Base.Order.ForwardOrdering}
    least_unused(sp::Int=0) = new(sp, sp, Base.Collections.PriorityQueue{Int, Int, Base.Order.ForwardOrdering}())
end

function unuse!(l::least_unused, val::Int)
    l.unused[val] = val
end

function use!(l::least_unused)
    if(isempty(l.unused))
        l.cur_pos += 1
        return l.cur_pos
    else
        return Base.Collections.dequeue!(l.unused)
    end
end

function use!(l::least_unused, n::Int)
    if n > l.cur_pos
        for k in (l.cur_pos+1):(n-1)
            unuse!(l, k)
        end
        l.cur_pos = n
    else
        l.unused[n] = -1
        Base.Collections.dequeue!(l.unused)
    end
end

function top(l::least_unused)
    if(isempty(l.unused))
        return l.cur_pos + 1
    else
        return Base.Collections.peek(l.unused)[1]
    end
end

##################################################################
# Interval tree to aggregate any function over an interval
##################################################################

type itreeNode{T<:Number}
    value::T
    pos::Int # position at which the value is obtained
    left::Int
    right::Int
    children::Vector{itreeNode{T}}
    itreeNode() = new(0, 0, 0, 0, Vector{itreeNode{T}}())
    itreeNode(v::T, p::Int, l::Int, r::Int) = new(v,p,l,r,Vector{itreeNode{T}}())
    itreeNode(v::T, p::Int, l::Int, r::Int, lc::itreeNode{T}, rc::itreeNode{T}) = new(v,p, l, r, [lc,rc])
    itreeNode(v::T, p::Int, l::Int, r::Int, ch::Vector{itreeNode{T}}) = new(v,p,l,r,ch)
end
itreeNode() = itreeNode{Int}()

function create_itree_closure{T}(v::Vector{T}, agg::Function, isleft::Function)
    function create_interval_tree(l::Int, r::Int)
        if l == r
            return itreeNode{T}(v[l], l, l, r)
        else
            mid_l = (l+r)>>1
            mid_r = mid_l + 1
            lhs = create_interval_tree(l, mid_l)
            rhs = create_interval_tree(mid_r, r)
            val = agg(lhs.value, rhs.value)
            left_side = isleft(lhs.value, rhs.value)
            pos = 0
            if left_side
                pos = lhs.pos
            else
                pos = rhs.pos
            end
            return itreeNode{T}(val, pos, l, r, lhs, rhs)
        end
    end
    return create_interval_tree
end

type itree{T}
    fn::Function
    isleft::Function
    x::Vector{T}
    root::itreeNode{T}
    default::T
    itree(v::Vector{T}, fn::Function, def::T) = new(fn, (x,y)->true, v, create_itree_closure(v,fn, (x,y)->true)(1, length(v)), def)
    itree(v::Vector{T}, fn::Function, leftfn::Function, def::T) = new(fn, leftfn, v, create_itree_closure(v,fn, leftfn)(1, length(v)), def)
end
itree(v::Vector{Int}, fn::Function, leftfn::Function, def::Int) = itree{Int}(v,fn,leftfn,def)
mintree(v::Vector{Int}) = itree{Int}(v, min, (x,y)->(x<=y), maximum(v)+1)
maxtree(v::Vector{Int}) = itree{Int}(v, max, (x,y)->(x>=y), minimum(v)-1)

# Query an itree on an interval
function query{T}(it::itree{T}, left::Int, right::Int)
    function qry_inner(cur::itreeNode{T}, l::Int, r::Int)
        if r < cur.left || l > cur.right
            return it.default
        elseif l <= cur.left && cur.right <= r
            return cur.value
        else
            mid_l = cur.children[1].right
            mid_r = cur.children[2].left
            if r <= mid_l
                return qry_inner(cur.children[1], l, r)
            elseif l >= mid_r
                return qry_inner(cur.children[2], l, r)
            else
                return it.fn(qry_inner(cur.children[1], l, min(r,mid_l)), qry_inner(cur.children[2], max(mid_r, l), r))
            end
        end
    end
    return qry_inner(it.root, left, right)
end

# Query an itree on an interval
function query_pos{T}(it::itree{T}, left::Int, right::Int)
    function qry_inner(cur::itreeNode{T}, l::Int, r::Int)
        if r < cur.left || l > cur.right
            return it.default
        elseif l <= cur.left && cur.right <= r
            return (cur.pos, cur.value)
        else
            mid_l = cur.children[1].right
            mid_r = cur.children[2].left
            if r <= mid_l
                return qry_inner(cur.children[1], l, r)
            elseif l >= mid_r
                return qry_inner(cur.children[2], l, r)
            else
                x = qry_inner(cur.children[1], l, min(r,mid_l))
                y = qry_inner(cur.children[2], max(mid_r, l), r)
                if it.isleft(x[2],y[2])
                    return (x[1], it.fn(x[2],y[2]))
                else
                    return (y[1], it.fn(x[2], y[2]))
                end
            end
        end
    end
    return qry_inner(it.root, left, right)[1]
end

##################################################################
# next smaller index algorithm in linear time
##################################################################

# Find the index of the next number in an array to be smaller than
# each element and to the right
function get_next_smaller(A::Vector{Int}, default_elt::Int=0)
    n = length(A)
    unmatched = Base.Collections.PriorityQueue{Int, Int, Base.Order.ForwardOrdering}()
    next_smaller = [default_elt for i in 1:n]

    for i in 1:n
        while !isempty(unmatched)
            (idx, val) = Base.Collections.peek(unmatched)
            val = -val
            if val > A[i]
                next_smaller[idx] = i
                Base.Collections.dequeue!(unmatched)
            else
                break
            end
        end
        unmatched[i] = -A[i]
    end
    return next_smaller
end

##################################################################
# Strongly connected components
##################################################################

# 1-indexed directed adjacency set graph
type directed_adjset_graph
    n::Int
    adj::Array{Set{Int},1}
    directed_adjset_graph(n_vert::Int) = new(n_vert, [Set{Int}() for i in 1:n_vert])
end

function add_edge!(G::directed_adjset_graph, u::Int, v::Int)
    if u == G.n+1 || v == G.n + 1
        push!(G.adj, Set{Int}())
        G.n += 1
    end

    push!(G.adj[u], v)
end

function is_edge(G::directed_adjset_graph, u::Int, v::Int)
    return (v in G.adj[u])
end

# Tarjan's algorithm to compute strongly connected components & condensation
function tarjan(G::directed_adjset_graph)
    index = 0
    dfs_stack   = Vector{Tuple}()
    ddone       = falses(G.n)
    ddescending = trues(G.n)
    S        = Vector{Int}()
    vindex   = zeros(Int, G.n)
    vlowlink = zeros(Int, G.n)
    visited  = falses(G.n)
    onStack  = falses(G.n)
    SCC_map  = zeros(Int, G.n)
    cur_SCC  = 0

    for vertex in 1:G.n
        if !visited[vertex]
            push!(dfs_stack, (vertex, vertex))
            while(!isempty(dfs_stack))
                v_prev,v = dfs_stack[end]
                if ddone[v]
                    pop!(dfs_stack)
                    continue
                end

                visited[v] = true

                # Set the depth index for v to the smallest unused index
                if ddescending[v]
                    ddescending[v] = false

                    vindex[v] = index
                    vlowlink[v] = index
                    index += 1
                    push!(S, v)
                    onStack[v] = true

                    for w in G.adj[v]
                        if !visited[w]
                            push!(dfs_stack, (v,w))
                        elseif onStack[w]
                            vlowlink[v] = min(vlowlink[v], vindex[w])
                        end
                    end
                else
                    v_prev, v = pop!(dfs_stack)
                    ddone[v] = true
                    vlowlink[v_prev] = min(vlowlink[v_prev], vlowlink[v])

                    # If v is a root node, pop the stack and generate an SCC
                    if vlowlink[v] == vindex[v]
                        cur_SCC += 1
                        while !isempty(S)
                            w = pop!(S)
                            onStack[w] = false
                            SCC_map[w] = cur_SCC
                            if w == v
                                break
                            end
                        end
                    end
                end
            end
        end
    end

    condensation = directed_adjset_graph(cur_SCC)
    for i in 1:G.n
        for j in G.adj[i]
            if SCC_map[i] != SCC_map[j]
               add_edge!(condensation, SCC_map[i], SCC_map[j])
           end
        end
    end

    return (SCC_map, condensation)
end

#################################################################
#          Connected components                                 #
#################################################################

# 1-indexed undirected adjacency set graph
type adjset_graph
    n::Int
    adj::Array{Set{Int},1}
    adjset_graph(n_vert::Int) = new(n_vert, [Set{Int}() for i in 1:n_vert])
end

function add_edge!(G::adjset_graph, u::Int, v::Int)
    if u == G.n+1 || v == G.n + 1
        push!(G.adj, Set{Int}())
        G.n += 1
    end

    push!(G.adj[u], v)
    push!(G.adj[v], u)
end

function is_edge(G::adjset_graph, u::Int, v::Int)
    return (v in G.adj[u])
end

immutable Edge
    u::Integer
    v::Integer
end

# Simple DFS (can be done much more efficiently with union-find)
function get_components(g::adjset_graph)
    components = zeros(Int, g.n)
    v0 = 1
    cpt = 1
    while v0 <= g.n
        if components[v0] == 0
            S = [v0]
            while !isempty(S)
                v = pop!(S)
                components[v] = cpt
                for w in g.adj[v]
                    if components[w] == 0
                        push!(S, w)
                    end
                end
            end
            cpt += 1
        end
        v0 += 1
    end
    return components
end

#################################################################
#          Biconnected components (adapted from LightGraphs.jl) #
#################################################################

type Biconnections
    low::Vector{Int}
    depth::Vector{Int}
    stack::Vector{Edge}
    biconnected_comps::Vector{Vector{Edge}}
    is_articulation::Vector{Bool}
    id::Int
end

function Biconnections(g::adjset_graph)
    n = g.n
    return Biconnections(zeros(Int, n), zeros(Int, n), Vector{Edge}(), Vector{Vector{Edge}}(), falses(n),0)
end

# Biconnected component dfs
function bc_dfs!(g::adjset_graph, state::Biconnections, u::Integer, v::Integer)
    children = 0
    state.id += 1
    state.depth[v] = state.id
    state.low[v] = state.depth[v]

    for w in g.adj[v]
        if state.depth[w] == 0
            children += 1
            push!(state.stack, Edge(min(v, w), max(v, w)))
            bc_dfs!(g, state, v, w)
            state.low[v] = min(state.low[v], state.low[w])

            if (u == v && children > 1) || (u != v && state.low[w] >= state.depth[v])
                state.is_articulation[v] = true
                e = Edge(0, 0)
                st = Vector{Edge}()
                while e != Edge(min(v, w),max(v, w))
                    e = pop!(state.stack)
                    push!(st, e)
                end
                push!(state.biconnected_comps, st)
            end

        elseif w != u && state.low[v] > state.depth[w]
            push!(state.stack, Edge(min(v, w), max(v, w)))
            state.low[v] = state.depth[w]
        end
    end
end

function do_bc_dfs(g::adjset_graph)
    state = Biconnections(g)
    for u in 1:g.n
        if state.depth[u] == 0
            bc_dfs!(g, state, u, u)
        end

        if !isempty(state.stack)
            push!(state.biconnected_comps, reverse(state.stack))
            empty!(state.stack)
        end
    end
    return state
end

function biconnected_components(g::adjset_graph)
    state = do_bc_dfs(g)
    return state.biconnected_comps
end

# Root a forest at the best possible roots, creating a parent array
# Root a forest at the best possible roots, creating a parent array
function root_forest(T::adjset_graph)
    vis_1 = falses(T.n)
    vis_2 = falses(T.n)
    vis_3 = falses(T.n)
    par = zeros(Int, T.n)
    dep = zeros(Int, T.n)
    v0 = 1
    while v0 <= T.n
        if par[v0] == 0
            v = v0

            # BFS 1: random node (v) to furthest point (w)
            vis_1[v] = true
            S = [v]
            w = v
            while !isempty(S)
                w = S[1]
                deleteat!(S, 1)
                for nbr in T.adj[w]
                    if !vis_1[nbr]
                        vis_1[nbr] = true
                        push!(S, nbr)
                    end
                end
            end

            # BFS 2: previous furthest point (v <- w_old) to furthest point from it (w)
            v = w
            S = [(v, 0)]
            vis_2[v] = true
            d = 0

            while !isempty(S)
                w, d = S[1]
                deleteat!(S,1)
                for nbr in T.adj[w]
                    if !vis_2[nbr]
                        vis_2[nbr] = true
                        push!(S, (nbr, d+1))
                    end
                end
            end

            # DFS: go from v to w & find half-way mark - the root of the current tree
            S = [v]
            npops = div(d+(d%2),2)

            root = v
            hst = []
            while !isempty(S)
                wcur = pop!(S)
                push!(hst, wcur)
                vis_3[wcur] = true

                if wcur == w
                    root = pop!(hst)
                    for idx in 1:npops
                        root = pop!(hst)
                    end
                    break
                end
                for nbr in T.adj[wcur]
                    if !vis_3[nbr]
                        vis_3[nbr] = true
                        push!(S, nbr)
                    end
                end
            end

            # Root the parent array at root with... yet another DFS!
            S = [root]
            par[root] = root # parent of root is itself by definition
            while !isempty(S)
                w = pop!(S)
                for nbr in T.adj[w]
                    if par[nbr] == 0
                        par[nbr] = w
                        dep[nbr] = dep[w] + 1
                        push!(S, nbr)
                    end
                end
            end
        end

        v0 += 1
    end
    return (par, dep)
end

# Get the block cut tree of an undirected graph
function block_cut_tree(g::adjset_graph)
    state = do_bc_dfs(g)
    a_pts = find(state.is_articulation)
    N0 = length(state.biconnected_comps)
    N1 = length(a_pts)
    N  = N0 + N1
    a_lkup = Dict{Int,Int}(a_pts[i] => N0+i for i in 1:N1)
    bctree = adjset_graph(N)
    v_lkup = zeros(Int,g.n)

    for cpt in 1:N0
        for e in state.biconnected_comps[cpt]
            if state.is_articulation[e.u]
                add_edge!(bctree, a_lkup[e.u], cpt)
                v_lkup[e.u] = a_lkup[e.u]
            else
                v_lkup[e.u] = cpt
            end
            if state.is_articulation[e.v]
                add_edge!(bctree, a_lkup[e.v], cpt)
                v_lkup[e.v] = a_lkup[e.v]
            else
                v_lkup[e.v] = cpt
            end
        end
    end

    return (v_lkup, bctree)
end

function scratch_bc()
    # sample code
    g = adjset_graph(9)
    add_edge!(g, 1,2)
    add_edge!(g, 1,3)
    add_edge!(g, 2,4)
    add_edge!(g, 3,4)
    add_edge!(g, 4,5)
    add_edge!(g, 4,6)
    add_edge!(g, 5,7)
    add_edge!(g, 6,7)
    add_edge!(g, 5,9)
    add_edge!(g, 7,8)
    add_edge!(g, 8,9)
    add_edge!(g, 7,10)
    add_edge!(g, 10,11)
    add_edge!(g, 10,12)

    biconnected_components(g)
end

#################################################################
#          Fast intersection of two sorted lists                #
#################################################################

function intersect_sorted(x::Vector{Int}, y::Vector{Int})
    if isempty(x) || isempty(y)
        return Vector{Int}
    end
    px = 1
    py = 1
    sol = Vector{Int}()
    while px <= length(x) && py <= length(y)
        while px <= length(x) && x[px] < y[py]
            px += 1
        end
        if px > length(x)
            break
        end

        while py <= length(y) && y[py] < x[px]
            py += 1
        end
        if py > length(y)
            break
        end

        while px <= length(x) && py <= length(y) && x[px] == y[py]
            push!(sol, x[px])
            px += 1
            py += 1
        end
    end
    return sol
end

##################################################################
# TRIE
##################################################################
type MinTrie
    value::Int
    children::Dict{Char,MinTrie}
end

max_stars_repo_name: DarioFanucchi/CompetitionCode

max_issues_repo_name: DarioFanucchi/CompetitionCode

max_forks_repo_name: DarioFanucchi/CompetitionCode


text:
type TraceNFXP{Tv, Tm} <: AbstractTrace
    Θ::Tm
    norm_ΔΘ::Tv
    ll::Tv
    Δll::Tv
    norm_g::Tv
end

TraceNFXP(K, n_var) = TraceNFXP(zeros(K, n_var),
                                zeros(K),
                                zeros(K),
                                zeros(K),
                                zeros(K))

type ConvergenceNFXP <: ConvergenceInfo
    # Outer loop
    flag::Bool
    outer::Int64
    # Inner loop
    iter::Vector{Int64}
    norm_sol::Float64
end
ConvergenceNFXP() = ConvergenceNFXP(false, 0, [0; 0], 0.)

outer_iterations(conv::ConvergenceNFXP) = conv.outer
inner_iterations(conv::ConvergenceNFXP) = conv.iter
contraction_iterations(conv::ConvergenceNFXP) = inner_iterations(conv)[1]
newton_iterations(conv::ConvergenceNFXP) = inner_iterations(conv)[2]

type NFXP{Tu<:AbstractUtility, Tm<:AbstractMatrix,
          S<:Stocas.AbstractState, Tv<:Stocas.AbstractValueFunction} <: EstimationMethod
    U::Tu
    X::S
    D::Data
    P::Vector{Vector{Float64}}
    logP::Vector{Vector{Float64}}
    Fᵁ::Tm
    V::Tv
    v::Vector{Vector{Float64}}
    dVu::Matrix{Float64}
    verbose::Bool
    K
end

NFXP(U, X, D; verbose = false, K = 50) = NFXP(U, X, D, [zeros(X.nX) for p = 1:length(U.U)], [zeros(X.nX) for p = 1:length(U.U)],
                    similar(X.F[1]), IntegratedValueFunction(X), [zeros(X.nX) for p = 1:length(U.U)], zeros(X.nX,length(U.Θ)), verbose, K)

function Base.display{M<:NFXP}(est_res::EstimationResults{M})
    @printf "Results of estimation\n"
    @printf " * Method: NFXP\n"
    @printf " * loglikelihood: %s\n" loglikelihood(est_res)
    if length(join(Optim.minimizer(est_res.trace), ",")) < 40
       @printf " * Estimate: [%s]\n" join(Optim.minimizer(est_res.trace), ",")
       @printf " * Std.err.: [%s]\n" join(stderr(est_res), ",")
   else
       @printf " * Estimate: [%s, ...]\n" join(Optim.minimizer(est_res.trace)[1:2], ",")
       @printf " * Std.err.: [%s, ...]\n" join(stderr(est_res)[1:2], ",")
   end
   @printf " Iterations\n"
   @printf " * Maximum likelihood: %d\n" Optim.iterations(est_res.trace)
   @printf " * Contractions: %d\n" est_res.conv.iter[1]
   @printf " * Newton steps: %d\n" est_res.conv.iter[2]
end

type NFXPCacheVars
    Θᵏ
    last_Θᵏ # Not Θᵏ⁺¹, but the last parameter vector that was used to solve the model
    score
end

function NFXPCacheVars(nvar, n)
         Θᵏ = zeros(nvar)
    last_Θᵏ = zeros(nvar)
    score   = zeros(n, nvar)
    NFXPCacheVars(Θᵏ, last_Θᵏ, score)
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 3c162d31-407a-4c50-a8a7-30153192f207
begin
	let
		using Pkg
		Pkg.activate(mktempdir())
		Pkg.Registry.update()
		Pkg.add("PlutoUI")
		Pkg.add("CairoMakie")
		Pkg.add("Statistics")
		Pkg.add("StatsBase")
		Pkg.add("ImageMorphology")
		Pkg.add("ImageFiltering")
		Pkg.add("CSV")
		Pkg.add("DataFrames")
		Pkg.add("GLM")
		Pkg.add(url="https://github.com/JuliaHealth/DICOM.jl")
		Pkg.add(url="https://github.com/Dale-Black/DICOMUtils.jl")
		Pkg.add(url="https://github.com/Dale-Black/PhantomSegmentation.jl")
		Pkg.add(url="https://github.com/Dale-Black/CalciumScoring.jl")
	end

	using PlutoUI
	using CairoMakie
	using Statistics
	using StatsBase: quantile!
	using ImageMorphology
	using ImageFiltering
	using CSV
	using DataFrames
	using GLM
	using DICOM
	using DICOMUtils
	using PhantomSegmentation
	using CalciumScoring
end

# ╔═╡ 0aeb13cf-6299-473a-908e-5b0be1b23c86
TableOfContents()

# ╔═╡ c4bba726-9f22-4f76-a4ae-966efd13afcf
md"""
## Load DICOMS

All you need to do is set `base_path` once and leave it. After that, the only thing that should change is the `VENDER`, once for every set, and the `SCAN_NUMBER`, once for each scan.
"""

# ╔═╡ 80ea14ec-0128-4899-bcae-5d75d6654a89
md"""
**Everything below should be automatic, just scroll through to visually inspect that things make sense**
"""

# ╔═╡ cb6af1f1-9e54-4172-9e65-bcf534840716
begin
	SCAN_NUMBER = 10
	VENDER = "Canon_Aquilion_One_Vision"
	BASE_PATH = "/Users/daleblack/Google Drive/Datasets/"
end

# ╔═╡ 47504ab4-a52f-418d-b9de-3468304c9578
root_path = string(BASE_PATH, VENDER)

# ╔═╡ d41f610c-928b-4a36-a7f8-ac50ff81f6db
dcm_path_list = dcm_list_builder(root_path)

# ╔═╡ 7efcfa5c-8d99-4e26-9538-efa1d2bfeea6
pth = dcm_path_list[SCAN_NUMBER]

# ╔═╡ 017121d9-a3bd-481a-ba89-7cd261f0ffe7
pth

# ╔═╡ 49fe69bc-7c19-4dc6-86df-ceaccbdd65f2
scan = basename(pth)

# ╔═╡ e95b1cbc-82a1-47db-8026-e2773f0d0cff
header, dcm_array, slice_thick_ori1 = dcm_reader(pth);

# ╔═╡ a4f1ebf0-6f99-42fe-96f3-765caf5f207d
md"""
## Helper Functions
"""

# ╔═╡ 8e874d2a-95bc-414b-91cd-1e6df8220c18
function collect_tuple(tuple_array)
	row_num = size(tuple_array)
	col_num = length(tuple_array[1])
	container = zeros(Int64, row_num..., col_num)
	for i in 1:length(tuple_array)
		container[i,:] = collect(tuple_array[i])
	end
	return container
end

# ╔═╡ 1fec01bc-445a-4961-b518-f25a0de737d3
function overlay_mask_bind(mask)
	indices = findall(x -> x == 1, mask)
	indices = Tuple.(indices)
	label_array = collect_tuple(indices)
	zs = unique(label_array[:,3])
	return PlutoUI.Slider(1:length(zs), default=3, show_value=true)
end

# ╔═╡ 4494460e-4ffc-4317-9a30-21ac8e005bea
function overlay_mask_plot(array, mask, var, title::AbstractString)
	indices = findall(x -> x == 1, mask)
	indices = Tuple.(indices)
	label_array = collect_tuple(indices)
	zs = unique(label_array[:,3])
	indices_lbl = findall(x -> x == zs[var], label_array[:,3])

	fig = Figure()
	ax = Makie.Axis(fig[1, 1])
	ax.title = title
	heatmap!(array[:, :, zs[var]], colormap=:grays)
	scatter!(label_array[:, 1][indices_lbl], label_array[:, 2][indices_lbl], markersize=0.6, color=:red)
	fig
end

# ╔═╡ 9d4357e3-5ee5-4bf9-8879-473a4da59812
md"""
## Segment Heart
"""

# ╔═╡ 4a947165-c565-4c6d-b63c-8cc714d0d039
masked_array, center_insert, mask = mask_heart(header, dcm_array, size(dcm_array, 3)÷2);

# ╔═╡ f1339241-169d-40c1-becc-8f230c12e835
@bind a PlutoUI.Slider(1:size(masked_array, 3), default=10, show_value=true)

# ╔═╡ 641043a6-6d00-4f20-96ad-4b588bda2893
heatmap(masked_array[:, :, a], colormap=:grays)

# ╔═╡ ede1723d-ac8b-43a5-a6d8-6a92fcbaa2fa
begin
	fig = Figure()

	ax = Makie.Axis(fig[1, 1])
	ax.title = "Raw DICOM Array"
	heatmap!(transpose(dcm_array[:, :, 25]), colormap=:grays)
	scatter!(center_insert[2]:center_insert[2]+1, center_insert[1]:center_insert[1]+1, markersize=10, color=:red)
	fig
end

# ╔═╡ 13671de8-9397-422f-81a4-57781e881178
begin
	fig2 = Figure()

	ax2 = Makie.Axis(fig2[1, 1])
	ax2.title = "Mask Array"
	heatmap!(transpose(mask), colormap=:grays)
	scatter!(center_insert[2]:center_insert[2]+1, center_insert[1]:center_insert[1]+1, markersize=10, color=:red)
	fig2
end

# ╔═╡ c07994fb-87c8-4a61-8dd9-7fe50b087b13
begin
	fig3 = Figure()

	ax3 = Makie.Axis(fig3[1, 1])
	ax3.title = "Masked DICOM Array"
	heatmap!(transpose(masked_array[:, :, 23]), colormap=:grays)
	scatter!(center_insert[2]:center_insert[2]+1, center_insert[1]:center_insert[1]+1, markersize=10, color=:red)
	fig3
end

# ╔═╡ 2988d4aa-a2f2-45a5-9e96-fb9c233e9430
md"""
## Segment Calcium Rod
"""

# ╔═╡ 06c46a2f-9f03-46ec-88ed-6219bd2a5f8d
calcium_image, slice_CCI, quality_slice, cal_rod_slice = mask_rod(masked_array, header);

# ╔═╡ 66c50557-49d0-4d9f-beae-404b6784f079
@bind c PlutoUI.Slider(1:size(calcium_image, 3), default=cal_rod_slice, show_value=true)

# ╔═╡ 8806c9d6-c416-47be-981e-a03d35d45602
heatmap(transpose(calcium_image[:, :, c]), colormap=:grays)

# ╔═╡ f30cd88e-c6ca-495d-8ef9-6f029182b5a0
md"""
## Segment Calcium Inserts
"""

# ╔═╡ 7d67fa1f-74de-4a18-9b0d-eb46fa98f420
angle_factor = -4

# ╔═╡ 3bb994fa-e4e8-4e98-bf1b-5e896f21e96e
mask_L_HD, mask_M_HD, mask_S_HD, mask_L_MD, mask_M_MD, mask_S_MD, mask_L_LD, mask_M_LD, mask_S_LD = mask_inserts(
		dcm_array, masked_array, header, slice_CCI, center_insert; angle_factor=angle_factor);

# ╔═╡ 0aece8c7-c8f7-49ca-b68e-b2fda92060c0
masks = mask_L_HD + mask_M_HD + mask_S_HD + mask_L_MD + mask_M_MD + mask_S_MD + mask_L_LD + mask_M_LD + mask_S_LD;

# ╔═╡ 5eb2cb0e-5353-4bd2-86ac-62f17be96f4c
heatmap(masks, colormap=:grays)

# ╔═╡ a488da37-5fd7-4d23-a958-d2fe3f1249b1
md"""
## Calibration Prep
"""

# ╔═╡ d7352bc6-2432-4926-bb70-c791a0050e6b
array_filtered = abs.(mapwindow(median, calcium_image[:, :, cal_rod_slice], (3, 3)));

# ╔═╡ e69626d2-a4eb-46f6-b3ca-491d427b7f70
bool_arr = array_filtered .> 0;

# ╔═╡ 88efbea6-138f-4885-9afb-543648e20152
bool_arr_erode = (((erode(erode(bool_arr)))));

# ╔═╡ a50ca0be-2613-498a-bcb8-5bdf10b4eeed
heatmap(bool_arr, colormap=:grays)

# ╔═╡ a2ea95bc-8840-4fbe-a29d-3618c629e9dc
heatmap(bool_arr_erode, colormap=:grays)

# ╔═╡ 29a2117d-0462-4aac-8ff0-060f1c6ba848
c_img = calcium_image[:, :, cal_rod_slice-1:cal_rod_slice+1];

# ╔═╡ 4c5b5a15-a456-4f67-977d-55a727f6ac7f
begin
	mask_cal_3D = Array{Bool}(undef, size(c_img))
	for z in 1:size(c_img, 3)
		mask_cal_3D[:, :, z] = bool_arr_erode
	end
end;

# ╔═╡ 0362be1a-514b-4133-8c79-b2f554d7ed9f
hist(c_img[mask_cal_3D])

# ╔═╡ af9b9f22-8a65-46f1-8e14-4b8a5191f764
cal_insert_mean = mean(c_img[mask_cal_3D])

# ╔═╡ aec152dd-ff56-49b5-be20-94d63874f624
# cal_insert_mean = quantile!(c_img[mask_cal_3D], 0.7)

# ╔═╡ 5176bb66-ad86-4c36-9ede-06403426a9bb
md"""
### Calibration Line
"""

# ╔═╡ 3d0415d5-5d6a-4ca7-9d79-8876e0618580
density_array_calc = [0, 200] # mg/cc

# ╔═╡ fe9c32a6-1580-4e7d-8d22-a4623b79db2e
intensity_array = [0, cal_insert_mean] # HU

# ╔═╡ 952a3a84-ce2d-4628-8c84-a2fc8f017a96
df_cal = DataFrame(:density => density_array_calc, :intensity => intensity_array)

# ╔═╡ 34ee7ab0-6d33-4d59-bb05-e2aa76a91f66
linearRegressor = lm(@formula(intensity ~ density), df_cal);

# ╔═╡ 620e74c0-7bca-4584-a4a4-7f95831270c0
linearFit = predict(linearRegressor)

# ╔═╡ 6142a256-8b23-495c-89ea-e26a99e93056
m = linearRegressor.model.pp.beta0[2]

# ╔═╡ 884ddd7d-7344-406d-a17f-6946c099fd91
b = linearRegressor.model.rr.mu[1]

# ╔═╡ 03e93f63-bfb4-4fdc-89ca-0a8a22c6658b
density(intensity) = (intensity - b) / m

# ╔═╡ 81be4377-f296-4bbe-86d8-fb98700f5863
intensity(ρ) = m*ρ + b

# ╔═╡ 67caa8b2-7fb7-49a0-b515-48a5da35f989
begin
	f = Figure()
	ax1 = Axis(f[1, 1])

	scatter!(density_array_calc, intensity_array)
	lines!(density_array_calc, linearFit, color = :red)
	ax1.title = "Calibration Line (Intensity vs Density)"
	ax1.ylabel = "Intensity (HU)"
	ax1.xlabel = "Density (mg/cm^3)"

	f
end

# ╔═╡ db84d9bb-e0da-4f81-88c8-c2a20c3ce006
md"""
# Score Large Inserts
"""

# ╔═╡ 707b8724-4f0c-43fd-b797-3badb3c617aa
arr = masked_array[:, :, slice_CCI-2:slice_CCI+2];

# ╔═╡ fe2ab4c6-0b7b-4d07-b84a-663d9d54aea5
single_arr = masked_array[:, :, slice_CCI];

# ╔═╡ 2d83f341-24e4-4f65-885b-e5852e1c5d96
md"""
## High Density
"""

# ╔═╡ 2d1f54af-2b00-49da-8c45-fd7cc3277234
begin
	mask_L_HD_3D = Array{Bool}(undef, size(arr))
	for z in 1:size(arr, 3)
		mask_L_HD_3D[:, :, z] = mask_L_HD
	end
end;

# ╔═╡ 4d08c2b9-b373-4fc4-ad4f-af7621f8e55f
md"""
#### Dilated mask
"""

# ╔═╡ dca4024a-de02-4f2a-a2eb-b458d8089c7d
dilated_mask_L_HD = dilate(dilate(mask_L_HD_3D));

# ╔═╡ 4da96d00-fc05-4b05-b41e-ff6527c2de36
@bind g2 overlay_mask_bind(dilated_mask_L_HD)

# ╔═╡ 258d3494-a618-4db6-84a2-189b24bb15db
overlay_mask_plot(arr, dilated_mask_L_HD, g2, "dilated mask")

# ╔═╡ 337109bb-42b2-4f6f-83f6-e0e98e8838c3
md"""
#### Ring (background) mask
"""

# ╔═╡ 6eece123-6547-4145-8c76-414a699823e6
ring_mask_L_HD = dilate(dilate(dilate(dilate(mask_L_HD_3D)))) - dilate(dilate(dilate(mask_L_HD_3D)));

# ╔═╡ 5faa2499-1017-4182-92bc-cb6f18401617
@bind g4 overlay_mask_bind(ring_mask_L_HD)

# ╔═╡ 6af0d533-6eb4-4129-8b12-a5eefe82b2c0
overlay_mask_plot(arr, ring_mask_L_HD, g4, "ring mask")

# ╔═╡ 45557e45-f916-43a9-a2c5-944b154ee264
begin
	single_ring_mask_L_HD = Bool.(ring_mask_L_HD[:, :, 3])
	s_bkg_L_HD = mean(single_arr[single_ring_mask_L_HD])
end

# ╔═╡ a3236070-89c2-4500-b5de-74dfe4aa4ba5
pixel_size = DICOMUtils.get_pixel_size(header)

# ╔═╡ 1a7c8177-6f86-465e-b0dd-2756ff0db034
S_Obj_HD = intensity(800)

# ╔═╡ 8305356e-9b8e-4ccf-83bd-2258b4bfa2d6
begin
	alg_L_HD = Integrated(arr[mask_L_HD_3D])
	ρ_hd = 0.8 # mg/mm^3
	mass_l_hd = score(s_bkg_L_HD, S_Obj_HD, pixel_size, ρ_hd, alg_L_HD)
end

# ╔═╡ 29d535c8-994d-42e1-ad34-48453a3f92f4
md"""
## Medium Density
"""

# ╔═╡ f8379cf2-ef15-42a3-a82f-63343ef7d9d9
begin
	mask_L_MD_3D = Array{Bool}(undef, size(arr))
	for z in 1:size(arr, 3)
		mask_L_MD_3D[:, :, z] = mask_L_MD
	end
end;

# ╔═╡ ec4c1348-6eac-4372-88b9-7323d20a80f9
md"""
#### Dilated mask
"""

# ╔═╡ f67d19c6-9496-480e-8fc2-35931434bb0b
dilated_mask_L_MD = dilate(dilate(mask_L_MD_3D));

# ╔═╡ 1e5bb04d-ecd8-4036-85e9-05ba522e65fd
@bind h2 overlay_mask_bind(dilated_mask_L_MD)

# ╔═╡ 5b01eec5-b967-4e52-ada2-d7e6f98e8743
overlay_mask_plot(arr, dilated_mask_L_MD, h2, "dilated mask")

# ╔═╡ bd50ca64-f6eb-4827-b8eb-a4186648d3c7
md"""
#### Ring (background) mask
"""

# ╔═╡ 7a920c7e-1859-4237-8499-b00fb530493b
ring_mask_L_MD = dilate(dilate(dilate(dilate(mask_L_MD_3D)))) - dilate(dilate(dilate(mask_L_MD_3D)));

# ╔═╡ 370c4549-69f5-4835-ae55-fc19e61c6c01
@bind h4 overlay_mask_bind(ring_mask_L_MD)

# ╔═╡ 835b7446-c0dc-4c31-a172-d52af2a81b89
overlay_mask_plot(arr, ring_mask_L_MD, h4, "ring mask")

# ╔═╡ d189b76f-d6ca-455a-9432-ba054377557a
begin
	single_ring_mask_L_MD = Bool.(ring_mask_L_MD[:, :, 3])
	s_bkg_L_MD = mean(single_arr[single_ring_mask_L_MD])
end

# ╔═╡ 16e21e05-b582-4986-b2e6-f7e40ea77f39
S_Obj_MD = intensity(400)

# ╔═╡ b7a1a6d7-82a6-40aa-afcb-5323a2642c56
begin
	alg_L_MD = Integrated(arr[mask_L_MD_3D])
	ρ_md = 0.4 # mg/mm^3
	mass_l_md = score(s_bkg_L_MD, S_Obj_MD, pixel_size, ρ_md, alg_L_MD)
end

# ╔═╡ a42e982c-16e5-45b3-a9a5-859375b72f05
md"""
## Low Density
"""

# ╔═╡ ebd8554c-8048-4739-a6bb-39cc73f7ee56
begin
	mask_L_LD_3D = Array{Bool}(undef, size(arr))
	for z in 1:size(arr, 3)
		mask_L_LD_3D[:, :, z] = mask_L_LD
	end
end;

# ╔═╡ 6d2174b3-16ce-4fe2-9ec1-9c945fa6629e
md"""
#### Dilated mask
"""

# ╔═╡ a7caa2fb-b7a3-4b72-bcb9-07ac9475451b
dilated_mask_L_LD = dilate(dilate(mask_L_LD_3D));

# ╔═╡ ba2c6db8-b246-45ba-b667-e7b1883894db
@bind i2 overlay_mask_bind(dilated_mask_L_LD)

# ╔═╡ a05fb79c-3782-460a-a4ee-54b943198e2e
overlay_mask_plot(arr, dilated_mask_L_LD, i2, "dilated mask")

# ╔═╡ fc0f8bd0-40dc-410b-a980-252179ecd3fc
md"""
#### Ring (background) mask
"""

# ╔═╡ 41313ede-e5d4-4df4-9881-a08c10ac227a
ring_mask_L_LD = dilate(dilate(dilate(dilate(mask_L_LD_3D)))) - dilate(dilate(dilate(mask_L_LD_3D)));

# ╔═╡ e9aa29a4-4480-44da-aeb8-1de14cf382b2
@bind i4 overlay_mask_bind(ring_mask_L_LD)

# ╔═╡ c3638503-8a36-4416-8842-88d4c5bda21c
overlay_mask_plot(arr, ring_mask_L_LD, i4, "ring mask")

# ╔═╡ d1e942c5-8ef5-443d-ab5e-083e0d81973e
begin
	single_ring_mask_L_LD = Bool.(ring_mask_L_LD[:, :, 3])
	s_bkg_L_LD = mean(single_arr[single_ring_mask_L_LD])
end

# ╔═╡ 54ba65d1-e7e8-40cb-b016-43686789e0e6
S_Obj_LD = intensity(200)

# ╔═╡ dc424553-ab65-4013-aa1e-0d8163ac0c6f
begin
	alg_L_LD = Integrated(arr[mask_L_LD_3D])
	ρ_ld = 0.2 # mg/mm^3
	mass_l_ld = score(s_bkg_L_LD, S_Obj_LD, pixel_size, ρ_ld, alg_L_LD)
end

# ╔═╡ 24b943da-9319-4360-9d41-992f41752bb8
md"""
# Score Medium Inserts
"""

# ╔═╡ ea1e7093-51c3-4058-ad9b-9a43e67b019f
md"""
## High Density
"""

# ╔═╡ 24636ba1-6ac6-42d9-9167-823cf31444a3
begin
	mask_M_HD_3D = Array{Bool}(undef, size(arr))
	for z in 1:size(arr, 3)
		mask_M_HD_3D[:, :, z] = mask_M_HD
	end
end;

# ╔═╡ fa5e0bef-3991-4b42-921d-60015da8000a
md"""
#### Dilated mask
"""

# ╔═╡ 21f76e13-8748-4eaf-a38e-4fd8a0e9d861
dilated_mask_M_HD = dilate(dilate(dilate(dilate(mask_M_HD_3D))));

# ╔═╡ 344026d1-d9b5-4dd0-bce5-5b668e9545cc
@bind j2 overlay_mask_bind(dilated_mask_M_HD)

# ╔═╡ f6e68cba-833a-467e-af86-658dcc5bf6c7
overlay_mask_plot(arr, dilated_mask_M_HD, j2, "dilated mask")

# ╔═╡ 33bf73da-96c1-4ea3-a6a2-c022c7b64bb9
md"""
#### Ring (background) mask
"""

# ╔═╡ 31d733e5-ca48-4047-9ba1-91b3e6655647
ring_mask_M_HD = dilate(dilate(dilate(dilate(dilate(mask_M_HD_3D))))) - dilate(dilate(dilate(dilate(mask_M_HD_3D))));

# ╔═╡ 894f3319-b750-490e-8e45-bc3710efb0ce
@bind j4 overlay_mask_bind(ring_mask_M_HD)

# ╔═╡ 82acb4b7-85ae-42c7-bf42-fe2ae1abdaf1
overlay_mask_plot(arr, ring_mask_M_HD, j4, "ring mask")

# ╔═╡ 3163aaaa-2463-4664-aa6c-6e44988aa337
begin
	single_ring_mask_M_HD = Bool.(ring_mask_M_HD[:, :, 3])
	s_bkg_M_HD = mean(single_arr[single_ring_mask_M_HD])
end

# ╔═╡ 8af7d835-c43b-40b7-a934-4f6ce7661535
begin
	alg_M_HD = Integrated(arr[mask_M_HD_3D])
	mass_m_hd = score(s_bkg_M_HD, S_Obj_HD, pixel_size, ρ_hd, alg_M_HD)
end

# ╔═╡ 18eae514-f9f8-4cd3-9626-0d702a8fda80
md"""
## Medium Density
"""

# ╔═╡ bbdb568d-c6d4-4eb7-b2d1-785a7df1ed6c
begin
	mask_M_MD_3D = Array{Bool}(undef, size(arr))
	for z in 1:size(arr, 3)
		mask_M_MD_3D[:, :, z] = mask_M_MD
	end
end;

# ╔═╡ 43432b7f-cb36-4fc7-9b2a-97920880004a
md"""
#### Dilated mask
"""

# ╔═╡ 5af97cce-0237-4f41-897a-688f9ed98652
dilated_mask_M_MD = dilate(dilate(dilate(dilate(mask_M_MD_3D))));

# ╔═╡ 5a320923-1966-4a68-b90d-af1da17e79ea
@bind k2 overlay_mask_bind(dilated_mask_M_MD)

# ╔═╡ 439dd6bb-ae38-4cd0-ad89-002c9d91a7c1
overlay_mask_plot(arr, dilated_mask_M_MD, k2, "dilated mask")

# ╔═╡ 63abca9e-2111-4799-be13-ab3f86d5b082
md"""
#### Ring (background) mask
"""

# ╔═╡ 7e7b5ad2-a944-4439-8361-0c88d85b9985
ring_mask_M_MD = dilate(dilate(dilate(dilate(dilate(dilate(mask_M_MD_3D)))))) - dilate(dilate(dilate(dilate(dilate(mask_M_MD_3D)))));

# ╔═╡ 495cc420-a24f-41d2-a0cb-116f778c0bff
@bind k4 overlay_mask_bind(ring_mask_M_MD)

# ╔═╡ 2e52ff34-db03-404c-ad98-c8f62df0d181
overlay_mask_plot(arr, ring_mask_M_MD, k4, "ring mask")

# ╔═╡ 0687ccfd-d58a-4fa4-aa00-815101997a75
begin
	single_ring_mask_M_MD = Bool.(ring_mask_M_MD[:, :, 3])
	s_bkg_M_MD = mean(single_arr[single_ring_mask_M_MD])
end

# ╔═╡ 980785a4-0c19-4208-89eb-2b2b4e85523d
begin
	alg_M_MD = Integrated(arr[mask_M_MD_3D])
	mass_m_md = score(s_bkg_M_MD, S_Obj_MD, pixel_size, ρ_md, alg_M_MD)
end

# ╔═╡ 48400ead-acd0-4c5b-8ff7-87ee2b07b9d1
md"""
## Low Density
"""

# ╔═╡ 7d68622e-6a35-4213-a74b-5c0736995579
begin
	mask_M_LD_3D = Array{Bool}(undef, size(arr))
	for z in 1:size(arr, 3)
		mask_M_LD_3D[:, :, z] = mask_M_LD
	end
end;

# ╔═╡ c8a9335b-c271-4097-be57-455119cd461e
md"""
#### Dilated mask
"""

# ╔═╡ 341c1d34-4b2d-4a10-ab04-da5d0dbd579f
dilated_mask_M_LD = dilate(dilate(dilate(dilate(dilate(mask_M_LD_3D)))));

# ╔═╡ 223efc65-64a0-4d27-b829-3288445967d4
@bind l2 overlay_mask_bind(dilated_mask_M_LD)

# ╔═╡ e2a2251d-9964-4fa3-becc-09d0b12cd7f0
overlay_mask_plot(arr, dilated_mask_M_LD, l2, "dilated mask")

# ╔═╡ f7c88984-fcaf-4173-8f8b-023979aa472f
md"""
#### Ring (background) mask
"""

# ╔═╡ 3f8f91db-d916-4bdc-a678-fdab41a9d05e
ring_mask_M_LD = dilate(dilate(dilate(dilate(dilate(dilate(mask_M_LD_3D)))))) - dilate(dilate(dilate(dilate(dilate(mask_M_LD_3D)))));

# ╔═╡ b03ba157-8b69-4616-8723-f0396bf184c6
@bind l4 overlay_mask_bind(ring_mask_M_LD)

# ╔═╡ f58e9581-8f7d-477e-8820-9c13e210b385
overlay_mask_plot(arr, ring_mask_M_LD, l4, "ring mask")

# ╔═╡ ceecd152-e1f6-497a-872b-c4de20345499
begin
	single_ring_mask_M_LD = Bool.(ring_mask_M_LD[:, :, 3])
	s_bkg_M_LD = mean(single_arr[single_ring_mask_M_LD])
end

# ╔═╡ 1a280788-84f9-4a3d-9084-b1ab9c2ba514
begin
	alg_M_LD = Integrated(arr[mask_M_LD_3D])
	mass_m_ld = score(s_bkg_M_LD, S_Obj_LD, pixel_size, ρ_ld, alg_M_LD)
end

# ╔═╡ 16d56ff5-238e-4916-bac6-56fe70164b54
md"""
# Score Small Inserts
"""

# ╔═╡ 37acb2c8-30e2-441d-90d2-884ea0fd2d49
md"""
## High Density
"""

# ╔═╡ d2d48580-177b-4f67-9bbf-022ce40f4771
begin
	mask_S_HD_3D = Array{Bool}(undef, size(arr))
	for z in 1:size(arr, 3)
		mask_S_HD_3D[:, :, z] = mask_S_HD
	end
end;

# ╔═╡ fdbf5185-6a1f-440b-aa81-81529e045b73
md"""
#### Dilated mask
"""

# ╔═╡ 62b45cc8-af83-48a1-bb33-73335a6081bf
dilated_mask_S_HD = dilate(dilate(dilate(dilate(dilate((mask_S_HD_3D))))));

# ╔═╡ b1a6f986-4d86-4c62-a662-552bdf9421f9
@bind m2 overlay_mask_bind(dilated_mask_S_HD)

# ╔═╡ de05e5fd-3972-41e9-bdaa-eba6631d8ecb
overlay_mask_plot(arr, dilated_mask_S_HD, m2, "dilated mask")

# ╔═╡ 339617d5-49c4-4bad-8c03-579438b55799
md"""
#### Ring (background) mask
"""

# ╔═╡ cb248de5-bcfd-4ce8-ba9b-c0c39067bbc3
ring_mask_S_HD = dilate(dilate(dilate(dilate(dilate(mask_S_HD_3D))))) - dilate(dilate(dilate(dilate(mask_S_HD_3D))));

# ╔═╡ b712ebf8-bf99-40c8-b861-19e6dfd968c0
@bind m4 overlay_mask_bind(ring_mask_S_HD)

# ╔═╡ d7bbcda2-37bf-43c2-9a4e-e6027d868d26
overlay_mask_plot(arr, ring_mask_S_HD, m4, "ring mask")

# ╔═╡ 10550a99-76b9-487e-922b-3239e7f811f4
begin
	single_ring_mask_S_HD = Bool.(ring_mask_S_HD[:, :, 3])
	s_bkg_S_HD = mean(single_arr[single_ring_mask_S_HD])
end

# ╔═╡ 933791e4-934d-40fd-8928-eccb68273059
begin
	alg_S_HD = Integrated(arr[mask_S_HD_3D])
	mass_s_hd = score(s_bkg_S_HD, S_Obj_HD, pixel_size, ρ_hd, alg_S_HD)
end

# ╔═╡ cdc09f53-4f83-46ff-8124-bdf54cea5097
md"""
## Medium Density
"""

# ╔═╡ a71cac4e-9ade-4ffc-a4f0-b00c9345221b
begin
	mask_S_MD_3D = Array{Bool}(undef, size(arr))
	for z in 1:size(arr, 3)
		mask_S_MD_3D[:, :, z] = mask_S_MD
	end
end;

# ╔═╡ d27098ee-43ae-4072-9843-11934f377df6
md"""
#### Dilated mask
"""

# ╔═╡ fbbd86c5-0acf-496a-b6ee-3dcb090435db
dilated_mask_S_MD = dilate(dilate(dilate(dilate(dilate(mask_S_MD_3D)))));

# ╔═╡ 6467393e-8590-45db-924b-0033e9ddb724
@bind n2 overlay_mask_bind(dilated_mask_S_MD)

# ╔═╡ a85ad33c-4d57-417d-b223-c3566f87c9ff
overlay_mask_plot(arr, dilated_mask_S_MD, n2, "dilated mask")

# ╔═╡ f1057c02-9060-4d44-a0d9-6b04245c9ed0
md"""
#### Ring (background) mask
"""

# ╔═╡ c32ec1cb-825b-477b-927d-7d72b1deb99d
ring_mask_S_MD = dilate(dilate(dilate(dilate(dilate(mask_S_MD_3D))))) - dilate(dilate(dilate(dilate(mask_S_MD_3D))));

# ╔═╡ 9b254664-72a6-4030-83b0-27db11c222ba
@bind n4 overlay_mask_bind(ring_mask_S_MD)

# ╔═╡ 60370f6d-51f3-436f-ab32-400de8d98f2b
overlay_mask_plot(arr, ring_mask_S_MD, n4, "ring mask")

# ╔═╡ d62d50a8-dd17-4d25-b353-fdfedb456ff0
begin
	single_ring_mask_S_MD = Bool.(ring_mask_S_MD[:, :, 3])
	s_bkg_S_MD = mean(single_arr[single_ring_mask_S_MD])
end

# ╔═╡ a77205bb-eab2-4a2a-90b7-2f4844f039de
begin
	alg_S_MD = Integrated(arr[mask_S_MD_3D])
	mass_s_md = score(s_bkg_S_MD, S_Obj_MD, pixel_size, ρ_md, alg_S_MD)
end

# ╔═╡ b3455451-651e-4f70-87f4-08f9db13d9b5
md"""
## Low Density
"""

# ╔═╡ 7d71e326-ab38-4ad3-9e20-abed47caf9e9
begin
	mask_S_LD_3D = Array{Bool}(undef, size(arr))
	for z in 1:size(arr, 3)
		mask_S_LD_3D[:, :, z] = mask_S_LD
	end
end;

# ╔═╡ 35b79339-bb2e-4987-a0f4-583a732c71fd
md"""
#### Dilated mask
"""

# ╔═╡ e0861340-3bef-4dd4-a21c-b0151d0f6d02
dilated_mask_S_LD = dilate(dilate(dilate(dilate(dilate(mask_S_LD_3D)))));

# ╔═╡ e35a7298-79fa-42b3-a2c7-faec6b5c9af4
@bind o2 overlay_mask_bind(dilated_mask_S_LD)

# ╔═╡ 86fac660-6ef4-4041-80d8-61eb6edc2a39
overlay_mask_plot(arr, dilated_mask_S_LD, o2, "dilated mask")

# ╔═╡ 005f8151-6b2b-4476-b92d-e63d0d216325
md"""
#### Ring (background) mask
"""

# ╔═╡ 8c379977-38dd-41e8-8956-05f48102cc67
ring_mask_S_LD = dilate(dilate(dilate(dilate(dilate(mask_S_LD_3D))))) - dilate(dilate(dilate(dilate(mask_S_LD_3D))));

# ╔═╡ b71e334a-8452-453c-a7a4-6d5f6eafd0bc
@bind o4 overlay_mask_bind(ring_mask_S_LD)

# ╔═╡ 13806613-6eca-4b9b-8216-894b25ff1b01
overlay_mask_plot(arr, ring_mask_S_LD, o4, "ring mask")

# ╔═╡ 9a8f53e1-4900-4b64-a09f-5a9b8fd5a0c5
begin
	single_ring_mask_S_LD = Bool.(ring_mask_S_LD[:, :, 3])
	s_bkg_S_LD = mean(single_arr[single_ring_mask_S_LD])
end

# ╔═╡ 05233816-9613-4cb6-90f0-bac0f01816c0
begin
	alg_S_LD = Integrated(arr[mask_S_LD_3D])
	mass_s_ld = score(s_bkg_S_LD, S_Obj_LD, pixel_size, ρ_ld, alg_S_LD)
end

# ╔═╡ b7f35164-9e08-49ea-aa40-f968fbbb92c2
md"""
# Results
"""

# ╔═╡ b5e14014-12bf-4004-bccc-2f3c759200b5
density_array = [0, 200, 400, 800]

# ╔═╡ 88edf2c7-3c5a-4b90-9c50-5d448a358ce8
inserts = [
	"Low Density",
	"Medium Density",
	"High Density"
]

# ╔═╡ 786eb037-e39e-41e1-9a63-6106800cb306
ground_truth_mass_large = [
	19.6,
	39.3,
	78.5
] # mg

# ╔═╡ fe51b3b3-31bc-454f-bc09-db5a7e6c6941
calculated_mass_large = [
	mass_l_ld,
	mass_l_md,
	mass_l_hd
]

# ╔═╡ 2867365f-f7fe-460c-b8fd-eb8640f7dc73
ground_truth_mass_medium = [
	4.2,
	8.5,
	17.0
]

# ╔═╡ 61cffd44-a5e8-4223-85ef-eb4f2b22de2b
calculated_mass_medium = [
	mass_m_ld,
	mass_m_md,
	mass_m_hd
]

# ╔═╡ 60324b4b-e880-4ea3-9831-fca775eb22ab
ground_truth_mass_small = [
	0.2,
	0.3,
	0.6
]

# ╔═╡ 7c4fe694-7343-4e87-aed4-8895802e34a9
calculated_mass_small = [
	mass_s_ld,
	mass_s_md,
	mass_s_hd
]

# ╔═╡ 9f12b7bd-ccc2-4108-9171-1902a0c75c30
df = DataFrame(
	scan = scan,
	inserts = inserts,
	ground_truth_mass_large = ground_truth_mass_large,
	calculated_mass_large = calculated_mass_large,
	ground_truth_mass_medium = ground_truth_mass_medium,
	calculated_mass_medium = calculated_mass_medium,
	ground_truth_mass_small = ground_truth_mass_small,
	calculated_mass_small = calculated_mass_small
)

# ╔═╡ fef2fcda-cc64-469e-9f17-946f06a3e4ae
begin
	fmass2 = Figure()
	axmass2 = Axis(fmass2[1, 1])

	scatter!(density_array[2:end], df[!, :ground_truth_mass_large], label="ground_truth_mass_large")
	scatter!(density_array[2:end], df[!, :calculated_mass_large], label="calculated_mass_large")

	axmass2.title = "Mass Measurements (Large)"
	axmass2.ylabel = "Mass (mg)"
	axmass2.xlabel = "Density (mg/cm^3)"

	xlims!(axmass2, 0, 850)
	ylims!(axmass2, 0, 100)

	fmass2[1, 2] = Legend(fmass2, axmass2, framevisible = false)

	fmass2
end

# ╔═╡ 4f9c5043-3485-45b8-b493-d62780134251
begin
	fmass3 = Figure()
	axmass3 = Axis(fmass3[1, 1])

	scatter!(density_array[2:end], df[!, :ground_truth_mass_medium], label="ground_truth_mass_medium")
	scatter!(density_array[2:end], df[!, :calculated_mass_medium], label="calculated_mass_medium")

	axmass3.title = "Mass Measurements (Medium)"
	axmass3.ylabel = "Mass (mg)"
	axmass3.xlabel = "Density (mg/cm^3)"

	xlims!(axmass3, 0, 850)
	ylims!(axmass3, 0, 25)

	fmass3[1, 2] = Legend(fmass3, axmass3, framevisible = false)

	fmass3
end

# ╔═╡ 197e0479-f15a-43f9-a7cf-5b5d58b5a901
begin
	fmass4 = Figure()
	axmass4 = Axis(fmass4[1, 1])

	scatter!(density_array[2:end], df[!, :ground_truth_mass_small], label="ground_truth_mass_small")
	scatter!(density_array[2:end], df[!, :calculated_mass_small], label="calculated_mass_small")

	axmass4.title = "Mass Measurements (Small)"
	axmass4.ylabel = "Mass (mg)"
	axmass4.xlabel = "Density (mg/cm^3)"

	xlims!(axmass4, 0, 850)
	ylims!(axmass4, 0, 1.5)

	fmass4[1, 2] = Legend(fmass4, axmass4, framevisible = false)

	fmass4
end

# ╔═╡ b030dcf1-fc55-4e32-be9b-6cba1d3456f0
md"""
### Save Results
"""

# ╔═╡ 5a9c7e99-0e6e-4162-8ff7-bff1935fd318
if ~isdir(string(cd(pwd, "..") , "/data/output/", VENDER))
	mkdir(string(cd(pwd, "..") , "/data/output/", VENDER))
end

# ╔═╡ 7976fbad-ae1f-40e8-961f-2ea68f752caa
output_path = string(cd(pwd, "..") , "/data/output/", VENDER, "/", scan, ".csv")

# ╔═╡ d8aeb68a-697d-4477-8b6c-14712f50a945
# CSV.write(output_path, df)

# ╔═╡ 3d6114a2-1356-4998-be56-800395262bc0
md"""
### Save full df
"""

# ╔═╡ cfca657a-da63-436d-b316-a4e9a63ae8e7
dfs = []

# ╔═╡ a2bf4e79-7f23-45f8-8261-8c10bcd7eb3c
push!(dfs, df)

# ╔═╡ b1a053ad-f18f-487d-8c38-01dbe09fd03a
if length(dfs) == 10
	global new_df = vcat(dfs[1:10]...)
	output_path_new = string(cd(pwd, "..") , "/data/output/", VENDER, "/full.csv")
	CSV.write(output_path_new, new_df)
end

# ╔═╡ Cell order:
# ╠═3c162d31-407a-4c50-a8a7-30153192f207
# ╠═0aeb13cf-6299-473a-908e-5b0be1b23c86
# ╟─c4bba726-9f22-4f76-a4ae-966efd13afcf
# ╟─80ea14ec-0128-4899-bcae-5d75d6654a89
# ╠═cb6af1f1-9e54-4172-9e65-bcf534840716
# ╠═47504ab4-a52f-418d-b9de-3468304c9578
# ╠═d41f610c-928b-4a36-a7f8-ac50ff81f6db
# ╠═7efcfa5c-8d99-4e26-9538-efa1d2bfeea6
# ╠═017121d9-a3bd-481a-ba89-7cd261f0ffe7
# ╠═49fe69bc-7c19-4dc6-86df-ceaccbdd65f2
# ╠═e95b1cbc-82a1-47db-8026-e2773f0d0cff
# ╟─a4f1ebf0-6f99-42fe-96f3-765caf5f207d
# ╟─8e874d2a-95bc-414b-91cd-1e6df8220c18
# ╟─1fec01bc-445a-4961-b518-f25a0de737d3
# ╟─4494460e-4ffc-4317-9a30-21ac8e005bea
# ╟─9d4357e3-5ee5-4bf9-8879-473a4da59812
# ╠═4a947165-c565-4c6d-b63c-8cc714d0d039
# ╟─f1339241-169d-40c1-becc-8f230c12e835
# ╠═641043a6-6d00-4f20-96ad-4b588bda2893
# ╠═ede1723d-ac8b-43a5-a6d8-6a92fcbaa2fa
# ╟─13671de8-9397-422f-81a4-57781e881178
# ╟─c07994fb-87c8-4a61-8dd9-7fe50b087b13
# ╟─2988d4aa-a2f2-45a5-9e96-fb9c233e9430
# ╠═06c46a2f-9f03-46ec-88ed-6219bd2a5f8d
# ╟─66c50557-49d0-4d9f-beae-404b6784f079
# ╠═8806c9d6-c416-47be-981e-a03d35d45602
# ╟─f30cd88e-c6ca-495d-8ef9-6f029182b5a0
# ╠═7d67fa1f-74de-4a18-9b0d-eb46fa98f420
# ╠═3bb994fa-e4e8-4e98-bf1b-5e896f21e96e
# ╠═0aece8c7-c8f7-49ca-b68e-b2fda92060c0
# ╠═5eb2cb0e-5353-4bd2-86ac-62f17be96f4c
# ╟─a488da37-5fd7-4d23-a958-d2fe3f1249b1
# ╠═d7352bc6-2432-4926-bb70-c791a0050e6b
# ╠═e69626d2-a4eb-46f6-b3ca-491d427b7f70
# ╠═88efbea6-138f-4885-9afb-543648e20152
# ╠═a50ca0be-2613-498a-bcb8-5bdf10b4eeed
# ╠═a2ea95bc-8840-4fbe-a29d-3618c629e9dc
# ╠═29a2117d-0462-4aac-8ff0-060f1c6ba848
# ╠═4c5b5a15-a456-4f67-977d-55a727f6ac7f
# ╠═0362be1a-514b-4133-8c79-b2f554d7ed9f
# ╠═af9b9f22-8a65-46f1-8e14-4b8a5191f764
# ╠═aec152dd-ff56-49b5-be20-94d63874f624
# ╟─5176bb66-ad86-4c36-9ede-06403426a9bb
# ╠═3d0415d5-5d6a-4ca7-9d79-8876e0618580
# ╠═fe9c32a6-1580-4e7d-8d22-a4623b79db2e
# ╠═952a3a84-ce2d-4628-8c84-a2fc8f017a96
# ╠═34ee7ab0-6d33-4d59-bb05-e2aa76a91f66
# ╠═620e74c0-7bca-4584-a4a4-7f95831270c0
# ╠═6142a256-8b23-495c-89ea-e26a99e93056
# ╠═884ddd7d-7344-406d-a17f-6946c099fd91
# ╠═03e93f63-bfb4-4fdc-89ca-0a8a22c6658b
# ╠═81be4377-f296-4bbe-86d8-fb98700f5863
# ╟─67caa8b2-7fb7-49a0-b515-48a5da35f989
# ╟─db84d9bb-e0da-4f81-88c8-c2a20c3ce006
# ╠═707b8724-4f0c-43fd-b797-3badb3c617aa
# ╠═fe2ab4c6-0b7b-4d07-b84a-663d9d54aea5
# ╟─2d83f341-24e4-4f65-885b-e5852e1c5d96
# ╠═2d1f54af-2b00-49da-8c45-fd7cc3277234
# ╟─4d08c2b9-b373-4fc4-ad4f-af7621f8e55f
# ╠═dca4024a-de02-4f2a-a2eb-b458d8089c7d
# ╟─4da96d00-fc05-4b05-b41e-ff6527c2de36
# ╠═258d3494-a618-4db6-84a2-189b24bb15db
# ╟─337109bb-42b2-4f6f-83f6-e0e98e8838c3
# ╠═6eece123-6547-4145-8c76-414a699823e6
# ╟─5faa2499-1017-4182-92bc-cb6f18401617
# ╠═6af0d533-6eb4-4129-8b12-a5eefe82b2c0
# ╠═45557e45-f916-43a9-a2c5-944b154ee264
# ╠═a3236070-89c2-4500-b5de-74dfe4aa4ba5
# ╠═1a7c8177-6f86-465e-b0dd-2756ff0db034
# ╠═8305356e-9b8e-4ccf-83bd-2258b4bfa2d6
# ╟─29d535c8-994d-42e1-ad34-48453a3f92f4
# ╠═f8379cf2-ef15-42a3-a82f-63343ef7d9d9
# ╟─ec4c1348-6eac-4372-88b9-7323d20a80f9
# ╠═f67d19c6-9496-480e-8fc2-35931434bb0b
# ╟─1e5bb04d-ecd8-4036-85e9-05ba522e65fd
# ╠═5b01eec5-b967-4e52-ada2-d7e6f98e8743
# ╟─bd50ca64-f6eb-4827-b8eb-a4186648d3c7
# ╠═7a920c7e-1859-4237-8499-b00fb530493b
# ╟─370c4549-69f5-4835-ae55-fc19e61c6c01
# ╠═835b7446-c0dc-4c31-a172-d52af2a81b89
# ╠═d189b76f-d6ca-455a-9432-ba054377557a
# ╠═16e21e05-b582-4986-b2e6-f7e40ea77f39
# ╠═b7a1a6d7-82a6-40aa-afcb-5323a2642c56
# ╟─a42e982c-16e5-45b3-a9a5-859375b72f05
# ╠═ebd8554c-8048-4739-a6bb-39cc73f7ee56
# ╟─6d2174b3-16ce-4fe2-9ec1-9c945fa6629e
# ╠═a7caa2fb-b7a3-4b72-bcb9-07ac9475451b
# ╟─ba2c6db8-b246-45ba-b667-e7b1883894db
# ╠═a05fb79c-3782-460a-a4ee-54b943198e2e
# ╟─fc0f8bd0-40dc-410b-a980-252179ecd3fc
# ╠═41313ede-e5d4-4df4-9881-a08c10ac227a
# ╟─e9aa29a4-4480-44da-aeb8-1de14cf382b2
# ╠═c3638503-8a36-4416-8842-88d4c5bda21c
# ╠═d1e942c5-8ef5-443d-ab5e-083e0d81973e
# ╠═54ba65d1-e7e8-40cb-b016-43686789e0e6
# ╠═dc424553-ab65-4013-aa1e-0d8163ac0c6f
# ╟─24b943da-9319-4360-9d41-992f41752bb8
# ╟─ea1e7093-51c3-4058-ad9b-9a43e67b019f
# ╠═24636ba1-6ac6-42d9-9167-823cf31444a3
# ╟─fa5e0bef-3991-4b42-921d-60015da8000a
# ╠═21f76e13-8748-4eaf-a38e-4fd8a0e9d861
# ╟─344026d1-d9b5-4dd0-bce5-5b668e9545cc
# ╠═f6e68cba-833a-467e-af86-658dcc5bf6c7
# ╟─33bf73da-96c1-4ea3-a6a2-c022c7b64bb9
# ╠═31d733e5-ca48-4047-9ba1-91b3e6655647
# ╟─894f3319-b750-490e-8e45-bc3710efb0ce
# ╠═82acb4b7-85ae-42c7-bf42-fe2ae1abdaf1
# ╠═3163aaaa-2463-4664-aa6c-6e44988aa337
# ╠═8af7d835-c43b-40b7-a934-4f6ce7661535
# ╟─18eae514-f9f8-4cd3-9626-0d702a8fda80
# ╠═bbdb568d-c6d4-4eb7-b2d1-785a7df1ed6c
# ╟─43432b7f-cb36-4fc7-9b2a-97920880004a
# ╠═5af97cce-0237-4f41-897a-688f9ed98652
# ╟─5a320923-1966-4a68-b90d-af1da17e79ea
# ╠═439dd6bb-ae38-4cd0-ad89-002c9d91a7c1
# ╟─63abca9e-2111-4799-be13-ab3f86d5b082
# ╠═7e7b5ad2-a944-4439-8361-0c88d85b9985
# ╟─495cc420-a24f-41d2-a0cb-116f778c0bff
# ╠═2e52ff34-db03-404c-ad98-c8f62df0d181
# ╠═0687ccfd-d58a-4fa4-aa00-815101997a75
# ╠═980785a4-0c19-4208-89eb-2b2b4e85523d
# ╟─48400ead-acd0-4c5b-8ff7-87ee2b07b9d1
# ╠═7d68622e-6a35-4213-a74b-5c0736995579
# ╟─c8a9335b-c271-4097-be57-455119cd461e
# ╠═341c1d34-4b2d-4a10-ab04-da5d0dbd579f
# ╟─223efc65-64a0-4d27-b829-3288445967d4
# ╠═e2a2251d-9964-4fa3-becc-09d0b12cd7f0
# ╟─f7c88984-fcaf-4173-8f8b-023979aa472f
# ╠═3f8f91db-d916-4bdc-a678-fdab41a9d05e
# ╟─b03ba157-8b69-4616-8723-f0396bf184c6
# ╠═f58e9581-8f7d-477e-8820-9c13e210b385
# ╠═ceecd152-e1f6-497a-872b-c4de20345499
# ╠═1a280788-84f9-4a3d-9084-b1ab9c2ba514
# ╟─16d56ff5-238e-4916-bac6-56fe70164b54
# ╟─37acb2c8-30e2-441d-90d2-884ea0fd2d49
# ╠═d2d48580-177b-4f67-9bbf-022ce40f4771
# ╟─fdbf5185-6a1f-440b-aa81-81529e045b73
# ╠═62b45cc8-af83-48a1-bb33-73335a6081bf
# ╟─b1a6f986-4d86-4c62-a662-552bdf9421f9
# ╠═de05e5fd-3972-41e9-bdaa-eba6631d8ecb
# ╟─339617d5-49c4-4bad-8c03-579438b55799
# ╠═cb248de5-bcfd-4ce8-ba9b-c0c39067bbc3
# ╟─b712ebf8-bf99-40c8-b861-19e6dfd968c0
# ╠═d7bbcda2-37bf-43c2-9a4e-e6027d868d26
# ╠═10550a99-76b9-487e-922b-3239e7f811f4
# ╠═933791e4-934d-40fd-8928-eccb68273059
# ╟─cdc09f53-4f83-46ff-8124-bdf54cea5097
# ╠═a71cac4e-9ade-4ffc-a4f0-b00c9345221b
# ╟─d27098ee-43ae-4072-9843-11934f377df6
# ╠═fbbd86c5-0acf-496a-b6ee-3dcb090435db
# ╟─6467393e-8590-45db-924b-0033e9ddb724
# ╠═a85ad33c-4d57-417d-b223-c3566f87c9ff
# ╟─f1057c02-9060-4d44-a0d9-6b04245c9ed0
# ╠═c32ec1cb-825b-477b-927d-7d72b1deb99d
# ╟─9b254664-72a6-4030-83b0-27db11c222ba
# ╠═60370f6d-51f3-436f-ab32-400de8d98f2b
# ╠═d62d50a8-dd17-4d25-b353-fdfedb456ff0
# ╠═a77205bb-eab2-4a2a-90b7-2f4844f039de
# ╟─b3455451-651e-4f70-87f4-08f9db13d9b5
# ╠═7d71e326-ab38-4ad3-9e20-abed47caf9e9
# ╟─35b79339-bb2e-4987-a0f4-583a732c71fd
# ╠═e0861340-3bef-4dd4-a21c-b0151d0f6d02
# ╟─e35a7298-79fa-42b3-a2c7-faec6b5c9af4
# ╠═86fac660-6ef4-4041-80d8-61eb6edc2a39
# ╟─005f8151-6b2b-4476-b92d-e63d0d216325
# ╠═8c379977-38dd-41e8-8956-05f48102cc67
# ╟─b71e334a-8452-453c-a7a4-6d5f6eafd0bc
# ╠═13806613-6eca-4b9b-8216-894b25ff1b01
# ╠═9a8f53e1-4900-4b64-a09f-5a9b8fd5a0c5
# ╠═05233816-9613-4cb6-90f0-bac0f01816c0
# ╟─b7f35164-9e08-49ea-aa40-f968fbbb92c2
# ╠═b5e14014-12bf-4004-bccc-2f3c759200b5
# ╠═88edf2c7-3c5a-4b90-9c50-5d448a358ce8
# ╠═786eb037-e39e-41e1-9a63-6106800cb306
# ╠═fe51b3b3-31bc-454f-bc09-db5a7e6c6941
# ╠═2867365f-f7fe-460c-b8fd-eb8640f7dc73
# ╠═61cffd44-a5e8-4223-85ef-eb4f2b22de2b
# ╠═60324b4b-e880-4ea3-9831-fca775eb22ab
# ╠═7c4fe694-7343-4e87-aed4-8895802e34a9
# ╠═9f12b7bd-ccc2-4108-9171-1902a0c75c30
# ╟─fef2fcda-cc64-469e-9f17-946f06a3e4ae
# ╟─4f9c5043-3485-45b8-b493-d62780134251
# ╟─197e0479-f15a-43f9-a7cf-5b5d58b5a901
# ╟─b030dcf1-fc55-4e32-be9b-6cba1d3456f0
# ╠═5a9c7e99-0e6e-4162-8ff7-bff1935fd318
# ╠═7976fbad-ae1f-40e8-961f-2ea68f752caa
# ╠═d8aeb68a-697d-4477-8b6c-14712f50a945
# ╟─3d6114a2-1356-4998-be56-800395262bc0
# ╠═cfca657a-da63-436d-b316-a4e9a63ae8e7
# ╠═a2bf4e79-7f23-45f8-8261-8c10bcd7eb3c
# ╠═b1a053ad-f18f-487d-8c38-01dbe09fd03a

max_stars_repo_name: MolloiLab/CAC-stanford-data

max_issues_repo_name: MolloiLab/CAC-stanford-data

max_forks_repo_name: MolloiLab/CAC-stanford-data