text:
using StructC14N
using Test

# Template is a NamedTuple
template = (xrange=NTuple{2, Number},
            yrange=NTuple{2, Number},
            title="Default string")

c = canonicalize(template, (xr=(1,2),))
@test c.xrange == (1, 2)
@test ismissing(c.yrange)
@test c.title == "Default string"

c = canonicalize(template, (xr=(1,2), yrange=(3, 4.), tit="Foo"))
@test c.xrange == (1, 2)
@test c.yrange == (3, 4.)
@test c.title == "Foo"

#  ...input is a Tuple
@test_throws AssertionError canonicalize(template, ((1,2)))
c = canonicalize(template, ((1,2), (3, 4.), "Foo"))
@test c.xrange == (1, 2)
@test c.yrange == (3, 4.)
@test c.title == "Foo"

#  ...inputs are given as keywords
c = canonicalize(template, xr=(1,2), yrange=(3, 4.), tit="Foo")
@test c.xrange == (1, 2)
@test c.yrange == (3, 4.)
@test c.title == "Foo"

# Template is a structure definition
struct AStruct
    xrange::NTuple{2, Number}
    yrange::NTuple{2, Number}
    title::Union{Missing, String}
end

@test_throws MethodError canonicalize(AStruct, (xr=(1,2), tit="Foo"))

c = canonicalize(AStruct, (xr=(1,2), yrang=(3, 4.)))
@test c.xrange == (1, 2)
@test c.yrange == (3, 4.)
@test ismissing(c.title)

c = canonicalize(AStruct, (xr=(1,2), yrang=(3, 4.), tit="Foo"))
@test c.xrange == (1, 2)
@test c.yrange == (3, 4.)
@test c.title == "Foo"

#  ...input is a Tuple
@test_throws AssertionError canonicalize(AStruct, ((1,2)))
c = canonicalize(AStruct, ((1,2), (3, 4.), "Foo"))
@test c.xrange == (1, 2)
@test c.yrange == (3, 4.)
@test c.title == "Foo"

#  ...inputs are given as keywords
c = canonicalize(AStruct, xr=(1,2), yrang=(3, 4.), tit="Foo")
@test c.xrange == (1, 2)
@test c.yrange == (3, 4.)
@test c.title == "Foo"

# Template is a structure instance
template = AStruct((0,0), (0.,0.), missing)
c = canonicalize(template, (xr=(1,2),))
@test c.xrange == (1, 2)
@test c.yrange == (0., 0.)
@test ismissing(c.title)

c = canonicalize(template, (yr=(1,2), xr=(3.3, 4.4), tit="Foo"))
@test c.yrange == (1, 2)
@test c.xrange == (3.3, 4.4)
@test c.title == "Foo"

#  ...input is a Tuple
@test_throws AssertionError canonicalize(template, ((1,2)))
c = canonicalize(template, ((1,2), (3, 4.), "Foo"))
@test c.xrange == (1, 2)
@test c.yrange == (3, 4.)
@test c.title == "Foo"

#  ...inputs are given as keywords
c = canonicalize(template, yr=(1,2), xr=(3.3, 4.4), tit="Foo")
@test c.yrange == (1, 2)
@test c.xrange == (3.3, 4.4)
@test c.title == "Foo"

configtemplate = (optStr=String,
                  optInt=Int,
                  optFloat=Float64)

configentry = "aa, 1, 2"
c = canonicalize(configtemplate, (split(configentry, ",")...,))
@test c.optStr == "aa"
@test c.optInt == 1
@test c.optFloat == 2.0

configentry = "optFloat=20, optStr=\"aaa\", optInt=10"
c = canonicalize(configtemplate, eval(Meta.parse("($configentry)")))
@test c.optStr == "aaa"
@test c.optInt == 10
@test c.optFloat == 20.0

function myparse(input)
    if input == "ten"
        return 10
    end
    return 1
end
configentry = "optFloat=20, optStr=\"aaa\", optInt=\"ten\""
c = canonicalize(configtemplate, eval(Meta.parse("($configentry)")), Dict(:optInt=>myparse))
@test c.optStr == "aaa"
@test c.optInt == 10
@test c.optFloat == 20.0


text:
# Parts of this code were taken / derived from Graphs.jl. See LICENSE for
# licensing details.

# A* shortest-path algorithm

function a_star_impl!(
    g::AbstractGraph,# the graph
    t::Integer, # the end vertex
    frontier,               # an initialized heap containing the active vertices
    colormap::Vector{Int},  # an (initialized) color-map to indicate status of vertices
    distmx::AbstractMatrix,
    heuristic::Function    # heuristic fn (under)estimating distance to target
    )

    while !isempty(frontier)
        (cost_so_far, path, u) = dequeue!(frontier)
        if u == t
            return path
        end

        for v in LightGraphs.out_neighbors(g, u)

            if colormap[v] < 2
                dist = distmx[u, v]
                colormap[v] = 1
                new_path = cat(1, path, Edge(u, v))
                path_cost = cost_so_far + dist
                enqueue!(frontier,
                (path_cost, new_path, v),
                path_cost + heuristic(v))
            end
        end
        colormap[u] = 2
    end
    Vector{Edge}()
end

"""
    a_star(g, s, t[, distmx][, heuristic])

Return a vector of edges comprising the shortest path between vertices `s` and `t`
using the [A* search algorithm](http://en.wikipedia.org/wiki/A%2A_search_algorithm).
An optional heuristic function and edge distance matrix may be supplied. If missing,
the distance matrix is set to [`DefaultDistance`](@ref) and the heuristic is set to
`n -> 0`.
"""
function a_star(
    g::AbstractGraph,  # the g

    s::Integer,                       # the start vertex
    t::Integer,                       # the end vertex
    distmx::AbstractMatrix{T} = weights(g),
    heuristic::Function = n -> 0
    ) where T
    # heuristic (under)estimating distance to target
    U = eltype(g)
    frontier = PriorityQueue(Tuple{T,Vector{Edge},U}, T)
    frontier[(zero(T), Vector{Edge}(), s)] = zero(T)
    colormap = zeros(Int, nv(g))
    colormap[s] = 1
    a_star_impl!(g, t, frontier, colormap, distmx, heuristic)
end

ext: jl


text:
#=

# 6: Surface tension simulation in 3D

WARNING: Takes very long time to compute (hours on cluster for 0.1 second)

=#

module drop

using Printf
using SmoothedParticles
using Parameters

#=
Declare constant parameters
=#

##geometrical
const dr = 3.7e-5          #average particle distance (decrease to make finer simulation)
const h = 3.0*dr           #size of kernel support
const rad = 1e-3
const deskw= 0.9h

##physical
const rho0 = 1000.   	   #fluid density
const m = rho0*dr^3        #particle mass
const mu = 0.1         #dynamic viscosity of water
const beta = 72e-3      #surface tension
const vol = dr^3
const g = -9.8*VECZ
const c = 10.0*max(sqrt(beta/rho0/dr), sqrt(4*norm(g)*rad))            #numerical speed of sound

##temporal
const dt = 0.3*dr/c
const t_end = 2e-5
#const t_end = 0.02
const dt_frame = max(t_end/50,dt)

#artificial
const s0 = dr*dr/100

const FLUID = 0.
const SOLID = 1.

@with_kw mutable struct Particle <: AbstractParticle
	x::RealVector #position
	v::RealVector = VEC0 #velocity
	a::RealVector = VEC0 #acceleration
	P::Float64 = 0. #pressure
	rho::Float64 = 0. #density
    rho0::Float64 = 0.
    n::RealVector = VEC0 #normal
    type::Float64
end

#=
Define geometry and make particles
=#

function make_system()
	grid = Grid(dr, :cubic)
    drop = Ball(0., 0., rad + h, rad)
    desk = Box(-2rad, -2rad, -deskw, 2rad, 2rad, 0.)
    dom = Box(-2rad, -2rad, -2deskw, 2rad, 2rad, 2.2rad)
    sys = ParticleSystem(Particle, dom, h)
	generate_particles!(sys, grid, drop, x -> Particle(x=x, type=FLUID))
    generate_particles!(sys, grid, desk, x -> Particle(x=x, type=SOLID))
	return sys
end

#=
Define particle interactions
=#

@inbounds function find_n!(p::Particle, q::Particle, r::Float64)
    p.n += 2*vol*vol*rDwendland3(h,r)*(p.x - q.x)
end

function reset_n!(p::Particle)
    p.n = VEC0
end

function normalize_n!(p::Particle)
    s = norm(p.n)
    p.n /= (s + s0)
end

@inbounds function find_rho!(p::Particle, q::Particle, r::Float64)
    p.rho += m*wendland3(h,r)
end

@inbounds function find_rho0!(p::Particle, q::Particle, r::Float64)
    p.rho0 += m*wendland3(h,r)
end

function find_pressure!(p::Particle)
	p.P = c^2*(p.rho - p.rho0)
end

@inbounds function internal_force!(p::Particle, q::Particle, r::Float64)
		ker = m*rDwendland3(h,r)
        #pressure
		p.a += -ker*(p.P/rho0^2 + q.P/rho0^2)*(p.x - q.x)
		#viscosity
        p.a += 2*ker*mu/rho0^2*(p.v - q.v)
        #surface tension
        p.a -= 2*beta/rho0^2*(
            (m*DDwendland3(h,r)-ker)*dot(p.x-q.x, p.n-q.n)*(p.x-q.x)/(r^2 + s0)
            +ker*(p.n-q.n)
        )
end

function reset_a!(p::Particle)
    p.a = VEC0
end

function reset_rho!(p::Particle)
    p.rho = 0.0
end

function move!(p::Particle)
	p.x += (p.type==FLUID)*dt*p.v
end

function accelerate!(p::Particle)
	p.v += (p.type==FLUID)*0.5*dt*(p.a+g)
end

