text:
export NearNeighborGraph

# NeighborIterator
## Hacky, more performant alternative to Iterators.flatten((neighborhood, (D[v] for D in extras if v in keys(D))))
## Currently, to ensure iterate returns a Union{Nothing,XXX}, if isempty(neighborhood) then no extras will be considered
struct NeighborIterator{T,X,N,D}
    v::X
    neighborhood::N
    extras::Vector{D}
end
function NeighborIterator{T}(v::X, neighborhood::N, extras::Vector{D}) where {T,X,N,D}
    NeighborIterator{T,X,N,D}(v, neighborhood, extras)
end
Base.IteratorSize(::NeighborIterator) = Base.SizeUnknown()
Base.eltype(::Type{<:NeighborIterator{T}}) where {T} = T
@inline function _start(itr::NeighborIterator)
    iter_result = iterate(itr.neighborhood)
    iter_result === nothing ? nothing : (iter_result..., 0)
end
@inline function Base.iterate(itr::NeighborIterator, state=_start(itr))
    state === nothing && return nothing
    if state[3] == 0
        iter_result = iterate(itr.neighborhood, state[2])
        iter_result === nothing ? (state[1], (state[1], state[2], 1)) : (state[1], (iter_result..., 0))
    else
        for i in state[3]:length(itr.extras)
            itr.v in keys(itr.extras[i]) && return (itr.extras[i][itr.v], (state[1], state[2], i + 1))
        end
        nothing
    end
end

struct NearNeighborGraph{T,SS<:SampleSet,BVP<:SteeringBVP,NC<:NeighborhoodCache,DS<:NearNeighborDataStructure,D}
    nodes::SS
    bvp::BVP
    edge_cache::NC
    NN_data_structure::DS
    init_neighbors::Vector{D}
    goal_neighbors::Vector{D}
end
function NearNeighborGraph{T}(nodes::SS, bvp::BVP, edge_cache::NC, NN_data_structure::DS,
                              init_neighbors::Vector{D}, goal_neighbors::Vector{D}) where {T,SS,BVP,NC,DS,D}
    NearNeighborGraph{T,SS,BVP,NC,DS,D}(nodes, bvp, edge_cache, NN_data_structure, init_neighbors, goal_neighbors)
end
function NearNeighborGraph(nodes::SampleSet,
                           bvp::SteeringBVP,
                           edge_cache::NeighborhoodCache=MemoizedNC(nodes, bvp),
                           NN_data_structure::NearNeighborDataStructure=NullNNDS();
                           include_controls::BoolVal=Val(includes_controls(edge_cache)))
    T = neighbor_info_type(nodes, bvp, include_controls=include_controls)
    D = Dict{indextype(nodes),T}
    init_neighbors = Vector{D}(undef, 1)
    goal_neighbors = Vector{D}(undef, 0)
    NearNeighborGraph{T}(nodes, bvp, edge_cache, NN_data_structure, init_neighbors, goal_neighbors)
end
Base.getindex(G::NearNeighborGraph, i) = G.nodes[i]
# Base.getindex(G::NearNeighborGraph, i, j) = SteeringEdge
function setinit!(G::NearNeighborGraph{NeighborInfo{X,D,U}}, x::State, r) where {X,D,U}
    G.nodes.init = x    # TODO: also possibly reset!(G.edge_cache); add r as a field probably (RRT = nothing)
    G.init_neighbors[1] = reverse_neighbor_dict(X(0), neighbors(G, X(0), r, dir=Val(:F)), NeighborInfo{X,D,U})
end
function setgoal!(G::NearNeighborGraph{NeighborInfo{X,D,U}}, f::Base.Callable, r; count=1) where {X,D,U}
    resize!(G.nodes.goal_samples, count)
    resize!(G.goal_neighbors, count)
    for i in 1:count
        G.nodes.goal_samples[i] = f()
        G.goal_neighbors[i] = reverse_neighbor_dict(X(-i), neighbors(G, X(-i), r, dir=Val(:B)), NeighborInfo{X,D,U})
    end
end
setgoal!(G::NearNeighborGraph, x, r) = setgoal!(G, () -> x, r, count=1)
neighbor_info_type(G::NearNeighborGraph{T}) where {T} = T
indextype(G::NearNeighborGraph) = indextype(G.nodes)
# costtype(G::NearNeighborGraph) = costtype(neighbor_info_type(G))
# controlstype(G::NearNeighborGraph) = controlstype(neighbor_info_type(G))
includes_controls(G::NearNeighborGraph) = includes_controls(neighbor_info_type(G))

function neighbors(G::NearNeighborGraph{T}, v::Int, r; dir::F_or_B=Val(:F)) where {T}
    neighborhood = neighbors(G.edge_cache, G.NN_data_structure, G.nodes, G.bvp, v, r,
                             dir=dir, include_controls=Val(includes_controls(T)))
    NeighborIterator{T}(v, neighborhood, dir === Val(:F) ? G.goal_neighbors : G.init_neighbors)
    # if dir === Val(:F)
    #     Iterators.flatten((neighborhood, (nbhd[v] for nbhd in G.goal_neighbors if v in keys(nbhd))))
    # else
    #     Iterators.flatten((neighborhood, (nbhd[v] for nbhd in G.init_neighbors if v in keys(nbhd))))
    # end
end
function neighbors(G::NearNeighborGraph{T,<:TiledSampleSet}, v::TiledIndex, r; dir::F_or_B=Val(:F)) where {T}
    on, off = v.on_lattice, v.off_lattice
    vz = zero_lattice_component(v)
    lattice_neighborhood = neighbors(G.edge_cache, G.NN_data_structure, G.nodes, G.bvp, vz, r,
                                     dir=dir, include_controls=Val(includes_controls(T)))
    neighborhood = (shift_lattice_component(n, on) for n in lattice_neighborhood)
    NeighborIterator{T}(v, neighborhood, dir === Val(:F) ? G.goal_neighbors : G.init_neighbors)
    # if dir === Val(:F)
    #     Iterators.flatten((neighborhood, (nbhd[v] for nbhd in G.goal_neighbors if v in keys(nbhd))))
    # else
    #     Iterators.flatten((neighborhood, (nbhd[v] for nbhd in G.init_neighbors if v in keys(nbhd))))
    # end
end
function neighbors(f::Function, G::NearNeighborGraph, v::SampleIndex, r; dir::F_or_B=Val(:F))
    Iterators.filter(f, neighbors(G, v, r, dir=dir))
end

max_stars_repo_name: pbpf/clean

max_issues_repo_name: pbpf/clean

max_forks_repo_name: pbpf/clean


text:
using KhepriRevit
using Test

@testset "KhepriRevit.jl" begin
    # Write your tests here.
end

ext: jl


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

# ╔═╡ 651419e4-4c13-4ddd-8e69-5cdf2b995a51
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

# ╔═╡ d9401636-eb09-45b2-9663-b8ab615531f9
TableOfContents()