function energy(p::Particle)::Float64
	kinetic = 0.5*m*dot(p.v, p.v)
	internal =  0.5*m*c^2*(p.rho - p.rho0)^2/rho0^2
    s = norm(p.n)
    tensile = beta*(s - s0*log(s/s0+1))
    potential = m*dot(p.x, -g)
    return kinetic + internal + tensile + potential
end

#=
Put everything into a time loop
=#

function verlet_step!(sys)
    apply!(sys, accelerate!)
    apply!(sys, move!)
    create_cell_list!(sys)
    apply!(sys, reset_rho!)
    apply!(sys, find_rho!, self = true)
    apply!(sys, reset_n!)
    apply!(sys, find_n!, self = true)
    apply!(sys, normalize_n!)
    apply!(sys, find_pressure!)
    apply!(sys, reset_a!)
    apply!(sys, internal_force!)
    apply!(sys, accelerate!)
end

function save_results!(out, sys, k)
    @printf("t = %.6e\n", k*dt)
    apply!(sys, reset_n!)
    apply!(sys, find_n!, self = true)
    save_frame!(out, sys, :v, :a, :P, :rho, :rho0, :type, :n)
end

function main()
    E0 = 0.0
	sys = make_system()
	out = new_pvd_file("results/drop")
    #initialization
    create_cell_list!(sys)
    apply!(sys, find_rho0!, self = true)
    apply!(sys, find_rho!, self = true)
    apply!(sys, find_pressure!)
    apply!(sys, find_n!)
    apply!(sys, normalize_n!)
    apply!(sys, internal_force!)
	for k in 0 : Int64(round(t_end/dt))
        verlet_step!(sys)
        if (k %  Int64(round(dt_frame/dt)) == 0)
            save_results!(out, sys, k)
            E = sum(energy, sys.particles)
            if k == 0
                E0 = E
            end
            E = E - E0
            @show E
            println("# of part. = ", length(sys.particles))
            println()
        end
	end
	save_pvd_file(out)
end ## function main

end ## module


text:
### VPRK Integrators ###
module VPRKIntegratorsTest

    export test_vprk_integrators

    using GeometricIntegrators.Config
    using GeometricIntegrators.Integrators
    using GeometricIntegrators.Solvers
    using GeometricIntegrators.Tableaus
    using GeometricIntegrators.TestProblems.LotkaVolterra
    using GeometricIntegrators.Utils
    using Test

    set_config(:nls_atol, 8eps())
    set_config(:nls_rtol, 2eps())
    set_config(:nls_solver, QuasiNewtonSolver)
    set_config(:jacobian_autodiff, true)

    using GeometricIntegrators.TestProblems.LotkaVolterra: Δt, nt

    iode = lotka_volterra_2d_iode()

    int  = IntegratorFIRK(lotka_volterra_2d_ode(), getTableauGLRK(8), Δt)
    sol  = integrate(int, nt)
    refx = sol.q[:,end]

    function test_vprk_integrators()
        @testset "$(rpad("VPRK integrators",80))" begin

            iint = Integrator(iode, TableauVPRK(:pglrk, 2, getCoefficientsGLRK(1), -1), Δt)
            isol = integrate(iint, nt)

            # println(rel_err(isol.q, refx))
            @test rel_err(isol.q, refx) < 2E-6

            iint = Integrator(iode, TableauVPRK(:pglrk, 4, getCoefficientsGLRK(2), +1), Δt)
            isol = integrate(iint, nt)

            # println(rel_err(isol.q, refx))
            @test rel_err(isol.q, refx) < 8E-7

            iint = Integrator(iode, TableauVPRK(:pglrk, 6, getCoefficientsGLRK(3), -1), Δt)
            isol = integrate(iint, nt)

            # println(rel_err(isol.q, refx))
            @test rel_err(isol.q, refx) < 4E-12

            vint = Integrator(iode, getTableauVPLobIIIA2(), Δt)
            vsol = integrate(vint, nt)

            # println(rel_err(isol.q, refx))
            @test rel_err(vsol.q, refx) < 4E-6

            vint = Integrator(iode, getTableauVPLobIIIA3(), Δt)
            vsol = integrate(vint, nt)

            # println(rel_err(vsol.q, refx))
            @test rel_err(vsol.q, refx) < 8E-7

            vint = Integrator(iode, getTableauVPLobIIIA4(), Δt)
            vsol = integrate(vint, nt)

            # println(rel_err(vsol.q, refx))
            @test rel_err(vsol.q, refx) < 3E-11

            vint = Integrator(iode, getTableauVPLobIIIB2(), Δt)
            vsol = integrate(vint, nt)

            # println(rel_err(vsol.q, refx))
            @test rel_err(vsol.q, refx) < 2E-6

            vint = Integrator(iode, getTableauVPLobIIIB3(), Δt)
            vsol = integrate(vint, nt)

            # println(rel_err(vsol.q, refx))
            @test rel_err(vsol.q, refx) < 8E-7

            vint = Integrator(iode, getTableauVPLobIIIB4(), Δt)
            vsol = integrate(vint, nt)

            # println(rel_err(vsol.q, refx))
            @test rel_err(vsol.q, refx) < 2E-11

            vint = IntegratorVPRKpStandard(iode, getTableauVPGLRK(1), Δt)
            isol = integrate(vint, nt)

            # println(rel_err(isol.q, refx))
            @test rel_err(isol.q, refx) < 1E-6

            vint = IntegratorVPRKpStandard(iode, getTableauVPGLRK(2), Δt)
            isol = integrate(vint, nt)

            # println(rel_err(isol.q, refx))
            @test rel_err(isol.q, refx) < 1E-11

            vint = IntegratorVPRKpStandard(iode, getTableauVPGLRK(3), Δt)
            isol = integrate(vint, nt)

            # println(rel_err(isol.q, refx))
            @test rel_err(isol.q, refx) < 5E-16

            vint = IntegratorVPRKpSymplectic(iode, getTableauVPGLRK(1), Δt)
            isol = integrate(vint, nt)

            # println(rel_err(isol.q, refx))
            @test rel_err(isol.q, refx) < 4E-6

            vint = IntegratorVPRKpSymplectic(iode, getTableauVPGLRK(2), Δt)
            isol = integrate(vint, nt)

            # println(rel_err(isol.q, refx))
            @test rel_err(isol.q, refx) < 1E-11

            vint = IntegratorVPRKpSymplectic(iode, getTableauVPGLRK(3), Δt)
            isol = integrate(vint, nt)

            # println(rel_err(isol.q, refx))
            @test rel_err(isol.q, refx) < 8E-12

            vint = IntegratorVPRKpSymmetric(iode, getTableauVPGLRK(1), Δt)
            isol = integrate(vint, nt)

            # println(rel_err(isol.q, refx))
            @test rel_err(isol.q, refx) < 1E-6

            vint = IntegratorVPRKpSymmetric(iode, getTableauVPGLRK(2), Δt)
            isol = integrate(vint, nt)

            # println(rel_err(isol.q, refx))
            @test rel_err(isol.q, refx) < 1E-11

            vint = IntegratorVPRKpSymmetric(iode, getTableauVPGLRK(3), Δt)
            isol = integrate(vint, nt)

            # println(rel_err(isol.q, refx))
            @test rel_err(isol.q, refx) < 4E-16

            vint = IntegratorVPRKpMidpoint(iode, getTableauVPGLRK(1), Δt)
            isol = integrate(vint, nt)

            # println(rel_err(isol.q, refx))
            @test rel_err(isol.q, refx) < 1E-6

            vint = IntegratorVPRKpMidpoint(iode, getTableauVPGLRK(2), Δt)
            isol = integrate(vint, nt)

            # println(rel_err(isol.q, refx))
            @test rel_err(isol.q, refx) < 1E-11

            vint = IntegratorVPRKpMidpoint(iode, getTableauVPGLRK(3), Δt)
            isol = integrate(vint, nt)

            # println(rel_err(isol.q, refx))
            @test rel_err(isol.q, refx) < 4E-16

            vintV1 = IntegratorVPRKpVariational(iode, getTableauVPGLRK(1), Δt)
            isolV1 = integrate(vintV1, nt)

            # println(rel_err(isolV1.q, refx))
            @test rel_err(isolV1.q, refx) < 8E-7

            vintV2 = IntegratorVPRKpVariational(iode, getTableauVPGLRK(2), Δt)
            isolV2 = integrate(vintV2, nt)

            # println(rel_err(isolV2.q, refx))
            @test rel_err(isolV2.q, refx) < 8E-8

            vintV3 = IntegratorVPRKpVariational(iode, getTableauVPGLRK(3), Δt)
            isolV3 = integrate(vintV3, nt)

            # println(rel_err(isolV3.q, refx))
            @test rel_err(isolV3.q, refx) < 1E-11

            vintQ1 = IntegratorVPRKpVariationalQ(iode, getTableauVPGLRK(1), Δt)
            isolQ1 = integrate(vintQ1, nt)

            # println(rel_err(isolQ1.q, refx))
            @test rel_err(isolQ1.q, refx) < 4E-5

            vintQ2 = IntegratorVPRKpVariationalQ(iode, getTableauVPGLRK(2), Δt)
            isolQ2 = integrate(vintQ2, nt)

            # println(rel_err(isolQ2.q, refx))
            @test rel_err(isolQ2.q, refx) < 2E-4

            vintQ3 = IntegratorVPRKpVariationalQ(iode, getTableauVPGLRK(3), Δt)
            isolQ3 = integrate(vintQ3, nt)

            # println(rel_err(isolQ3.q, refx))
            @test rel_err(isolQ3.q, refx) < 1E-8

            vintP1 = IntegratorVPRKpVariationalP(iode, getTableauVPGLRK(1), Δt)
            isolP1 = integrate(vintP1, nt)

            # println(rel_err(isolP1.q, refx))
            @test rel_err(isolP1.q, refx) < 8E-7

            vintP2 = IntegratorVPRKpVariationalP(iode, getTableauVPGLRK(2), Δt)
            isolP2 = integrate(vintP2, nt)

            # println(rel_err(isolP2.q, refx))
            @test rel_err(isolP2.q, refx) < 8E-8

            vintP3 = IntegratorVPRKpVariationalP(iode, getTableauVPGLRK(3), Δt)
            isolP3 = integrate(vintP3, nt)

            # println(rel_err(isolP3.q, refx))
            @test rel_err(isolP3.q, refx) < 1E-11

            @test rel_err(isolV1.q, isolP1.q[:,end]) == 0
            @test rel_err(isolV2.q, isolP2.q[:,end]) == 0
            @test rel_err(isolV3.q, isolP3.q[:,end]) == 0
        end
    end
end

using .VPRKIntegratorsTest
test_vprk_integrators()


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

# ╔═╡ 14b6c3e4-48dd-11ec-336c-6918bf024f43
begin
	let
		using Pkg
		Pkg.activate(mktempdir())
		Pkg.Registry.update()
		Pkg.add("Revise")
		Pkg.add("PlutoUI")
		Pkg.add("ImageFiltering")
		Pkg.add("Images")
		Pkg.add("ImageMorphology")
		Pkg.add("ImageSegmentation")
		Pkg.add("ImageComponentAnalysis")
		Pkg.add("DataFrames")
		Pkg.add("Statistics")
		Pkg.add("CairoMakie")
		Pkg.add("DataStructures")
		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/Phantoms.jl")
	end

	using Revise
	using PlutoUI
	using ImageFiltering
	using Images
	using ImageMorphology
	using ImageSegmentation
	using ImageComponentAnalysis
	using DataFrames
	using Statistics
	using CairoMakie
	using DataStructures
	using DICOM
	using DICOMUtils
	using Phantoms
end

# ╔═╡ 9506ae45-f63f-4fd6-8467-0597cbde7006
TableOfContents()

# ╔═╡ 482d8669-c750-450b-aae8-03afc2a90007
md"""
## `dcm_list_builder`
"""

# ╔═╡ f53772f8-d9f0-4562-896d-6a6112cd8735
function dcm_list_builder(path)
    dcm_path_list = []
    for (dirpath, dirnames, filenames) in walkdir(path, topdown=true)
        if (dirpath in dcm_path_list) == false
            for filename in filenames
                try
                    tmp_str = string(dirpath, "/", filename)
                    ds = dcm_parse(tmp_str)
                    if (dirpath in dcm_path_list) == false
                        push!(dcm_path_list, dirpath)
					end
				catch
                    nothing
				end
			end
        else
                nothing
		end
	end
    return dcm_path_list
end

# ╔═╡ 3f313b73-2785-4299-a1e9-bc5d4e1ad422
md"""
## `dcm_reader`
"""

# ╔═╡ 202cef5d-85f1-4673-9e6c-52ba80e44de0
function dcm_reader(dcm_path)
    dcm_files = []
    for (dirpath, dirnames, filenames) in walkdir(dcm_path, topdown=false)
        for filename in filenames
            try
                if (filename == "DIRFILE") == false
                    dcm_file = string(dirpath, "/", filename)
                    dcm_parse(dcm_file)
                    push!(dcm_files, dcm_file)
				end
			catch
				nothing
			end
		end
	end

    read_RefDs = true
	local RefDs
    while read_RefDs
        for index in range(1, length(dcm_files))
            try
                RefDs = dcm_parse(dcm_files[index])
                read_RefDs = false
                break
			catch
                nothing
			end
		end
	end

	header = RefDs.meta
	slice_thick_ori = header[(0x0018, 0x0050)]
	rows, cols = Int(header[(0x0028, 0x0010)]), Int(header[(0x0028, 0x0011)])

    ConstPixelDims = (rows, cols, length(dcm_files))
    dcm_array = zeros(ConstPixelDims...)

    instances = []
    for filenameDCM in dcm_files
        try
            ds = dcm_parse(filenameDCM)
			head = ds.meta
			InstanceNumber = head[(0x0020, 0x0013)]
            push!(instances, InstanceNumber)
		catch
            nothing
		end
	end

    sort!(instances)

    index = 0
    for filenameDCM in dcm_files
        try
            ds = dcm_parse(filenameDCM)
			head = ds.meta
			InstanceNumber = head[(0x0020, 0x0013)]
			index = findall(x -> x==InstanceNumber, instances)
			pixel_array = head[(0x7fe0, 0x0010)]
            dcm_array[:, :, index] = pixel_array
            index += 1
		catch
            nothing
		end
	end

    RescaleSlope = header[(0x0028, 0x1053)]
	RescaleIntercept = header[(0x0028, 0x1052)]
    dcm_array = dcm_array .* RescaleSlope .+ RescaleIntercept
    return RefDs.meta, dcm_array, slice_thick_ori
end

# ╔═╡ 39688933-f68f-47b4-a9dc-aff25bdc6123
md"""
## Load DICOMs
"""

# ╔═╡ 0883a26a-54ba-4670-a19b-a6a0a23f5aee
root_path = "/Users/daleblack/Google Drive/Datasets/Canon_Aquilion_One_Vision"

# ╔═╡ 59976bc3-271e-4418-874c-b553b59e730a
dcm_path_list = dcm_list_builder(root_path)

# ╔═╡ f6a29f28-b787-43d6-8b6b-5e163cbbbcec
dcm_reader(dcm_path_list[1]);

# ╔═╡ 8652d37a-7dc2-4123-b20f-87a8981baa26
dcm_path_list[2]

# ╔═╡ 9b02d4cc-e9b3-47aa-8d57-e0befd64a5fc
header, dcm_array, slice_thick_ori = dcm_reader(dcm_path_list[2]);

# ╔═╡ 2185c4e4-a518-4133-8174-662e0dd744eb
heatmap(transpose(dcm_array[:, :, 8]), colormap=:grays)

# ╔═╡ d9bbaeb8-b027-4a58-9ddb-0b015aa867de
maximum(dcm_array)

# ╔═╡ 7d7b1fd5-e244-4010-ac61-6f4e1925f216
minimum(dcm_array)

# ╔═╡ 237cfd63-3479-468d-954f-2ebd05376fca
md"""
# Whole heart mask
"""

# ╔═╡ 1b755e07-926a-42bd-871e-a693bb037c81
md"""
## `find_circle`
"""

# ╔═╡ d2d1d02c-e592-4c73-afa2-80f89dd5534e
function find_circle(point_1, point_2, point_3)
    x1, y1 = point_1
    x2, y2 = point_2
    x3, y3 = point_3

    x12 = x1 - x2
    x13 = x1 - x3
    y12 = y1 - y2
    y13 = y1 - y3
    y31 = y3 - y1
    y21 = y2 - y1
    x31 = x3 - x1
    x21 = x2 - x1

    sx13 = x1^2 - x3^2
    sy13 = y1^2 - y3^2
    sx21 = x2^2 - x1^2
    sy21 = y2^2 - y1^2

    f = (((sx13) * (x12) + (sy13) * (x12) + (sx21) * (x13) + (sy21) * (x13)) ÷ (2 * ((y31) * (x12) - (y21) * (x13))))

    g = (((sx13) * (y12) + (sy13) * (y12) + (sx21) * (y13) + (sy21) * (y13)) ÷ (2 * ((x31) * (y12) - (x21) * (y13))))

    # eqn of circle be x^2 + y^2 + 2*g*x + 2*f*y + c = 0 where center is (h = -g, k = -f)
    center_insert = [-g, -f]

    return center_insert
end

# ╔═╡ 2d16b7cb-584e-4ff8-84a1-4e5670bd49d3
find_circle([309, 309], [312, 200], [155, 155])

# ╔═╡ a69293d3-6093-4b87-abc6-11f1ec9398c3
md"""
## `mask_heart`
"""

# ╔═╡ dfcda44d-5b2c-4a50-9fb2-44cb70b5488e
"""
    mask_heart(
        header;
        array_used=nothing,
        radius_val=95,
        slice_used_center=nothing
        )
Given a QRM Phantom with a heart insert, this function will create a mask of the whole heart for image-processing purposes.
"""
function mask_heart(header, array_used, slice_used_center; radius_val=95)
    pixel_size = Phantoms.get_pixel_size(header)

    radius = (radius_val / 2) / pixel_size[1]
    central_image = copy(array_used[:, :, slice_used_center])
    central_image = Int.(central_image .< -200)
    kern = Int.(round(5 / pixel_size[1]))
    if kern % 2 == 0
        kern += 1
    end
    central_image = mapwindow(median, central_image, (kern, kern))
    center = [size(central_image, 1) ÷ 2, size(central_image, 2) ÷ 2]
    a = copy(central_image)
    local point_1
    for index in 1:(size(central_image, 2) ÷ 2)
        if (
            central_image[center[1] + index, center[2] + index] == 1 &&
            central_image[center[1] + index, center[2] + index + 5] == 1
        )
            point_1 = [center[1] + index, center[2] + index]
            break
        else
            a[center[1] + index, center[2] + index] = 2
        end
    end

    local point_2
    for index in 1:(size(central_image, 2) ÷ 2)
        if (
            central_image[center[1] + index, center[2] - index] == 1 &&
            central_image[center[1] + index, center[2] - index - 5] == 1
        )
            point_2 = [center[1] + index, center[2] - index]
            break
        else
            a[center[1] + index, center[2] - index] = 2
        end
    end

    local point_3
    for index in 1:(size(central_image, 2) ÷ 2)
        if (
            central_image[center[1] - index, center[2] - index] == 1 &&
            central_image[center[1] - index, center[2] - index - 5] == 1
        )
            point_3 = [center[1] - index, center[2] - index]
            break
        else
            a[center[1] - index, center[2] - index] = 2
        end
    end

    center_insert = find_circle(point_1, point_2, point_3)
    rows, cols = Int(header[(0x0028, 0x0010)]), Int(header[(0x0028, 0x0011)])
    Y, X = collect(1:rows), collect(1:cols)'
    dist_from_center = @. sqrt((X - center_insert[2])^2 + (Y - center_insert[1])^2)

    mask = dist_from_center .<= radius[1]
    masked_array = zeros(size(array_used))
    for index in 1:size(array_used, 3)
        masked_array[:, :, index] = array_used[:, :, index] .* mask
    end

    return masked_array, center_insert, mask