# ╔═╡ ceb200df-c17b-476f-ae95-a13e572a243e
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.
"""

# ╔═╡ 35c2524a-28a4-49eb-9fe2-e4232084af1a
begin
	SCAN_NUMBER = 1
	VENDER = "Canon_Aquilion_One_Vision"
	BASE_PATH = "/Users/daleblack/Google Drive/Datasets/"
end

# ╔═╡ 7c5f2ec6-26bb-4a2a-9dd1-a7f3b0de10b0
md"""
**Everything below should be automatic, just scroll through to visually inspect that things make sense**
"""

# ╔═╡ 7dab92d6-13c3-4c1e-8b51-ad80a709a80b
root_path = string(BASE_PATH, VENDER)

# ╔═╡ 9bcfc64a-170b-461c-baa1-ac5c1ab73579
dcm_path_list = dcm_list_builder(root_path)

# ╔═╡ 041192c8-875e-43f3-881e-dc993955c91f
pth = dcm_path_list[SCAN_NUMBER]

# ╔═╡ 57cd84bb-9307-44eb-87c1-9a663e10e5c2
scan = basename(pth)

# ╔═╡ 807cfdd5-379a-4c65-ab15-534f1e28b65c
header, dcm_array, slice_thick_ori1 = dcm_reader(pth);

# ╔═╡ 948b9234-f845-48f5-a246-18cbc961ad67
md"""
## Helper Functions
"""

# ╔═╡ baa55458-c8e2-4852-b486-ef8c854f77bf
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

# ╔═╡ 7434ae0d-06f0-44fc-b9e2-198fb9afb97c
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

# ╔═╡ 7325a8d4-6e0a-4087-ad64-d0874fb26aad
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=1, color=:red)
	fig
end

# ╔═╡ 63c84140-6510-4186-a0bf-dab85d19e2cc
function create_mask(array, mask)
	@assert size(array) == size(mask)
	idxs = findall(x -> x == true, mask)
	overlayed_mask = zeros(size(array))
	for idx in idxs
		overlayed_mask[idx] = array[idx]
	end
	return overlayed_mask
end

# ╔═╡ a9516b63-80c0-42ff-9bac-a9ad311fbe9c
md"""
## Segment Heart
"""

# ╔═╡ fb4542a8-3406-48bb-b255-ebde03adb951
masked_array, center_insert, mask = mask_heart(header, dcm_array, size(dcm_array, 3)÷2);

# ╔═╡ 835490d8-332c-400a-a2f7-2dcff791d3d6
@bind a PlutoUI.Slider(1:size(masked_array, 3), default=10, show_value=true)

# ╔═╡ a2172056-1f1c-40c9-ac6f-816c80dd09ef
heatmap(masked_array[:, :, a], colormap=:grays)

# ╔═╡ 38042368-0be5-4c64-a4e6-c596201cdd2e
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

# ╔═╡ a467ddd4-e2d3-445e-808b-298fea9b927c
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

# ╔═╡ dfb18ee0-7aa7-43c3-bf7e-914b2ee4b729
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

# ╔═╡ 65b0a9dc-efc2-4eb7-9995-1b8bd8113a63
md"""
## Segment Calcium Rod
"""

# ╔═╡ 8ae2afdc-ae49-49e6-ba59-d20373907ba7
calcium_image, slice_CCI, quality_slice, cal_rod_slice = mask_rod(masked_array, header);

# ╔═╡ d635ed56-ef46-428a-ac7e-f2fc5c10e4d6
@bind c PlutoUI.Slider(1:size(calcium_image, 3), default=cal_rod_slice, show_value=true)

# ╔═╡ e4052629-fe1c-4b9f-8d74-fe7dce34de06
heatmap(transpose(calcium_image[:, :, c]), colormap=:grays)

# ╔═╡ 630ec973-b16d-46b3-bb3b-229d407087b1
cal_rod_slice

# ╔═╡ bdb61513-e02a-43df-be3e-164f415e061b
md"""
## Segment Calcium Inserts
"""

# ╔═╡ 5ec27948-098e-4acf-961a-0bc36307f807
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
);

# ╔═╡ 480f30c9-b0cb-4f49-8767-717d6213902e
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;

# ╔═╡ 02cd946e-98a3-4043-b482-7a19dcd1f157
heatmap(masks, colormap=:grays)

# ╔═╡ 84f78e6d-c73f-4e74-bfb2-1de22990d1c8
md"""
## Mass cal factor
"""

# ╔═╡ fe8c8d86-816d-46f8-9068-663c4a45ce97
begin
	output = calc_output(masked_array, header, slice_CCI, 130, trues(3, 3))
    insert_centers = calc_centers(dcm_array, output, header, center_insert, slice_CCI)
end

# ╔═╡ 0604f5cf-a731-405e-9c3c-a9f918352f96
center_large_LD = insert_centers[:Large_LD]

# ╔═╡ fd78b7f0-a83e-4365-a770-16483010376a
rows, cols = Int(header[tag"Rows"]), Int(header[tag"Columns"])

# ╔═╡ 8e1ca1c1-ea27-4f0e-89cd-b466053e9485
pixel_size = DICOMUtils.get_pixel_size(header)

# ╔═╡ 3c53f529-9fe6-47e9-86f9-41de9ee92b87
function mass_calibration(
    dcm_array, center_large_LD, center, cal_rod_slice, rows, cols, spacing
)
    center_LD = center_large_LD
    dist_x = abs(center_LD[1] - center[1])
    dist_y = abs(center_LD[2] - center[2])

    if dist_x == 0
        mass_center_x = center[1]
        if center_LD[2] > center[2]
            mass_center_y = round(center[1] - round(23 / spacing[1], RoundUp), RoundUp)
        else
            mass_center_y = round(center[1] + round(23 / spacing[1], RoundUp), RoundUp)
        end
    elseif dist_y == 0
        mass_center_y = center[2]
        if center_LD[1] > center[1]
            mass_center_x = round(center[1] - round(23 / spacing[1], RoundUp), RoundUp)
        else
            mass_center_x = round(center[0] + round(23 / spacing[1], RoundUp), RoundUp)
        end

    else
        mass_angle = atan(dist_y / dist_x)
        dist_x = (23 / spacing[1]) * cos(mass_angle)
        dist_y = (23 / spacing[1]) * sin(mass_angle)

        if (center_LD[1] < center[1] && center_LD[2] < center[2])
            mass_center_x = round(center[1] + dist_x, RoundUp)
            mass_center_y = round(center[2] + dist_y, RoundUp)
        elseif (center_LD[1] < center[1] && center_LD[2] > center[2])
            mass_center_x = round(center[1] + dist_x, RoundUp)
            mass_center_y = round(center[2] - dist_y, RoundUp)
        elseif (center_LD[1] > center[1] && center_LD[2] < center[2])
            mass_center_x = round(center[1] - dist_x, RoundUp)
            mass_center_y = round(center[2] + dist_y, RoundUp)
        elseif (center_LD[1] > center[1] && center_LD[2] > center[2])
            mass_center_x = round(center[1] - dist_x, RoundUp)
            mass_center_y = round(center[2] - dist_y, RoundUp)
        end
    end

    mass_cal_center = [mass_center_y, mass_center_x]
    x_distance = abs(center[1] - mass_cal_center[2])
    angled_distance = sqrt(
        (center[1] - mass_cal_center[2])^2 + (center[2] - mass_cal_center[2])^2
    )
    angle_0_200HA = acos(x_distance / angled_distance) * 180 / π
    mask_0HU = PhantomSegmentation.create_circular_mask(
        cols, rows, mass_cal_center, Int(round(6.9 / spacing[1]))
    )
	@info size(mask_0HU)
    masked_0HU = mask_0HU .* dcm_array[:, :, cal_rod_slice]
    nonzero_count = length(findall(x -> x != 0, masked_0HU))
    mean_0HU = sum(masked_0HU) / nonzero_count
    std_0HU = 0

    for voxel in vec(masked_0HU)
        if voxel != 0
            std_0HU += (voxel - mean_0HU)^2
        end
    end

    std_0HU = sqrt(std_0HU / nonzero_count)
    mask_200HU = PhantomSegmentation.create_circular_mask(
        cols, rows, (center[2], center[1]), Int(round(6.9 / spacing[1]))
    )
    masked_200HU = mask_200HU .* dcm_array[:, :, cal_rod_slice]
    nonzero_count_200HU = length(findall(x -> x != 0, masked_200HU))
    mean_200HU = sum(masked_200HU) / nonzero_count_200HU
    mass_cal_factor = 0.2 / (mean_200HU - mean_0HU)
    water_rod_metrics = mean_0HU, std_0HU

    return mass_cal_factor, angle_0_200HA, water_rod_metrics
end

# ╔═╡ 64d8f848-cc5b-4b51-90bd-aacca098593b
mass_cal_factor, angle_0_200HA, water_rod_metrics = mass_calibration(masked_array, center_large_LD, center_insert, cal_rod_slice, rows, cols, pixel_size)

# ╔═╡ 11e23023-d88b-4509-a562-51d5255c6b76
md"""
# Score Large Inserts
"""

# ╔═╡ 9c14d354-9f0d-45d5-8dc2-98a719b0465c
arr = masked_array[:, :, slice_CCI-2:slice_CCI+2];

# ╔═╡ 5ee01db8-2a9d-491d-8305-adff24c6dd09
single_arr = masked_array[:, :, slice_CCI];

# ╔═╡ c0f2a223-d737-4823-80c0-76801ae26745
md"""
## High Density
"""

# ╔═╡ c986c5a2-dffb-4200-af4b-fde0e94c5361
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;

# ╔═╡ cc598eba-d013-4984-aa5a-6f94f8f03f19
md"""
#### Dilated mask
"""

# ╔═╡ 630c79f5-16c3-4ce7-b197-16348783c2a0
dilated_mask_L_HD = dilate(dilate(mask_L_HD_3D));

# ╔═╡ cf80857d-f58d-46a8-b7fc-f44dfbd8d2a6
@bind g2 overlay_mask_bind(dilated_mask_L_HD)

# ╔═╡ 5e82dcc1-431f-4ba7-9adf-df62beaa0ca9
overlay_mask_plot(arr, dilated_mask_L_HD, g2, "dilated mask")

# ╔═╡ 53b8bada-351b-42d1-8807-79ed375b8186
alg = Agatston()

# ╔═╡ dc41660b-6620-4cd6-859b-b5fc3acc7063
overlayed_mask_l_hd = create_mask(arr, dilated_mask_L_HD);

# ╔═╡ 4cd5d9bd-1b56-4cd4-b46c-168cb3da6b4d
agat_l_hd, mass_l_hd = score(overlayed_mask_l_hd, pixel_size, mass_cal_factor, alg)

# ╔═╡ c164f2e6-edea-45f6-8cde-4da1ba1820a0
md"""
## Medium Density
"""

# ╔═╡ c1ba29be-fad6-40f5-a828-122ea002ee74
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;

# ╔═╡ 4b2fe0ac-f265-4a29-9ef0-b452e5a66d22
md"""
#### Dilated mask
"""

# ╔═╡ 4493b19e-bfb7-42da-9d5d-457b30707f14
dilated_mask_L_MD = dilate(dilate(mask_L_MD_3D));

# ╔═╡ 52242114-a86b-43ff-a10b-580805e48b27
@bind h2 overlay_mask_bind(dilated_mask_L_MD)

# ╔═╡ 3e627949-fe0f-4f3a-8f94-9b8b8eb04736
overlay_mask_plot(arr, dilated_mask_L_MD, h2, "dilated mask")

# ╔═╡ 9049f063-bcc7-4ddf-90db-adf813b2e645
overlayed_mask_l_md = create_mask(arr, dilated_mask_L_MD);

# ╔═╡ a6810274-a229-4e81-ba76-86a116c2cec1
agat_l_md, mass_l_md = score(overlayed_mask_l_md, pixel_size, mass_cal_factor, alg)

# ╔═╡ 03d7f139-8283-42fc-b93c-cf52a2f1a171
md"""
## Low Density
"""

# ╔═╡ 04511af5-6cba-4d2e-83c5-423347e9d12a
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;

# ╔═╡ b2a7843a-08e1-4fb7-95c7-4c30c66b2403
md"""
#### Dilated mask
"""

# ╔═╡ 149d69f2-2280-4925-8c73-47499a6cbf01
dilated_mask_L_LD = dilate(dilate(mask_L_LD_3D));

# ╔═╡ 797cab84-8ea9-4ec0-8ce9-13e0de6371ba
@bind i2 overlay_mask_bind(dilated_mask_L_LD)

# ╔═╡ d1cd6b6f-fd66-4255-aa30-eb9a2ba08809
overlay_mask_plot(arr, dilated_mask_L_LD, i2, "dilated mask")

# ╔═╡ 063add05-86dc-44a6-8b7a-0eba4f04da3f
overlayed_mask_l_ld = create_mask(arr, dilated_mask_L_LD);

# ╔═╡ 9c8e037d-f157-4b19-9e91-88ad2a219e0e
agat_l_ld, mass_l_ld = score(overlayed_mask_l_ld, pixel_size, mass_cal_factor, alg)

# ╔═╡ 4633f872-8ef1-42fd-baa8-3b7a32583dd7
md"""
# Score Medium Inserts
"""

# ╔═╡ d443e864-7eaa-4f4a-b2b9-642c8f7bf165
md"""
## High Density
"""

# ╔═╡ 8e96f01a-5cd3-4dde-a926-d81f851900fa
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;

# ╔═╡ 7155e9f1-73a2-4d37-a7f3-92d1676f8058
md"""
#### Dilated mask
"""

# ╔═╡ 91820ed4-6a31-4132-92f3-b138fdebc482
dilated_mask_M_HD = dilate(dilate(dilate(dilate(mask_M_HD_3D))));

# ╔═╡ 5fe3c5f4-ed73-4a35-b71a-dd8829b59da4
@bind j2 overlay_mask_bind(dilated_mask_M_HD)

# ╔═╡ 52c7be1a-2daa-4b7c-89c4-d72463535aa9
overlay_mask_plot(arr, dilated_mask_M_HD, j2, "dilated mask")

# ╔═╡ 74636d4e-aaef-4fa0-8c55-5afc49418397
overlayed_mask_m_hd = create_mask(arr, dilated_mask_M_HD);

# ╔═╡ 3fefa81c-89d4-416e-be01-174a8a543dd6
agat_m_hd, mass_m_hd = score(overlayed_mask_m_hd, pixel_size, mass_cal_factor, alg)

# ╔═╡ f7dda15e-b22d-4a33-9efa-da167abe6538
md"""
## Medium Density
"""

# ╔═╡ 1bea5ec6-a17a-4640-89ea-c6666499b2fc
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;

# ╔═╡ da9714c0-fead-4343-a50c-a8a2c7f1d9f4
md"""
#### Dilated mask
"""

# ╔═╡ d602f6af-0130-4a51-be9e-465627f6730d
dilated_mask_M_MD = dilate(dilate(dilate(dilate(mask_M_MD_3D))));

# ╔═╡ e36c2ef8-260b-418f-b8a8-52ebf731ebe3
@bind k2 overlay_mask_bind(dilated_mask_M_MD)

# ╔═╡ efbd8bf7-ce73-463b-a4ca-820f7213b3ee
overlay_mask_plot(arr, dilated_mask_M_MD, k2, "dilated mask")

# ╔═╡ 5e154bdf-84b1-4c97-bc71-1c632ce8d5bd
overlayed_mask_m_md = create_mask(arr, dilated_mask_M_MD);

# ╔═╡ 7c9ca207-83b0-48e9-982f-f2408a956d9c
agat_m_md, mass_m_md = score(overlayed_mask_m_md, pixel_size, mass_cal_factor, alg)

# ╔═╡ f8e48d57-c674-4915-8c79-589bde4cc677
md"""
## Low Density
"""

# ╔═╡ c6bbd448-0e32-4eeb-b7d4-f9718ec18304
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;

# ╔═╡ 79fea025-1886-4929-8d96-d9f0505c432d
md"""
#### Dilated mask
"""

# ╔═╡ 28ba6334-8adb-4a55-a48d-19f8e4f5f728
dilated_mask_M_LD = dilate(dilate(dilate(dilate(dilate(mask_M_LD_3D)))));

# ╔═╡ 81ed8294-fa65-4164-86a8-1c219ab321ef
@bind l2 overlay_mask_bind(dilated_mask_M_LD)

# ╔═╡ 2c0861f2-7d97-4b75-a889-c66e27b09f13
overlay_mask_plot(arr, dilated_mask_M_LD, l2, "dilated mask")

# ╔═╡ 7bfd57f2-2027-49cb-abab-8ef11409645c
overlayed_mask_m_ld = create_mask(arr, dilated_mask_M_LD);

# ╔═╡ f4f21b51-6b56-4b97-812e-b9b095c22505
agat_m_ld, mass_m_ld = score(overlayed_mask_m_ld, pixel_size, mass_cal_factor, alg)

# ╔═╡ f02905f4-f15f-4876-b26e-47d70b9c488e
md"""
# Score Small Inserts
"""

# ╔═╡ 8df6e736-67af-4fdf-8eeb-18680b703d8a
md"""
## High Density
"""

# ╔═╡ d9695ec4-ad31-47ae-9baa-aa411ac46e1c
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;

# ╔═╡ e574c7b9-09c6-4fc5-820c-78ceaf5a3c80
md"""
#### Dilated mask
"""

# ╔═╡ f0abacd4-eb6b-494d-ac93-342c1e29bd91
dilated_mask_S_HD = dilate(dilate(dilate(dilate(dilate((mask_S_HD_3D))))));

# ╔═╡ 5b43337b-1ad0-4640-b285-1cd0bd032410
@bind m2 overlay_mask_bind(dilated_mask_S_HD)

# ╔═╡ 5072be8f-34ed-46e7-aad3-9315525f7232
overlay_mask_plot(arr, dilated_mask_S_HD, m2, "dilated mask")

# ╔═╡ 9bc5292c-d669-46b4-8827-8ff950e44c6c
overlayed_mask_s_hd = create_mask(arr, dilated_mask_S_HD);

# ╔═╡ e9e619f8-5b94-45e9-80d3-cdf8faa4523a
agat_s_hd, mass_s_hd = score(overlayed_mask_s_hd, pixel_size, mass_cal_factor, alg)

# ╔═╡ 49d9fb49-5a12-4b14-b759-5cb6263eafed
md"""
## Medium Density
"""

# ╔═╡ 551eb59d-79d5-4098-93ff-c624fee56119
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;

# ╔═╡ 31803d25-bc2f-44b3-84b2-88291d2138ae
md"""
#### Dilated mask
"""

# ╔═╡ c91a36be-99bc-495a-9969-954918164375
dilated_mask_S_MD = dilate(dilate(dilate(dilate(dilate(mask_S_MD_3D)))));

# ╔═╡ b4af4774-916b-4712-bdf8-811c289d316c
@bind n2 overlay_mask_bind(dilated_mask_S_MD)

# ╔═╡ e66946b7-fcd4-4cfd-9500-ca0b0ececf26
overlay_mask_plot(arr, dilated_mask_S_MD, n2, "dilated mask")

# ╔═╡ fa8ce4b9-898a-4706-97c1-9398b15779ea
overlayed_mask_s_md = create_mask(arr, dilated_mask_S_MD);

# ╔═╡ 70581d09-b3eb-4c79-bc83-89945b22881e
agat_s_md, mass_s_md = score(overlayed_mask_s_md, pixel_size, mass_cal_factor, alg)

# ╔═╡ 16585130-17c0-46bd-9b35-c55dad77819f
md"""
## Low Density
"""

# ╔═╡ 98b22d21-4a54-4b2b-9ca7-7f32c0442195
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;

# ╔═╡ e45981da-de12-4484-9a07-f6ca3f1c9cf1
md"""
#### Dilated mask
"""

# ╔═╡ 2cff2e76-a209-45d1-8970-bb7b101b91c0
dilated_mask_S_LD = dilate(dilate(dilate(dilate(dilate(mask_S_LD_3D)))));

# ╔═╡ 3495f22c-8883-48ef-81d2-fcc670b654be
@bind o2 overlay_mask_bind(dilated_mask_S_LD)

# ╔═╡ 8856ad47-91ac-4f59-89a5-26e7106f16c7
overlay_mask_plot(arr, dilated_mask_S_LD, o2, "dilated mask")

# ╔═╡ c3b274e2-7bd5-440e-a977-c29e1c20b5ce
overlayed_mask_s_ld = create_mask(arr, dilated_mask_S_LD);

# ╔═╡ cb920bb8-7c28-48d5-af49-8f94caa61f94
# agat_s_ld = score(overlayed_mask_s_ld, pixel_size, alg)

# ╔═╡ 996d788e-9d3e-46a6-9532-59002db192d0
agat_s_ld, mass_s_ld = score(overlayed_mask_s_ld, pixel_size, mass_cal_factor, alg)

# ╔═╡ be6ac424-cd70-4dda-99ee-6c08c131a56b
md"""
# Results
"""

# ╔═╡ 3d021782-5514-4266-bfc5-a42b4ab5ec66
density_array = [0, 200, 400, 800]

# ╔═╡ cdf5731b-091f-47d9-a992-8a3c5ef3f878
inserts = [
	"Low Density",
	"Medium Density",
	"High Density"
]

# ╔═╡ fa458085-3a0c-4d8d-b3d3-9a3cb6f4b36f
md"""
## Agatston
"""

# ╔═╡ 59474596-ba9c-470a-8007-c656623b07fb
calculated_agat_large = [
	agat_l_ld,
	agat_l_md,
	agat_l_hd
]

# ╔═╡ b3909180-04f4-471a-b1e3-457c5431bc19
calculated_agat_medium = [
	agat_m_ld,
	agat_m_md,
	agat_m_hd
]

# ╔═╡ 3dfc2f6e-7ba2-4482-9240-8a6800fa4a2d
calculated_agat_small = [
	agat_s_ld,
	agat_s_md,
	agat_s_hd
]

# ╔═╡ 6652c87e-44dc-4a4b-b291-c5978c97a616
md"""
## Mass
"""

# ╔═╡ 9968b6e4-e852-4766-a30a-43d23be32b7b
ground_truth_mass_large = [
	19.6,
	39.3,
	78.5
] # mg

# ╔═╡ b8bf59f2-24fa-43fe-8260-e85ed3028b93
calculated_mass_large = [
	mass_l_ld,
	mass_l_md,
	mass_l_hd
]

# ╔═╡ a48c296d-1e1b-4ecd-89d9-645577bc877d
ground_truth_mass_medium = [
	4.2,
	8.5,
	17.0
]

# ╔═╡ 84b7fdd1-76cc-4556-bf99-b93b1d93d555
calculated_mass_medium = [
	mass_m_ld,
	mass_m_md,
	mass_m_hd
]

# ╔═╡ da4af805-291d-4a52-b3c8-f1f8d98bc1d8
ground_truth_mass_small = [
	0.2,
	0.3,
	0.6
]

# ╔═╡ 749a5fe6-d063-4728-ab2b-530e694a9f17
calculated_mass_small = [
	mass_s_ld,
	mass_s_md,
	mass_s_hd
]

# ╔═╡ 35175bc9-1fec-4bd2-b9a8-381b69996f4a
df = DataFrame(
	inserts = inserts,
	calculated_agat_large = calculated_agat_large,
	calculated_agat_medium = calculated_agat_medium,
	calculated_agat_small = calculated_agat_small,
	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,
	mass_cal_factor = mass_cal_factor
)

# ╔═╡ bcacf403-5859-4b30-b8f5-0106ca01cd2c
begin
	fmass22 = Figure()
	axmass22 = Axis(fmass22[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")

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

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

	fmass22[1, 2] = Legend(fmass22, axmass22, framevisible = false)

	fmass22
end

# ╔═╡ b83fc531-4867-4056-bd54-af8eeab50496
begin
	fmass32 = Figure()
	axmass32 = Axis(fmass32[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")

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

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

	fmass32[1, 2] = Legend(fmass32, axmass32, framevisible = false)

	fmass32
end

# ╔═╡ 3d20bc0e-9156-41d9-93b3-d044ce50b0ef
begin
	fmass42 = Figure()
	axmass42 = Axis(fmass42[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")

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

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

	fmass42[1, 2] = Legend(fmass42, axmass42, framevisible = false)

	fmass42
end

# ╔═╡ c9e8865a-9143-4a5d-8b0c-aae193168d0d
md"""
### Save Results
"""

# ╔═╡ 2629fd24-e10e-4c42-a388-9f4b680bf550
if ~isdir(string(cd(pwd, "..") , "/data/output/", VENDER, "_agatston"))
	mkdir(string(cd(pwd, "..") , "/data/output/", VENDER, "_agatston",))
end

# ╔═╡ 618cf15a-4ccc-48a3-b122-e346d99d0635
output_path = string(cd(pwd, "..") , "/data/output/", VENDER, "_agatston", "/", scan, ".csv")

# ╔═╡ 62d446bf-6e50-4a6b-aa3c-f8e97e99490f
CSV.write(output_path, df)

# ╔═╡ Cell order:
# ╠═651419e4-4c13-4ddd-8e69-5cdf2b995a51
# ╠═d9401636-eb09-45b2-9663-b8ab615531f9
# ╟─ceb200df-c17b-476f-ae95-a13e572a243e
# ╠═35c2524a-28a4-49eb-9fe2-e4232084af1a
# ╟─7c5f2ec6-26bb-4a2a-9dd1-a7f3b0de10b0
# ╠═7dab92d6-13c3-4c1e-8b51-ad80a709a80b
# ╠═9bcfc64a-170b-461c-baa1-ac5c1ab73579
# ╠═041192c8-875e-43f3-881e-dc993955c91f
# ╠═57cd84bb-9307-44eb-87c1-9a663e10e5c2
# ╠═807cfdd5-379a-4c65-ab15-534f1e28b65c
# ╟─948b9234-f845-48f5-a246-18cbc961ad67
# ╟─baa55458-c8e2-4852-b486-ef8c854f77bf
# ╟─7434ae0d-06f0-44fc-b9e2-198fb9afb97c
# ╟─7325a8d4-6e0a-4087-ad64-d0874fb26aad
# ╟─63c84140-6510-4186-a0bf-dab85d19e2cc
# ╟─a9516b63-80c0-42ff-9bac-a9ad311fbe9c
# ╠═fb4542a8-3406-48bb-b255-ebde03adb951
# ╟─835490d8-332c-400a-a2f7-2dcff791d3d6
# ╟─a2172056-1f1c-40c9-ac6f-816c80dd09ef
# ╟─38042368-0be5-4c64-a4e6-c596201cdd2e
# ╟─a467ddd4-e2d3-445e-808b-298fea9b927c
# ╟─dfb18ee0-7aa7-43c3-bf7e-914b2ee4b729
# ╟─65b0a9dc-efc2-4eb7-9995-1b8bd8113a63
# ╠═8ae2afdc-ae49-49e6-ba59-d20373907ba7
# ╟─d635ed56-ef46-428a-ac7e-f2fc5c10e4d6
# ╠═e4052629-fe1c-4b9f-8d74-fe7dce34de06
# ╠═630ec973-b16d-46b3-bb3b-229d407087b1
# ╟─bdb61513-e02a-43df-be3e-164f415e061b
# ╠═5ec27948-098e-4acf-961a-0bc36307f807
# ╠═480f30c9-b0cb-4f49-8767-717d6213902e
# ╠═02cd946e-98a3-4043-b482-7a19dcd1f157
# ╟─84f78e6d-c73f-4e74-bfb2-1de22990d1c8
# ╠═fe8c8d86-816d-46f8-9068-663c4a45ce97
# ╠═0604f5cf-a731-405e-9c3c-a9f918352f96
# ╠═fd78b7f0-a83e-4365-a770-16483010376a
# ╠═8e1ca1c1-ea27-4f0e-89cd-b466053e9485
# ╠═3c53f529-9fe6-47e9-86f9-41de9ee92b87
# ╠═64d8f848-cc5b-4b51-90bd-aacca098593b
# ╟─11e23023-d88b-4509-a562-51d5255c6b76
# ╠═9c14d354-9f0d-45d5-8dc2-98a719b0465c
# ╠═5ee01db8-2a9d-491d-8305-adff24c6dd09
# ╟─c0f2a223-d737-4823-80c0-76801ae26745
# ╠═c986c5a2-dffb-4200-af4b-fde0e94c5361
# ╟─cc598eba-d013-4984-aa5a-6f94f8f03f19
# ╠═630c79f5-16c3-4ce7-b197-16348783c2a0
# ╟─cf80857d-f58d-46a8-b7fc-f44dfbd8d2a6
# ╠═5e82dcc1-431f-4ba7-9adf-df62beaa0ca9
# ╠═53b8bada-351b-42d1-8807-79ed375b8186
# ╠═dc41660b-6620-4cd6-859b-b5fc3acc7063
# ╠═4cd5d9bd-1b56-4cd4-b46c-168cb3da6b4d
# ╟─c164f2e6-edea-45f6-8cde-4da1ba1820a0
# ╠═c1ba29be-fad6-40f5-a828-122ea002ee74
# ╟─4b2fe0ac-f265-4a29-9ef0-b452e5a66d22
# ╠═4493b19e-bfb7-42da-9d5d-457b30707f14
# ╟─52242114-a86b-43ff-a10b-580805e48b27
# ╠═3e627949-fe0f-4f3a-8f94-9b8b8eb04736
# ╠═9049f063-bcc7-4ddf-90db-adf813b2e645
# ╠═a6810274-a229-4e81-ba76-86a116c2cec1
# ╟─03d7f139-8283-42fc-b93c-cf52a2f1a171
# ╠═04511af5-6cba-4d2e-83c5-423347e9d12a
# ╟─b2a7843a-08e1-4fb7-95c7-4c30c66b2403
# ╠═149d69f2-2280-4925-8c73-47499a6cbf01
# ╟─797cab84-8ea9-4ec0-8ce9-13e0de6371ba
# ╠═d1cd6b6f-fd66-4255-aa30-eb9a2ba08809
# ╠═063add05-86dc-44a6-8b7a-0eba4f04da3f
# ╠═9c8e037d-f157-4b19-9e91-88ad2a219e0e
# ╟─4633f872-8ef1-42fd-baa8-3b7a32583dd7
# ╟─d443e864-7eaa-4f4a-b2b9-642c8f7bf165
# ╠═8e96f01a-5cd3-4dde-a926-d81f851900fa
# ╟─7155e9f1-73a2-4d37-a7f3-92d1676f8058
# ╠═91820ed4-6a31-4132-92f3-b138fdebc482
# ╟─5fe3c5f4-ed73-4a35-b71a-dd8829b59da4
# ╠═52c7be1a-2daa-4b7c-89c4-d72463535aa9
# ╠═74636d4e-aaef-4fa0-8c55-5afc49418397
# ╠═3fefa81c-89d4-416e-be01-174a8a543dd6
# ╟─f7dda15e-b22d-4a33-9efa-da167abe6538
# ╠═1bea5ec6-a17a-4640-89ea-c6666499b2fc
# ╟─da9714c0-fead-4343-a50c-a8a2c7f1d9f4
# ╠═d602f6af-0130-4a51-be9e-465627f6730d
# ╟─e36c2ef8-260b-418f-b8a8-52ebf731ebe3
# ╠═efbd8bf7-ce73-463b-a4ca-820f7213b3ee
# ╠═5e154bdf-84b1-4c97-bc71-1c632ce8d5bd
# ╠═7c9ca207-83b0-48e9-982f-f2408a956d9c
# ╟─f8e48d57-c674-4915-8c79-589bde4cc677
# ╠═c6bbd448-0e32-4eeb-b7d4-f9718ec18304
# ╟─79fea025-1886-4929-8d96-d9f0505c432d
# ╠═28ba6334-8adb-4a55-a48d-19f8e4f5f728
# ╟─81ed8294-fa65-4164-86a8-1c219ab321ef
# ╠═2c0861f2-7d97-4b75-a889-c66e27b09f13
# ╠═7bfd57f2-2027-49cb-abab-8ef11409645c
# ╠═f4f21b51-6b56-4b97-812e-b9b095c22505
# ╟─f02905f4-f15f-4876-b26e-47d70b9c488e
# ╟─8df6e736-67af-4fdf-8eeb-18680b703d8a
# ╠═d9695ec4-ad31-47ae-9baa-aa411ac46e1c
# ╟─e574c7b9-09c6-4fc5-820c-78ceaf5a3c80
# ╠═f0abacd4-eb6b-494d-ac93-342c1e29bd91
# ╟─5b43337b-1ad0-4640-b285-1cd0bd032410
# ╠═5072be8f-34ed-46e7-aad3-9315525f7232
# ╠═9bc5292c-d669-46b4-8827-8ff950e44c6c
# ╠═e9e619f8-5b94-45e9-80d3-cdf8faa4523a
# ╟─49d9fb49-5a12-4b14-b759-5cb6263eafed
# ╠═551eb59d-79d5-4098-93ff-c624fee56119
# ╟─31803d25-bc2f-44b3-84b2-88291d2138ae
# ╠═c91a36be-99bc-495a-9969-954918164375
# ╟─b4af4774-916b-4712-bdf8-811c289d316c
# ╠═e66946b7-fcd4-4cfd-9500-ca0b0ececf26
# ╠═fa8ce4b9-898a-4706-97c1-9398b15779ea
# ╠═70581d09-b3eb-4c79-bc83-89945b22881e
# ╟─16585130-17c0-46bd-9b35-c55dad77819f
# ╠═98b22d21-4a54-4b2b-9ca7-7f32c0442195
# ╟─e45981da-de12-4484-9a07-f6ca3f1c9cf1
# ╠═2cff2e76-a209-45d1-8970-bb7b101b91c0
# ╟─3495f22c-8883-48ef-81d2-fcc670b654be
# ╠═8856ad47-91ac-4f59-89a5-26e7106f16c7
# ╠═c3b274e2-7bd5-440e-a977-c29e1c20b5ce
# ╠═cb920bb8-7c28-48d5-af49-8f94caa61f94
# ╠═996d788e-9d3e-46a6-9532-59002db192d0
# ╟─be6ac424-cd70-4dda-99ee-6c08c131a56b
# ╠═3d021782-5514-4266-bfc5-a42b4ab5ec66
# ╠═cdf5731b-091f-47d9-a992-8a3c5ef3f878
# ╟─fa458085-3a0c-4d8d-b3d3-9a3cb6f4b36f
# ╠═59474596-ba9c-470a-8007-c656623b07fb
# ╠═b3909180-04f4-471a-b1e3-457c5431bc19
# ╠═3dfc2f6e-7ba2-4482-9240-8a6800fa4a2d
# ╠═35175bc9-1fec-4bd2-b9a8-381b69996f4a
# ╟─6652c87e-44dc-4a4b-b291-c5978c97a616
# ╠═9968b6e4-e852-4766-a30a-43d23be32b7b
# ╠═b8bf59f2-24fa-43fe-8260-e85ed3028b93
# ╠═a48c296d-1e1b-4ecd-89d9-645577bc877d
# ╠═84b7fdd1-76cc-4556-bf99-b93b1d93d555
# ╠═da4af805-291d-4a52-b3c8-f1f8d98bc1d8
# ╠═749a5fe6-d063-4728-ab2b-530e694a9f17
# ╟─bcacf403-5859-4b30-b8f5-0106ca01cd2c
# ╟─b83fc531-4867-4056-bd54-af8eeab50496
# ╟─3d20bc0e-9156-41d9-93b3-d044ce50b0ef
# ╟─c9e8865a-9143-4a5d-8b0c-aae193168d0d
# ╠═2629fd24-e10e-4c42-a388-9f4b680bf550
# ╠═618cf15a-4ccc-48a3-b122-e346d99d0635
# ╠═62d446bf-6e50-4a6b-aa3c-f8e97e99490f

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:
using ADCME
using ADSeismic
using PyPlot
using DelimitedFiles
using JLD2

scale = 201

param = AcousticPropagatorParams(NX=scale, NY=scale,
    NSTEP=1000, DELTAT=1e-4,  DELTAX=1.0, DELTAY=1.0,
    PropagatorKernel = 0, Rcoef = 1e-8)

rc = Ricker(param, 30.0, 200.0, 1e6)
srci = [div(param.NX,2)]
srcj = [div(param.NY,5)]
srcv = reshape(rc, :, 1)
src = AcousticSource(srci, srcj, srcv)

layers = ones(param.NX+2, param.NY+2)
n_piece = div(param.NX + 1, 3) + 1
for k = 1:3
    i_interval = (k-1)*n_piece+1:min(k*n_piece, param.NX+2)
    layers[:, i_interval] .= 0.5 + (k-1)*0.25
end

C = placeholder(3300*layers)
model = AcousticPropagatorSolver(param, src, C)

sess = Session(); init(sess)
u = run(sess, model.u)

p = visualize_wavefield(u, param)
saveanim(p, "wavefield.gif")

visualize_model(3300*layers, param)
savefig("model.png")

close("all")
plot(LinRange(0, 1, 1000), run(sess, rc))
xlabel("Time")
savefig("ricker.png")

@save "data.jld2" u


text:
using GenServers
using SafeTestsets, Test

@safetestset "GenServers.jl" begin
    include("../examples/stack.jl")
    using .Stack

    gs = start([1,2,3])
    @test info(gs) == [1,2,3]
    push(gs, 4)
    @test info(gs) == [1,2,3,4]
    @test pop(gs) == 4
    @test info(gs) == [1,2,3]
end


text:
using OctopusArm
using Test

@testset "OctopusArm.jl" begin
    # Write your tests here.
end


text:
# Parameter learning by Stochastic Mini-Batch Expectation-Maximimzation

"""
Learn weights using the Expectation Maximization algorithm.
"""
mutable struct SEM <: ParameterLearner
spn::SumProductNetwork
layers::Vector{Vector{Int}}
previous::SumProductNetwork # for performing temporary updates and applying momentum
diff::Vector{Float64} # to store derivatives
values::Vector{Float64} # to store logprobabilities
# dataset::AbstractMatrix
score::Float64     # score (loglikelihood)
prevscore::Float64 # for checking convergence
tolerance::Float64 # tolerance for convergence criterion
steps::Integer   # number of learning steps (epochs)
minimumvariance::Float64 # minimum variance for Gaussian leaves
SEM(spn::SumProductNetwork) = new(spn, layers(spn), deepcopy(spn), Array{Float64}(undef,length(spn)), Array{Float64}(undef,length(spn)), NaN, NaN, 1e-4, 0, 0.5)
end

"""
Verifies if algorithm has converged.