end

# ╔═╡ 57730e9f-05b2-4e3a-bff3-a3a6403a66e0
masked_array, center_insert, mask = mask_heart(header, dcm_array, size(dcm_array, 3)÷2);

# ╔═╡ 19e0f1b3-86c3-4ffd-aa39-dd7e1e509312
center_insert

# ╔═╡ 55665585-8076-4e23-a989-71fc9909a57b
begin
	fig = Figure()

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

# ╔═╡ 78ffb833-27c5-4df5-9989-47d57aca47fe
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

# ╔═╡ 1ac4e0fd-79fb-41d1-b2c7-8b5e011c80e2
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

# ╔═╡ 5147ce8d-2d44-443b-a6a2-edba166356a8
@bind a2 PlutoUI.Slider(1:size(masked_array, 3), default=10, show_value=true)

# ╔═╡ 004f4a1e-9d11-44ef-8622-69c311edb2cb
heatmap(transpose(masked_array[:, :, a2]), colormap=:grays)

# ╔═╡ da57bd3c-0d8b-453c-9c4f-ea846850cbec
md"""
# Calcium rod mask
"""

# ╔═╡ 11e0c502-5b1f-4e60-bb43-a8c04a4e335b
md"""
## `get_calcium_slices`
"""

# ╔═╡ c4c6841a-88f4-4ca9-9353-0e9daa6d4dff
"""
	get_calcium_slices(dcm_array, header; calcium_threshold=130)

Returns the slices that contain the calcium vessel inserts `slice_dict`
and the slices that contain the calcium calibration rod insert `large_index`.

The calcium rod slice `large_index` usually omits the first and last slice
of the rod. So to account for all of the slices containing the calcium rod,
one would want a range like so: `(large_index - 1):(large_index + 1)`
"""
function get_calcium_slices(dcm_array, header; calcium_threshold=130)
    array = copy(dcm_array)
    array = Int.(array .> (1.1 * calcium_threshold))

    pixel_size = Phantoms.get_pixel_size(header)
    CCI_5mm_num_pixels = Int(round(π * (5 / 2)^2 / pixel_size[1]^2))
    cal_rod_num_pixels = Int(round(π * (20 / 2)^2 / pixel_size[1]^2))

    kern = Int.(round(5 / pixel_size[1]))
    if kern % 2 == 0
        kern += 1
    end

    slice_dict = Dict()
    large_index = []
    cal_rod_dict = Dict()
    for idx in 1:size(array, 3)
        array_filtered = mapwindow(median, array[:, :, idx], (kern, kern))
        components = ImageComponentAnalysis.label_components(array_filtered)
        a1 = analyze_components(components, BasicMeasurement(; area=true, perimeter=true))
        a2 = analyze_components(components, BoundingBox(; box_area=true))
        df = leftjoin(a1, a2; on=:l)

        count_5mm = 0
        count = 0
        for row in eachrow(df)
            count += 1
            df_area = Int(round(row[:area]))

            r1_1 = Int(round(CCI_5mm_num_pixels * 0.6))
            r1_2 = Int(round(CCI_5mm_num_pixels * 1.5))
            r2_1 = Int(round(cal_rod_num_pixels * 0.7))
            r2_2 = Int(round(cal_rod_num_pixels * 1.3))

            if df_area in r1_1:r1_2
                count_5mm += 1
            elseif df_area in r2_1:r2_2
                indices = row[:box_indices]
                x_point = ((indices[1][end] - indices[1][1]) ÷ 2) + indices[1][1]
                y_point = ((indices[2][end] - indices[2][1]) ÷ 2) + indices[2][1]
                cal_rod_dict[count] = [x_point, y_point]
            end
        end

        if count_5mm > 0 && count_5mm < 4
            slice_dict[idx] = count_5mm
        end

        poppable_keys = []
        for key in cal_rod_dict
            start_coordinate = [key[2][1], key[2][2]]

            x_right = 0
            while array_filtered[start_coordinate[1], start_coordinate[2] + x_right] == 1
                x_right += 1
            end

            x_left = 0
            while array_filtered[start_coordinate[1], start_coordinate[2] - x_left] == 1
                x_left += 1
            end

            y_top = 0
            while array_filtered[start_coordinate[1] + y_top, start_coordinate[2]] == 1
                y_top += 1
            end

            y_bottom = 0
            while array_filtered[start_coordinate[1] - y_bottom, start_coordinate[2]] == 1
                y_bottom += 1
            end

            x_dist = x_right + x_left
            y_dist = y_top + y_bottom

            range1 = round(0.7 * y_dist):round(1.2 * y_dist)
            if ((x_dist in range1) == false) ||
                ((round(0.7 * y_dist) == 0) && (round(1.2 * y_dist) == 0))
                push!(poppable_keys, key)
            else
                nothing
            end
        end

        for key in poppable_keys
            pop!(cal_rod_dict)
        end

        if length(cal_rod_dict) == 0
            nothing
        else
            append!(large_index, idx)
        end
    end
    return slice_dict, large_index
end

# ╔═╡ 462fd944-d4ce-4d7b-82ad-40342daff354
slice_dict, large_index = get_calcium_slices(masked_array, header);

# ╔═╡ 8959750c-291c-4c2a-a6a1-218f8ed82ab1
slice_dict, large_index

# ╔═╡ 68c16f9e-8520-41ef-912a-9e61acf2ad1d
md"""
## `get_calcium_center_slices`
"""

# ╔═╡ f49fe2a8-8429-4fd1-98e3-6e75189a3ac2
"""
	get_calcium_center_slices(dcm_array, slice_dict, large_index)

Returns the slices that contain the calcium vessel inserts `slice_dict`
and the center of the calcium calibration rod slice `flipped_index`
"""
function get_calcium_center_slices(dcm_array, slice_dict, large_index)
    flipped_index = Int(round(median(large_index)))
    edge_index = []
    if flipped_index < (size(dcm_array, 3) / 2)
        flipped = -1
        for element in large_index
            if element > (size(dcm_array, 3) / 2)
                append!(edge_index, element)
            end
        end
        if length(edge_index) == 0
            nothing
        else
            for index_edge in minimum(edge_index):size(dcm_array, 3)
                try
                    delete!(slice_dict, index_edge)
                catch
                    nothing
                end
            end
            for element2 in edge_index
                deleteat!(large_index, findall(x -> x == element2, large_index))
            end
        end

        for element in 1:maximum(large_index)
            try
                delete!(slice_dict, element)
            catch
                nothing
            end
        end
    else
        flipped = 1
        for element in large_index
            if element < (size(dcm_array, 3) / 2)
                append!(edge_index, element)
            end
        end
        if length(edge_index) == 0
            nothing
        else
            for index_edge in 1:maximum(edge_index)
                try
                    delete!(slice_dict, index_edge)
                catch
                    nothing
                end
            end
            for element2 in edge_index
                deleteat!(large_index, findall(x -> x == element2, large_index))
            end
        end
        for element in minimum(large_index):size(dcm_array, 3)
            try
                delete!(slice_dict, element)
            catch
                nothing
            end
        end
    end
    return slice_dict, flipped, flipped_index
end

# ╔═╡ 494163b1-cb52-45a9-815e-991633018cd5
slice_dict2, flipped, flipped_index = get_calcium_center_slices(masked_array, slice_dict, large_index)

# ╔═╡ 880db52b-e39c-4578-a6de-ce4a6afa038e
md"""
## `poppable_keys`
"""

# ╔═╡ 6613bd97-7d53-4b8b-a69c-12e040f4aa36
"""
    poppable_keys(flipped, flipped_index, header, slice_dict)

No idea what the function is
"""
function poppable_keys(flipped, flipped_index, header, slice_dict)
    SliceThickness = header[(0x0018, 0x0050)]
    poppable_keys = []
    if flipped == -1
        for key in slice_dict
            if key[1] > (flipped_index + (55 / SliceThickness))
                append!(poppable_keys, key)
            elseif flipped == 1
                for key in slice_dict
                    if key[1] < (flipped_index - (55 / SliceThickness))
                        append!(poppable_keys, key)
                    end
                end
            end
        end
    end
    for key in poppable_keys
        pop!(slice_dict)
    end
    return slice_dict
end

# ╔═╡ bff84e32-1013-44f2-9eda-eb5f6bac4c5a
pop_keys = poppable_keys(flipped, flipped_index, header, slice_dict)

# ╔═╡ 451dd3c6-f740-47df-a5b9-85749caf04fb
md"""
## `compute_CCI`
"""

# ╔═╡ 1ac3d601-0057-475e-a6ac-9522618f7bb9
function compute_CCI(dcm_array, header, slice_dict, flipped; calcium_threshold=130)
	SliceThickness = header[(0x0018,0x0050)]
	max_key, _ = maximum(zip(values(slice_dict), keys(slice_dict)))
    max_keys = []
    for key in slice_dict
        if key[2] == max_key
            append!(max_keys, key[1])
		end
	end
    slice_CCI = Int(floor(median(max_keys)))

    array = copy(dcm_array)
    array = Int.(array .> calcium_threshold)

    calcium_image = array .* dcm_array
    quality_slice = Int.(round(slice_CCI - flipped * (20 / SliceThickness)))

    cal_rod_slice = slice_CCI + (flipped * Int(30 / SliceThickness))

    return calcium_image, slice_CCI, quality_slice, cal_rod_slice
end

# ╔═╡ 5d407e49-fe05-4748-be36-b232b5ce3ccb
calcium_image1, slice_CCI1, quality_slice1, cal_rod_slice1 = compute_CCI(masked_array, header, slice_dict, flipped);