Convergence is defined as absolute difference of score.
Requires at least 2 steps of optimization.
"""
converged(learner::SEM) = learner.steps < 2 ? false : abs(learner.score - learner.prevscore) < learner.tolerance

# TODO: take filename os CSV File object as input and iterate over file to decreae memory footprint
"""
Improvement step for learning weights using the Expectation Maximization algorithm. Returns improvement in negative loglikelihood.

spn[i].weights[j] = spn[i].weights[k] * backpropagate(spn)[i]/sum(spn[i].weights[k] * backpropagate(spn)[i] for k=1:length(spn[i].children))

## Arguments

- `learner`: EMParamLearner struct
- `data`: Data Matrix
- `learningrate`: learning inertia, used to avoid forgetting in minibatch learning [default: 1.0 correspond to no inertia, i.e., previous weights are discarded after update]
- `smoothing`: weight smoothing factor (= pseudo expected count) [default: 1e-4]
- `minimumvariance`: minimum variance for Gaussian leaves [default: learner.minimumvariance]
"""
function update(learner::SEM, Data::AbstractMatrix, learningrate::Float64 = 1.0, smoothing::Float64 = 1e-4, minimumvariance::Float64 = learner.minimumvariance)

    numrows, numcols = size(Data)

    spn_p = learner.spn      # current parameters
    spn_n = learner.previous # updated parameters

    # # @assert numcols == spn._numvars "Number of columns should match number of variables in network."
    # m, n = size(spn._weights)
    # weights = nonzeros(spn._weights)
    # childrens = rowvals(spn._weights)
    # # oldweights = fill(τ, length(weights))
    # newweights = fill(τ, length(weights))
    # maxweights = similar(newweights)
    score = 0.0 # data loglikelihood
    sumnodes = filter(i -> isa(spn_p[i], SumNode), 1:length(spn_p))
    # gaussiannodes = filter(i -> isa(spn[i],GaussianDistribution), 1:length(spn))
    # if length(gaussiannodes) > 0
    #     means = Dict{Integer,Float64}(i => 0.0 for i in gaussiannodes)
    #     squares = Dict{Integer,Float64}(i => 0.0 for i in gaussiannodes)
    #     denon = Dict{Integer,Float64}(i => 0.0 for i in gaussiannodes)
    # end
    #diff = zeros(Float64, length(spn))
    diff = learner.diff
    values = learner.values
    #values = similar(diff)
    # TODO beter exploit multithreading
    # Compute expected weights
    for t = 1:numrows
        datum = view(Data, t, :)
        #lv = logpdf!(values,spn,datum) # propagate input Data[i,:]
        lv = plogpdf!(values,spn_p,learner.layers,datum) # parallelized version
        @assert isfinite(lv) "logvalue of datum $t is not finite: $lv"
        score += lv
        #TODO: implement multithreaded version of backpropagate
        backpropagate!(diff, spn_p, values) # backpropagate derivatives
        Threads.@threads for i in sumnodes # update each node in parallel
            @inbounds for (k,j) in enumerate(spn_p[i].children)
                # @assert isfinite(diff[i]) "derivative of node $i is not finite: $(diff[i])"
                # @assert !isnan(values[j]) "value of node $j is NaN: $(values[j])"
                if isfinite(values[j])
                    δ = spn_p[i].weights[k]*diff[i]*exp(values[j]-lv) # improvement
                    @assert isfinite(δ) "improvement to weight ($i,$j):$(spn_p[i].weights[k]) is not finite: $δ, $(diff[i]), $(values[j]), $(exp(values[j]-lv))"
                else
                    δ = 0.0
                end
                spn_n[i].weights[k] = ((t-1)/t)*spn_n[i].weights[k] + δ/t # running average for improved precision
            end
    #         # for j in children(spn,i) #spn._backward[i]
    #         #     oldweights[j,i] += spn._weights[j,i]*diff[i]*exp(values[j]-lv)
    #         # end
    #         for k in nzrange(spn._weights,i)
    #             j = childrens[k]
    #             #oldweights[k] += weights[k]*diff[i]*exp(values[j]-lv)
    #             Δ = log(weights[k]) + log(diff[i]) + values[j] - lv
    #             if t == 1
    #                 maxweights[k] = Δ
    #                 newweights[k] = 1.0
    #             else
    #                 if Δ > maxweights[k]
    #                     newweights[k] = exp(log(newweights[k])+maxweights[k]-Δ)+1.0
    #                     maxweights[k] = Δ
    #                 elseif isfinite(Δ) && isfinite(maxweights[k])
    #                     newweights[k] += exp(Δ-maxweights[k])
    #                 end
    #                 @assert isfinite(newweights[k]) "Infinite weight: $(newweights[k])"
    #             end
    #         end
        end
    #     for i in gaussiannodes
    #         α = diff[i]*exp(values[i]-lv)
    #         denon[i] += α
    #         means[i] += α*datum[spn[i].scope]
    #         squares[i] += α*datum[spn[i].scope]^2
    #     end
    end
    # # add regularizer to avoid degenerate distributions
    # Do update
    # TODO: implement momentum acceleration (nesterov acceleration, Adam, etc)
    # newweights =  log.(newweights) .+ maxweights
    @inbounds Threads.@threads for i in sumnodes
        spn_n[i].weights .+= smoothing/length(spn_n[i].weights) # smoothing factor to prevent degenerate probabilities
        spn_n[i].weights .*= learningrate/sum(spn_n[i].weights) # normalize weights
        # online update: θ[t+1] = (1-η)*θ[t] + η*update(θ[t])
        spn_n[i].weights .+= (1.0-learningrate)*spn_p[i].weights
        # cache[i].weights .*= learningrate/sum(cache[i].weights) # normalize weights
        # spn[i].weights .*= 1.0-learningrate # apply update with inertia strenght given by learning rate
        # @assert sum(spn[i].weights) ≈ 1.0-learningrate "Unnormalized weight vector at node $i: $(sum(spn[i].weights)) | $(spn[i].weights)"
        # spn[i].weights .+= cache[i].weights
        @assert sum(spn_n[i].weights) ≈ 1.0 "Unnormalized weight vector at node $i: $(sum(spn_n[i].weights)) | $(spn_n[i].weights) | $(spn_p[i].weights)"
        # for (k,j) in enumerate(spn[i].children)
        #     # spn[i].weights .*= cache[i].weights
            # Z = sum(spn[i].weights)
            # spn[i].weights ./= Z
            # normexp!(spn[i].weights, prev[i].weights, τ) # if weights are in log
        # end
        #     chval = nzrange(spn._weights,i)
    #     normexp!(view(newweights,chval), view(weights,chval), τ)
    #     if !isfinite(sum(weights[chval])) # some numerical problem occurred, set weights to uniform
    #         weights[chval] .= 1/length(chval)
    #     end
    #     # @assert sum(weights[chval]) ≈ 1.0 "Unormalized weight vector: $(sum(weights[chval])) | $(weights[chval])"
    #     # add regularizer to avoid degenerate distributions
    #     # Z = sum(oldweights[:,i]) + τ*length(children(spn,i))
    #     # for j in  children(spn,i) #spn._backward[i]
    #     #    spn._weights[j,i] = (oldweights[j,i]+τ)/Z
    #     #    # if isnan(Z)
    #     #    #    @warn "Not a number for weight $i -> $j"
    #     #    # end
    #     # end
    end
    #println(spn)

    # for i in gaussiannodes
    #     spn[i].mean = means[i]/denon[i]
    #     spn[i].variance = squares[i]/denon[i] - (spn[i].mean)^2
    #     if spn[i].variance < minimumvariance
    #         spn[i].variance = minimumvariance
    #     end
    # end
    learner.previous = spn_p
    learner.spn = spn_n
    learner.steps += 1
    learner.prevscore = learner.score
    learner.score = -score/numrows
    return learner.prevscore-learner.score
end


text:
using PlutoSliderServer
using FromFile
using Test

toml_path = tempname()

write(toml_path, PlutoSliderServer.sample_config_toml_file)

result = PlutoSliderServer.get_configuration(toml_path; Export_binder_url="yayy")

@test result.Export.binder_url == "yayy"
@test result.Pluto.compiler.threads == 1


text:
# developed with Julia 1.1.1
#
# EMS simulation with a dummy controller

using EMSx

include("arguments.jl")

args = parse_commandline()

const controller = EMSx.DummyController()

EMSx.simulate_sites(controller,
    joinpath(args["save"], "dummy"),
    args["price"],
    args["metadata"],
    args["test"])


text:
# Test KNITRO with MathProgBase unit-tests
#
using MathProgBase, KNITRO

include(joinpath(dirname(pathof(MathProgBase)),"..", "test", "nlp.jl"))
convexnlptest(KnitroSolver())
rosenbrocktest(KnitroSolver(opttol=1e-7))
nlptest(KnitroSolver())
nlptest_nohessian(KnitroSolver())


text:
"""
Differentiation
"""
function Δ(X,k,dt)
    return (X[:,k] - X[:, k-1])/dt
end

"""
Linspace
"""
function linspace(a,b,n)
    return [collect(range(a, stop=b, length=n))]
end

function linearpath(point_start, point_end, TotalPoint)

    initpath = point_start*ones(1,TotalPoint )

    for i in 1:size(initpath)[2]
        step = (point_end - point_start)*i/size(initpath)[2]
        initpath[:,i] = initpath[:,i] + step
    end
    return initpath
end


text:
# --- Utilities

valueof(::Val{x}) where x = x

ieltype(T) =
    if Base.IteratorEltype(T) isa Base.HasEltype
        eltype(T)
    else
        Any
    end

ieltype(bc::Broadcasted) = Base.promote_op(bc.f, eltype.(bc.args)...)

avaltype(x) = avaltype(typeof(x))
avaltype(T::Type) = valtype(T)
avaltype(T::Type{<:Union{AbstractArray,AbstractString}}) = eltype(T)

prefixed_type_name(@nospecialize x) =
    sprint(show, typeof(x), context = :module => Base)
# `:module => Base` to enforce that the type is prefixed even when
# `typeof(x)` is imported in `Main`.  See:
# https://github.com/JuliaLang/julia/pull/29466

const DenseSubVector{T} =
    SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int64}}, true}

const _non_executable_transducer_msg = """
Output type of the transducer is inferred to be a `Union{}`.  This
probably means one or more of the composed transducers throw.
"""

@inline _poptail(xs) = _poptail_impl(xs...)
@inline _poptail_impl(a) = (), a
@inline function _poptail_impl(a, xs...)
    head, tail = _poptail_impl(xs...)
    return (a, head...), tail
end

Base.@pure nthtype(::Val{1}, ::Type{T}) where {T <: Tuple} =
    Base.tuple_type_head(T)

Base.@pure nthtype(::Val{n}, ::Type{T}) where {n, T <: Tuple} =
    nthtype(Val(n - 1), Base.tuple_type_tail(T))

_cljapiurl(name) =
    "https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/$name"
_cljref(name) =
    "[`$name` in Clojure]($(_cljapiurl(name)))"
_thx_clj(name) =
    "This API is modeled after $(_cljref(name))."

# https://github.com/JuliaLang/julia/pull/30575
const _true_str = sprint(show, true; context=:typeinfo => Bool)
const _false_str = sprint(show, false; context=:typeinfo => Bool)


text:
### Abstract parent types for Adaptive samplers
abstract AbstractAdaptiveSampler <: AbstractSampler
abstract AbstractAdaptiveSamplerState <: AbstractSamplerState{BaseSample}

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

### Type holding the parameters for an adaptive Metropolis sampler (i.e with a normal density the covariance of which is adapted to the landscape)
immutable AdaptiveNormalSampler{T<:AbstractFloat} <: AbstractAdaptiveSampler
    initialscalefactor::T
    nparas::Int
    AdaptiveNormalSampler(s::T,n::Int) = new(s,n)
end

###Factory functions
@inline _sampler(::Type{Val{:adaptive}},s::AbstractFloat,n::Integer) = AdaptiveNormalSampler{typeof(s)}(s,n)

###Size functions
@inline numparas(sampler::AdaptiveNormalSampler) = sampler.nparas

###Name functions (used in println for generic AbstractSampler)
@inline samplername(sampler::AdaptiveNormalSampler) = "Adaptive Normal Metropolis"

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

### Type holding the running mean and covariance for the adaptive sampler with normal density
type NormalRunningState{T<:AbstractFloat}
    iteration::Int
    runningmean::Vector{T}
    runningcov::Matrix{T}
    scale::T
    tempdiff::Vector{T}
    densitycov::Matrix{T}
    NormalRunningState(i::Int,m::Vector{T},c::Matrix{T},s::T) = new(i,m,c,s,similar(m),similar(c))
end

### Factory function
@inline function _runningstate{T<:AbstractFloat}(::Type{Val{:normal}},n::Int,s::T)
    NormalRunningState{T}(0,zeros(T,n),zeros(T,n,n),s)
end

@inline _covupdatefactor{T<:AbstractFloat}(i::Int,::Type{T}) = i==1?zero(T):T((i-2)/(i-1))

@inline function _updatetempdiff!{T<:AbstractFloat}(d::Vector{T},m::Vector{T},v::Vector{T})
    @simd for i=1:length(v)
        @inbounds d[i] = v[i] - m[i]
    end
    d
end

@inline function _updaterunningcov!{T<:AbstractFloat}(c::Matrix{T},d::Vector{T},ci::Int)
    f = _covupdatefactor(ci,T)
    nrows = length(d)
    for j=1:size(c,2)
        @simd for i=1:nrows
            @inbounds c[i,j] = f*c[i,j] + d[i]*d[j]/ci
        end
    end
    c
end

@inline function _updaterunningmean!{T<:AbstractFloat}(m::Vector{T},v::Vector{T},ci::Int)
    f::T = T(ci-1)
    @simd for i=1:length(m)
        @inbounds m[i] = (f*m[i] + v[i])/ci
    end
    m
end

@inline function _updatedensitycov!{T<:AbstractFloat}(d::Matrix{T},c::Matrix{T},s::T)
    f = s*s
    copy!(d,c)
    @simd for i=1:size(d,1)
        @inbounds d[i,i] += f #sic, this comes from Ben's code - check with him
    end
    d
end

### Update running mean and covariance
### The formula for running covariance can be derived from the recurrence formulas
### for running covariance, derived by Knuth (see www.johndcook.com/blog/standard_deviation)
@inline function updaterunning!(r::NormalRunningState,v::Vector)
    r.iteration += 1
    _updatetempdiff!(r.tempdiff,r.runningmean,v)
    _updaterunningcov!(r.runningcov,r.tempdiff,r.iteration)
    _updaterunningmean!(r.runningmean,v,r.iteration)
    _updatedensitycov!(r.densitycov,r.runningcov,r.scale)
    r
end

function show(io::IO,r::NormalRunningState)
    println(io,"NormalRunningState with fields:")
    println(io,"  iteration: ",r.iteration)
    println(io,"  runningmean: ",r.runningmean)
    println(io,"  runningcov: ",r.runningcov)
    println(io,"  scale: ",r.scale)
    nothing
end

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

### Type holding the state of the Markov Chain for an Adaptive Metropolis Sampler (with a normal density)
type AdaptiveNormalSamplerState{T<:AbstractFloat} <: AbstractAdaptiveSamplerState
    density::NormalDensity
    from::BaseSample
    proposals::BaseSample
    acceptance::Vector
    #element to keep track of running mean and covariance
    #this is only needed for the main state, not the auxiliary states, hence wrapping it into a "nullable"
    runningstate::Nullable{NormalRunningState{T}}
    AdaptiveNormalSamplerState(n::NormalDensity,f::BaseSample,p::BaseSample,a::Vector,r::Nullable{NormalRunningState{T}}) = new(n,f,p,a,r)
end

###Factory functions
@inline function _samplerstate{N<:Number,T<:AbstractFloat}(s::AdaptiveNormalSampler,nsamples::Integer,::Type{N},::Type{T},auxiliary::Bool)
    nparas = numparas(s)
    f = _samples(Val{:base},nparas,1,N,T)
    p = _samples(Val{:base},nparas,nsamples,N,T)
    n = _density(Val{:normal},s.initialscalefactor*s.initialscalefactor*eye(N,nparas))
    r = auxiliary?Nullable{NormalRunningState{N}}():Nullable(_runningstate(Val{:normal},nparas,s.initialscalefactor))
    AdaptiveNormalSamplerState{N}(n,f,p,zeros(T,nsamples),r)
end

###Name functions (used in println for AbstractSamplerState)
@inline samplerstatename(state::AdaptiveNormalSamplerState) = "Adaptive Normal Metropolis"

@inline function updaterunning!(state::AdaptiveNormalSamplerState)
    if !isnull(state.runningstate)
        r = get(state.runningstate)
        updaterunning!(r,state.from.values)
        state.density = update!(state.density,r.densitycov)
    end
    state
end

###update the sampler state, called from SMHRunner
@inline function prepare!(state::AdaptiveNormalSamplerState,updatefrom::Bool=false)
    if updatefrom
        #set the current point around which we are sampling
        setfrom!(state,state.proposals)
    end
    updaterunning!(state)
end

###update the auxiliary samplerstate with information from the indicator samplerstate
@inline function prepareauxiliary!(indicatorstate::AdaptiveNormalSamplerState,auxiliarystate::AdaptiveNormalSamplerState)
    #set the auxiliary point
    setfrom!(auxiliarystate,indicatorstate.proposals)
    #set the proposal density of the auxiliary state to that of the indicatorstate
    auxiliarystate.density = indicatorstate.density
    #return the auxiliary state
    auxiliarystate
end

@inline prepareindicator!(indicatorstate::AdaptiveNormalSamplerState) = updaterunning!(indicatorstate)

###update the indicator state when starting an iteration of GMHRunner
@inline function prepareindicator!(indicatorstate::AdaptiveNormalSamplerState,auxiliarystate::AdaptiveNormalSamplerState,i::Int)
    #set the current point around which we are sampling
    setfrom!(indicatorstate,auxiliarystate.proposals,i)
    #update the running mean and coveriance
    updaterunning!(indicatorstate)
end

###propose new samples from the current point
@inline function propose!(state::AdaptiveNormalSamplerState)
    #propose n values using the covariance around mean 0.0
    propose!(state.density,state.proposals.values)
    #add the current from to the proposals
    offset!(state.proposals,state.from.values)
    #return the proposals
    state
end

###for a symmetric proposal density, no adjustment of the acceptance rate is needed
function acceptance!(state::AdaptiveNormalSamplerState)
    fl = state.from.loglikelihood[1] + state.from.logprior[1]
    for i=1:numsamples(state.proposals)
        @inbounds state.acceptance[i] = state.proposals.loglikelihood[i] + state.proposals.logprior[i] - fl
    end
    state.acceptance
end

@inline function tune!(state::AdaptiveNormalSamplerState,scalefactor::AbstractFloat)
    if !isnull(state.runningstate) #only the main state needs updating
        get(state.runningstate).scale *= scalefactor
    end
    state
end

function getsamplerstatevars(state::AdaptiveNormalSamplerState)
    Dict("density"=>state.density,"runningstate"=>state.runningstate)
end


text:
using Pluto
server_options = Pluto.Configuration.ServerOptions(; host="0.0.0.0", port=1234, launch_browser=false)
security_options = Pluto.Configuration.SecurityOptions(; require_secret_for_access=false)
Pluto.run(Pluto.Configuration.Options(; server=server_options, security=security_options))


text:
function update_chainrules_primitives!(;force=false)
    @info "update_chainrules_primitives!() is deprecated, you can safely remove this call"
end


text:
# Unstable Bickley jet in GeophysicalFlows

ENV["GKSwstype"] = "nul"
using Printf
using Random
using Plots

using FFTW: rfft, irfft
using CUDA
using JLD2

using FourierFlows
using GeophysicalFlows: TwoDNavierStokes

include("Bickley.jl")

using .Bickley

function run(; dev = CPU(), nx = 128, dt = 1e-2)

    nsubs = ceil(Int, 2 / dt) # number of steps between each plot
    nsteps = 100 * nsubs # total number of steps

    problem = TwoDNavierStokes.Problem(dev; nx = nx,
                                            Lx = 4π,
                                            ny = nx,
                                            Ly = 4π,
                                            dt = dt,
                                            stepper = "FilteredRK4")

    ϵ = 0.1
    k = 0.5
    ℓ = 0.5

    DeviceArray = dev isa CPU ? Array : CuArray

    grid = problem.grid
    x = repeat(reshape(grid.x, grid.nx, 1), 1, grid.ny)
    y = repeat(reshape(grid.y, 1, grid.ny), grid.nx, 1)

    ψᵢ = DeviceArray(@. Bickley.Ψ(y, grid.Ly) + ϵ * Bickley.ψ̃(x, y, k, ℓ))

    ψᵢh = rfft(ψᵢ)
    CUDA.@allowscalar ψᵢh[1, 1] = 0

    ψᵢ = irfft(ψᵢh, grid.nx)

    @. problem.vars.zetah = - grid.Krsq * ψᵢh

    ζᵢ = irfft(problem.vars.zetah, problem.grid.nx)

    TwoDNavierStokes.set_zeta!(problem, ζᵢ)
    TwoDNavierStokes.set_c!(problem, DeviceArray(sin.(y/2)))

    # ## Output

    # We choose folder for outputing `.jld2` files and snapshots (`.png` files).
    filename = "geophysical_flows_unstable_bickley_jet_Nh$nx.jld2"

    isfile(filename) && rm(filename, force=true)

    output = Output(problem, filename,
        (:ζh, p -> Array(p.vars.zetah)),
        (:ζ,  p -> Array(irfft(p.vars.zetah, p.grid.nx))),
        (:ch, p -> Array(p.vars.ch)),
        (:c,  p -> Array(irfft(p.vars.ch, p.grid.nx))),
    )

    saveproblem(output)

    start_time = time_ns()

    e₀ = TwoDNavierStokes.energy(problem)
    χ₀ = TwoDNavierStokes.tracer_variance(problem)

    for j = 0:Int(nsteps/nsubs)

        stepforward!(problem, nsubs)
        TwoDNavierStokes.updatevars!(problem)
        saveoutput(output)

        e = TwoDNavierStokes.energy(problem)
        χ = TwoDNavierStokes.tracer_variance(problem)

        log = @sprintf("Step: %04d, t: %d, Δe/e₀: %.2e, Δχ/χ₀: %.2e, wall_time: %.2f seconds",
                       problem.clock.step,
                       problem.clock.t,
                       e / e₀ - 1,
                       χ / χ₀ - 1,
                       (time_ns() - start_time) * 1e-9,
                      )
        println(log)

    end

    run_time = (time_ns() - start_time) * 1e-9

    return output.path, run_time, nsteps
end

function visualize(filename)

    file = jldopen(filename)

    iterations = parse.(Int, keys(file["snapshots/t"]))

    x = file["grid/x"]
    y = file["grid/y"]

    animation = @animate for (i, iter) in enumerate(iterations)

        @info "Plotting frame $i of $(length(iterations))"

        t = file["snapshots/t/$iter"]
        ζ = file["snapshots/ζ/$iter"]
        c = file["snapshots/c/$iter"]

        kwargs = Dict(
                      :aspectratio => 1,
                      :linewidth => 0,
                      :ticks => nothing,
                      :colorbar => :none,
                      :clims => (-1, 1),
                      :xlims => (-2π, 2π),
                      :ylims => (-2π, 2π),
                     )

        ζ_plot = heatmap(x, y, clamp.(ζ, -1, 1)'; color=:balance, kwargs...)
        c_plot = heatmap(x, y, clamp.(c, -1, 1)'; color=:thermal, kwargs...)

        ω_title = @sprintf("ζ at t = %.1f", t)
        c_title = @sprintf("c at t = %.1f", t)

        plot(ζ_plot, c_plot, title = [ω_title c_title], size=(4000, 2000))
    end

    close(file)

    name = filename[1:end-5]

    gif(animation, name * ".gif", fps=8)

    return nothing
end

cfl = 0.1
for nx in (32,
           64, 128, 256, 512,
           1024, 2048)
           #32, 64, 128, 256, 512)

    @show dt = cfl * 4π / nx
    filename, run_time, nsteps = run(dev=GPU(), nx=nx, dt=dt)

    DOF = nx^2
    @info "Results for nx = $nx !"
    @show cost = run_time / (nsteps * DOF)

    filename = "geophysical_flows_unstable_bickley_jet_Nh$nx.jld2"
    visualize(filename)
end


text:
module FIX

using DataStructures
using DandelionWebSockets

global const TAGS_INT_STRING = Dict{Int64, String}()
global const TAGS_STRING_INT = Dict{String, Int64}()

import Base: start, next, done, length, collect, close

abstract type AbstractMessageHandler <: DandelionWebSockets.WebSocketHandler end
export AbstractMessageHandler, FIXClient, send_message, start, close
export onFIXMessage

function onFIXMessage(this::AbstractMessageHandler, x::Any)
    T = typeof(this)
    X = typeof(x)
    throw(ErrorException("Method `onFIXMessage` is not implemented by $T for argument type $X"))
end

function __init__()
    global TAGS_INT_STRING
    global TAGS_STRING_INT

    fid = open(joinpath(@__DIR__, "../etc/tags.csv"), "r");
    line_number = 0
    while !eof(fid)
        line_number += 1
        line = readline(fid);
        data = split(line, ",");
        if length(data) != 2
            close(fid)
            throw(ErrorException("Invalid data in 'etc/tags.csv' file on line $line_number: $line"))
        end
        tag = parse(Int64, String(data[1]))
        val = String(data[2])

        TAGS_INT_STRING[tag] = val
        TAGS_STRING_INT[val] = tag
    end
    close(fid)

    return nothing
end

mutable struct FIXClientTasks
    read::Nullable{Task}
    function FIXClientTasks()
        return new(Nullable{Task}())
    end
end

mutable struct Container{T}
    data::T
end

include("parse.jl")
include("management.jl")

struct FIXClient{T <: IO, H <: AbstractMessageHandler}
    stream::T
    handler::H
    delimiter::Char
    m_head::Dict{Int64, String}
    m_tasks::FIXClientTasks
    m_messages::FIXMessageManagement
    m_lock::ReentrantLock
    #optimization
    m_intmap::Dict{Int64, String}
    function FIXClient(stream::T,
                        handler::H,
                        header::Dict{Int64, String},
                        ratelimit::RateLimit;
                        delimiter::Char = Char(1)) where {T <: IO, H <: AbstractMessageHandler}
        m_intmap = Dict{Int64, String}()
        [m_intmap[id] = string(id) for id = 1 : 9999]
        return new{T, H}(stream,
                        handler,
                        delimiter,
                        header,
                        FIXClientTasks(),
                        FIXMessageManagement(ratelimit),
                        ReentrantLock(),
                        m_intmap)
    end
end

checksum(this::String)::Int64 = sum([Int(x) for x in this]) % 256
fixjoin(this::OrderedDict{Int64, String}, delimiter::Char)::String = join([string(k) * "=" * v for (k, v) in this], delimiter) * delimiter

function fixjoin(this::OrderedDict{Int64, String}, delimiter::Char, IntMap::Dict{Int64, String})::String
    join([IntMap[k] * "=" * v for (k, v) in this], delimiter) * delimiter
end

function fixmessage(this::FIXClient, msg::Dict{Int64, String})::OrderedDict{Int64, String}
    ordered = OrderedDict{Int64, String}()

    #header
    ordered[8] = this.m_head[8]
    ordered[9] = ""
    ordered[35] = msg[35] #message type
    ordered[49] = this.m_head[49] #SenderCompID
    ordered[56] = this.m_head[56] #TargetCompID
    ordered[34] = getNextOutgoingMsgSeqNum(this)
    ordered[52] = ""
    #body
    body_length = 0
    for (k, v) in msg
        if k != 8 && k != 9 && k != 10
            ordered[k] = v
        end
    end

    for (k, v) in ordered
        if k != 8 && k != 9
            body_length += length(string(k)) + 1 + length(v) + 1 #tag=value|
        end
    end
    ordered[9] = string(body_length)

    #tail
    msg = fixjoin(ordered, this.delimiter, this.m_intmap)
    c = checksum(msg)
    c_str = string(c)
    while length(c_str) < 3
        c_str = '0' * c_str
    end
    ordered[10] = c_str

    #done
    return ordered
end

function send_message_fake(this::FIXClient, msg::Dict{Int64, String})
    lock(this.m_lock)

    msg = fixmessage(this, msg)
    msg_str = fixjoin(msg, this.delimiter, this.m_intmap)
    # write(this.stream, msg_str)
    onSent(this.m_messages, msg)

    unlock(this.m_lock)

    return (msg, msg_str)
end

function send_message(this::FIXClient, msg::Dict{Int64, String})
    lock(this.m_lock)

    msg = fixmessage(this, msg)
    msg_str = fixjoin(msg, this.delimiter, this.m_intmap)
    write(this.stream, msg_str)
    onSent(this.m_messages, msg)

    unlock(this.m_lock)

    return (msg, msg_str)
end

function start(this::FIXClient)
    this.m_tasks.read = Nullable( @async begin
        while true
            incoming = readavailable(this.stream)
            if isempty(incoming)
                @printf("[%ls] EMPTY FIX MESSAGE\n", now())
                break
            end

            for (_, msg) in fixparse(incoming)
                onGet(this, msg)
                onFIXMessage(this.handler, msg)
            end
        end
        @printf("[%ls] FIX: read task done\n", now())
    end
    )
    return this.m_tasks
end

function close(this::FIXClient)
    close(this.stream)
end

function onGet(this::FIXClient, msg::DICTMSG)
    onGet(this.m_messages, msg)
end

function getOpenOrders(this::FIXClient)::Vector{Dict{Int64, String}}
    return getOpenOrders(this.m_messages)
end

function getNextOutgoingMsgSeqNum(this::FIXClient)
    return getNextOutgoingMsgSeqNum(this.m_messages)
end

function getPosition(this::FIXClient, instrument::String)
    return getPosition(this.m_messages, instrument)
end

function getPositions(this::FIXClient)
    return getPositions(this.m_messages)
end

function numMsgsLeftToSend(this::FIXClient)
    return numleft(this.m_messages.outgoing.ratelimit, now())
end

end


text:
module Primeau_2x2x2
#=
This circulation and grid was created by François Primeau (fprimeau@uci.edu)
It was taken from his Notebook at
https://github.com/fprimeau/BIOGEOCHEM_TEACHING/blob/master/Intro2TransportOperators.ipynb

The simple box model we consider is embeded in a 2×2×2 "shoebox".
It has 5 wet boxes and 3 dry boxes.
See image at
https://github.com/fprimeau/BIOGEOCHEM_TEACHING/blob/master/boxmodel.png
=#

using LinearAlgebra, SparseArrays
using Unitful
using Reexport
@reexport using OceanGrids            # To store the grid
using ..CirculationGeneration
CG = CirculationGeneration

function build_wet3D()
    wet3D = trues(2, 2, 2)
    wet3D[[4,7,8]] .= false # land points
    return wet3D
end

function build_grid()
    elon = [0,180,360] * u"°"
    elat = [-90, 0, 90] * u"°"
    edepth = [0, 200, 3700] * u"m"
    return OceanGrid(elon, elat, edepth, build_wet3D())
end

function build_T(grid)
    # From Archer et al. [2000]
    v3D = array_of_volumes(grid)
    nb = length(v3D)

    # Antarctic Circumpoloar Current 1 -> 3 -> 1
    ACC = 100e6u"m^3/s"
    T  = CG.T_advection(ACC, [1, 3, 1], v3D, nb)
    # Meridional Overturning Circulation 1 -> 2 -> 6 -> 5 -> 1
    MOC = 15e6u"m^3/s"
    T += CG.T_advection(MOC, [1, 2, 6, 5, 1], v3D, nb)
    # vertical mixing at "high northern latitudes" 2 <-> 6
    MIX = 10e6u"m^3/s"
    T += CG.T_diffusion(MIX, 2, 6, v3D, nb)

    # Only keep wet points
    iwet = findall(vec(grid.wet3D))

    return T[iwet, iwet]
end

"""
    load

Returns `grd` and `T` (in that order).
"""
function load()
    print("Creating François Primeau's 2x2x2 model")
    grd = build_grid()
    T = build_T(grd)
    println(" ✔")
    return grd, ustrip.(T)
end

end

export Primeau_2x2x2


text:
using Cxx
cxx"""
#include <iostream>
"""
function foo()
    for i = 1:10
        icxx"""
            std::cout << "foo\n";
        """
    end
end
foo()

function bar()
    for i = 1:10
        icxx"""
           std::cout << $i << "\n";
        """
   end
end
bar()

function baz()
    for i = 1:10
        x = icxx" $i + 1;"
        @assert x == i + 1
    end
end
baz()

function foobar()
    for i = 1:10
        x = icxx"""
            if ($i > 5)
                return 2;
            else
                return 1;
        """
        @assert x == (i > 5 ? 2 : 1)
    end
end
foobar()


text:
# reminder, to type H*, do H\^+
#cd("/Users/Frank/.julia/dev/Bridge")
outdir="output/bolus/"

using Bridge, StaticArrays, Distributions
using Bridge:logpdfnormal
using Test, Statistics, Random, LinearAlgebra
using Bridge.Models
using DelimitedFiles
using DataFrames
using CSV
using RCall

sk = 0 # skipped in evaluating loglikelihoodν

νHparam = true
simlongpath = true
obs_scheme =["full","firstcomponent"][2]

# settings in case of νH - parametrisation
ϵ = 10^(-3)
Σdiagel = 10^(-4)

# settings sampler
iterations = 20000
skip_it = 400# 1000
subsamples = 0:skip_it:iterations

ρ = 0.0#95

L = @SMatrix [.5 .5 ]
#L = @SMatrix [1. 0. ; 0.0 1.0]

m, d = size(L)
Σ = SMatrix{m,m}(Σdiagel*I)
dp = 2 # dprime

################## specify target process
struct Diffusion <: ContinuousTimeProcess{ℝ{2}}
    α::Float64
    β::Float64
    λ::Float64
    μ::Float64
    σ1::Float64
    σ2::Float64
end

# pdfemg(x, μ, σ, λ) = λ/2*exp(λ/2*(2μ + λ*σ^2 - 2x)).*erfc((μ + λ*σ^2 - x)/(sqrt(2)*σ))
# dose(t, c) = pdfemg(t, c...)*c[3]
# dose(t, c) = 1. *(t > c)

Bridge.b(t, x, P::Diffusion) = ℝ{2}(P.α*dose(t) -(P.λ + P.β)*x[1] + P.μ*x[2],  P.λ*x[1] -P.μ*x[2])  # reminder mu = k-lambda
Bridge.σ(t, x, P::Diffusion) = @SMatrix [P.σ1 0.0 ;0.0  P.σ1]
Bridge.constdiff(::Diffusion) = true

################## specify auxiliary process
struct DiffusionAux <: ContinuousTimeProcess{ℝ{2}}
    α::Float64
    β::Float64
    λ::Float64
    μ::Float64
    σ1::Float64
    σ2::Float64
end

#Random.seed!(42)
Bridge.B(t, P::DiffusionAux) = @SMatrix [ -P.λ - P.β P.μ ;  P.λ  -P.μ]
Bridge.β(t, P::DiffusionAux) = ℝ{2}(P.α*dose(t),0.0)
Bridge.σ(t, P::DiffusionAux) = @SMatrix [P.σ1 0.0 ;0.0   P.σ2]
Bridge.constdiff(::DiffusionAux) = true
Bridge.b(t, x, P::DiffusionAux) = Bridge.B(t,P) * x + Bridge.β(t,P)
Bridge.a(t, P::DiffusionAux) = Bridge.outer(Bridge.σ(t, P))
DiffusionAux(P::Diffusion) = DiffusionAux(P.α, P.β,P.λ,P.μ,P.σ1,P.σ2)

dose(t) = 2*(t/2)/(1+(t/2)^2)
FT = 70.; VB = 20.; PS = VE = 15.; HE = 0.4; DT = 2.4 # Favetto-Samson
DT = 1. # out choice
#Ptrue = Diffusion(FT/(1-HE), FT/(VB*(1-HE)), PS/(VB*(1-HE)),  PS/VE, sqrt(2),0.2)
Ptrue = Diffusion(FT/(1-HE), FT/(VB*(1-HE)), PS/(VB*(1-HE)), 1.5* PS/VE, sqrt(2),0.2)

################## simulate discrete time data
if simlongpath
    # Simulate one long path
    # Random.seed!(2)
    x0 = ℝ{2}(0.0, 0.0)
    #x0 = ℝ{2}(-8.0, 1.0)
    T_long = 12.0#10.0
    dt = 0.0001
    tt_long = 0.:dt:T_long
    W_long = sample(tt_long, Wiener{ℝ{dp}}())
    X_long = solve(Euler(), x0, W_long, Ptrue)
    # Extract partial observations
    lt = length(tt_long)
    # obsnum = 10
    # if obsnum > 2
    #     obsind = sort(sample(2:lt-1,obsnum-2,replace=false))
    #     pushfirst!(obsind,1)
    #     push!(obsind,lt)
    # elseif obsnum==2
    #     obsind = [1, lt]
    # else
    #     error("provide valid number of observations ")
    # end
    obsnum = 15
    if obsnum > 2
        obsind = 1:(lt÷obsnum):lt
        obsnum = length(obsind)
    elseif obsnum==2
         obsind = [1, lt]
    else
         error("provide valid number of observations ")
     end
    _pairs(collection) = Base.Generator(=>, keys(collection), values(collection))
    V_ = SamplePath(collect(_pairs(X_long))[obsind])
    V = SamplePath(V_.tt, map(y -> (L*y)[1:m] .+ (cholesky(Σ)).U' * randn(m), V_.yy))
end

obsnum = length(V)
segnum = obsnum-1

longpath = [Any[tt_long[j], d, X_long.yy[j][d]] for d in 1:2, j in 1:5:length(X_long) ][:]
obs = [Any[V.tt[j], dind, V.yy[j][dind]] for dind in 1:1, j in 1:length(V) ][:]
obsDf = DataFrame(time=map(x->x[1],obs), component = map(x->x[2],obs),value=map(x->x[3],obs) )
longpathDf = DataFrame(time=map(x->x[1],longpath), component = map(x->x[2],longpath),value=map(x->x[3],longpath) )

################## define gpudate function
function gpupdate(ν::SVector, H⁺::SMatrix, Σ, L, v)
    if all(diag(H⁺) .== Inf)
        H⁺_ = SMatrix(inv(L' * inv(Σ) * L))
        V_ = (L' * inv(Σ) * L)\(L' * inv(Σ) *  v)
        return V_, H⁺_
    else
        Z = I - H⁺*L'*inv(Σ + L*H⁺*L')*L
        return SVector(Z*H⁺*L'*inv(Σ)*v + Z*ν), Z*H⁺
    end
end

################## initialisation
tX = ℝ{d}
tW = ℝ{dp} #Float64 #ℝ{dp}# #
typeX = SamplePath{tX}
typeW = SamplePath{tW}
typeQ = Bridge.PartialBridgeνH{ℝ{d},Diffusion,DiffusionAux,ℝ{d},SArray{Tuple{d,d},Float64,2,d^2}}
XX = Vector{typeX}(undef, segnum)
XXtemp = Vector{typeX}(undef, segnum)
WW = Vector{typeW}(undef, segnum)
Q = Vector{typeQ}(undef, segnum)
Qᵒ = Vector{typeQ}(undef, segnum) # needed when parameter estimation is done

βinit = 8.0#Ptrue.β#3.0 ##10.0 #0.7 * Ptrue.μ
σ1init = 0.5
P = Diffusion(Ptrue.α, βinit,Ptrue.λ,Ptrue.μ,σ1init,Ptrue.σ2)
Pᵒ = deepcopy(P)
Pt = DiffusionAux(P)
Ptᵒ = DiffusionAux(P)
#H⁺i = Vector{typeof(Hend⁺)}(undef, segnum)

dtimp = 0.001  # mesh width for imputed paths
τ(t, T0, Tend) = T0 +(t-T0)*(2-(t-T0)/(Tend-T0))

# solve backward recursion on [0,T]
νend = SVector{d}(zeros(d))
Hend⁺ = SMatrix{d,d}(I/ϵ)
νend, Hend⁺ = gpupdate(νend, Hend⁺, Σ, L, V.yy[end])
Hrightmost⁺ = deepcopy(Hend⁺)
νrightmost = deepcopy(νend)

for i in segnum:-1:1
    # update on interval (t[i-1],t[i])
    tt_ = τ.(V.tt[i]:dtimp:V.tt[i+1],V.tt[i],V.tt[i+1])
    #tt_ = V.tt[i]:dtimp:V.tt[i+1],V.tt[i],V.tt[i+1]
    XX[i] = Bridge.samplepath(tt_, zero(tX)) # initialise
    WW[i] = Bridge.samplepath(tt_, zero(tW)) # initialise
#    H⁺i[i] = Hend⁺
    Q[i], νend, Hend⁺ = Bridge.partialbridgeνH(tt_, P, Pt, νend, Hend⁺)
    Qᵒ[i] = Q[i]
    global νend, Hend⁺ = gpupdate(νend, Hend⁺, Σ, L, V.yy[i])
end
#H⁺i[1] = Hend⁺

#elapsed_time= @elapsed begin

# simulate forward initial path
xstart = νend #+ √(Hend⁺) * randn(d)  # note that this is really ν(0)
for i in 1:segnum
    tt = Q[i].tt
    WW[i] = sample(tt,Wiener{tW}())
    Bridge.solve!(Euler(), XX[i], xstart, WW[i], Q[i])
    xstart = XX[i].yy[end] # starting point for next segment
end

################## visualisation of generated data, discrete time data and initial path
XXinit = vcat(XX...)
longpath = [Any[tt_long[j], D, X_long.yy[j][D]] for D in 1:d, j in 1:5:length(X_long) ][:]
initpath = [Any[XXinit.tt[j], D, XXinit.yy[j][D]] for D in 1:d, j in 1:1:length(XXinit) ][:]
obs = [Any[V.tt[j], dind, V.yy[j][dind]] for dind in 1:m, j in 1:length(V) ][:]
obsDf = DataFrame(time=map(x->x[1],obs), component = map(x->x[2],obs),value=map(x->x[3],obs) )
longpathDf = DataFrame(time=map(x->x[1],longpath), component = map(x->x[2],longpath),value=map(x->x[3],longpath) )
initpathDf = DataFrame(time=map(x->x[1],initpath), component = map(x->x[2],initpath),value=map(x->x[3],initpath) )

@rput obsDf
@rput longpathDf
@rput initpathDf
R"""
library(ggplot2)
library(tidyverse)