# ╔═╡ de1ac95e-2399-45c1-b337-dcfe800394a2
slice_CCI1, quality_slice1, cal_rod_slice1

# ╔═╡ a8f6ece5-9f59-4af1-bcf7-7dafa839bae2
@bind a PlutoUI.Slider(1:size(calcium_image1, 3), default=10, show_value=true)

# ╔═╡ 267d1d6a-5c35-480e-8425-f000f7f31c7f
begin
	fig4 = Figure()

	ax4 = Makie.Axis(fig4[1, 1])
	heatmap!(transpose(calcium_image1[:, :, a]), colormap=:grays)
	ax5 = Makie.Axis(fig4[1,2])
	heatmap!(transpose(dcm_array[:, :, a]), colormap=:grays)
	fig4
end

# ╔═╡ 9a77474e-a853-4369-a6f3-a9e76a683a27
md"""
## `mask_rod`
"""

# ╔═╡ a9f499f8-ca32-4310-82d7-f20b790dc04c
"""
	mask_rod(dcm_array, header; calcium_threshold=130)
"""
function mask_rod(dcm_array, header; calcium_threshold=130)
    slice_dict, large_index = get_calcium_slices(
        dcm_array, header; calcium_threshold=calcium_threshold
    )
    slice_dict, flipped, flipped_index = get_calcium_center_slices(
        dcm_array, slice_dict, large_index
    )
    slice_dict = poppable_keys(flipped, flipped_index, header, slice_dict)
    calcium_image, slice_CCI, quality_slice, cal_rod_slice = compute_CCI(
        dcm_array, header, slice_dict, flipped; calcium_threshold=calcium_threshold
    )
    return calcium_image, slice_CCI, quality_slice, cal_rod_slice
end

# ╔═╡ 4eece882-8416-4330-b95a-1cfaeac89aac
calcium_image, slice_CCI, quality_slice, cal_rod_slice = mask_rod(masked_array, header);

# ╔═╡ d6cb5578-e924-4e78-be0b-82d8866fbee8
slice_CCI

# ╔═╡ 214d485e-df7e-4a45-b52f-10f854d839a1
@bind b PlutoUI.Slider(1:size(calcium_image, 3), default=10, show_value=true)

# ╔═╡ 03622243-d2af-4d1b-a6e1-fd50670e48a8
heatmap(transpose(calcium_image[:, :, b]), colormap=:grays)

# ╔═╡ 93861c7a-4086-4404-842e-ca91c2ad6094
heatmap(transpose(calcium_image[:, :, slice_CCI]), colormap=:grays)

# ╔═╡ ee208b64-3cf3-428b-a5f6-74793a82fd93
md"""
# Calcium inserts mask
"""

# ╔═╡ 5930b9c5-2f85-465b-abe0-223580142d49
md"""
## `angle_calc`
"""

# ╔═╡ 480bd663-d822-4417-b7f0-6679d26a58b4
function angle_calc(side1, side2)
    #Calculate angle between two sides of rectangular triangle
    if side1 == 0
        angle = 0
	elseif side2 == 0
        angle = π / 2
    else
        angle = atan(side1 / side2)
	end

    return angle
end

# ╔═╡ f498c292-367d-4441-873e-29bec99af2a6
angle_calc(4, 3)

# ╔═╡ 56ce9966-939e-48cd-9c70-2e441bfaea58
md"""
## `create_circular_mask`
"""

# ╔═╡ 228dd95c-fdab-4a5a-a079-e14238ff5e83
function create_circular_mask(h, w, center_circle, radius_circle)
	Y, X = collect(1:h), collect(1:w)'
    dist_from_center = sqrt.((X .- center_circle[1]).^2 .+ (Y .- center_circle[2]).^2)

    mask = dist_from_center .<= radius_circle

    return mask
end

# ╔═╡ 0f572e01-a690-4699-8a8c-6684c27085bb
mask1 = create_circular_mask(40, 40, [20, 20], 1);

# ╔═╡ ce5b9c0a-4f25-4b8f-bf16-17d05b7bfff4
heatmap(mask1, colormap=:grays)

# ╔═╡ 5a8963f1-3ec8-4dbe-b567-df2648219b7f
md"""
## `calc_output`
"""

# ╔═╡ 92fca992-b36d-48f3-803e-f92e6463c291
"""
	calc_output(dcm_array, CCI_slice, calcium_threshold=130, comp_connect=trues(3, 3))

Calculate the output of a dcm_array
"""
function calc_output(
    dcm_array, header, CCI_slice, calcium_threshold=130, comp_connect=trues(3, 3)
)
    # Actual scoring for CCI insert
    # First step is to remove slices without calcium from arrays
    PixelSpacing = Phantoms.get_pixel_size(header)
    SliceThickness = header[(0x0018, 0x0050)]
    CCI_min = Int((CCI_slice - round(5 / SliceThickness, RoundUp)) - 1)
    CCI_max = Int((CCI_slice + round(5 / SliceThickness, RoundUp)) + 1)
    central_CCI = Int(round(CCI_max - CCI_min) / 2)

    if CCI_min < 0
        CCI_min = 0
    end
    if CCI_max > size(dcm_array, 3)
        CCI_max = size(dcm_array, 3)
    end

    CCI_array = copy(dcm_array[:, :, CCI_min:CCI_max])

    image_kernel = Int(round(3 / PixelSpacing[1]))
    if image_kernel % 2 == 0
        image_kernel += 1
    end

    CCI_array_binary = copy(CCI_array)
    CCI_array_binary = Int.(CCI_array_binary .> 1.0 * calcium_threshold)
    inp =
        CCI_array_binary[:, :, central_CCI - 1] +
        CCI_array_binary[:, :, central_CCI] +
        CCI_array_binary[:, :, central_CCI + 1]
    components = ImageComponentAnalysis.label_components(inp, comp_connect)
    a1 = analyze_components(components, BasicMeasurement(; area=true, perimeter=true))
    a2 = analyze_components(components, BoundingBox(; box_area=true))
    df = leftjoin(a1, a2; on=:l)
    centroids = []
    for row in eachrow(df)
        indices = row[:box_indices]
        x_point = ((indices[1][end] - indices[1][1]) ÷ 2) + indices[1][1]
        y_point = ((indices[2][end] - indices[2][1]) ÷ 2) + indices[2][1]
        push!(centroids, (x_point, y_point))
    end

    centroids = deleteat!(centroids, 1)

    i1 = mapwindow(
        median, CCI_array_binary[:, :, central_CCI - 1], (image_kernel, image_kernel)
    )
    i2 = mapwindow(
        median, CCI_array_binary[:, :, central_CCI], (image_kernel, image_kernel)
    )
    i3 = mapwindow(
        median, CCI_array_binary[:, :, central_CCI + 1], (image_kernel, image_kernel)
    )

    image_for_center = i1 + i2 + i3

    components2 = ImageComponentAnalysis.label_components(image_for_center, comp_connect)
    components2 = Int.(components2 .> 0)
    components2 = ImageComponentAnalysis.label_components(components2, comp_connect)

    b1 = analyze_components(components2, BasicMeasurement(; area=true, perimeter=true))
    b2 = analyze_components(components2, BoundingBox(; box_area=true))
    df2 = leftjoin(b1, b2; on=:l)
    centroids2 = []
    for row in eachrow(df2)
        indices = row[:box_indices]
        x_point = ((indices[1][end] - indices[1][1]) ÷ 2) + indices[1][1]
        y_point = ((indices[2][end] - indices[2][1]) ÷ 2) + indices[2][1]
        push!(centroids2, (y_point, x_point))
    end

    output = length(unique(components2)), components2, df2, centroids2
    return output
end

# ╔═╡ fb320178-0d33-4402-b075-7bf7541f94e4
output = calc_output(masked_array, header, slice_CCI1);

# ╔═╡ 68dae231-9907-44b2-afca-a5dac57998bb
output

# ╔═╡ 14f46c1f-62fa-45ff-9c78-b01404c466da
heatmap(transpose(output[2]))

# ╔═╡ feca7b76-fc10-4aca-adff-4cd323dabc74
md"""
## `center_points`
"""

# ╔═╡ ae9ab557-553c-474a-8865-37807b9270fb
"""
	center_points(output, header, tmp_center, CCI_slice)

Function ...
"""
function center_points(dcm_array, output, header, tmp_center, CCI_slice)
    PixelSpacing = Phantoms.get_pixel_size(header)
    rows, cols = Int(header[(0x0028, 0x0010)]), Int(header[(0x0028, 0x0011)])
    sizes = []
    for row in eachrow(output[3])
        area = row[:area]
        append!(sizes, area)
    end

    centroids = output[4]
    largest = Dict()
    for index in 1:length(centroids)
        x = centroids[index][1]
        y = centroids[index][2]
        dist_loc = sqrt((tmp_center[2] - x)^2 + (tmp_center[1] - y)^2)
        dist_loc *= PixelSpacing[1]
        if dist_loc > 31
            largest[index] = [round(y), round(x)]
        else
            nothing
        end
    end

    max_dict = Dict()
    radius = round(2.5 / PixelSpacing[1], RoundUp)
    for key in largest
        tmp_arr = create_circular_mask(rows, cols, (key[2][2], key[2][1]), radius)
        tmp_arr = @. abs(tmp_arr * dcm_array[:, :, CCI_slice]) +
            abs(tmp_arr * dcm_array[:, :, CCI_slice - 1]) +
            abs(tmp_arr * dcm_array[:, :, CCI_slice + 1])
        tmp_arr = @. ifelse(tmp_arr == 0, missing, tmp_arr)
        max_dict[key[1]] = median(skipmissing(tmp_arr))
    end
    large1_index, large1_key = maximum(zip(values(max_dict), keys(max_dict)))
    pop!(max_dict, large1_key)
    large2_index, large2_key = maximum(zip(values(max_dict), keys(max_dict)))
    pop!(max_dict, large2_key)
    large3_index, large3_key = maximum(zip(values(max_dict), keys(max_dict)))

    center1 = largest[large1_key]
    center2 = largest[large2_key]
    center3 = largest[large3_key]

    center = find_circle(center1, center2, center3)
    return center, center1, center2, center3
end

# ╔═╡ 39ce8635-28a7-4eee-875e-325e85d142d5
center_insert

# ╔═╡ 42e86a48-1438-4edd-ac2d-f307d534e165
center, center1, center2, center3 = center_points(dcm_array, output, header, center_insert, slice_CCI1)

# ╔═╡ 0f73d3a7-c6fa-49ee-8fdf-357b0fe94ec8
md"""
## `calc_centers`
"""

# ╔═╡ 9ac2f9ed-92bb-4540-add5-422df2a3f6da
"""
	calc_centers(dcm_array, output, header, tmp_center, CCI_slice)

Function ...
"""
function calc_centers(dcm_array, output, header, tmp_center, CCI_slice)
    PixelSpacing = Phantoms.get_pixel_size(header)
    center, center1, center2, center3 = center_points(
        dcm_array, output, header, tmp_center, CCI_slice
    )
    centers = Dict()
    for size_index4 in (center1, center2, center3)
        center_index = size_index4
        side_x = abs(center[1] - center_index[1])
        side_y = abs(center[2] - center_index[2])

        angle = angle_calc(side_x, side_y)
        if (center_index[1] < center[1] && center_index[2] < center[2])
            medium_calc = [
                center_index[1] + (12.5 / PixelSpacing[1]) * sin(angle),
                (center_index[2] + (12.5 / PixelSpacing[2]) * cos(angle)),
            ]
            low_calc = [
                center_index[1] + (25 / PixelSpacing[1]) * sin(angle),
                (center_index[2] + (25 / PixelSpacing[2]) * cos(angle)),
            ]
        elseif (center_index[1] < center[1] && center_index[2] > center[2])
            medium_calc = [
                center_index[1] + (12.5 / PixelSpacing[1]) * sin(angle),
                (center_index[2] - (12.5 / PixelSpacing[2]) * cos(angle)),
            ]
            low_calc = [
                center_index[1] + (25 / PixelSpacing[1]) * sin(angle),
                (center_index[2] - (25 / PixelSpacing[2]) * cos(angle)),
            ]
        elseif (center_index[1] > center[1] && center_index[2] < center[2])
            medium_calc = [
                center_index[1] - (12.5 / PixelSpacing[1]) * sin(angle),
                (center_index[2] + (12.5 / PixelSpacing[2]) * cos(angle)),
            ]
            low_calc = [
                center_index[1] - (25 / PixelSpacing[1]) * sin(angle),
                (center_index[2] + (25 / PixelSpacing[2]) * cos(angle)),
            ]
        elseif (center_index[1] > center[1] && center_index[2] > center[2])
            medium_calc = [
                center_index[1] - (12.5 / PixelSpacing[1]) * sin(angle),
                (center_index[2] - (12.5 / PixelSpacing[2]) * cos(angle)),
            ]
            low_calc = [
                center_index[1] - (25 / PixelSpacing[1]) * sin(angle),
                (center_index[2] - (25 / PixelSpacing[2]) * cos(angle)),
            ]
        elseif (side_x == 0 && center_index[2] < center[2])
            medium_calc = [center_index[1], center_index[2] + (12.5 / PixelSpacing[2])]
            low_calc = [center_index[1], center_index[2] + (25 / PixelSpacing[2])]
        elseif (side_x == 0 && center_index[2] > center[2])
            medium_calc = [center_index[1], center_index[2] - (12.5 / PixelSpacing[2])]
            low_calc = [center_index[1], center_index[2] - (25 / PixelSpacing[2])]
        elseif (center_index[1] > center[1] && side_y == 0)
            medium_calc = [center_index[1] - (12.5 / PixelSpacing[1]), center_index[2]]
            low_calc = [center_index[1] - (25 / PixelSpacing[1]), center_index[2]]
        elseif (center_index[1] > center[1] && side_y == 0)
            medium_calc = [center_index[1] + (12.5 / PixelSpacing[1]), center_index[2]]
            low_calc = [(center_index[1] + (25 / PixelSpacing[1])), center_index[1]]
        else
            error("unknown angle")
        end

        if size_index4 == center1
            centers[:Large_HD] = Int.(round.(center_index))
            centers[:Medium_HD] = Int.(round.(medium_calc))
            centers[:Small_HD] = Int.(round.(low_calc))

        elseif size_index4 == center2
            centers[:Large_MD] = Int.(round.(center_index))
            centers[:Medium_MD] = Int.(round.(medium_calc))
            centers[:Small_MD] = Int.(round.(low_calc))

        elseif size_index4 == center3
            centers[:Large_LD] = Int.(round.(center_index))
            centers[:Medium_LD] = Int.(round.(medium_calc))
            centers[:Small_LD] = Int.(round.(low_calc))

        else
            nothing
        end
    end
    return centers
end

# ╔═╡ 533a77b8-d541-4277-8123-cbfe416fd263
dict = calc_centers(dcm_array, output, header, center_insert, slice_CCI1)

# ╔═╡ 4c8ed7b6-392c-438b-b6d4-946f001c7980
dict[:Large_HD], dict[:Medium_HD], dict[:Small_HD], dict[:Large_MD], dict[:Medium_MD], dict[:Small_MD]

# ╔═╡ 891d81f2-f30e-4365-b0ba-b0cebae9f60f
md"""
## `mask_inserts`
"""

# ╔═╡ efad1a40-b466-4e8e-a4ae-ab3c14efc3f6
"""
	mask_inserts(
		dcm_array, masked_array, header, CCI_slice, center_insert;
		calcium_threshold=130, comp_connect=trues(3, 3)
		)

Function ...
"""
function mask_inserts(
    dcm_array,
    masked_array,
    header,
    CCI_slice,
    center_insert;
    calcium_threshold=130,
    comp_connect=trues(3, 3),
)
    output = calc_output(masked_array, header, CCI_slice, calcium_threshold, comp_connect)
    insert_centers = calc_centers(dcm_array, output, header, center_insert, CCI_slice)

    PixelSpacing = Phantoms.get_pixel_size(header)
    rows, cols = Int(header[tag"Rows"]), Int(header[tag"Columns"])

    mask_L_HD = create_circular_mask(
        cols,
		rows,
		insert_centers[:Large_HD],
		(round(5 / PixelSpacing[1], RoundUp) / 2) + 1
	)
    mask_L_MD = create_circular_mask(
        cols,
		rows,
		insert_centers[:Large_MD],
		(round(5 / PixelSpacing[1], RoundUp) / 2) + 1
    )
    mask_L_LD = create_circular_mask(
        cols,
		rows,
		insert_centers[:Large_LD],
		(round(5 / PixelSpacing[1], RoundUp) / 2) + 1
    )
    mask_M_HD = create_circular_mask(
        cols,
        rows,
        insert_centers[:Medium_HD],
        (round(3 / PixelSpacing[1], RoundUp) / 2) + 1,
    )
    mask_M_MD = create_circular_mask(
        cols,
        rows,
        insert_centers[:Medium_MD],
        (round(3 / PixelSpacing[1], RoundUp) / 2) + 1,
    )
    mask_M_LD = create_circular_mask(
        cols,
        rows,
        insert_centers[:Medium_LD],
        (round(3 / PixelSpacing[1], RoundUp) / 2) + 1,
    )
    mask_S_HD = create_circular_mask(
        cols,
		rows,
		insert_centers[:Small_HD],
		(round(1 / PixelSpacing[1], RoundUp) / 2) + 1
    )
    mask_S_MD = create_circular_mask(
        cols,
		rows,
		insert_centers[:Small_MD],
		(round(1 / PixelSpacing[1], RoundUp) / 2) + 1
    )
    mask_S_LD = create_circular_mask(
        cols,
		rows,
		insert_centers[:Small_LD],
		(round(1 / PixelSpacing[1], RoundUp) / 2) + 1
    )

    return transpose(mask_L_HD), transpose(mask_M_HD), transpose(mask_S_HD), transpose(mask_L_MD), transpose(mask_M_MD), transpose(mask_S_MD), transpose(mask_L_LD), transpose(mask_M_LD), transpose(mask_S_LD)
end

# ╔═╡ 9cf4a86b-d41f-4814-b14e-0854176e73d4
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);

# ╔═╡ 492574dc-5674-462d-9df1-a77e0f445994
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;

# ╔═╡ 33b38675-7434-46e5-808c-ad4b0cbb0408
heatmap(transpose(mask_L_HD), colormap=:grays)

# ╔═╡ 4ff3a9d0-9681-4d49-a165-6787245a2a9f
@bind a3 PlutoUI.Slider(1:size(masked_array, 3), default=26, show_value=true)

# ╔═╡ 105d699f-1d2c-4291-b5c2-3b7f18d81012
heatmap(transpose(masked_array[:, :, a3]), colormap=:grays)

# ╔═╡ 6e7c6d1a-f25d-46f9-88c7-85400d724736
md"""
## Visualize
"""

# ╔═╡ 3cdf5136-6315-45b2-9f7a-b8868f6436d4
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

# ╔═╡ bcdb7fe9-b699-48ae-a260-f7d4636a3ec7
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=25, show_value=true)
end

# ╔═╡ 0d7b6f09-badd-45fc-936f-ff516c622909
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.5, color=:red)
	fig
end

# ╔═╡ ee5cd3f6-0549-48fa-95dd-ab4f88308a31
begin
	masks_3D = Array{Bool}(undef, size(dcm_array))
	for z in 1:size(dcm_array, 3)
		masks_3D[:, :, z] = masks
	end
end;

# ╔═╡ 5a799514-6d20-43f3-adf9-7be531495108
@bind v1 overlay_mask_bind(masks_3D)

# ╔═╡ 3d39f7b6-1f73-4783-9ee0-b770d7b35988
overlay_mask_plot(dcm_array, masks_3D, v1, "masks overlayed")