initpathDf$component <- as.factor(initpathDf$component)
longpathDf$component <- as.factor(longpathDf$component)
p <- ggplot() +
  ylab("") + geom_path(aes(x=time,y=value,group=component),data=initpathDf)+
   geom_point(aes(x=time,y=value),data=obsDf,colour="red")+
  geom_path(aes(x=time,y=value,colour=component),data=longpathDf)+theme_minimal()+
  theme(legend.position="bottom")
  #facet_wrap(~component,ncol=1,scales='free_y')
"""

XXᵒ = deepcopy(XX)
XXtemp = deepcopy(XX)
WWᵒ = deepcopy(WW)

# save some of the paths
XXsave = Any[]
if 0 in subsamples
    push!(XXsave, deepcopy(XX))
end

acc = 0
accparams = 0
mhsteps = 0
mhstepsparams = 0

Hzero⁺ = SMatrix{d,d}(0.1*I)

param(P) = [P.β P.σ1]
logπ(P) = logpdf(Gamma(1,100),P.β) + logpdf(Gamma(1,100),P.σ1)
logq(P, Pᵒ) = 0.0
function propose(σ, P)
#    Diffusion(P.α,P.β,P.λ,P.μ + σ * randn(),P.σ1,P.σ2)
    Diffusion(P.α,P.β + σ * randn(), P.λ,P.μ,P.σ1 + σ * randn(),P.σ2)
end

################## MH-algorithm

C = [param(P)]
for iter in 1:iterations
    finished = false
    klow = 1
    ind = segnum:-1:1
    kup = obsnum

    xstart = XX[1].yy[1]
    xstartᵒ = xstart
    updateparams = rand(Bool)#false#
    while !finished
        if updateparams # update all segments
            ind = (segnum==2) ? (1:1) : (segnum:-1:1)
            kup = obsnum
        else
            segnum_update = sample(1:obsnum-klow)   # number of segments to update (could also do min(2,obsnum-klow))
            kup = klow + segnum_update  # update on interval [t(klow), t(kup)]
            ind = (segnum_update==1) ? (1:1) : ((kup-1):-1:klow) # indices of segments to update
        end
        #hasbegin = ind[end]==1
        hasend = ind[1]==segnum

        # initialise νend, Hend⁺, νendᵒ, Hend⁺ᵒ
        νend  = hasend ? νrightmost :  XX[ind[1]].yy[end]  # initialise on rightmost segment
        Hend⁺ = hasend ? Hrightmost⁺ : Hzero⁺ # initialise on rightmost segment
        νendᵒ = deepcopy(νend)
        Hend⁺ᵒ = deepcopy(Hend⁺)

        # propose new parameter value
        Pᵒ = updateparams ? propose(.02, P) : P
        Ptᵒ = DiffusionAux(Pᵒ)

        # compute guiding term
        for i in ind
            tt = Q[i].tt
            Q[i], νend, Hend⁺ = Bridge.partialbridgeνH(tt, P, Pt, νend, Hend⁺)
            if updateparams
                Qᵒ[i], νendᵒ, Hend⁺ᵒ = Bridge.partialbridgeνH(tt, Pᵒ, Ptᵒ, νendᵒ, Hend⁺ᵒ)
            else
                Qᵒ[i], νendᵒ, Hend⁺ᵒ = Q[i], νend, Hend⁺
            end
            if !(i==last(ind))  # on the left endpoint we don't need to do gupdate
                νend, Hend⁺ = gpupdate(νend, Hend⁺, Σ, L, V.yy[i])
                if updateparams
                    νendᵒ, Hend⁺ᵒ = gpupdate(νendᵒ, Hend⁺ᵒ, Σ, L, V.yy[i])
                else
                    νendᵒ, Hend⁺ᵒ = νend, Hend⁺
                end
            end
        end

        ### to get the starting point updated correctly, we need to omit the very last step of gpupdate

        global diffll = 0.0  # difference in loglikehood

        # simulate guided proposal
        for i in reverse(ind)
            #tt = Q[i].tt
            if !updateparams
                sample!(WWᵒ[i], Wiener{ℝ{dp}}())
                WWᵒ[i].yy .= ρ * WW[i].yy + sqrt(1-ρ^2) * WWᵒ[i].yy
            else
                WWᵒ[i].yy .= WW[i].yy
            end
            if i==1  # starting point
                xstart = XX[1].yy[1]
                if updateparams
                    u = randn()
                    xstartᵒ = rand(Bool) ? (xstart .+ 0.1 *  ℝ{2}(u, -u)) : xstart  # with prob 0.5 propose joint update on starting point and parameter
                else # propose new starting point
                    u = randn()
                    xstartᵒ = xstart .+ 0.1 *  ℝ{2}(u, -u)
                end
                diffll += logpdfnormal(xstartᵒ-νendᵒ, Bridge.symmetrize(Hend⁺ᵒ))-logpdfnormal(xstart-νend, Bridge.symmetrize(Hend⁺))               # plus possibly log q(X0|X0o) = log q(X0o|X0)
                # need to put  it here, as xstart and xstartᵒ change later on
            else
                xstart = XXtemp[i-1].yy[end]
                xstartᵒ = XXᵒ[i-1].yy[end]
            end
            solve!(Euler(), XXtemp[i], xstart, WW[i], Q[i])
            solve!(Euler(), XXᵒ[i], xstartᵒ, WWᵒ[i], Qᵒ[i])
        end

        # update loglikelihood
        for i in ind
            diffll += llikelihood(LeftRule(), XXᵒ[i],  Qᵒ[i]) - llikelihood(LeftRule(), XXtemp[i],  Q[i])
        end
        if updateparams # this term is still not correct
            diffll +=  (V.tt[end]-V.tt[1]) *(tr(Bridge.B(0.0, Ptᵒ)) - tr(Bridge.B(0.0,Pt)))        + logπ(Pᵒ) - logπ(P)
        end

        println("ind:  ",ind,"  updateparms= ",updateparams,"  diff_ll:   ",round(diffll, digits=3))
        # MH step
        if log(rand()) <= diffll
            print("✓")
            for i in ind
                 XX[i], XXᵒ[i] = XXᵒ[i], XX[i]
                 WW[i], WWᵒ[i] = WWᵒ[i], WW[i]
            end
            P = Pᵒ
            Pt = Ptᵒ

            if updateparams
                accparams += 1
                push!(C, param(P))
            else
                acc += 1
            end
        end
        klow = kup # adjust counter
        finished = (klow==obsnum) # all segments have been updated (this is correct when a single sweep is done for parameter update)
        mhsteps += !updateparams
        mhstepsparams += updateparams
    end
    println("iteration finished")
    if iter in subsamples
        push!(XXsave, deepcopy(XX))
    end
end

@info "Done."*"\x7"^6

ave_acc_perc = 100*round(acc/mhsteps, digits=2)
ave_acc_percparams = 100*round(accparams/mhstepsparams, digits=2)
println("Average acceptance percentage: ",ave_acc_perc,"\n")
println("Average acceptance percentage params: ",ave_acc_percparams,"\n")

trueval = param(Ptrue)
est = DataFrame(iteration = 1:length(C), beta= map(x->x[1],C), sigma1=map(x->x[2],C))

@rput est
#@rput P
@rput trueval

R"""
plot.ts(est$beta,ylim=c(0,10)); abline(h=trueval[1],col='red')
lines(est$sigma1,col='blue'); abline(h=trueval[2],col='red')
"""

#error("STOP HERE")

limp = length(vcat(XXsave[1]...))
iterates = [Any[s,  vcat(XXsave[i]...).tt[j], dind, vcat(XXsave[i]...).yy[j][dind]] for dind in 1:d, j in 1:10:limp, (i,s) in enumerate(subsamples) ][:]
iteratesaverage = [Any[s,  vcat(XXsave[i]...).tt[j], mean(vcat(XXsave[i]...).yy[j])] for j in 1:10:limp, (i,s) in enumerate(subsamples) ][:]

write2csv = true
if write2csv
    # write long path to csv file
    f = open(outdir*"truepath.csv","w")
    headl = "time, component, value \n"
    write(f, headl)
    writedlm(f, longpath, ',')
    close(f)
    # write mcmc iterates to csv file
    fn = outdir*"iterates.csv"
    f = open(fn,"w")
    headl = "iteration, time, component, value \n"
    write(f, headl)
    writedlm(f,iterates,",")
    close(f)
    # also average of iterates
    fn = outdir*"iteratesaverage.csv"
    f = open(fn,"w")
    headl = "iteration, time,  value \n"
    write(f, headl)
    writedlm(f,iteratesaverage,",")
    close(f)
    # write observations to csv file
    fn = outdir*"observations.csv"
    f = open(fn,"w")
    headl = "time, component, value \n"
    write(f, headl)
    writedlm(f,obs,",")
    close(f)

#    CSV.write("parestimates.csv",est)
    fn = outdir*"parestimates.csv"
    f = open(fn,"w")
    headl = "iterate, beta, sigma1 \n"
    write(f,headl)
    truevals = [0 param(Ptrue)]
    estim = hcat(1:length(C), first.(C),last.(C))
    writedlm(f,vcat(truevals,estim),",")
    # iterate 0 corresponds to the true values
    close(f)
end

# plotting the results
iteratesDf = DataFrame(iteration = map(x->x[1],iterates), time=map(x->x[2],iterates), component = map(x->x[3],iterates),value=map(x->x[4],iterates) )
iteratesaverageDf = DataFrame(iteration = map(x->x[1],iteratesaverage), time=map(x->x[2],iteratesaverage), value = map(x->x[3],iteratesaverage) )
@rput iteratesDf
@rput iteratesaverageDf

R"""
library(ggplot2)
library(tidyverse)
iteratesDfsub <- iteratesDf #%>% filter(iteration > 5)
p <- ggplot() +
  geom_path(mapping=aes(x=time,y=value,colour=iteration,group=iteration),data=iteratesDfsub) +
  scale_colour_gradient(low='green',high='blue')+
   ylab("") + geom_point(aes(x=time,y=value),data=obsDf,colour="red")+
   geom_path(aes(x=time,y=value),data=longpathDf,colour="yellow")+
   facet_wrap(~component,ncol=1,scales='free_y') +
theme_minimal()

show(p)
"""

R"""
ggplot() +
  geom_path(mapping=aes(x=time,y=value,colour=iteration,group=iteration),data=iteratesaverageDf) +
  scale_colour_gradient(low='green',high='blue')+
   ylab("") + geom_point(aes(x=time,y=value),data=obsDf,colour="red")+
   #geom_path(aes(x=time,y=value),data=longpathDf,colour="yellow")+
   theme_minimal()
"""

writeinfo = true
if writeinfo
    # write info to txt file
    fn = outdir*"info.txt"
    f = open(fn,"w")
    #write(f, "Elapsed time: ", string(elapsed_time),"\n")
    #write(f, "Choice of observation schemes: ",obs_scheme,"\n")
    #write(f, "Easy conditioning (means going up to 1 for the rough component instead of 2): ",string(easy_conditioning),"\n")
    write(f, "Number of iterations: ",string(iterations),"\n")
    write(f, "Skip every ",string(skip_it)," iterations, when saving to csv","\n\n")
    write(f, "Starting point: ",string(x0),"\n")
    write(f, "End time T: ", string(V.tt[end]),"\n")
    #write(f, "Endpoint v: ",string(v),"\n")
    write(f, "Noise Sigma: ",string(Σ),"\n")
    write(f, "Regularisation parameter epsilon", string(ϵ),"\n")
    write(f, "L: ",string(L),"\n\n")
    write(f, "Mesh width: ",string(dtimp),"\n")
    write(f, "rho (Crank-Nicholsen parameter: ",string(ρ),"\n")
    write(f, "Average acceptance percentage: ",string(ave_acc_perc),"\n\n")
    write(f, "Backward type parametrisation in terms of nu and H? ",string(νHparam),"\n")
    close(f)
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

# using filters from paper
# Allauzen et al. Filters for Efficient Composition of Weighted Finite-State Transducers
# main differences here are that we use ⟂ = 0 and base 1 index for the states
abstract type CompositionFilter end

export Trivial
"""
`Trivial`