# ╔═╡ 6d5f62e8-fdda-42df-a253-9b55d1395a22
slice_CCI

# ╔═╡ Cell order:
# ╠═14b6c3e4-48dd-11ec-336c-6918bf024f43
# ╠═9506ae45-f63f-4fd6-8467-0597cbde7006
# ╟─482d8669-c750-450b-aae8-03afc2a90007
# ╠═f53772f8-d9f0-4562-896d-6a6112cd8735
# ╟─3f313b73-2785-4299-a1e9-bc5d4e1ad422
# ╠═202cef5d-85f1-4673-9e6c-52ba80e44de0
# ╠═f6a29f28-b787-43d6-8b6b-5e163cbbbcec
# ╟─39688933-f68f-47b4-a9dc-aff25bdc6123
# ╠═0883a26a-54ba-4670-a19b-a6a0a23f5aee
# ╠═59976bc3-271e-4418-874c-b553b59e730a
# ╠═8652d37a-7dc2-4123-b20f-87a8981baa26
# ╠═9b02d4cc-e9b3-47aa-8d57-e0befd64a5fc
# ╠═2185c4e4-a518-4133-8174-662e0dd744eb
# ╠═d9bbaeb8-b027-4a58-9ddb-0b015aa867de
# ╠═7d7b1fd5-e244-4010-ac61-6f4e1925f216
# ╟─237cfd63-3479-468d-954f-2ebd05376fca
# ╟─1b755e07-926a-42bd-871e-a693bb037c81
# ╠═d2d1d02c-e592-4c73-afa2-80f89dd5534e
# ╠═2d16b7cb-584e-4ff8-84a1-4e5670bd49d3
# ╟─a69293d3-6093-4b87-abc6-11f1ec9398c3
# ╠═dfcda44d-5b2c-4a50-9fb2-44cb70b5488e
# ╠═57730e9f-05b2-4e3a-bff3-a3a6403a66e0
# ╠═19e0f1b3-86c3-4ffd-aa39-dd7e1e509312
# ╟─55665585-8076-4e23-a989-71fc9909a57b
# ╟─78ffb833-27c5-4df5-9989-47d57aca47fe
# ╟─1ac4e0fd-79fb-41d1-b2c7-8b5e011c80e2
# ╟─5147ce8d-2d44-443b-a6a2-edba166356a8
# ╠═004f4a1e-9d11-44ef-8622-69c311edb2cb
# ╟─da57bd3c-0d8b-453c-9c4f-ea846850cbec
# ╟─11e0c502-5b1f-4e60-bb43-a8c04a4e335b
# ╠═c4c6841a-88f4-4ca9-9353-0e9daa6d4dff
# ╠═462fd944-d4ce-4d7b-82ad-40342daff354
# ╠═8959750c-291c-4c2a-a6a1-218f8ed82ab1
# ╟─68c16f9e-8520-41ef-912a-9e61acf2ad1d
# ╠═f49fe2a8-8429-4fd1-98e3-6e75189a3ac2
# ╠═494163b1-cb52-45a9-815e-991633018cd5
# ╟─880db52b-e39c-4578-a6de-ce4a6afa038e
# ╠═6613bd97-7d53-4b8b-a69c-12e040f4aa36
# ╠═bff84e32-1013-44f2-9eda-eb5f6bac4c5a
# ╟─451dd3c6-f740-47df-a5b9-85749caf04fb
# ╠═1ac3d601-0057-475e-a6ac-9522618f7bb9
# ╠═5d407e49-fe05-4748-be36-b232b5ce3ccb
# ╠═de1ac95e-2399-45c1-b337-dcfe800394a2
# ╟─a8f6ece5-9f59-4af1-bcf7-7dafa839bae2
# ╠═267d1d6a-5c35-480e-8425-f000f7f31c7f
# ╟─9a77474e-a853-4369-a6f3-a9e76a683a27
# ╠═a9f499f8-ca32-4310-82d7-f20b790dc04c
# ╠═4eece882-8416-4330-b95a-1cfaeac89aac
# ╠═d6cb5578-e924-4e78-be0b-82d8866fbee8
# ╟─214d485e-df7e-4a45-b52f-10f854d839a1
# ╠═03622243-d2af-4d1b-a6e1-fd50670e48a8
# ╠═93861c7a-4086-4404-842e-ca91c2ad6094
# ╟─ee208b64-3cf3-428b-a5f6-74793a82fd93
# ╟─5930b9c5-2f85-465b-abe0-223580142d49
# ╠═480bd663-d822-4417-b7f0-6679d26a58b4
# ╠═f498c292-367d-4441-873e-29bec99af2a6
# ╟─56ce9966-939e-48cd-9c70-2e441bfaea58
# ╠═228dd95c-fdab-4a5a-a079-e14238ff5e83
# ╠═0f572e01-a690-4699-8a8c-6684c27085bb
# ╠═ce5b9c0a-4f25-4b8f-bf16-17d05b7bfff4
# ╟─5a8963f1-3ec8-4dbe-b567-df2648219b7f
# ╠═92fca992-b36d-48f3-803e-f92e6463c291
# ╠═fb320178-0d33-4402-b075-7bf7541f94e4
# ╠═68dae231-9907-44b2-afca-a5dac57998bb
# ╠═14f46c1f-62fa-45ff-9c78-b01404c466da
# ╟─feca7b76-fc10-4aca-adff-4cd323dabc74
# ╠═ae9ab557-553c-474a-8865-37807b9270fb
# ╠═39ce8635-28a7-4eee-875e-325e85d142d5
# ╠═42e86a48-1438-4edd-ac2d-f307d534e165
# ╟─0f73d3a7-c6fa-49ee-8fdf-357b0fe94ec8
# ╠═9ac2f9ed-92bb-4540-add5-422df2a3f6da
# ╠═533a77b8-d541-4277-8123-cbfe416fd263
# ╠═4c8ed7b6-392c-438b-b6d4-946f001c7980
# ╟─891d81f2-f30e-4365-b0ba-b0cebae9f60f
# ╠═efad1a40-b466-4e8e-a4ae-ab3c14efc3f6
# ╠═9cf4a86b-d41f-4814-b14e-0854176e73d4
# ╠═492574dc-5674-462d-9df1-a77e0f445994
# ╠═33b38675-7434-46e5-808c-ad4b0cbb0408
# ╟─4ff3a9d0-9681-4d49-a165-6787245a2a9f
# ╠═105d699f-1d2c-4291-b5c2-3b7f18d81012
# ╟─6e7c6d1a-f25d-46f9-88c7-85400d724736
# ╠═3cdf5136-6315-45b2-9f7a-b8868f6436d4
# ╠═bcdb7fe9-b699-48ae-a260-f7d4636a3ec7
# ╠═0d7b6f09-badd-45fc-936f-ff516c622909
# ╠═ee5cd3f6-0549-48fa-95dd-ab4f88308a31
# ╟─5a799514-6d20-43f3-adf9-7be531495108
# ╠═3d39f7b6-1f73-4783-9ee0-b770d7b35988
# ╠═6d5f62e8-fdda-42df-a253-9b55d1395a22

max_stars_repo_name: MolloiLab/CAC-stanford-data

max_issues_repo_name: MolloiLab/CAC-stanford-data

max_forks_repo_name: MolloiLab/CAC-stanford-data


text:
greater(g::TimeOrderedGF) = g.G
lesser(g::TimeOrderedGF) = g.L
retarded(g::TimeOrderedGF) = g.R
advanced(g::TimeOrderedGF) = adjoint(retarded(g))

greater(g::TimeOrderedConvolution) = (retarded(g.A) .* transpose(g.dts)) * greater(g.B) + greater(g.A) * (g.dts .* advanced(g.B)) |> skew_hermitify!
lesser(g::TimeOrderedConvolution) = (retarded(g.A) .* transpose(g.dts)) * lesser(g.B) + lesser(g.A) * (g.dts .* advanced(g.B)) |> skew_hermitify!


text:
__precompile__()

module Equations

    export Equation, ODE, PODE, SODE, DAE, PDAE, do_nothing

    include("equations/equations.jl")

end


text:
function save_noise!(W::NoiseWrapper)

end

function interpolate!(W::NoiseWrapper,t)
  W.source(t)
end

function interpolate!(out1,out2,W::NoiseWrapper,t)
  W.source(out1,out2,t)
end

function calculate_step!(W::NoiseWrapper,dt)
  if isinplace(W)
    W(W.dW,W.dZ,W.curt+dt)
    W.dW .-= W.curW
    if W.Z != nothing
      W.dZ .-= W.curZ
    end
  else
    new_W, new_Z = W(W.curt+dt)
    W.dW = new_W - W.curW
    if W.Z != nothing
      W.dZ = new_Z - W.curZ
    end
  end
  W.dt = dt
end

function accept_step!(W::NoiseWrapper,dt,setup_next=true)
  if isinplace(W)
    W.curW .+= W.dW
  else
    W.curW += W.dW
  end
  W.curt += W.dt
  push!(W.W,copy(W.curW))
  push!(W.t,copy(W.curt))
  if W.Z != nothing
    if isinplace(W)
      W.curZ .+= W.dZ
    else
      W.curZ += W.dZ
    end

    push!(W.Z,copy(W.curZ))
  end

  W.dt = dt #dtpropose
  if setup_next
    calculate_step!(W,dt)
  end
end

function reject_step!(W::NoiseWrapper,dtnew)
  calculate_step!(W::NoiseWrapper,dtnew)
end

function setup_next_step!(W::NoiseWrapper)
  calculate_step!(W,W.dt)
end


text:
using Base.Test

# @testset "Deep-Autoencoder with *cross entropy* and *BetaCV* loss" begin

using Flux
using Flux.Tracker
using Flux: crossentropy, throttle

using FluxMO.Clustering: knn_clustering
using FluxMO.Validation: betacv, betacv_pairwise