Simplest composition filter, can be used for epsilon-free WFSTs.
"""
struct Trivial{W,D} <: CompositionFilter
  i3::D
  q_max::D
  rho::Vector{W}
  ϵL::D
  Trivial{W,D}(ϵL::D) where {W,D} = new{W,D}(1,2,[one(W),one(W)],ϵL)
end

function phi(filter::Trivial, e1, e2, q3)
  ϵL = filter.ϵL
  o_e1, i_e2 = get_olabel(e1), get_ilabel(e2)
  if (o_e1 == i_e2) ||
    (o_e1 == ϵL && i_e2 == 0) ||
    (o_e1 == 0 && i_e2 == ϵL)
    if iseps(get_ilabel(e1)) && iseps(get_olabel(e2))
      # avoids self eps:eps loops
      return e1, e2, 0
    else
      return e1, e2, 1
    end
  else
    return e1, e2, 0
  end
end

export EpsMatch
"""
`EpsMatch`

Can be used for WFSTs containing epsilon labels.
Avoids redundant epsilon paths and giving priority to those that match epsilon labels.

```julia
julia> A = txt2fst("0 1 a a 3.0
       1 2 b <eps> 77.0
       2 3 c <eps> 33.0
       3 4 d d 44.0
       4 90.0
       ",Dict("a"=>1,"b"=>2,"c"=>3,"d"=>4,"e"=>5))
WFST #states: 5, #arcs: 4, #isym: 5, #osym: 5
|1/0.0f0|
a:a/3.0f0 → (2)
(2)
b:ϵ/77.0f0 → (3)
(3)
c:ϵ/33.0f0 → (4)
(4)
d:d/44.0f0 → (5)
((5/90.0f0))

julia> B = txt2fst("0 1 a d 15.0
       1 2 <eps> e 19.0
       2 3 d a 12.0
       3 55.0
       ",Dict("a"=>1,"b"=>2,"c"=>3,"d"=>4,"e"=>5))
WFST #states: 4, #arcs: 3, #isym: 5, #osym: 5
|1/0.0f0|
a:d/15.0f0 → (2)
(2)
ϵ:e/19.0f0 → (3)
(3)
d:a/12.0f0 → (4)
((4/55.0f0))

julia> C = ∘(A,B;filter=EpsMatch)
WFST #states: 5, #arcs: 4, #isym: 5, #osym: 5
|1/0.0f0|
a:d/18.0f0 → (2)
(2)
b:e/96.0f0 → (3)
(3)
c:ϵ/33.0f0 → (4)
(4)
d:a/56.0f0 → 5)
((5/145.0f0))

```
"""
struct EpsMatch{W,D} <: CompositionFilter
  i3::D
  q_max::D
  rho::Vector{W}
  ϵL::D
  EpsMatch{W,D}(ϵL::D) where {W,D} =
  new{W,D}(1,3,[one(W),one(W),one(W)],ϵL)
end

function phi(filter::EpsMatch, e1, e2, q3)
  ϵL = filter.ϵL
  o_e1, i_e2 = get_olabel(e1), get_ilabel(e2)
  if (o_e1 == i_e2) &&
    !(iseps(o_e1) || o_e1 == ϵL) &&
    !(iseps(i_e2) || i_e2 == ϵL)
    return e1, e2, 1
  elseif (o_e1 == 0) && (i_e2 == 0) && (q3 == 1)
    return e1, e2, 1
  elseif (o_e1 == ϵL) && (i_e2 == 0) && (q3 != 3)
    return e1, e2, 2
  elseif (o_e1 == 0) && (i_e2 == ϵL) && (q3 != 2)
    return e1, e2, 3
  else
    return e1, e2, 0
  end
end

export EpsSeq
"""
`EpsSeq`

Can be used for WFSTs containing epsilon labels.
Avoids redundant epsilon paths.
Gives priority to epsilon paths consisting of epsilon-arcs in `A` followed by epsilon-arcs in `B`.

```julia
julia> A = txt2fst("0 1 a a 3.0
       1 2 b <eps> 77.0
       2 3 c <eps> 33.0
       3 4 d d 44.0
       4 90.0
       ",Dict("a"=>1,"b"=>2,"c"=>3,"d"=>4,"e"=>5))
WFST #states: 5, #arcs: 4, #isym: 5, #osym: 5
|1/0.0f0|
a:a/3.0f0 → (2)
(2)
b:ϵ/77.0f0 → (3)
(3)
c:ϵ/33.0f0 → (4)
(4)
d:d/44.0f0 → (5)
((5/90.0f0))

julia> B = txt2fst("0 1 a d 15.0
       1 2 <eps> e 19.0
       2 3 d a 12.0
       3 55.0
       ",Dict("a"=>1,"b"=>2,"c"=>3,"d"=>4,"e"=>5))
WFST #states: 4, #arcs: 3, #isym: 5, #osym: 5
|1/0.0f0|
a:d/15.0f0 → (2)
(2)
ϵ:e/19.0f0 → (3)
(3)
d:a/12.0f0 → (4)
((4/55.0f0))

julia> C = ∘(A,B;filter=EpsSeq)
WFST #states: 6, #arcs: 5, #isym: 5, #osym: 5
|1/0.0f0|
a:d/18.0f0 → (2)
(2)
b:ϵ/77.0f0 → (3)
(3)
c:ϵ/33.0f0 → (4)
(4)
ϵ:e/19.0f0 → (5)
(5)
d:a/56.0f0 → (6)
((6/145.0f0))
```
"""
struct EpsSeq{W,D} <: CompositionFilter
  i3::D
  q_max::D
  rho::Vector{W}
  ϵL::D
  EpsSeq{W,D}(ϵL::D) where {W,D} =
  new{W,D}(1,2,[one(W),one(W)],ϵL::D)
end

function phi(filter::EpsSeq, e1, e2, q3)
  ϵL = filter.ϵL
  o_e1, i_e2 = get_olabel(e1), get_ilabel(e2)
  if (o_e1 == i_e2) &&
    !(iseps(o_e1) || o_e1 == ϵL) &&
    !(iseps(i_e2) || i_e2 == ϵL)
    return e1, e2, 1
  elseif (o_e1 == 0) && (i_e2 == ϵL) && (q3 == 1)
    return e1, e2, 1
  elseif (o_e1 == ϵL) && (i_e2 == 0)
    return e1, e2, 2
  else
    return e1, e2, 0
  end
end


text:
# Squares building

W,(WW,EW,FW) = Lar.cuboid([1,1], true)
square = (W,EW)

squares = Lar.Struct([  square,  Lar.t(-1, -1),  Lar.s(3, 3),  square  ])

V,EV = Lar.struct2lar(squares)

model = CAGD.Model(V);
CAGD.addModelCells!(model, 1, EV);

# Lines adding
lines = CAGD.Model([
    -2.0  1.0  0.0  3.0
     0.0  3.0  3.0  0.0
]);
CAGD.addModelCells!(lines, 1, [[1, 2], [3, 4]]);

# Merge the models
CAGD.uniteModels!(model, lines)

# Planar Arrangment

A_model = CAGD.planar_arrangement(model, sparse(Int8[0,0,0,0,1,1,1,1,0,0]));

A_model.G
Matrix(A_model.T[1])
Matrix(A_model.T[2])

# Display Models

if todisplay
    VV = [[k] for k=1:size(A_model.G ,2)];
    triangulated_faces = CAGD.triangulate2D(A_model);
    FVs = convert(Array{Lar.Cells}, triangulated_faces);

    GL.VIEW(GL.GLExplode(A_model.G, FVs, 1.0, 1.0, 1.0, 99, 1));
    GL.VIEW(GL.GLExplode(A_model.G, FVs, 1.2, 1.2, 1.2, 99, 1));
end;


text:
# These are the varints used in blocks! There are other varints used in the
# txindex data base!

# TODO: there should be a type VarInt and a read function

# TODO: little/big endian
"""
    read_varint(io::IO)::UInt64

Read a varint from `io`, find documentation
(here)[https://en.bitcoin.it/wiki/Protocol_documentation]
"""
function read_varint(io::IO)::UInt64
    i = read(io, UInt8)
    if     i <  0xFD UInt64(i)
    elseif i == 0xFD UInt64(read(io, UInt16))
    elseif i == 0xFE UInt64(read(io, UInt32))
    elseif i == 0xFF UInt64(read(io, UInt64))
    end
end

read_varint(x::BCIterator) = read_varint(x.io)

# We have to assume that the varint was always saved in the most space saving
# way for the following to work, no idea if this is guaranteed:
#
# TODO: this could return a tuple, this would avoid allocations but make this
# type unstable
"""
    to_varint(x::T)::Vector{UInt8}

Take an unsigned type and convert it into an varint
(here)[https://en.bitcoin.it/wiki/Protocol_documentation]

"""
function to_varint(x::T) where T
    if     x < T(0xFD)        UInt8[ UInt8(x) ]
    elseif x < T(0xFFFF)      UInt8[ 0xFD, to_byte_tuple(UInt16(x))... ]
    elseif x < T(0xFFFF_FFFF) UInt8[ 0xFE, to_byte_tuple(UInt32(x))... ]
    else                      UInt8[ 0xFF, to_byte_tuple(UInt64(x))... ]
    end
end


text:
module ResNet9

using Flux
using Flux: convfilter, Zeros
using Base.Iterators: takewhile

export resnet9

conv(chs) = Conv((3, 3), chs, identity, pad = (1, 1), bias = false)

function convblock(chs::Pair; pool = false)
    _, out = chs
    cv = conv(chs)
    bn = BatchNorm(out, relu)
    pool && return cv, bn, MaxPool((2, 2))
    cv, bn
end

convblock(ch::Int) = convblock(ch=>ch)

resblock(ch) = SkipConnection(Chain(convblock(ch)..., convblock(ch)...), +)

function classifier(ch, nc, drp)
    fl = [GlobalMaxPool(), Flux.flatten]
    fc = []
    if drp > 0
         push!(fc, Dropout(drp))
    end
    (nc > 2) ? push!(fc, Dense(ch, nc), softmax) : push!(fc, Dense(ch, 1, sigmoid))
    fl..., fc...
end

"""
    resnet9(;inchannels, nclasses, basewidth = 64, expansion = 2)
Build a ResNet9 network.
* inchannels: number of input channels.
* nclasses: number of classes, if `nclasses` = 2 the network ends with a fully
  connected layer with a sigmoid activation function, otherwise it ends with a
  fully connected layer followed by a softmax function.
* basewidth: base number of channels.
* expansion: factor of channels expansion.
"""
function resnet9(;inchannels, nclasses, dropout = 0, basewidth = 64, expansion = 2) # groupnorm
    ch1 = basewidth
    ch2 = expansion * ch1
    ch3 = expansion * ch2
    ch4 = expansion * ch3
    Chain(convblock(inchannels=>ch1)...,                      # Input layer
          convblock(ch1=>ch2, pool = true)..., resblock(ch2), # Layer 1
          convblock(ch2=>ch3, pool = true)...,                # Layer 2
          convblock(ch3=>ch4, pool = true)..., resblock(ch4), # Layer 3
          classifier(ch4, nclasses, dropout)...)              # Classifier
end

"""
    resnet9(model::Chain; inchannels, nclasses, dropout)
Build a ResNet9 model from an existing (pre-trained) model and change the first
or last layers according to `inchannels`, `nclasses` and `dropout`.
"""
function resnet9(model::Chain;
                 inchannels = inchan(model),
                 nclasses = nclass(model),
                 dropout = dpout(model))
    ich = inchan(model)
    if inchannels != ich
        l = model[1]
        iout = size(l.weight, ndims(l.weight))
        model = Chain(conv(inchannels=>iout), model[2:end]...)
    end
    if nclasses != nclass(model)
        fc = fclayer(model)
        fcin = size(fc.weight, 2)
        i = (findall(l -> isa(l, GlobalMaxPool), model) |> first) - 1 #takewhile(l -> !isa(l, GlobalMaxPool), model) |> collect
        back = model[1:i]
        front = classifier(fcin, nclasses, dropout)
        model = Chain(back..., front...)
    elseif dropout != dpout(model)
        i = (findall(l -> isa(l, GlobalMaxPool), model) |> first) + 1
        j = findall(l -> isa(l, Dense), model) |> first
        back = model[1:i]
        front = model[j:end]
        front = (dropout > 0) ? [Dropout(dropout), front...] : front
        model = Chain(back..., front...)
    end
    model
end

inchan(m) = size(m[1].weight, ndims(m[1].weight) - 1)

fclayer(m) = filter(l -> isa(l, Dense), m.layers) |> first

nclass(m) = fclayer(m) |> l -> size(l.weight, 1)

dpout(m) = m[isa.(m, Dropout)] |> c -> (isempty(c.layers)) ? 0 : c.layers[1].p

# """
#     resnet9(model::Chain; nclasses, dropout)
# Build a ResNet9 model from an existing (pre-trained) model and change the first
# layer according to `inchannels`.
# """
# function resnet9(model::Chain; inchannels)
#     ch = nchout(model[1:2])
#     Chain(convblock(inchannels=>ch)..., model[3:end]...)
# end

# function modelarch(model::Chain)
#     l = model[1]
#     ich = size(l.weight, ndims(l.weight) - 1)
#     fcin, fcout = 0, 0
#     for l ∈ model
#         if isa(l, Dense)
#             fcin = size(l.weight, 2)
#             fcout = size(l.weight, 1)
#         end
#     end
#     ich, fcin, fcout
# end

# function ncl(m::Chain)
#     # for l ∈
#     #     takewhile(l -> !isa(l, Dense), m)

#     # size(d.weight, 1)
# end

# nchin(l) = 0
# nchin(l::Conv) = size(l.weight, ndims(l.weight) - 1)
# nchin(l::BatchNorm) = l.chs
# nchin(l::SkipConnection) = nchin(l.layers)
# function nchin(c::Chain)
#     nch = 0
#     for l ∈ c
#         nch = nchin(l)
#         if nch > 0
#             break
#         end
#     end
#     if nch == 0
#         error("Failed to find number of channels.") |> throw
#     end
#     nch
# end

# nchout(l) = 0
# nchout(l::Conv) = size(l.weight, ndims(l.weight))
# nchout(l::BatchNorm) = l.chs
# nchout(l::SkipConnection) = nchout(l.layers)
# function nchout(c::Chain)
#     layers = [c...] |> reverse
#     nch = 0
#     for l ∈ layers
#         nch = nchout(l)
#         if nch > 0
#             break
#         end
#     end
#     if nch == 0
#         error("Failed to find number of channels.") |> throw
#     end
#     nch
# end

end