using Plots
gr()

function train(seed::Int; mode = :with_betacv)

    srand(seed)

    N = 1000    # samples
    M = 100     # data size
    L = 2       # latent space size

    # generate correlated data
    D = map(_->abs.(cumsum(rand(M)./round(Int,N/2))),1:N)
    # cor(D[1],D[2])

    X = deepcopy(D[1:end-1])
    Y = deepcopy(D[2:end])

    a = sin
    s = 1
    m = Chain(
        Dense(    M, 100*s, a),
        Dense(100*s,  10*s, a),
        Dense( 10*s,     L, a),
        Dense(    L,  10*s, a),
        Dense( 10*s, 100*s, a),
        Dense(100*s,     M, sigmoid),
    )
    embd_layer = 3

    BCV_take = 20

    epochs = 1

    kNN = 20

    randomize = false

    function supervised_betacv(i::Int)

        # Embed all samples from X into the latent space of size L
        embds_tracked = map(x-> m[1:embd_layer](x), X)
        # 999-element Array{TrackedArray{…,Array{Float64,1}},1}:
        # param([0.63072, -0.046672])
        # ⋮

        # untracked embedded points as matrix
        embds  = hcat(map(ta->ta.data, embds_tracked)...)
        # 2×999 Array{Float64,2}:

        # Cluster the current untracked embeddings to generate
        # a supervised scenario.
        # This clould also be done by DBSCAN, OPTICS, K-Means ...
        clustering, _ = knn_clustering(embds, k = kNN)
        # 221-element Array{Array{Int64,1},1}:
        # Array{Int64,1}[[350, 837, 466, 364, 600, 964, 271, 976, 1, 804], …

        # take only n first clusters to reduce backtracking...
        take = 1:BCV_take

        # take only n random clusters to reduce backtracking...
        if randomize
            cn = length(clustering)
            take = unique(rand(1:cn, min(BCV_take,cn)))
        end

        # Obtain tracked embedded values from the clustering
        embds_clustered_tracked = map(c->map(i->embds_tracked[i],c), clustering[take])
        # 221-element Array{Array{TrackedArray{…,Array{Float64,1}},1},1}:
        # TrackedArray{…,Array{Float64,1}}[param([0.311645, 0.0285734]), param([0.309032, 0.019196]), …]
        # ⋮

        # Validate the clustering via BetaCV measure (small is good).
        # This is done with the tracked values, as it should influence
        # the model weights and biases towards optimizing this measure.
        bcv = @time betacv(embds_clustered_tracked)
        # Flux.Tracker.TrackedReal{Float64}

        # untracked betacv calculates fine.
        # embds_clustered = map(c->map(i->embds[:,i],c), clustering)
        # bcv = betacv(embds_clustered)

        # embds_clustered_tracked_matrix = map(c->hcat(map(i->embds_tracked[i],c)), clustering[take])

        # try classic betacv with for loops on tracked matrix
        # bcvm = @time betacv(embds_clustered_tracked_matrix)

        # try matrix optimized betacv for better backpropagation
        # bcvp = @time betacv_pairwise(embds_clustered_tracked_matrix)
        # ERROR: MethodError: no method matching one(::Type{TrackedArray{…,Array{Float64,1}}})
        # Closest candidates are:
        # one(::Type{Measures.Length{:mm,Float64}}) at /home/sebastian/.julia/v0.6/Plots/src/layouts.jl:13
        # one(::Type{Measures.Length{:pct,Float64}}) at /home/sebastian/.julia/v0.6/Plots/src/layouts.jl:31
        # one(::BitArray{2}) at bitarray.jl:427
        # ...
        # Stacktrace:
        # [1] result_type at /home/sebastian/.julia/v0.6/Distances/src/metrics.jl:194 [inlined]
        # [2] pairwise(::Distances.Euclidean, ::Array{TrackedArray{…,Array{Float64,1}},2}, ::Array{TrackedArray{…,Array{Float64,1}},2}) at /home/sebastian/.julia/v0.6/Distances/src/generic.jl:120
        # [3] intra_cluster_weights_pairwise(::Array{Array{TrackedArray{…,Array{Float64,1}},2},1}) at /home/sebastian/develop/julia/dev/FluxMO/src/betacv.jl:48
        # [4] betacv_pairwise(::Array{Array{TrackedArray{…,Array{Float64,1}},2},1}) at /home/sebastian/develop/julia/dev/FluxMO/src/betacv.jl:93
        # [5] macro expansion at ./util.jl:237 [inlined]
        # [6] (::#supervised_betacv#11{Array{Array{Float64,1},1},Flux.Chain,Int64})(::Int64) at ./none:63
        # [7] (::#loss#18{Symbol,Flux.Chain,#supervised_betacv#11{Array{Array{Float64,1},1},Flux.Chain,Int64}})(::Array{Float64,1}, ::Array{Float64,1}, ::Int64) at ./none:87
        # [8] #train!#130(::Flux.#throttled#14, ::Function, ::Function, ::Base.Iterators.Zip{Array{Array{Float64,1},1},Base.Iterators.Zip2{Array{Array{Float64,1},1},Array{Int64,1}}}, ::Flux.Optimise.##71#75) at /home/sebastian/.julia/v0.6/Flux/src/optimise/train.jl:39
        # [9] (::Flux.Optimise.#kw##train!)(::Array{Any,1}, ::Flux.Optimise.#train!, ::Function, ::Base.Iterators.Zip{Array{Array{Float64,1},1},Base.Iterators.Zip2{Array{Array{Float64,1},1},Array{Int64,1}}}, ::Function) at ./<missing>:0
        # [10] #train#1(::Symbol, ::Function, ::Int64) at ./none:110
        # [11] train() at ./none:3
        # [12] eval(::Module, ::Any) at ./boot.jl:235

        bcv
    end

    last_bcv = param(0.0)
    last_ce = param(0.0)

    function loss(x, y, i = -1)

        # unsupervised
        ŷ = m(x)
        prev_ce = deepcopy(last_ce.tracker.data)
        last_ce = crossentropy(ŷ, y)

        # supervised

        # NOTE:
        # This is where I like to apply the second optimization...
        # How to backtrack this function efficently?
        #   * Do I have to treat it as second model with an supervised setting?
        #   * or can I apply the `betacv` as intendet to the existing model?

        if mode == :with_betacv
            # apply after a fraction of all seen samples
            # do not apply when called manually from the callback
            if i != -1 && i % 100 == 0
            # if i != -1 && i % floor(Int, (N-1)/3) == 0
                println("betacv() after $i samples $kNN $BCV_take \t")
                prev_bcv = deepcopy(last_bcv.tracker.data)

                last_bcv = supervised_betacv(i)
                println("crossentropy: ", last_ce,  "\tdiff: ", last_ce.tracker.data  - prev_ce)
                println("betacv:       ", last_bcv, "\tdiff: ", last_bcv.tracker.data - prev_bcv)
            end
        end

        # NOTE: Does this even help to reduce the Stacktrace?
        Flux.truncate!(m)
        # Flux.reset!(m)

        # may fail for `tanh` as activation function...
        if isnan(last_ce)
            # @show x, ŷ, y, i
            @show last_ce, prev_ce
            @show i
            @show ŷ
        end

        # optimize both metrics
        last_ce + last_bcv
    end

    opt = Flux.ADAM(params(m))

    function callback()
        # sqrt(N) random samples for approx. loss
        ns = rand(1:N-1, round(Int, sqrt(N)))
        ls = [loss(X[n], Y[n]) for n in ns]
        println("Train: loss:  ", sum(ls)/length(ls), "\tstd:  ", std(ls))

        # one embedding example
        embd = m[1:embd_layer](X[10]).data
        n = rand(1:N-1)
        embd_r = m[1:embd_layer](X[n]).data
        println("Embedding:    ", embd," (10)   \t",embd_r," ($n)")

        println()
    end

    for epoch in 1:epochs
        info("Epoch: $epoch (Seed: $seed)")
        Flux.train!(loss, zip(X, Y, collect(1:length(X))), opt, cb=throttle(callback,3))
    end

    return X,Y,m,last_ce.tracker.data, [seed,N,M,L,embd_layer,a,BCV_take,epochs,kNN,randomize]
end

seed = rand(1:10000)
# seed = 8651  # fails after 500 samples...

X_bcv,Y_bcv,model_bcv,ce_bcv,opts = train(seed)
data_bcv = deepcopy(hcat(map(x->model_bcv[1:3](x).data, X_bcv)...))
clustering, _ = knn_clustering(data_bcv)
data_clustered = map(c->hcat(map(i->data_bcv[:,i],c)), clustering)
bcv_bcv = betacv(data_clustered)

ce_bcv_s = @sprintf "%1.2e" ce_bcv
bcv_bcv_s = @sprintf "%1.5f" bcv_bcv

X_ce,Y_ce,model_ce,ce_ce,_ = train(seed, mode = :other)
data_ce = deepcopy(hcat(map(x->model_ce[1:3](x).data, X_ce)...))
clustering, _ = knn_clustering(data_ce)
data_clustered = map(c->hcat(map(i->data_ce[:,i],c)), clustering)
bcv_ce = betacv(data_clustered)
ce_ce_s = @sprintf "%1.2e" ce_ce
bcv_ce_s = @sprintf "%1.5f" bcv_ce

p = plot(
    scatter(data_bcv[1,:], data_bcv[2,:],
    label=["embedded (ce+bcv)"],
    title="crossentropy: $ce_bcv_s\nbetacv: $bcv_bcv_s"),

    scatter(data_ce[1,:], data_ce[2,:],
    label=["embedded (ce)"],
    title="crossentropy: $ce_ce_s\nbetacv: $bcv_ce_s"),
)

savefig(p, string(now(), "_", opts, "_", ".png"))
# end

max_stars_repo_name: swiesend/FluxMO

max_issues_repo_name: swiesend/FluxMO

max_forks_repo_name: swiesend/FluxMO