text:
"""
The `apply2list` function applies the function `fun`, the first argument, for
each element on `list` but avoiding `getfield` calls inside the loop.
The second argument of the macro is the `PairwiseListMatrix` that will be
iterated. `fun` should take two arguments; the first is the `list` field of
the `PairwiseListMatrix`, and the second is the index over that list at a
given iteration.

```jldoctest
julia> using PairwiseListMatrices

julia> plm = PairwiseListMatrix([1,2,3], false)
3×3 PairwiseListMatrix{Int64,false,Array{Int64,1}}:
 0  1  2
 1  0  3
 2  3  0

julia> apply2list((list, k) -> println(list[k]), plm)
1
2
3

```
"""
function apply2list(fun, plm)
    list = plm.list
    for k in 1:length(list)
        fun(list, k)
    end
end

"""
The `apply2diag` function applies the function `fun`, the first argument,
for each element on the `diag` field of a`PairwiseListMatrix{T,false,VT}`
but avoiding `getfield` calls inside the loop. The second argument of the
macro is the `PairwiseListMatrix` that will be iterated. `fun` should take
two arguments; the first is the `diag` field of the `PairwiseListMatrix`,
and the second is the index over that vector of diagonal elements at a
given iteration.

```jldoctest
julia> using PairwiseListMatrices

julia> plm = PairwiseListMatrix([1,2,3], false)
3×3 PairwiseListMatrix{Int64,false,Array{Int64,1}}:
 0  1  2
 1  0  3
 2  3  0

julia>apply2diag((diag, k) -> diag[k] += 10k, plm)

julia> plm
3×3 PairwiseListMatrix{Int64,false,Array{Int64,1}}:
 10   1   2
  1  20   3
  2   3  30

```
"""
function apply2diag(fun, plm::PairwiseListMatrix{T,false,VT}) where {T, VT}
    diag = plm.diag
    for k in 1:length(diag)
        fun(diag, k)
    end
end

"""
The `apply2upper` function applies the `fun` function over the upper triangular
part of the `PairwiseListMatrix`. Set `use_diag` to `true` if the diagonal
needs to be included in the iteration. The function should take four
arguments, first the `list` and diag fields of the `PairwiseListMatrix`,
the second is the index `k` over that `list`, and the third and fourth are
the `i` and `j` indexes for the position `k` in the upper triangular
part of the matrix.

```jldoctest
julia> using PairwiseListMatrices

julia> plm = PairwiseListMatrix([1,2,3], true)
2×2 PairwiseListMatrix{Int64,true,Array{Int64,1}}:
 1  2
 2  3

julia> mat = zeros(Int, 2, 2)
2×2 Array{Int64,2}:
 0  0
 0  0

julia> apply2upper((list, k, i, j) -> mat[i, j] = list[k], plm; use_diag = true)

julia> mat
2×2 Array{Int64,2}:
 1  2
 0  3

```
"""
function apply2upper(fun, plm::PairwiseListMatrix{T,false,VT}; use_diag::Bool=false) where {T, VT}
    N = plm.nelements
    if use_diag
        k = 0
        diag = plm.diag
        list = plm.list
        for i in 1:N
            for j in i:N
                if i != j
                    k += 1
                    fun(list, k, i, j)
                else
                    let list = diag, k = i
                        fun(list, k, i, j)
                    end
                end
            end
        end
    else
        k = 0
        list = plm.list
        for i in 1:(N-1)
            for j in (i+1):N
                k += 1
                fun(list, k, i, j)
            end
        end
    end
end

function apply2upper(fun, plm::PairwiseListMatrix{T,true,VT}; use_diag::Bool=false) where {T, VT}
    N = plm.nelements
    if hasdiagonal(plm)
        if use_diag
            k = 0
            list = plm.list
            for i in 1:N
                for j in i:N
                    k += 1
                    fun(list, k, i, j)
                end
            end
        else
            k = 0
            list = plm.list
            for i in 1:N
                for j in i:N
                    k += 1
                    if i != j
                        fun(list, k, i, j)
                    end
                end
            end
        end
    end
end


text:
include("test_exports.jl")
include("test_static_clip.jl")
include("test_plane_tet_intersection.jl")
include("test_poly_eight.jl")
include("test_utility.jl")
include("test_quadrature.jl")

ext: jl


text:
using DoubleFloats, AccurateArithmetic
using Test

using Random
phi = Base.MathConstants.golden

include("basic.jl")
include("concrete_accuracy.jl")
include("function_accuracy.jl")
include("arithmetic.jl")


text:
@testset "Unindexed table basics" begin
    va = [1,   2,   3]
    vb = [2.0, 4.0, 6.0]
    t = Table{(:a, :b), NamedTuple{(:a,:b), Tuple{Int, Float64}}, Tuple{Vector{Int}, Vector{Float64}}, Tuple{}}((va, vb), ())

    @test Table(a = va, b = vb) === t
    @test Table((a = va, b = vb)) === t
    @test Table((a = va, b = vb), ()) === t
    @test Table{(:a, :b)}(b = vb, a = va) === t
    @test Table{(:a, :b)}((b = vb, a = va)) === t
    @test Table{(:a, :b)}((b = vb, a = va), ()) === t

    @test colnames(t) === (:a, :b)
    @test columns(t) === (a = va, b = vb)
    @test MinimumViableTables.getdata(t) === (va, vb)
    @test t.a === va
    @test t.b === vb
    @test_throws ErrorException t.c
    @test length(t) === 3
    @test size(t) === (3,)
    @test axes(t) === (Base.OneTo(3),)
    @test keys(t) === Base.OneTo(3)

    @test @inferred(t[1]) === (a = 1, b = 2.0)
    @test @inferred(t[2]) === (a = 2, b = 4.0)
    @test @inferred(t[3]) === (a = 3, b = 6.0)

    @test @inferred(first(t)) === (a = 1, b = 2.0)

    @test @inferred(similar(t)) isa typeof(t)
    @test @inferred(copy(t)) == t
    @test !(copy(t) === t)

    @test project(t, (:a,)) === Table(a = va)
    @test project(t, (:b,)) === Table(b = vb)
    @test project(t, (:b, :a)) === Table(b = vb, a = va)

    @test @inferred(Project((:a,))(t)) === Table(a = va)
    @test @inferred(Project((:b,))(t)) === Table(b = vb)
    @test @inferred(Project((:b, :a))(t)) === Table(b = vb, a = va)
end


text:
################################################################################
# MIT License
#
# Copyright (c) 2021 Hajime Nakagami<nakagami@gmail.com>
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
################################################################################

mutable struct Connection <: DBInterface.Connection
    wp::Firebird.WireProtocol
    transaction::Transaction

    function Connection(host::String, username::String, password::String, db_name::String, opts::Dict)
        username = uppercase(username)
        port = UInt16(haskey(opts, :port) ? opts[:port] : 3050)
        wire_crypt = Bool(haskey(opts, :wire_crypt) ? opts[:wire_crypt] : true)
        create_new = Bool(haskey(opts, :create_new) ? opts[:create_new] : false)
        page_size = Int32(haskey(opts, :page_size) ? opts[:page_size] : 4096)
        timezone = haskey(opts, :timezone) ? opts[:timezone] : ""

        wp = WireProtocol(host, username, password, port)
        client_public, client_secret = get_client_seed()
        _op_connect(wp, db_name, username, password, wire_crypt, client_public)
        parse_connect_response(wp, username, password, wire_crypt, client_public, client_secret)
        if create_new
            _op_create(wp, db_name, username, password, timezone, page_size)
        else
            _op_attach(wp, db_name, username, password, timezone)
        end
        db_handle, _, _ = _op_response(wp)
        if db_handle < 0
            throw(DomainError("connection error"))
        end
        wp.db_handle = db_handle
        transaction = Transaction(wp)
        conn = new(wp, transaction)
    end

end

function isopen(conn::Connection)::Bool
    isopen(conn.wp.chan)
end

"""
    DBInterface.connect(Firebird.Connection, host::String, user::String, passwd::String; db::String, db_name::String)

Connect to a Firebird database
"""
DBInterface.connect(::Type{Connection}, host::String, user::String, password::String, db_name::String; kwargs...) =
    Connection(host, user, password, db_name, Dict(kwargs))

"""
    DBInterface.close!(Firebird.Connection)

Close a `Firebird.Connection` opened by `DBInterface.connect`.
"""
function DBInterface.close!(conn::Connection)
    close!(conn.transaction)
    close!(conn.wp.chan)
end


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

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 el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 7e46e0b6-04c6-11eb-18f7-bf3077b15a72
begin
	using Pkg
	Pkg.activate(".")
end

# ╔═╡ 2b54d586-0474-11eb-22bb-73ff01024985
begin
	using Revise
	using InteractiveUtils
	using PlutoUI
	using RvSpectMLBase
	using RvSpectML
	using EchelleInstruments
	using DataFrames
	using Query
	using Statistics
	using Plots
	using RvSpectMLPlots
end

# ╔═╡ a4d84168-048d-11eb-1c8f-0fa3136a1e0f
size(order_ccfs1,2)

# ╔═╡ f06adcc4-0480-11eb-08ba-5dc0b9d3062b
md"Mask width factor: $(@bind mask_scale_factor Slider(0.1:0.1:5; default=2, show_value=true))"

# ╔═╡ 8bf3f9bc-0479-11eb-144e-39bcfb215074
md"CCF Residuals? $(@bind plt_ccf_resid CheckBox(default=false))
CCF Timeseries? $(@bind plt_ccf_timeseries CheckBox(default=false))
CCF Orders? $(@bind plt_ccf_orders CheckBox(default=false))
CCF Order Timeseries? $(@bind plt_ccf_order_timeseries CheckBox(default=false))"

# ╔═╡ e0f9e742-048c-11eb-171f-4f1f3dea0f21
  begin
	if plt_ccf_orders
	local plt = plot()

	for ord in 1:size(order_ccfs1,2)
		if !plt_ccf_resid
			plot!(plt,v_grid_order_ccfs1, order_ccfs1[:,ord,1]./maximum(order_ccfs1[:,ord,1],dims=1) .+ 0.00*ord, label=:none)
			plot!(plt,v_grid_order_ccfs2, order_ccfs2[:,ord,1]./maximum(order_ccfs2[:,ord,1],dims=1) .+ 0.00*ord, label=:none)
		else
				plot!(plt,v_grid_order_ccfs1, order_ccfs1[:,ord,1]./maximum(order_ccfs1[:,ord,1],dims=1) .-  order_ccfs2[:,ord,1]./maximum(order_ccfs2[:,ord,1],dims=1) .+ 0.00*ord, label=:none)
		end
	end
    #local colorscale = cgrad(:balance)
    #plt = heatmap(v_grid_order_ccfs,collect(1:size(order_ccfs,2)),zvals', c=colorscale, clims=(-maximum(abs.(zvals)),maximum(abs.(zvals))) )
    #add_time_gap_lines(plt,order_list_timeseries.times)
    xlabel!("v (m/s)")
    ylabel!("Order #")
    title!("CCF(v,t)-<CCF>(v) vs time")
    plt
	end
end

# ╔═╡ 4156a452-0476-11eb-2a93-6f6ba585ab3f
md"### Setup environment"

# ╔═╡ 37855416-0477-11eb-315c-4b4bca632ae5
begin
	expres_data_path = joinpath(homedir(),"Data/EXPRES/inputs")
	target_subdir = "101501"
 	fits_target_str = "101501"
 	output_dir = "examples/output/"
    df_files = EchelleInstruments.make_manifest(expres_data_path, target_subdir, EXPRES );
end;

# ╔═╡ 3ca4f348-0477-11eb-35a9-4ba6addfd1b4
begin
	# Configuration from param.in
	linelist_for_ccf_filename = "G8.espresso.mas"
	ccf_mid_velocity = -5e3
	max_spectra_to_use = 50
	df_files_use = df_files |>
	    @filter( _.target == fits_target_str ) |>
	    @take(max_spectra_to_use) |>
	    DataFrame;
end;

# ╔═╡ 8b7b2098-0486-11eb-2abe-c530be84ba4e
"MJD = $(df_files_use[obs_highlight,:bjd])  airmass = $(df_files_use[obs_highlight,:airmass])   epoch = $( df_files_use[obs_highlight,:expres_epoch])"

# ╔═╡ 4ad9e038-0477-11eb-2363-0f093d35f3ca
md"#### Read data files & extract each order"

# ╔═╡ 334533ae-0478-11eb-2c9e-07eedf379844
begin
	pipeline_plan = PipelinePlan()
	dont_make_plot!(pipeline_plan, :movie)
end;

# ╔═╡ 4d93dcfc-0477-11eb-2619-c10f7fd1b909
begin
	spectra = map(row->EXPRES.read_data(row), eachrow(df_files_use))
	dont_need_to!(pipeline_plan,:read_spectra)
end;

# ╔═╡ abdf1064-0487-11eb-1f99-eb6563e505d7
md"Select Plot Obs# $(@bind obs_idx1 Slider(1:length(spectra),default=1,show_value=true))    Select Ref Obs# $(@bind obs_idx2 Slider(1:length(spectra),default=2,show_value=true))"

# ╔═╡ 5e945f0e-0490-11eb-1745-65666efed48b
obs_idx2

# ╔═╡ eb708bb0-0476-11eb-0ff3-bd45dc247e5a
function extract_orders(data::AbstractVector{ST}; verbose::Bool = false, orders_to_use = RvSpectMLBase.orders_to_use_default(first(data).inst) ) where { ST<:AbstractSpectra }
      order_list_timeseries = make_order_list_timeseries(data, orders_to_use=orders_to_use)
      order_list_timeseries = filter_bad_chunks(order_list_timeseries,verbose=true)
end;

# ╔═╡ 29638118-0475-11eb-19db-bdcac9172074
begin
	olt1 = make_order_list_timeseries(spectra, orders_to_use=[obs_idx1])
    #olt1 = filter_bad_chunks(olt1,verbose=true)
	dont_need_to!(pipeline_plan,:extract_orders)
end

# ╔═╡ ff7bab4c-04be-11eb-3de1-27f14eef3b91
begin
	olt2 = make_order_list_timeseries(spectra, orders_to_use=[obs_idx2])
    #olt1 = filter_bad_chunks(olt1,verbose=true)
end

# ╔═╡ 07c981e6-04bf-11eb-3c02-b9d8d6018719
begin
	order_list_timeseries = make_order_list_timeseries(spectra, orders_to_use=[obs_idx1])
    #order_list_timeseries = filter_bad_chunks(order_list_timeseries,verbose=true)
end

# ╔═╡ b8eba050-04bf-11eb-16b7-358930a15ffd
RvSpectMLBase.make_orders_into_chunks(spectra[1], first(spectra).inst, orders_to_use = [5])

# ╔═╡ 059cc5ca-0478-11eb-3f29-6b50cb83c78e
line_list_df = prepare_line_list(linelist_for_ccf_filename, spectra, pipeline_plan,  v_center_to_avoid_tellurics=ccf_mid_velocity, Δv_to_avoid_tellurics = 30e3, recalc=true);

# ╔═╡ 08c3f4d2-047d-11eb-0b79-e1f9035e3df7
md"#### Actual calculations"

# ╔═╡ 0bae871e-0484-11eb-0b0e-234b276e02f3
begin
	lambda_min = minimum(line_list_df.lambda)
	lambda_max = maximum(line_list_df.lambda)
end;

# ╔═╡ 3c588554-0482-11eb-2fff-cd4645f903f5
@bind mask_lambda_min  Slider(lambda_min:10:lambda_max; default=lambda_min, show_value=true)

# ╔═╡ d93361c2-0483-11eb-3da2-db2a8ff59b90
@bind mask_lambda_max Slider(lambda_min:10:lambda_max; default=lambda_max, show_value=true)

# ╔═╡ 17d6a8b4-0482-11eb-330b-cd9ed27943bf
line_list_use_df = @from i in line_list_df begin
    @where (mask_lambda_min <= i.lambda <= mask_lambda_max)
	#@where (mask_weight_min <= i.weight <= mask_weight_max)
    @select i
    @collect DataFrame
end;

# ╔═╡ 30571030-04c7-11eb-1516-dbd5589c12b9
function ccf_orders(order_list_timeseries::AbstractChunkListTimeseries,  line_list_df::DataFrame; verbose::Bool = false, range_no_mask_change::Real=30e3, ccf_mid_velocity::Real=0.0, mask_scale_factor::Real=1, mask_type::Symbol = :tophat )
      if mask_type == :tophat
         mask_shape = EchelleCCFs.TopHatCCFMask(order_list_timeseries.inst, scale_factor=mask_scale_factor)
      elseif mask_type == :gaussian
         mask_shape = EchelleCCFs.GaussianCCFMask(order_list_timeseries.inst, σ_scale_factor=mask_scale_factor)
      elseif mask_type == :supergaussian
         mask_shape = EchelleCCFs.SuperGaussianCCFMask(order_list_timeseries.inst, σ_scale_factor=mask_scale_factor)
      elseif mask_type == :halfcos
         mask_shape = EchelleCCFs.CosCCFMask(order_list_timeseries.inst, scale_factor=mask_scale_factor)
      else
         @error("Requested mask shape (" * string(mask_type) * " not avaliable.")
      end

      line_list = EchelleCCFs.BasicLineList(line_list_df.lambda, line_list_df.weight)
      ccf_plan = EchelleCCFs.BasicCCFPlan(mask_shape = mask_shape, line_list=line_list, midpoint=ccf_mid_velocity, range_no_mask_change=range_no_mask_change)
      v_grid = EchelleCCFs.calc_ccf_v_grid(ccf_plan)
      order_ccfs = EchelleCCFs.calc_order_ccf_chunklist_timeseries(order_list_timeseries, ccf_plan)
	return order_ccfs
end

# ╔═╡ f1def888-047b-11eb-1b72-63189bf2650a
(order_ccfs, v_grid_order_ccfs) = ccf_orders(order_list_timeseries, line_list_df);

# ╔═╡ e22f9f90-04c7-11eb-2519-a74e022bcd38
order_ccfs

# ╔═╡ 85b04c2e-047c-11eb-330a-e321c5501b45
  begin
	if plt_ccf_orders
	local plt = plot()
	#local zvals = order_ccfs[:,:,obs_highlight]./maximum(order_ccfs[:,:,obs_highlight],dims=1).-0*plt_ccf_resid*mean(order_ccfs[:,:,obs_highlight]./maximum(order_ccfs[:,:,obs_highlight],dims=1),dims=2)
	local zvals = order_ccfs[:,:,obs_idx1]./maximum(order_ccfs[:,:,idx1],dims=1).- plt_ccf_resid.*reshape(mean(order_ccfs[:,:,obs_idx2]./maximum(order_ccfs[:,:,obs_idx2],dims=1),dims=3),size(order_ccfs)[1:2])
    local colorscale = cgrad(:balance)
    plt = heatmap(v_grid_order_ccfs,collect(1:size(order_ccfs,2)),zvals', c=colorscale, clims=(-maximum(abs.(zvals)),maximum(abs.(zvals))) )
    add_time_gap_lines(plt,order_list_timeseries.times)
    xlabel!("v (m/s)")
    ylabel!("Order #")
    title!("CCF(v,t)-<CCF>(v) vs time")
    plt
	end
end

# ╔═╡ 5496cbaa-0476-11eb-0637-4d3d72087b10
md"### Pluto utils"

# ╔═╡ dd0b46c2-0472-11eb-08ef-b9b5984cf08c
@bind pluto_width Radio(["800px"=>"small", "1200px"=>"medium", "2400px"=>"large", "max-width"=>"max"], default="max-width")

# ╔═╡ d2a8c018-047e-11eb-2f3f-83b241d73c3d
begin
	mask_shape_keys = ["Gaussian", "Tophat", "Half Cos"]
	mask_shape_vals = [:gaussian, :tophat,:halfcos]
	mask_shape_dict = Dict(zip(mask_shape_keys,mask_shape_vals))
end;

# ╔═╡ 82a0f3b4-047f-11eb-2fb0-45fe12555aed
md"Mask shape: $(@bind mask_shape Radio(map(x->Pair(x[1],x[2]),zip(mask_shape_keys,mask_shape_vals)),default=\"Tophat\") )"

# ╔═╡ c623d82e-048c-11eb-3bde-43cb4234f60f
begin
	(ccfs1, v_grid1) = ccf_total(olt1, line_list_use_df, pipeline_plan,  mask_type=mask_shape_dict[mask_shape], mask_scale_factor=mask_scale_factor, ccf_mid_velocity=ccf_mid_velocity, recalc=true,  calc_ccf_var=false)
	#(order_ccfs1, v_grid_order_ccfs1) = ccf_orders(olt1, line_list_df,mask_type=mask_shape_dict[mask_shape], mask_scale_factor=mask_scale_factor, ccf_mid_velocity=ccf_mid_velocity)
end;

# ╔═╡ eeea64e8-0488-11eb-029e-8d0c711f4715
begin
	local plt = plot()
	plot!(plt, v_grid1,ccfs1./maximum(ccfs1,dims=1) .- 0.0 .* mean(ccfs1./maximum(ccfs1,dims=1),dims=2),labels=:none)
	#plot!(plt, v_grid_order_ccfs1, ccfs1./maximum(ccfs1,dims=1).-plt_ccf_resid.*ccfs2./maximum(ccfs2,dims=1) ,labels=:none)
	#plot!(plt, v_grid_order_ccfs2, ccfs2./maximum(ccfs2,dims=1).-plt_ccf_resid.*ccfs2./maximum(ccfs2,dims=1) ,labels=:none)
	#scatter!(plt, v_grid2, ccfs1[:,obs_idx1]./maximum(ccfs1[:,obs_idx1]) .- plt_ccf_resid*mean(ccfs./maximum(ccfs,dims=1),dims=2),markersize=1.2,labels=:none)
	xlabel!("RV (m/s)")
	plt_ccf_resid ? ylabel!("CCF") : ylabel!("CCF")
	plt
end

# ╔═╡ 75c284dc-0489-11eb-34ff-d90b6385687a
ccfs1./maximum(ccfs1,dims=1)

# ╔═╡ 44c3e344-048e-11eb-2bac-cd263f61abae
begin
	(ccfs2, v_grid2) = ccf_total(olt2, line_list_use_df, pipeline_plan,  mask_type=mask_shape_dict[mask_shape], mask_scale_factor=mask_scale_factor, ccf_mid_velocity=ccf_mid_velocity, recalc=true,  calc_ccf_var=false)
	#(order_ccfs2, v_grid_order_ccfs2) = ccf_orders(olt2, line_list_df, pipeline_plan,mask_type=mask_shape_dict[mask_shape], mask_scale_factor=mask_scale_factor, ccf_mid_velocity=ccf_mid_velocity,recalc=true)
	#ccfs2 = reshape(sum(order_ccfs2,dims=3),(size(order_ccfs2,1),size(order_ccfs2,2)))
	#v_grid2 = v_grid_order_ccfs2
	#ccfs2
end;

# ╔═╡ ec169d3a-048e-11eb-178a-75dc8b1b3d95
 extrema(ccfs1.-ccfs2)

# ╔═╡ 5328c0d2-0478-11eb-357c-e750df71d6ff
begin
	ENV["JULIA_NUM_THREADS"] = 4
	(ccfs, v_grid) = ccf_total(order_list_timeseries, line_list_use_df, pipeline_plan,  mask_type=mask_shape_dict[mask_shape], mask_scale_factor=mask_scale_factor, ccf_mid_velocity=ccf_mid_velocity, recalc=true,  calc_ccf_var=false)
	#((ccfs, ccf_vars), v_grid) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan,  mask_type=mask_shape_dict[mask_shape], mask_scale_factor=mask_scale_factor, ccf_mid_velocity=ccf_mid_velocity, recalc=true,  calc_ccf_var=true)
end;

# ╔═╡ 7f3e45ca-0478-11eb-2c41-4db39b42331e
begin
	if plt_ccf_timeseries
		local plt = plot()
	#plot!(plt, v_grid,ccfs./maximum(ccfs,dims=1) .- mean(ccfs./maximum(ccfs,dims=1),2),labels=:none)
	plot!(plt, v_grid, ccfs./maximum(ccfs,dims=1) .- plt_ccf_resid*mean(ccfs./maximum(ccfs,dims=1),dims=2),labels=:none)
	scatter!(plt, v_grid, ccfs[:,obs_highlight]./maximum(ccfs[:,obs_highlight]) .- plt_ccf_resid*mean(ccfs./maximum(ccfs,dims=1),dims=2),markersize=1.2,labels=:none)
	xlabel!("RV (m/s)")
	plt_ccf_resid ? ylabel!("CCF") : ylabel!("CCF")
	plt
	end
end

# ╔═╡ 78f39af8-0487-11eb-0948-bf2a888bd088
begin
	if plt_ccf_orders
	local plt = plot()
	#plot!(plt, v_grid,ccfs./maximum(ccfs,dims=1) .- mean(ccfs./maximum(ccfs,dims=1),2),labels=:none)
	plot!(plt, v_grid, ccfs./maximum(ccfs,dims=1) .- plt_ccf_resid*mean(ccfs./maximum(ccfs,dims=1),dims=2),labels=:none)
	scatter!(plt, v_grid, ccfs[:,obs_highlight]./maximum(ccfs[:,obs_highlight]) .- plt_ccf_resid*mean(ccfs./maximum(ccfs,dims=1),dims=2),markersize=1.2,labels=:none)
	xlabel!("RV (m/s)")
	plt_ccf_resid ? ylabel!("CCF") : ylabel!("CCF")
	plt
	end
end

# ╔═╡ 237fc58e-047c-11eb-0853-05ab5af54793
  begin
	if plt_ccf_timeseries
	local plt = plot()
	local zvals = ccfs./maximum(ccfs,dims=1).-plt_ccf_resid*mean(ccfs./maximum(ccfs,dims=1),dims=2)
    local colorscale = cgrad(:balance)
    plt = heatmap(v_grid,collect(1:size(ccfs,2)),zvals', c=colorscale, clims=(-maximum(abs.(zvals)),maximum(abs.(zvals))) )
    add_time_gap_lines(plt,order_list_timeseries.times)
    xlabel!("v (m/s)")
    ylabel!("Observation #")
    title!("CCF(v,t)-<CCF>(v) vs time")
    plt
	end
end

# ╔═╡ 2a705298-047a-11eb-3967-b3207243a5c3
begin
	mean_norm_ccf_total = mean(ccfs./maximum(ccfs,dims=1),dims=2)
	ccf1_norm = maximum(ccfs[:,obs_idx1])
end;

# ╔═╡ b2c72e50-0484-11eb-1d15-b953a0640a69
@bind mask_weight_min  Slider(0.0:0.05:0.95; default=0.0, show_value=true)

# ╔═╡ b27dcd6e-0484-11eb-0ab6-e73b9b4e80ec
@bind mask_weight_max Slider(0.05:0.05:1.0; default=1.0, show_value=true)

# ╔═╡ a3737596-0473-11eb-0902-258a8ea5a515
HTML("<style> main { max-width:" * pluto_width * "; } </style> ")

# ╔═╡ 7ae40ec2-0475-11eb-3191-f75c848caaa0
HTML("""<style> input[type="range"] { width:250px; } </style>""")

# ╔═╡ 38a87a16-0475-11eb-1e64-b99500c57983
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ Cell order:
# ╠═c623d82e-048c-11eb-3bde-43cb4234f60f
# ╠═44c3e344-048e-11eb-2bac-cd263f61abae
# ╠═ec169d3a-048e-11eb-178a-75dc8b1b3d95
# ╠═5e945f0e-0490-11eb-1745-65666efed48b
# ╠═eeea64e8-0488-11eb-029e-8d0c711f4715
# ╠═e0f9e742-048c-11eb-171f-4f1f3dea0f21
# ╠═a4d84168-048d-11eb-1c8f-0fa3136a1e0f
# ╠═abdf1064-0487-11eb-1f99-eb6563e505d7
# ╠═f06adcc4-0480-11eb-08ba-5dc0b9d3062b
# ╠═82a0f3b4-047f-11eb-2fb0-45fe12555aed
# ╠═8bf3f9bc-0479-11eb-144e-39bcfb215074
# ╠═7f3e45ca-0478-11eb-2c41-4db39b42331e
# ╟─78f39af8-0487-11eb-0948-bf2a888bd088
# ╠═75c284dc-0489-11eb-34ff-d90b6385687a
# ╟─8b7b2098-0486-11eb-2abe-c530be84ba4e
# ╠═3c588554-0482-11eb-2fff-cd4645f903f5
# ╠═d93361c2-0483-11eb-3da2-db2a8ff59b90
# ╠═237fc58e-047c-11eb-0853-05ab5af54793
# ╟─4156a452-0476-11eb-2a93-6f6ba585ab3f
# ╠═7e46e0b6-04c6-11eb-18f7-bf3077b15a72
# ╠═2b54d586-0474-11eb-22bb-73ff01024985
# ╠═37855416-0477-11eb-315c-4b4bca632ae5
# ╠═3ca4f348-0477-11eb-35a9-4ba6addfd1b4
# ╟─4ad9e038-0477-11eb-2363-0f093d35f3ca
# ╠═334533ae-0478-11eb-2c9e-07eedf379844
# ╠═4d93dcfc-0477-11eb-2619-c10f7fd1b909
# ╠═eb708bb0-0476-11eb-0ff3-bd45dc247e5a
# ╠═29638118-0475-11eb-19db-bdcac9172074
# ╠═ff7bab4c-04be-11eb-3de1-27f14eef3b91
# ╠═07c981e6-04bf-11eb-3c02-b9d8d6018719
# ╠═b8eba050-04bf-11eb-16b7-358930a15ffd
# ╠═059cc5ca-0478-11eb-3f29-6b50cb83c78e
# ╟─17d6a8b4-0482-11eb-330b-cd9ed27943bf
# ╟─08c3f4d2-047d-11eb-0b79-e1f9035e3df7
# ╠═5328c0d2-0478-11eb-357c-e750df71d6ff
# ╠═0bae871e-0484-11eb-0b0e-234b276e02f3
# ╠═2a705298-047a-11eb-3967-b3207243a5c3
# ╠═30571030-04c7-11eb-1516-dbd5589c12b9
# ╠═e22f9f90-04c7-11eb-2519-a74e022bcd38
# ╠═f1def888-047b-11eb-1b72-63189bf2650a
# ╠═85b04c2e-047c-11eb-330a-e321c5501b45
# ╟─5496cbaa-0476-11eb-0637-4d3d72087b10
# ╠═dd0b46c2-0472-11eb-08ef-b9b5984cf08c
# ╠═d2a8c018-047e-11eb-2f3f-83b241d73c3d
# ╟─b2c72e50-0484-11eb-1d15-b953a0640a69
# ╟─b27dcd6e-0484-11eb-0ab6-e73b9b4e80ec
# ╟─a3737596-0473-11eb-0902-258a8ea5a515
# ╠═7ae40ec2-0475-11eb-3191-f75c848caaa0
# ╟─38a87a16-0475-11eb-1e64-b99500c57983


text:
# Map-MCMC for mixtures.

# This code assumes the following functions have been defined:
#   likelihood(x[j],theta)
#   log_v_prior(v)
#   log_theta_prior(theta)
#   theta_prop(theta)
#   log_theta_prop(theta,thetap)
#   new_thetas(m)

function sampler(x,n_samples,m,c,sigma,zeta)
    # x = data (array of datapoints)
    # n_samples = # of MCMC iterations
    # m = maximum # of mixture components
    # c = cutoff point for weight map
    # sigma = scale of weight proposals
    # zeta = power to raise likelihood to

    n = length(x)

    # initialize state
    v = zeros(m); v[1] = 2*c; v[2:end] = c/2 # latent weights
    q = max(v-c,0) # mapped weights (unnormalized)
    s = sum(q)
    theta = new_thetas(m)

    # initialize vars for computing likelihood
    L = [likelihood(x[j],theta[i]) for i=1:m, j=1:n]
    M = vec(q'*L)  # mixture density with unnormalized weights
    ll = sum(log(M)) - n*log(s)  # log-lik
    Mp = zeros(n)  # Mp and Lp will hold proposed values
    Lp = zeros(n)

    # record-keeping
    k_r = zeros(n_samples)
    m_r = zeros(n_samples,m)
    l_r = zeros(n_samples,m)
    v_r = zeros(n_samples,m)
    nta = 0  # number of theta proposals accepted
    nva = 0  # number of v proposals accepted

    # draw samples
    for iter = 1:n_samples
        # update parameters with Metropolis-Hastings moves
        for i = 1:m
            thetap = theta_prop(theta[i])
            llp = -n*log(s)
            for j = 1:n
                Lp[j] = likelihood(x[j],thetap)
                Mp[j] = max(M[j] + q[i]*(Lp[j] - L[i,j]), 0)
                # Note: max(.,0) prevents negative values due to roundoff error.
                llp += log(Mp[j])
            end
            # compute acceptance probability
            top = log_theta_prior(thetap) + zeta*llp + log_theta_prop(thetap,theta[i])
            bot = log_theta_prior(theta[i]) + zeta*ll + log_theta_prop(theta[i],thetap)
            p_accept = min(1, exp(top-bot))
            # accept or reject
            if rand() < p_accept
                theta[i] = copy(thetap)
                for j = 1:n; L[i,j] = Lp[j]; end
                M,Mp = Mp,M
                ll = llp
                nta += 1
            end
        end

        # update weights with Metropolis-Hastings moves
        for i = 1:m
            vp = v[i]*exp(randn()*sigma)
            qp = max(vp-c,0)
            sp = s + qp - q[i]
            if sp > 0
                llp = -n*log(sp)
                for j = 1:n
                    Mp[j] = max(M[j] + (qp - q[i])*L[i,j], 0)
                    llp += log(Mp[j])
                end
                # compute acceptance probability
                top = log_v_prior(vp) + zeta*llp - log(v[i])
                bot = log_v_prior(v[i]) + zeta*ll - log(vp)
                p_accept = min(1, exp(top-bot))
                # accept or reject
                if rand() < p_accept
                    v[i] = vp
                    q[i] = qp
                    s = sp
                    M,Mp = Mp,M
                    ll = llp
                    nva += 1
                end
            end
        end

        # record
        k_r[iter] = sum(q.>0)  # number of active components
        for i = 1:m
            m_r[iter,i] = theta[i][1]  # means
            l_r[iter,i] = theta[i][2]  # precisions
            v_r[iter,i] = v[i]  # latent weights
        end

        # visualize mixture components during sampling
        if false
            figure(10); clf(); hold(true)
            xs = linspace(-5,5,1000)
            for i = 1:m
                if q[i]>0
                    plot(xs,(q[i]/s)*normpdf(xs,theta[i][1],theta[i][2]))
                end
            end
            draw()
        end
    end
    art = nta/(m*n_samples)  # MH acceptance rate for theta proposals
    arv = nva/(m*n_samples)  # MH acceptance rate for v proposals
    return (q/s),theta,k_r,v_r,art,arv,m_r,l_r
end

max_stars_repo_name: jwmi/CoarsenedPosterior

max_issues_repo_name: jwmi/CoarsenedPosterior

max_forks_repo_name: jwmi/CoarsenedPosterior


text:
mutable struct PlaneWaveMW{T,P}
  direction::P
  polarisation::P
  wavenumber::T
  amplitude::T
end

function PlaneWaveMW(d,p,k,a = 1)
    T = promote_type(eltype(d), eltype(p), typeof(k), typeof(a))
    P = similar_type(typeof(d), T)
    PlaneWaveMW{T,P}(d,p,k,a)
end

"""
    planewavemw3d(;direction, polarization, wavenumber[, amplitude=1])

Create a plane wave solution to Maxwell's equations.
"""
planewavemw3d(;
    direction    = error("missing arguement `direction`"),
    polarization = error("missing arguement `polarization`"),
    wavenumber   = error("missing arguement `wavenumber`"),
    amplitude    = 1,
    ) = PlaneWaveMW(direction, polarization, wavenumber, amplitude)

function (e::PlaneWaveMW)(x)
  k = e.wavenumber
  d = e.direction
  u = e.polarisation
  a = e.amplitude
  a * exp(-im * k * dot(d, x)) * u
end

function curl(field::PlaneWaveMW)
  k = field.wavenumber
  d = field.direction
  u = field.polarisation
  a = field.amplitude
  v = d × u
  b = -a * im * k
  PlaneWaveMW(d, v, k, b)
end

*(a::Number, e::PlaneWaveMW) = PlaneWaveMW(e.direction, e.polarisation, e.wavenumber, a*e.amplitude)

mutable struct CrossTraceMW{F} <: Functional
  field::F
end

mutable struct TangTraceMW{F} <: Functional
  field::F
end

cross(::NormalVector, p::Function) = CrossTraceMW(p)
cross(::NormalVector, p::PlaneWaveMW) = CrossTraceMW(p)
cross(t::CrossTraceMW, ::NormalVector) = TangTraceMW(t.field)

function (ϕ::CrossTraceMW)(p)
  F = ϕ.field
  x = cartesian(p)
  n = normal(p)
  return n × F(x)
end

function (ϕ::TangTraceMW)(p)
  F = ϕ.field
  x = cartesian(p)
  n = normal(p)
  return (n × F(x)) × n
end

integrand(::TangTraceMW, gx, ϕx) = gx[1] ⋅ ϕx
integrand(::CrossTraceMW, test_vals, field_val) = test_vals[1] ⋅ field_val

struct NDotTrace{F} <: Functional
  field::F
end

(ϕ::NDotTrace)(p) = dot(normal(p), ϕ.field(cartesian(p)))
integrand(::NDotTrace, g, ϕ) = dot(g.value, ϕ)
LinearAlgebra.dot(::NormalVector, f) = NDotTrace(f)

mutable struct CurlGreen{T,U,V}
  wavenumber::T
  source::U
  position::V
end

function (f::CurlGreen)(x)
  γ = im * f.wavenumber
  r = x-f.position
  R = norm(r)
  g = exp(-γ*R)/(4π*R)
  j = f.source
  return -γ/R * (1 + 1/(γ*R)) * g * (r × j)
end

mutable struct CurlCurlGreen{T,U,V}
  wavenumber::T
  source::U
  position::V
end

cross(::NormalVector, p::CurlGreen) = CrossTraceMW(p)
cross(::NormalVector, p::CurlCurlGreen) = CrossTraceMW(p)

function (f::CurlCurlGreen)(x)
  γ = im * f.wavenumber
  r = x - f.position
  R = norm(r)
  g = exp(-γ*R)/(4π*R)
  j = f.source
  return (-γ^2/R^2 * (r × j) × r + (1/R^2 + γ/R)/R^2 * (3r * dot(r,j) - R^2 * j)) * g
end

curl(f::CurlGreen) = CurlCurlGreen(f.wavenumber, f.source, f.position)
function curl(f::CurlCurlGreen)
  κ = f.wavenumber
  j = κ^2 * f.source
  x = f.position
  return CurlGreen(κ, j, x)
end

Base.:*(a::Number, f::CurlGreen) = CurlGreen(f.wavenumber, a*f.source, f.position)
Base.:*(a::Number, f::CurlCurlGreen) = CurlCurlGreen(f.wavenumber, a*f.source, f.position)


text:
D2s = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("C_2^1", [-1.0 -1.2246467991473532e-16 0.0; 1.2246467991473532e-16 -1.0 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2(2)", [1.0 -2.4492935982947064e-16 1.4997597826618576e-32; -2.4492935982947064e-16 -1.0 1.2246467991473532e-16; -1.4997597826618576e-32 -1.2246467991473532e-16 -1.0]),
    Symmetry.Symel("C_2(1)", [1.0 -4.898587196589413e-16 -2.999519565323715e-32; -4.898587196589413e-16 -1.0 -1.2246467991473532e-16; 2.999519565323715e-32 1.2246467991473532e-16 -1.0])]

D3s = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("C_3^1", [-0.4999999999999998 -0.8660254037844387 0.0; 0.8660254037844387 -0.4999999999999998 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_3^2", [-0.5000000000000003 0.8660254037844384 0.0; -0.8660254037844384 -0.5000000000000003 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2(2)", [-0.5000000000000002 -0.8660254037844387 1.0605752387249073e-16; -0.8660254037844387 0.5000000000000009 6.123233995736765e-17; -1.0605752387249073e-16 -6.123233995736765e-17 -1.0]),
    Symmetry.Symel("C_2(3)", [-0.4999999999999991 0.8660254037844394 -1.0605752387249068e-16; 0.8660254037844394 0.4999999999999998 6.123233995736771e-17; 1.0605752387249068e-16 -6.123233995736771e-17 -1.0]),
    Symmetry.Symel("C_2(1)", [1.0 -1.2212453270876724e-15 -7.477970903957902e-32; -1.2212453270876724e-15 -1.0 -1.2246467991473532e-16; 7.477970903957902e-32 1.2246467991473532e-16 -1.0])]

D4s = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("C_4^1", [6.123233995736766e-17 -1.0 0.0; 1.0 6.123233995736766e-17 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2^1", [-1.0 -1.2246467991473532e-16 0.0; 1.2246467991473532e-16 -1.0 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_4^3", [-1.8369701987210297e-16 1.0 0.0; -1.0 -1.8369701987210297e-16 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2(2)", [-1.0 1.2246467991473532e-16 1.2246467991473532e-16; 1.2246467991473532e-16 1.0 -7.498798913309288e-33; -1.2246467991473532e-16 7.498798913309288e-33 -1.0]),
    Symmetry.Symel("C_2(3)", [1.0 -2.4492935982947064e-16 1.4997597826618576e-32; -2.4492935982947064e-16 -1.0 1.2246467991473532e-16; -1.4997597826618576e-32 -1.2246467991473532e-16 -1.0]),
    Symmetry.Symel("C_2(4)", [-1.0 3.6739403974420594e-16 -1.2246467991473532e-16; 3.6739403974420594e-16 1.0 2.2496396739927864e-32; 1.2246467991473532e-16 -2.2496396739927864e-32 -1.0]),
    Symmetry.Symel("C_2(1)", [1.0 -4.898587196589413e-16 -2.999519565323715e-32; -4.898587196589413e-16 -1.0 -1.2246467991473532e-16; 2.999519565323715e-32 1.2246467991473532e-16 -1.0])]

D5s = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("C_5^1", [0.30901699437494745 -0.9510565162951535 0.0; 0.9510565162951535 0.30901699437494745 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_5^2", [-0.8090169943749473 -0.5877852522924731 0.0; 0.5877852522924731 -0.8090169943749473 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_5^3", [-0.8090169943749475 0.587785252292473 0.0; -0.587785252292473 -0.8090169943749475 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_5^4", [0.30901699437494734 0.9510565162951535 0.0; -0.9510565162951535 0.30901699437494734 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2(2)", [-0.8090169943749473 0.5877852522924734 1.1647083184890926e-16; 0.5877852522924734 0.8090169943749481 -3.784366730434151e-17; -1.1647083184890926e-16 3.784366730434151e-17 -1.0]),
    Symmetry.Symel("C_2(3)", [0.3090169943749479 -0.951056516295154 7.198293278059968e-17; -0.951056516295154 -0.3090169943749471 9.907600726170918e-17; -7.198293278059968e-17 -9.907600726170918e-17 -1.0]),
    Symmetry.Symel("C_2(4)", [0.30901699437494834 0.9510565162951539 -7.198293278059966e-17; 0.9510565162951539 -0.30901699437494745 9.90760072617092e-17; 7.198293278059966e-17 -9.90760072617092e-17 -1.0]),
    Symmetry.Symel("C_2(5)", [-0.8090169943749475 -0.5877852522924731 -1.1647083184890926e-16; -0.5877852522924731 0.8090169943749481 -3.78436673043415e-17; 1.1647083184890926e-16 3.78436673043415e-17 -1.0]),
    Symmetry.Symel("C_2(1)", [1.0 -2.2204460492503136e-16 -1.3596310734468913e-32; -2.2204460492503136e-16 -1.0 -1.2246467991473532e-16; 1.3596310734468913e-32 1.2246467991473532e-16 -1.0])]

D6s = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("C_6^1", [0.5000000000000001 -0.8660254037844386 0.0; 0.8660254037844386 0.5000000000000001 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_3^1", [-0.4999999999999998 -0.8660254037844388 0.0; 0.8660254037844388 -0.4999999999999998 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2^1", [-1.0 -3.885780586188048e-16 0.0; 3.885780586188048e-16 -1.0 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_3^2", [-0.5000000000000006 0.8660254037844385 0.0; -0.8660254037844385 -0.5000000000000006 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_6^5", [0.49999999999999944 0.8660254037844392 0.0; -0.8660254037844392 0.49999999999999944 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2(2)", [-0.4999999999999998 0.8660254037844388 1.0605752387249068e-16; 0.8660254037844388 0.4999999999999998 -6.123233995736767e-17; -1.0605752387249068e-16 6.123233995736767e-17 -1.0]),
    Symmetry.Symel("C_2(3)", [-0.5000000000000004 -0.8660254037844385 1.0605752387249072e-16; -0.8660254037844385 0.5000000000000007 6.123233995736764e-17; -1.0605752387249072e-16 -6.123233995736764e-17 -1.0]),
    Symmetry.Symel("C_2(4)", [1.0 -7.771561172376096e-16 4.7587087570641186e-32; -7.771561172376096e-16 -1.0 1.2246467991473532e-16; -4.7587087570641186e-32 -1.2246467991473532e-16 -1.0]),
    Symmetry.Symel("C_2(5)", [-0.4999999999999989 0.8660254037844395 -1.0605752387249067e-16; 0.8660254037844395 0.49999999999999933 6.123233995736773e-17; 1.0605752387249067e-16 -6.123233995736773e-17 -1.0]),
    Symmetry.Symel("C_2(6)", [-0.5000000000000013 -0.8660254037844378 -1.0605752387249073e-16; -0.8660254037844378 0.5000000000000009 -6.123233995736757e-17; 1.0605752387249073e-16 6.123233995736757e-17 -1.0]),
    Symmetry.Symel("C_2(1)", [1.0 -1.6653345369377344e-15 -1.019723305085168e-31; -1.6653345369377344e-15 -1.0 -1.2246467991473532e-16; 1.019723305085168e-31 1.2246467991473532e-16 -1.0])]

D2cn = ["A", "B1", "B2", "B3"]
D2ct = C2vct
D3cn = ["A1", "A2", "E"]
D3ct = C3vct
D4cn = ["A1", "A2", "B1", "B2", "E"]
D4ct = C4vct
D5cn = ["A1", "A2", "E1", "E2"]
D5ct = C5vct
D6cn = ["A1", "A2", "B1", "B2", "E1", "E2"]
D6ct = C6vct


text:
#=
Macro for taking normal Julia-like functions of Pseudo Assembly and emitting
a vector of AsapInstructions.

This will probably be a multi-pass operation, where first the expression tree
of the function is converted into AsapInstructions with reckless abandon, and
then several processing steps to check if any assembly instructions were
provided incorrectly, if the source/destinations are invalid, etc.
=#

macro asap4asm(expr)
    return esc(convertcode(expr))
end

# Convenience functions for converting symbols to either all caps or all
# lowercase.
symbol_uc(x::Symbol) = Symbol(uppercase(string(x)))
symbol_lc(x::Symbol) = Symbol(lowercase(string(x)))

# Need to make this mutable since we will potentially modifying and changing
# some instructions before separating the instructions and labels.
mutable struct InstructionLabelTarget
    instruction   :: AsapInstruction
    label         :: Union{Symbol,Nothing}
    branch_target :: Union{Symbol,Nothing}
end

function convertcode(expr::Expr)
    # This macro does the following:
    # 1. Extracts the body of the passed function using MacroTools' excellent
    #    "splitdef" function.
    # 2. Drop in an empty Vector{InstructionLabelTarget} at the beginning of
    #    the function.
    # 3. Walk through the body of the function in order, recording the labels
    #    it sees and converting pseudo-assembly function calls into
    #    InstructionLabelTarget and appending them to the vector.
    # 4. Replace each return by passing the Vector{InstructionLabelTarget} with
    #    a call to an "assembler" function that will analyze the code to resolve
    #    labels and RPT instructions before returning a final program as a
    #    Vector{AsapInstruction}

    # 1 - split function.
    split_expr = splitdef(expr)
    function_body = split_expr[:body]

    # 2 - generate a symbol for the intermediate program and insert it into the
    # beginning of the function body.
    program_vector_sym = gensym("program")
    startvector!(function_body, program_vector_sym)

    # 3 - Replace all instructions with "push!" to the program_vector.
    function_body = replace_instructions(function_body, program_vector_sym)

    # 4 - Replace "return" statements with calls to the assembler.
    # Insert a return statement at the end as well as a default catch.
    function_body = replace_returns(function_body, program_vector_sym)

    # Rebuild the function and return the modified source code.
    split_expr[:body] = function_body
    return MacroTools.combinedef(split_expr)
end

# Put the expression
#
# program = InstructionLabelTarget[]
#
# At the beginning of the function body - expept that "program" is some
# complicated thing given by "gensym" so we don't have to worry about it
# clobbering anything in the function.
function startvector!(function_body :: Expr, program_vector_sym)
    pushfirst!(function_body.args, :($program_vector_sym = InstructionLabelTarget[]))
end

# Container for holding kwargs to give to the Keyword constructor of
# AsapInstruction for a given instruction.
#
# Since I changed from using Symbols to encode Source/Destination locations
# to an Enum, the branch_target can no longer be stored in the instruction
# itself. Thus, we have to record branch targets as a separate field in this
# container.
mutable struct KwargHolder
    kwargs          :: Vector{Pair{Symbol,Any}}
    branch_target   :: Union{Symbol,Nothing}
end
KwargHolder() = KwargHolder(Pair{Symbol,Any}[], nothing)

function replace_instructions(function_body :: Expr, program_vector_sym)
    # Record the last seen label. Attach this to all found assembly instructions.
    #
    # MacroTools.postwalk will encounter expressions in program order, so it is
    # sufficient to keep track of the most recently seen label.
    #
    # When assembling, the label will be resolved to the first instruction that
    # has that label.
    label :: Union{Symbol,Nothing} = nothing

    # Walk through the expression - replacing pseudo-assembly with parsed
    # instruction types.
    new_body = MacroTools.postwalk(function_body) do ex
        if @capture(ex, op_(args__))
            # Check if the name of the function being called is a Assembly
            # instruction. If not - leave it alone.
            if isop(op)
                # Parse the arguments for this operation. Expect to recieve a
                # vector of Pair{Symbol,Any} that represent the arguments to
                # be passed to an AsapInstructionKeyword constructor.
                holder = getkwargs(op, args)

                # Check if branch target is in here. If so, extract it and
                # record it.
                kwargs = [Expr(:kw, p[1], p[2]) for p in holder.kwargs]
                branch_target = holder.branch_target

                # Create an expression pushing the pair of this label and
                # a AsapInstruction to the program vector.
                new_ex = :(push!(
                     $program_vector_sym,
                     InstructionLabelTarget(
                        AsapInstructionKeyword(;$(kwargs...)),
                        $(Meta.quot(label)),
                        $(Meta.quot(branch_target)),
                     )))

                # replace the expression
                return new_ex
            end

        # Record macro calls
        # The item returned by the @capture will be either a QuoteNode or
        # a Symbol. For consistency sake, always turn it into a Symbol.
        #
        # It seems that code captured from the REPL will be lowered to a Symbol
        # while code lowered from a file yields a QuoteNode. Not entirely sure
        # what is going on there.
        elseif @capture(ex, @label label_quotenode_or_symbol_)
            label = Symbol(label_quotenode_or_symbol)
            # Replace the label with nothing - removing it from the code.
            return nothing
        end

        # Default fallback - return expression unchanged.
        return ex
    end

    return new_body
end

# Helper macro.
macro pairpush!(kwargs, syms...)
    symmaps = [:($(Meta.quot(sym)) => $sym) for sym in syms]
    return esc(:(push!($kwargs, $(symmaps...))))
end

function getkwargs(op, args)
    # Need to get a collection of keyword arguments to give to the caller for
    # this operation and arguments.
    holder = KwargHolder()
    push!(holder.kwargs, :op => :(AsapSim.$op))

    # Get the InstructionDefinition for this operation.
    def = getdef(op)

    # Hendle special instrucitons first.
    if op == :BR || op == :BRL
        # Get the destination label
        destination_label = first(args)

        # Add the branch target to the front of the vector. Upper level routine
        # will extract it.
        holder.branch_target = destination_label

        # Record if this is a return.
        if destination_label == :back
            push!(holder.kwargs, :isreturn => true)
        end

        # Remove the destination argument from the vector to avoid using it
        # more than once.
        popfirst!(args)

        # (TODO) - this is temporary. Need to really first perform a check to
        # see if the caller is providing branch options. Otherwise, need to be
        # a little more precise with how this is handled
        push!(holder.kwargs, :src1 => Loc(IMMEDIATE, make_branch_mask(args)))

        # Get options for this instruction.
        append!(holder.kwargs, getoptions(args))

    # General instruction decoding.
    else
        # If this instruction definition has a destination - extract it and
        # remove the destination from the argument list.
        if def.dest
            dest = extract_loc(first(args))
            @pairpush! holder.kwargs dest
            popfirst!(args)
        end

        # Parse out src1 and src2
        if def.src1
            src1 = extract_loc(first(args))
            @pairpush! holder.kwargs src1
            popfirst!(args)
        end
        if def.src2
            src2 = extract_loc(first(args))
            @pairpush! holder.kwargs src2
        end

        # Get options for this instruction.
        append!(holder.kwargs, getoptions(args))
    end

    # Mark the op type and special attributes.
    push!(holder.kwargs, :optype => def.optype)
    if def.signed
        push!(holder.kwargs, :signed => true)
    end
    if def.saturate
        push!(holder.kwargs, :saturate => true)
    end
    return holder
end

extract_loc(val::Integer) = :(Loc(AsapSim.IMMEDIATE, $val))
function extract_loc(sym::Symbol)
    # Convert the symbol to all uppercase letters for matching with the
    # enumerated source/destination locations.
    sym_upper = symbol_uc(sym)

    # If this symbol is a valid source or destination keyword, return the
    # constructor for a Loc based on the SrcDest enum.
    #
    # Prefix it with "AsapSim" to avoid scoping issues.
    if sym_upper in SourceDestSymbols
        return :(Loc(AsapSim.$sym_upper))

    # Otherwise, check if this is a common keyword. If so, pass the symbol
    # directly to the constructor
    #
    # i.e., extract_loc(:ibuf0) == :(Loc(:ibuf0))
    elseif sym in keys(SourceAliases)
        return :(Loc($(Meta.quot(sym))))

    # Otherwise, assume that this is some call-time construct and don't modify
    # the argument at all.
    else
        return :(Loc($sym))
    end
end

function extract_loc(expr::Expr)
    # Check if this expression is a :ref expression (like "dmem[0]"), if so
    # return the two args (i.e. :dmem and 0)
    if @capture(expr, sym_[ind_])
        sym_upper = symbol_uc(sym)
        if sym_upper in SourceDestSymbols
            return :(Loc(AsapSim.$sym_upper, $ind))
        else
            return :(Loc($expr))
        end

    # Otherwise, assume this expression is some computation for an immediate
    # and pass the whole expression.
    else
        return :(Loc($expr))
    end
end

function replace_returns(function_body, program_vector_sym; return_fn = :assemble)
    # Create the expression to replace the returns with.
    new_ex = :(return $return_fn($program_vector_sym))

    # Add an empty return at the end of the function body.
    # This will get replaced with the above expression during the pass below.
    push!(function_body.args, :(return))

    # Postwalk again, replace returns.
    new_body = MacroTools.postwalk(function_body) do ex
        if @capture(ex, return args__)
            # TODO: Maybe give error if there are arguments other than "nothing"
            return new_ex
        end

        # Default return
        return ex
    end

    return new_body
end

# Function for constructing the immediate mask for the branch instruction.
function make_branch_mask(args)
    # Iterate through the arguments in the order the bits are set for convenience.
    possible_args = (:N, :Z, :C, :O, :F0, :F1, :PI, :PO, :OB, :MI, :MO, :U, :neg)

    mask = 0

    for (bit, arg) in enumerate(possible_args)
        indices = findall(x -> x == arg, args)
        # Check if this argument exists. If so, set the corresponding bit of the
        # mask.
        if length(indices) > 0
            mask |= 1 << (bit - 1)
            # Delete this entry in the args.
            deleteat!(args, indices)
        end
    end

    # Check if the mask is still zero, if so, set the unconditional bit of the
    # mask.
    if mask == 0
        mask |= 1 << 11
    end
    return mask
end

function make_stall_mask(args)
    # Similar to branch mask.
    possible_args = (
        :ibuf0,
        :ibuf1,
        :packet_in,
        :memory_in,
        :obuf_full,
        :packet_out,
        :memory_out,
    )

    for (bit, arg) in enumerate(possible_args)
        indices = find(x -> x == arg, args)
        if length(indices) > 0
            mask |= 1 << (bit - 1)
            deleteat!(args, indices)
        end
    end

    # No default bit needs to be set for the stall mask.
    return mask
end


text:
@testset "190.reverse-bits.jl" begin
    @test reverse_bit(0b00000010100101000001111010011100) ==
          0b00111001011110000010100101000000
    @test reverse_bit(0b11111111111111111111111111111101) ==
          0b10111111111111111111111111111111
    @test reverse_bit(0b00000010100101000001111010011100) ==
          0b00111001011110000010100101000000
    @test reverse_bit(0b11111111111111111111111111111101) ==
          0b10111111111111111111111111111111
end


text:
"""
     module JCMolecule
The module required for determination of molecular coordinate-based properties
(such as bond lengths, bond angles, and dihedral angles). Import this module
into the script when you wish to determine such information. Note that this
module is not strictly necessary for every calculation.
"""
module JCMolecule

Base.include(@__MODULE__,"MoleculeAnalysis.jl")

"""
     run(input_info::Dict{String,Dict{String,Any}})
Execute the JuliaChem molecular coordinate analysis functions.

One input variable is required:
1. input_info = Information gathered from the input file.

No variables are output.

Thus, proper use of the Molecule.run() function would look like this:

```
Molecule.run(input_info)
```
"""
function run(mol::Molecule)
    comm=MPI.COMM_WORLD

    if (MPI.Comm_rank(comm) == 0)
        println("---------------------------------------------------------------------------------")
        println("                       ========================================          ")
        println("                             MOLECULAR COORDINATE ANALYSIS               ")
        println("                       ========================================          ")
        println("")
    end

    #== print coordinates ==#
    print_xyz(mol)

    #== compute bond lengths ==#
    #bond_lengths = analyze_bond_lengths(mol)

    #== compute bond angles ==#
    #bond_angles = analyze_bond_angles(mol,bond_lengths)

    if (MPI.Comm_rank(comm) == 0)
        println("                       ========================================          ")
        println("                                END COORDINATE ANALYSIS                  ")
        println("                       ========================================          ")
    end
end
export run

end


text:
include("ModuleA.jl")

import MyModule: foo, MyType
#equivalent to
import MyModule.foo, MyModule.MyType

mytype = MyType(2)
@show foo(mytype)
@show MyType

@show MyModule.foo(mytype)
@show MyModule.MyType
@show MyModule.bar(2)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
#
# wrappedarrays.jl --
#
# Management of object wrapped into Julia arrays.
#
#------------------------------------------------------------------------------
#
# This file is part of InterProcessCommunication.jl released under the MIT
# "expat" license.
#
# Copyright (C) 2016-2019, Éric Thiébaut
# (https://github.com/emmt/InterProcessCommunication.jl).
#

"""
```julia
WrappedArray(mem, [T [, dims...]]; offset=0)
```

yields a Julia array whose elements are stored in the "memory" object `mem`.
Argument `T` is the data type of the elements of the returned array and
argument(s) `dims` specify the dimensions of the array.  If `dims` is omitted
the result is a vector of maximal length (accounting for the offset and the
size of the `mem` object).  If `T` is omitted, `UInt8` is assumed.

Keyword `offset` may be used to specify the address (in bytes) relative to
`pointer(mem)` where is stored the first element of the array.

The size of the memory provided by `mem` must be sufficient to store all
elements (accounting for the offset) and the alignment of the elements in
memory must be a multiple of `Base.datatype_alignment(T)`.

Another possibility is:

```julia
WrappedArray(mem, dec)
```

where `mem` is the "memory" object and `dec` is a function in charge of
decoding the array type and layout given the memory object.  The decoder is
applied to the memory object as follow:

```julia
dec(mem) -> T, dims, offset
```

which must yield the data type `T` of the array elements, the dimensions `dims`
of the array and the offset of the first element relative to `pointer(mem)`.

## Restrictions

The `mem` object must extend the methods `pointer(mem)` and `sizeof(mem)` which
must respectively yield the base address of the memory provided by `mem` and the
number of available bytes.  Furthermore, this memory is assumed to be available
at least until object `mem` is reclaimed by the garbage collector.

## Shared Memory Arrays

```julia
WrappedArray(id, T, dims; perms=0o600, volatile=true)
```

creates a new wrapped array whose elements (and a header) are stored in shared
memory identified by `id` (see [`SharedMemory`](@ref) for a description of `id`
and for keywords).  To retrieve this array in another process, just do:

```julia
WrappedArray(id; readonly=false)
```

## See Also

[`SharedMemory`](@ref).

"""
function WrappedArray(mem::M, ::Type{T} = UInt8;
                      offset::Integer = 0) where {M,T}
    ptr, siz = _check_wrapped_array_arguments(mem, T, offset)
    siz ≥ sizeof(T) ||
        throw_argument_error("insufficient memory for at least one element")
    number = div(siz, sizeof(T))
    return WrappedArray{T,1,M}(ptr, (number,), mem)
end

WrappedArray(mem::M, ::Type{T}, dims::Integer...; kwds...) where {M,T} =
    WrappedArray(mem, T, dims; kwds...)

function WrappedArray(mem::M, ::Type{T}, dims::NTuple{N,Integer};
                      offset::Integer = 0) where {T,N,M}
    ptr, siz = _check_wrapped_array_arguments(mem, T, offset)
    number = checkdims(dims)
    siz ≥ sizeof(T)*number ||
        throw_argument_error("insufficient memory for array")
    return WrappedArray{T,N,M}(ptr, dims, mem)
end

function WrappedArray(mem, dec::Function)
    T, dims, offset = dec(mem)
    isa(T, DataType) || error("`dec(mem)[1]` must be a data type")
    isa(dims, Tuple{Vararg{Integer}}) ||
        error("`dec(mem)[2]` must be a tuple of dimensions")
    isa(offset, Integer) || error("`dec(mem)[3]` must be an integer")
    return WrappedArray(mem, T, dims; offset = offset)
end

function WrappedArray(id::Union{AbstractString,ShmId,Key},
                      ::Type{T}, dims::Vararg{Integer,N}; kwds...) where {T,N}
    return WrappedArray(id, T, convert(NTuple{N,Int}, dims); kwds...)
end

function WrappedArray(id::Union{AbstractString,ShmId,Key},
                      ::Type{T}, dims::NTuple{N,Integer}; kwds...) where {T,N}
    return WrappedArray(id, T, convert(NTuple{N,Int}, dims); kwds...)
end

function WrappedArray(id::Union{AbstractString,ShmId,Key},
                      ::Type{T}, dims::NTuple{N,Int};
                      kwds...) where {T,N}
    num = checkdims(dims)
    off = _wrapped_array_header_size(N)
    siz = off + sizeof(T)*num
    mem = SharedMemory(id, siz; kwds...)
    write(mem, WrappedArrayHeader, T, dims)
    return WrappedArray(mem, T, dims; offset=off)
end

function WrappedArray(id::Union{AbstractString,ShmId,Key}; kwds...)
    mem = SharedMemory(id; kwds...)
    T, dims, off = read(mem, WrappedArrayHeader)
    return WrappedArray(mem, T, dims; offset=off)
end

function _check_wrapped_array_arguments(mem::M, ::Type{T},
                                        offset::Integer) where {M,T}
    offset ≥ 0 || throw_argument_error("offset must be nonnegative")
    isbitstype(T) || throw_argument_error("illegal element type (", T, ")")
    ptr, len = get_memory_parameters(mem)
    align = Base.datatype_alignment(T)
    addr = ptr + offset
    rem(convert(Int, addr), align) == 0 ||
        throw_argument_error("base address must be a multiple of ",
                             align, " bytes")
    return (convert(Ptr{T}, addr),
            convert(Int, len) - convert(Int, offset))
end

# FIXME: push!, pop!, append!, resize! cannot be extended for WrappedVectors
# unless it is possible to query the size of the memory object, in fact many
# things are doable if the address and size of memory object can be retrieved.
# However, push!, append!, resize!, ... would require to rewrap the buffer.

# The following methods come for free (with no performance penalties) because a
# WrappedArray is a subtype of DenseArray:
#
#    Base.eltype, Base.elsize, Base.ndims, Base.first, Base.endof,
#    Base.eachindex, ...

# FIXME: extend Base.view?

@inline Base.parent(obj::WrappedArray) = obj.arr

Base.length(obj::WrappedArray) = length(parent(obj))
Base.sizeof(obj::WrappedArray) = sizeof(parent(obj))
Base.size(obj::WrappedArray) = size(parent(obj))
Base.size(obj::WrappedArray, d::Integer) = size(parent(obj), d)
Base.axes(obj::WrappedArray) = axes(parent(obj))
Base.axes(obj::WrappedArray, d::Integer) = axes(parent(obj), d)
@inline Base.axes1(obj::WrappedArray) = Base.axes1(parent(obj))
Base.strides(obj::WrappedArray) = strides(parent(obj))
Base.stride(obj::WrappedArray, d::Integer) = stride(parent(obj), d)
Base.elsize(::Type{WrappedArray{T,N,M}}) where {T,N,M} = elsize(Array{T,N})

Base.IndexStyle(::Type{<:WrappedArray}) = Base.IndexLinear()

@inline Base.getindex(A::WrappedArray, i::Int) = begin
    @boundscheck checkbounds(A, i)
    @inbounds val = parent(A)[i]
    val
end

@inline Base.setindex!(A::WrappedArray, val, i::Int) = begin
    @boundscheck checkbounds(A, i)
    @inbounds parent(A)[i] = val
    A
end

@inline Base.checkbounds(::Type{Bool}, A::WrappedArray, i::Int) =
    (i % UInt) - 1 < length(A)

Base.copy(obj::WrappedArray) = copy(parent(obj))

Base.copyto!(dest::WrappedArray, src::AbstractArray) =
    (copyto!(dest.arr, src); dest)

Base.reinterpret(::Type{T}, obj::WrappedArray) where {T} =
    reinterpret(T, parent(obj))

Base.reshape(obj::WrappedArray, dims::Tuple{Vararg{Int}}) =
    reshape(parent(obj), dims)

# Extend `Base.unsafe_convert` for `ccall`.  Note that this also make `pointer`
# applicable and that the 2 following definitions are needed to avoid
# ambiguities and cover all cases.
unsafe_convert(::Type{Ptr{T}}, obj::WrappedArray{T}) where {T} =
    unsafe_convert(Ptr{T}, parent(obj))
unsafe_convert(::Type{Ptr{S}}, obj::WrappedArray{T}) where {S,T} =
    unsafe_convert(Ptr{S}, parent(obj))

# Make a wrapped array iterable:
@static if isdefined(Base, :iterate) # VERSION ≥ v"0.7-alpha"
    Base.iterate(iter::WrappedArray) = iterate(iter.arr)
    Base.iterate(iter::WrappedArray, state) = iterate(iter.arr, state)
    Base.IteratorSize(iter::WrappedArray) = Base.IteratorSize(iter.arr)
    Base.IteratorEltype(iter::WrappedArray) = Base.IteratorEltype(iter.arr)
else
    Base.start(iter::WrappedArray) = start(iter.arr)
    Base.next(iter::WrappedArray, state) = next(iter.arr, state)
    Base.done(iter::WrappedArray, state) = done(iter.arr, state)
    Base.iteratorsize(iter::WrappedArray) = Base.iteratorsize(iter.arr)
    Base.iteratoreltype(iter::WrappedArray) = Base.iteratoreltype(iter.arr)
end

#------------------------------------------------------------------------------
# WRAPPED ARRAYS WITH HEADER

# Magic number
const _WA_MAGIC = UInt32(0x57412D31) # "WA-1" = Wrapped Array version 1

# We require at least 64 bytes (512 bits) alignment for the first element of
# the array (to warrant that SIMD vectors are correctly aligned), this value is
# equal to the constant JL_CACHE_BYTE_ALIGNMENT used for the elements of Julia
# arrays.
const _WA_ALIGN = 64

const _WA_TYPES = (( 1, Int8,       "signed 8-bit integer"),
                   ( 2, UInt8,      "unsigned 8-bit integer"),
                   ( 3, Int16,      "signed 16-bit integer"),
                   ( 4, UInt16,     "unsigned 16-bit integer"),
                   ( 5, Int32,      "signed 32-bit integer"),
                   ( 6, UInt32,     "unsigned 32-bit integer"),
                   ( 7, Int64,      "signed 64-bit integer"),
                   ( 8, UInt64,     "unsigned 64-bit integer"),
                   ( 9, Float32,    "32-bit floating-point"),
                   (10, Float64,    "64-bit floating-point"),
                   (11, ComplexF32, "64-bit complex"),
                   (12, ComplexF64, "128-bit complex"))

const _WA_ETYPES = DataType[T for (i, T, str) in _WA_TYPES]
const _WA_DESCRS = String[str for (i, T, str) in _WA_TYPES]
const _WA_IDENTS = Dict(T => i for (i, T, str) in _WA_TYPES)

"""
```julia
WrappedArrayHeader(T, N)
```

yields a structure `WrappedArrayHeader` instanciated for an array with `N`
dimensions and whose element type is `T`.

## Possible Usages

```julia
read(src, WrappedArrayHeader) -> T, dims, off
```

reads wrapped array header in `src` and yields the element type `T`, dimensions
`dims` and offset `off` of the first array element relative to the the base of
`src`.

```julia
write(dst, WrappedArrayHeader, T, dims)
```

writes wrapped array header in `dst` for element type `T` and dimensions
`dims`.

```julia
WrappedArray(mem, x -> read(x, WrappedArrayHeader)) -> arr
```

retrieves a wrapped array `arr` whose header and elements are stored in the
memory provided by object `mem`.

"""
function WrappedArrayHeader(::Type{T}, N::Int) where {T}
    N ≥ 1 || throw_argument_error("illegal number of dimensions (", N, ")")
    haskey(_WA_IDENTS, T) || throw_argument_error("unsupported data type (",
                                                  T, ")")
    off = _wrapped_array_header_size(N)
    return WrappedArrayHeader(_WA_MAGIC, _WA_IDENTS[T], N, off)
end

_wrapped_array_header_size(N::Integer) =
    roundup(sizeof(WrappedArrayHeader) + N*sizeof(Int64), _WA_ALIGN)

WrappedArrayHeader(::Type{T}, N::Integer) where {T} =
    WrappedArrayHeader(T, convert(Int, N))

function Base.write(mem, ::Type{WrappedArrayHeader},
                    ::Type{T}, dims::Vararg{<:Integer,N}) where {T,N}
    write(mem, WrappedArrayHeader, T, convert(NTuple{N,Int}, dims))
end

function Base.write(mem, ::Type{WrappedArrayHeader},
                    ::Type{T}, dims::NTuple{N,<:Integer}) where {T,N}
    write(mem, WrappedArrayHeader, T, convert(NTuple{N,Int}, dims))
end

function Base.write(mem, ::Type{WrappedArrayHeader},
                    ::Type{T}, dims::NTuple{N,Int}) where {T,N}
    # Check arguments, then write header and dimensions.
    hdr = WrappedArrayHeader(T, N)
    off = hdr.offset
    ptr, siz = get_memory_parameters(mem)
    num = checkdims(dims)
    siz ≥ off + sizeof(T)*num ||
        throw_argument_error("insufficient size of memory block")
    unsafe_store!(convert(Ptr{WrappedArrayHeader}, ptr), hdr)
    addr = convert(Ptr{Int64}, ptr + sizeof(WrappedArrayHeader))
    for i in 1:N
        unsafe_store!(addr, dims[i], i)
    end
end

function Base.read(mem, ::Type{WrappedArrayHeader})
    ptr, siz = get_memory_parameters(mem)
    siz ≥ sizeof(WrappedArrayHeader) ||
        throw_argument_error("insufficient size of memory block for header")
    hdr = unsafe_load(convert(Ptr{WrappedArrayHeader}, ptr))
    hdr.magic == _WA_MAGIC ||
        throw_error_exception("invalid magic number (0x",
                              string(hdr.magic, base=16), ")")
    etype = Int(hdr.etype)
    1 ≤ etype ≤ length(_WA_ETYPES) ||
        throw_error_exception("invalid element type identifier (", etype, ")")
    ndims = Int(hdr.ndims)
    1 ≤ ndims ||
        throw_error_exception("invalid number of dimensions (", ndims, ")")
    off = Int(hdr.offset)
    off == _wrapped_array_header_size(ndims) ||
        throw_error_exception("invalid offset (", off, ")")
    addr = convert(Ptr{Int64}, ptr + sizeof(WrappedArrayHeader))
    dims = ntuple(i -> convert(Int, unsafe_load(addr, i)), ndims)
    num = checkdims(dims)
    T = _WA_ETYPES[etype]
    siz ≥ off + sizeof(T)*num ||
        throw_error_exception("insufficient size of memory block (", Int(siz),
                              " < ", Int(off + sizeof(T)*num), ")")
    return T, dims, off
end


text:
module TripletEmbeddings

using Printf
using Random
using Logging
using Distances
using StatsFuns
using Statistics: mean
using SparseArrays
using LinearAlgebra
using InvertedIndices

export Embedding, ndims, nitems, V, Gram, L,
       Triplet, Triplets, checktriplets, sampletriplet, sampletriplets,
       LabeledTriplet, LabeledTriplets,
       STE, ρSTE, tSTE, kernel,
       fit!, gradient, tgradient, tgradient!, cost, tcost,
       procrustes, procrustes!, mse, apply,
       show

import Base: size, getindex, length, iterate, ndims, broadcast, -, *, show
import Distances: pairwise

include("embedding.jl")
include("triplets.jl")
include("utilities.jl")
include("losses.jl")
include("STE.jl")
include("rhoSTE.jl")
include("tSTE.jl")
include("fit.jl")
include("procrustes.jl")

end # module


text:
using Markdown

"""
    TerminalLogger(stream=stderr, min_level=$ProgressLevel; meta_formatter=default_metafmt,
                   show_limited=true, right_justify=0)

Logger with formatting optimized for interactive readability in a text console
(for example, the Julia REPL). This is an enhanced version of the terminal
logger `Logging.ConsoleLogger` which comes installed with Julia by default.

Log levels less than `min_level` are filtered out.

Message formatting can be controlled by setting keyword arguments:

* `meta_formatter` is a function which takes the log event metadata
  `(level, _module, group, id, file, line)` and returns a color (as would be
  passed to printstyled), prefix and suffix for the log message.  The
  default is to prefix with the log level and a suffix containing the module,
  file and line location.
* `show_limited` limits the printing of large data structures to something
  which can fit on the screen by setting the `:limit` `IOContext` key during
  formatting.
* `right_justify` is the integer column which log metadata is right justified
  at. The default is zero (metadata goes on its own line).
"""
struct TerminalLogger <: AbstractLogger
    stream::IO
    min_level::LogLevel
    meta_formatter
    show_limited::Bool
    right_justify::Int
    message_limits::Dict{Any,Int}
    sticky_messages::StickyMessages
    bartrees::Vector{Node{ProgressBar}}
end
function TerminalLogger(stream::IO=stderr, min_level=ProgressLevel;
                        meta_formatter=default_metafmt, show_limited=true,
                        right_justify=0)
    TerminalLogger(
        stream,
        min_level,
        meta_formatter,
        show_limited,
        right_justify,
        Dict{Any,Int}(),
        StickyMessages(stream),
        Union{}[],
    )
end

shouldlog(logger::TerminalLogger, level, _module, group, id) =
    get(logger.message_limits, id, 1) > 0

min_enabled_level(logger::TerminalLogger) = logger.min_level

function default_logcolor(level)
    level < Info  ? :blue :
    level < Warn  ? Base.info_color()  :
    level < Error ? Base.warn_color()  :
                    Base.error_color()
end

function default_metafmt(level, _module, group, id, file, line)
    color = default_logcolor(level)
    prefix = (level == Warn ? "Warning" : string(level))*':'
    suffix = ""
    Info <= level < Warn && return color, prefix, suffix
    _module !== nothing && (suffix *= "$(_module)")
    if file !== nothing
        _module !== nothing && (suffix *= " ")
        suffix *= file
        if line !== nothing
            suffix *= ":$(isa(line, UnitRange) ? "$(first(line))-$(last(line))" : line)"
        end
    end
    !isempty(suffix) && (suffix = "@ " * suffix)
    return color, prefix, suffix
end

# Length of a string as it will appear in the terminal (after ANSI color codes
# are removed)
function termlength(str)
    N = 0
    in_esc = false
    for c in str
        if in_esc
            if c == 'm'
                in_esc = false
            end
        else
            if c == '\e'
                in_esc = true
            else
                N += 1
            end
        end
    end
    return N
end

function format_message(message, prefix_width, io_context)
    formatted = sprint(show, MIME"text/plain"(), message, context=io_context)
    msglines = split(chomp(formatted), '\n')
    if length(msglines) > 1
        # For multi-line messages it's possible that the message was carefully
        # formatted with vertical alignemnt. Therefore we play it safe by
        # prepending a blank line.
        pushfirst!(msglines, SubString(""))
    end
    msglines
end

function format_message(message::AbstractString, prefix_width, io_context)
    # For strings, use Markdown to do the formatting. The markdown renderer
    # isn't very composable with other text formatting so this is quite hacky.
    message = Markdown.parse(message)
    prepend_prefix = !isempty(message.content) &&
                     message.content[1] isa Markdown.Paragraph
    if prepend_prefix
        # Hack: We prepend the prefix here to allow the markdown renderer to be
        # aware of the indenting which will result from prepending the prefix.
        # Without this we will get many issues of broken vertical alignment.
        # Avoid collisions: using placeholder from unicode private use area
        placeholder = '\uF8FF'^prefix_width
        pushfirst!(message.content[1].content, placeholder)
    end
    formatted = sprint(show, MIME"text/plain"(), message, context=io_context)
    msglines = split(chomp(formatted), '\n')
    # Hack': strip left margin which can't be configured in Markdown
    # terminal rendering.
    msglines = [startswith(s, "  ") ? s[3:end] : s for s in msglines]
    if prepend_prefix
        # Hack'': Now remove the prefix from the rendered markdown.
        msglines[1] = replace(msglines[1], placeholder=>""; count=1)
    elseif !isempty(msglines[1])
        pushfirst!(msglines, SubString(""))
    end
    msglines
end

# Formatting of values in key value pairs
function showvalue(io, key, msg)
    if key === :exception && msg isa Vector && length(msg) > 1 && msg[1] isa Tuple{Exception,Any}
        if VERSION >= v"1.2"
            # Ugly code path to support passing exception=Base.catch_stack() We
            # don't have a useful "Exception Stack" type to look for here until
            # https://github.com/JuliaLang/julia/pull/29901 gets unstuck.
            Base.show_exception_stack(io, msg)
        else
            # v1.0 and 1.1 don't have Base.show_exception_stack
            show(io, "text/plain", msg)
        end
    else
        show(io, "text/plain", msg)
    end
end
function showvalue(io, key, e::Tuple{Exception,Any})
    ex,bt = e
    showerror(io, ex, bt; backtrace = bt!==nothing)
end
showvalue(io, key, ex::Exception) = showerror(io, ex)

# Generate a text representation of all key value pairs, split into lines with
# per-line indentation as an integer.
function format_key_value_pairs(kwargs, io_context)
    msglines = Tuple{Int,String}[]
    valbuf = IOBuffer()
    dsize = displaysize(io_context)
    rows_per_value = max(1, dsize[1]÷(length(kwargs)+1))
    valio = IOContext(valbuf, IOContext(io_context, :displaysize=>(rows_per_value,dsize[2]-3)))
    for (key,val) in kwargs
        showvalue(valio, key, val)
        vallines = split(String(take!(valbuf)), '\n')
        if length(vallines) == 1
            push!(msglines, (2,SubString("$key = $(vallines[1])")))
        else
            push!(msglines, (2,SubString("$key =")))
            append!(msglines, ((3,line) for line in vallines))
        end
    end
    msglines
end

function findbar(bartree, id)
    if !(bartree isa AbstractArray)
        bartree.data.id === id && return bartree
    end
    for node in bartree
        found = findbar(node, id)
        found === nothing || return found
    end
    return nothing
end

function foldtree(op, acc, tree)
    for node in tree
        acc = foldtree(op, acc, node)
    end
    if !(tree isa AbstractArray)
        acc = op(acc, tree)
    end
    return acc
end

const BAR_MESSAGE_ID = gensym(:BAR_MESSAGE_ID)

function handle_progress(logger, progress, kwargs)
    node = findbar(logger.bartrees, progress.id)
    if node === nothing
        # Don't do anything when it's already done:
        (progress.done || something(progress.fraction, 0.0) >= 1) && return

        parentnode = findbar(logger.bartrees, progress.parentid)
        bar = ProgressBar(
            fraction = progress.fraction,
            name = progress.name,
            id = progress.id,
            parentid = progress.parentid,
        )
        if parentnode === nothing
            node = Node(bar)
            pushfirst!(logger.bartrees, node)
        else
            bar.level = parentnode.data.level + 1
            node = addchild(parentnode, bar)
        end
    else
        bar = node.data
        set_fraction!(bar, progress.fraction)
        if progress.name != ""
            bar.name = progress.name
        end
        node.data = bar
    end
    if progress.done
        if isroot(node)
            deleteat!(logger.bartrees, findfirst(x -> x === node, logger.bartrees))
        else
            prunebranch!(node)
        end
    end

    bartxt = sprint(context = :displaysize => displaysize(logger.stream)) do io
        foldtree(true, logger.bartrees) do isfirst, node
            isfirst || println(io)
            printprogress(io, node.data)
            false  # next `isfirst`
        end
    end
    if isempty(bartxt)
        pop!(logger.sticky_messages, BAR_MESSAGE_ID)
    else
        bartxt = sprint(context = logger.stream) do io
            printstyled(io, bartxt; color=:green)
        end
        push!(logger.sticky_messages, BAR_MESSAGE_ID => bartxt)
    end

    # "Flushing" non-sticky message should be done after the sticky
    # message is re-drawn:
    if progress.done
        ensure_done!(bar)
        donetxt = sprint(context = :displaysize => displaysize(logger.stream)) do io
            printprogress(io, bar)
        end
        printstyled(logger.stream, donetxt; color=:light_black)
        println(logger.stream)
    end
end

function handle_message(logger::TerminalLogger, level, message, _module, group, id,
                        filepath, line; maxlog=nothing,
                        sticky=nothing, kwargs...)
    if maxlog !== nothing && maxlog isa Integer
        remaining = get!(logger.message_limits, id, maxlog)
        logger.message_limits[id] = remaining - 1
        remaining > 0 || return
    end

    progress = asprogress(level, message, _module, group, id, filepath, line; kwargs...)
    if progress !== nothing
        handle_progress(logger, progress, kwargs)
        return
    end

    color,prefix,suffix = logger.meta_formatter(level, _module, group, id, filepath, line)

    # Generate a text representation of the message
    dsize = displaysize(logger.stream)
    context = IOContext(logger.stream, :displaysize=>(dsize[1],dsize[2]-2))
    msglines = format_message(message, textwidth(prefix), context)
    # Add indentation level
    msglines = [(0,l) for l in msglines]
    if !isempty(kwargs)
        ctx = logger.show_limited ? IOContext(context, :limit=>true) : context
        append!(msglines, format_key_value_pairs(kwargs, ctx))
    end

    # Format lines as text with appropriate indentation and with a box
    # decoration on the left.
    minsuffixpad = 2
    buf = IOBuffer()
    iob = IOContext(buf, logger.stream)
    nonpadwidth = 2 + (isempty(prefix) || length(msglines) > 1 ? 0 : length(prefix)+1) +
                  msglines[end][1] + termlength(msglines[end][2]) +
                  (isempty(suffix) ? 0 : length(suffix)+minsuffixpad)
    justify_width = min(logger.right_justify, dsize[2])
    if nonpadwidth > justify_width && !isempty(suffix)
        push!(msglines, (0,SubString("")))
        minsuffixpad = 0
        nonpadwidth = 2 + length(suffix)
    end
    for (i,(indent,msg)) in enumerate(msglines)
        boxstr = length(msglines) == 1 ? "[ " :
                 i == 1                ? "┌ " :
                 i < length(msglines)  ? "│ " :
                                         "└ "
        printstyled(iob, boxstr, bold=true, color=color)
        if i == 1 && !isempty(prefix)
            printstyled(iob, prefix, " ", bold=true, color=color)
        end
        print(iob, " "^indent, msg)
        if i == length(msglines) && !isempty(suffix)
            npad = max(0, justify_width - nonpadwidth) + minsuffixpad
            print(iob, " "^npad)
            printstyled(iob, suffix, color=:light_black)
        end
        println(iob)
    end

    msg = take!(buf)
    if sticky !== nothing
        # Ensure we see the last message, even if it's :done
        push!(logger.sticky_messages, id=>String(msg))
        if sticky == :done
            pop!(logger.sticky_messages, id)
        end
    else
        write(logger.stream, msg)
    end
    nothing
end


text:
import JSON

json(next, r::Resource, req, id) = begin
    try
        req[:body] = JSON.parse(req[:body] |> String)
    catch
    end
    res = next(req, id)
    isa(res, Union{Dict, Vector}) || return res
    res = res |> JSON.json |> Response
    res.headers["Content-Type"] = "application/json"
    res
end

cors(r::Resource, req, id, res) = begin
    ACAO = "Access-Control-Allow-Origin"
    ACRH = "Access-Control-Request-Headers"
    ACAH = "Access-Control-Allow-Headers"
    ACRM = "Access-Control-Request-Method"
    ACAM = "Access-Control-Allow-Methods"
    res.headers[ACAO] = "*"
    if haskey(req[:headers], ACRH)
        res.headers[ACAH] = req[:headers][ACRH]
    end
    if haskey(req[:headers], ACRM)
        res.headers[ACAM] = req[:headers][ACRM]
    end
end

cors(next::Function, r::Resource, req, id) = begin
    res = next(req, id)
    cors(r, req, id, res)
    res
end


text:
using Parameters
include("common.jl")
# ATP-dependent K channel (KATP) parameters
@with_kw struct KATPParams
    k_o = 5.4
    # Temperature adjustment
    FT = 1.0266
    # Default Michaelis constant for sodium
    KM_NA = 25.9
    # Default Michaelis constant for sodium
    KM_MG = 0.65 * sqrt(k_o + 5.0)
    # electrical distance of magnesium
    δMG = 0.32
    # electrical distance of sodium
    δNA = 0.35
    # Channel density (/cm²)
    σ = 6.0E8
    # Max channel open probability
    P₀ = 0.91
    # Unitary conductance of KATP channel
    γ = 35.375E-9 * (k_o / 5.4)^2.4
    # Mg factor of KATP channel
    IMAX = σ * P₀ * FT * γ
end

#=
ATP-inhibited K current from Gauthier et. al, 2012, based on Ferrero et al. 1996
Scalar version
=#
function ikatp(vm, eK, na_i, atp_i, adp_i, mg_i, KM_MG, KM_NA, δMG, δNA, IMAX)
    f_mg = _mm(KM_MG, mg_i * _ra(2 * δMG * vm))
    f_na = _hills(KM_NA, na_i * _ra(δNA * vm), 2)
    km_atp = 35.8e-3 + 17.9 * 1000^(-0.44) * adp_i^0.56
    h = 1.9 + 0.74 * exp(-adp_i)
    f_atp = _hills(km_atp, atp_i, h)
    iKatp = IMAX * f_mg * f_na * f_atp * (vm - eK)
end

function ikatp(vm, eK, na_i, atp_i, adp_i, mg_i, pKATP::KATPParams)
    @unpack KM_MG, KM_NA, δMG, δNA, IMAX = pKATP
    ikatp(vm, eK, na_i, atp_i, adp_i, mg_i, KM_MG, KM_NA, δMG, δNA, IMAX)
end


text:
module ZeroOrOne

export ZeroOne
"""
Zero or some (but we don't know how one)
"""
struct ZeroOne <: Number
    some::Bool
end

Base.promote_rule(::Type{ZeroOne}, ::Type{N}) where {N<:Number} = N

Base.:(==)(a::ZeroOne, b::ZeroOne) = a.some == b.some
Base.:(<)(a::ZeroOne, b::ZeroOne) = a.some < b.some
Base.hash(z::ZeroOne, h::UInt) = hash(0xcfe01124, hash(a.some, h))

Base.zero(::ZeroOne) = ZeroOne(false)
Base.one(::ZeroOne) = ZeroOne(true)
Base.Bool(a::ZeroOne) = a.some
Base.Int8(a::ZeroOne) = Int8(Bool(a))
Base.Int(a::ZeroOne) = Int(Bool(a))

# ZeroOne(b::Bool)
ZeroOne(b::Integer) = ZeroOne(b ≠ 0)

Base.:+(a::ZeroOne) = a
Base.:+(a::ZeroOne, b::ZeroOne) = ZeroOne(a.some | b.some)
Base.:*(a::ZeroOne, b::ZeroOne) = ZeroOne(a.some & b.some)
Base.abs(a::ZeroOne) = a
Base.abs2(a::ZeroOne) = a
Base.max(a::ZeroOne, b::ZeroOne) = a ≥ b ? a : b
Base.min(a::ZeroOne, b::ZeroOne) = a ≤ b ? a : b

Base.:-(a::ZeroOne) = -Int(a)
Base.:-(a::ZeroOne, b::ZeroOne) = Int(a) - Int(b)

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

export One
"""
One (always has the same value)
"""
struct One <: Number end

Base.promote_rule(::Type{One}, ::Type{N}) where {N<:Number} = N

Base.:(==)(::One, ::One) = true
Base.:(<)(::One, ::One) = false
Base.hash(::One, h::UInt) = hash(0xcfe01124, h)

Base.zero(::One) = zero(ZeroOne)
Base.one(::One) = One()
Base.Bool(::One) = true
Base.Int8(::One) = Int8(Bool(One()))
Base.Int(::One) = Int(Bool(One()))

One(b::Bool) = (@assert b; One())
One(b::Integer) = (@assert b == 1; One())

Base.:+(::One) = One()
Base.:+(::One, ::One) = One()
Base.:*(::One, ::One) = One()
Base.abs(::One) = One()
Base.abs2(::One) = One()
Base.max(::One, ::One) = One()
Base.min(::One, ::One) = One()

Base.:-(::One) = -1
Base.:-(::One, ::One) = 0

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

Base.promote_rule(::Type{ZeroOne}, ::Type{One}) = ZeroOne

ZeroOne(::One) = ZeroOne(Bool(One()))
One(a::ZeroOne) = One(Bool(a))

end


text:
"""
    TeXExpr(head::Symbol, args::Vector)

A TeX expression, represented by a head and any number of arguments.

See the documentation to see the valid combinations of head and arguments,
and the meaning of the latter.

Fields
======
    - head::Symbol
    - args::Vector
"""
struct TeXExpr
    head::Symbol
    args::Vector

    TeXExpr(head, args::Vector) = new(head, args)
end

TeXExpr(head) = TeXExpr(head, [])
TeXExpr(head, args...) = TeXExpr(head, collect(args))
TeXExpr(head, arg) = TeXExpr(head, [arg])

function Base.Char(texexpr::TeXExpr)
    if texexpr.head in [:char, :symbol, :digit]
        return texexpr.args[1]
    end

    throw(ArgumentError("cannot convert TeXExpr with head $(texexpr.head) to Char."))
end

"""
    manual_texexpr(data)

Convenience function to manually create a TeXExpr from tuples or latex strings.

# Example
julia> manual_texexpr((:expr, 'a', (:spaced, '+'), raw"\\omega^2"))
TeXExpr :expr
├─ 'a'
├─ TeXExpr :spaced
│  └─ '+'
└─ TeXExpr :decorated
   ├─ TeXExpr :symbol
   │  ├─ 'ω'
   │  └─ "\\omega"
   ├─ nothing
   └─ '2'
```
"""
function manual_texexpr(tuple::Tuple)
    head = tuple[1]
    args = []

    if head in [:char, :digit, :symbol]
        return TeXExpr(head, tuple[2])
    end

    for arg in tuple[2:end]
        push!(args, manual_texexpr(arg))
    end

    return TeXExpr(head, args)
end

function manual_texexpr(str::LaTeXString)
    expr = texparse(str[2:end])
    if length(expr.args) == 1
        return first(expr.args)
    else
        return TeXExpr(:group, expr.args)
    end
end

manual_texexpr(char::Char) = TeXExpr(:char, char)
manual_texexpr(any) = any

head(texexpr::TeXExpr) = texexpr.head
head(::Char) = :char

Base.copy(texexpr::TeXExpr) = TeXExpr(texexpr.head, deepcopy(texexpr.args))

AbstractTrees.children(texexpr::TeXExpr) = texexpr.args
AbstractTrees.printnode(io::IO, texexpr::TeXExpr) = print(io, "TeXExpr :$(texexpr.head)")

function Base.show(io::IO, texexpr::TeXExpr)
    print_tree(io, texexpr, maxdepth=10)
end

function Base.:(==)(tex1::TeXExpr, tex2::TeXExpr)
    childs1 = children(tex1)
    childs2 = children(tex2)

    length(childs1) != length(childs2) && return false

    return all(childs1 .== childs2)
end


text:
# This file is a part of GERDADeepLearning.jl, licensed under the MIT License (MIT).

export FileKey
immutable FileKey #datafile name
  run::Integer
  event_type::Symbol # :cal, :phy
  fileid::AbstractString
  id::AbstractString

  function FileKey(id::AbstractString)
    id = strip(id)
    if id[1:6] != "gerda-"
      throw(ArgumentError("$id is not a valid file identifier"))
    end
    run = parse(Int, id[10:13])
    event_type = Symbol(id[32:34])
    fileid = id[15:30]
    return new(run, event_type, fileid, id)
  end
end

function Base.string(key::FileKey)
  return key.id
end

export KeyList
immutable KeyList
  name::String #.txt name
  entries::Vector{FileKey} #datafiles
end

export path #returns path of actual root datafile
function path(base_dir::AbstractString, key::FileKey, tier)
  subfolder = (tier == :tier1 || tier == :tier2 || tier == :tierW) ? "ged" : "all"
  str_tier = string(tier)
  str_run = "run"*lpad(key.run, 4, 0) # pad e.g. 89 to 0089
  str_event_type = string(key.event_type)
  filename = "gerda-$str_run-$(key.fileid)-$str_event_type-$subfolder-$str_tier.root"
  return joinpath(base_dir, str_tier, subfolder, str_event_type, str_run, filename)
end

export parse_keylist #parses whats in the txt
function parse_keylist(keylist_file::AbstractString, name::String) #path to file and filename
  entries = [FileKey(id) for id in open(readlines, keylist_file)] # gets actual data file names
  return KeyList(name, entries) # a .txt file with its entries
end

export name
function listname(keylist::KeyList)
  return keylist.name
end

function Base.merge(keylists::KeyList...)
  return KeyList("", vcat([kl.entries for kl in keylists]...))
end

function Base.length(keylist::KeyList)
  return length(keylist.entries)
end

function Base.start(keylist::KeyList)
  return start(keylist.entries)
end

function Base.done(keylist::KeyList, state)
  return done(keylist.entries, state)
end

function Base.next(keylist::KeyList, state)
  return next(keylist.entries, state)
end

phase2_detectors = ["GD91A", "GD35B", "GD02B", "GD00B", "GD61A", "GD89B", "GD02D", "GD91C", # string1
		"ANG5", "RG1", "ANG3", # string2
		"GD02A", "GD32B", "GD32A", "GD32C", "GD89C", "GD61C",  "GD76B", "GD00C", # string3
		"GD35C", "GD76C", "GD89D", "GD00D", "GD79C", "GD35A", "GD91B", "GD61B", # string4
		"ANG2", "RG2", "ANG4", # string 5
		"GD00A", "GD02C", "GD79B", "GD91D", "GD32D", "GD89A", "ANG1", # string6
		"GTF112", "GTF32", "GTF45_2" # string7
    ]

export get_detector_index
function get_detector_index(detector_name::String)
  return find(n->n==detector_name, phase2_detectors)[1] - 1
end

export get_detector_name
function get_detector_name(detector_index::Integer)
  return phase2_detectors[detector_index+1]
end

export parse_detectors
function parse_detectors(list::Vector)
  return [(isa(det,Integer) ? phase2_detectors[det+1] : det) for det in list]
end

BEGes_GERDA_II() = ["GD91A", "GD35B", "GD02B", "GD00B", "GD61A", "GD89B", "GD02D", "GD91C", "GD02A", "GD32B", "GD32A", "GD32C", "GD89C", "GD61C", "GD76B", "GD00C", "GD35C", "GD76C", "GD89D", "GD00D", "GD79C", "GD35A", "GD91B", "GD61B", "GD00A", "GD02C", "GD79B", "GD91D", "GD32D", "GD89A"]

Coax_GERDA_II() = ["ANG1", "ANG2", "ANG3", "ANG4", "ANG5", "RG1", "RG2"]

Natural_GERDA_II() = ["GTF112", "GTF32", "GTF45_2"]

Used_GERDA_II() = ["GD91A", "GD35B", "GD02B", "GD00B", "GD61A", "ANG5", "RG1", "ANG3", "GD02A", "GD32B", "GD32A", "GD32C", "GD89C", "GD61C", "GD76B", "GD00C", "GD35C", "GD76C", "GD89D", "GD00D", "GD79C", "GD35A", "GD91B", "GD61B", "ANG2", "RG2", "ANG4", "GD00A", "GD02C", "GD79B", "GD91D", "GD32D", "GD89A", "ANG1"]


text:
#ENV["CPLEX_STUDIO_BINARIES"]="/home/software/cplex/128/cplex/bin/x86-64_linux"
using DelimitedFiles, LinearAlgebra
include("src/core.jl")

function getProblemDataBeasley(pathName)
    testData=readdlm(pathName)
    n=testData[1,1]
    μ=testData[2:n+1, 1]
    μ=convert(Array{Float64}, μ)
    σ=testData[2:n+1, 2]
    Σ=zeros(n,n)
    for iter=n+2:size(testData,1)
        # This loop inserts diagonal elements twice. This is ok as benign, and a "fix" would require more lines of code.
        Σ[testData[iter, 1], testData[iter, 2]]=σ[testData[iter, 1]]*σ[testData[iter, 2]]*testData[iter, 3]
        Σ[testData[iter, 2], testData[iter, 1]]=σ[testData[iter, 1]]*σ[testData[iter, 2]]*testData[iter, 3]
    end
    X=sqrt(Σ)
    Y=(X*X')\(X*μ)
    d=X'*Y-μ
    d=round.(d, digits=8) #avoids round-off errors when Σ is full rank
    A=zeros(1,size(X,2))
    l=zeros(size(A,1))
    u=zeros(size(A,1))
    min_investment=zeros(n)
    f=size(X, 2)
    m=size(A, 1)
    k=5
    γ=(100/sqrt(n))*ones(n)

    return SparsePortfolioData(μ, Y, d, X, A, l, u, k, n, m, f, min_investment, γ)
end

for portfolioNumber=1:5
    pathName = "/Users/Ryancw2/Dropbox (MIT)/Sparse Portfolios Project/OptimalMarkowitzPortfolios.jl/data/port"
    pathName=pathName*string(portfolioNumber)*".txt"
    thePortfolioData=getProblemDataBeasley(pathName)

    sparse_Portfolio=cutting_planes_portfolios(thePortfolioData; useSOCPLB=false, useHeuristic=false);
    @show @test(abs.(sum(sparse_Portfolio.w)-1.0)<1e-4)
    println("Variance is:", sparse_Portfolio.portfolioVariance)
    println("Return is:", sparse_Portfolio.expectedReturn)
    println("Reg term is:", 0.5*sum(sparse_Portfolio.w[i]^2/thePortfolioData.γ[sparse_Portfolio.indices[i]] for i=1:size(sparse_Portfolio.indices,1)))
    open("results/Table3ResultsOutput.csv","a") do fp
    println(fp, "Instance, cardinality, regularizer, method, solvetime, nodecount, iters")
    println(fp, "beasley"*string(portfolioNumber), ",", thePortfolioData.k, ",",
                     "100/sqrt(n)", ",", "OA", ",",  sparse_Portfolio.solveTime, ",", sparse_Portfolio.nodeCount, ",", sparse_Portfolio.numIters)
    end

    # # Test against Big-M for correctness
    sparse_Portfolio_Cplex_bigM=cplex_raw_bigM(thePortfolioData, 300.0, 1e-8, 1)
    @show @test(abs.(sum(sparse_Portfolio_Cplex_bigM.w)-1.0)<1e-4)
    println("Variance is:", sparse_Portfolio_Cplex_bigM.portfolioVariance)
    println("Return is:", sparse_Portfolio_Cplex_bigM.expectedReturn)
    println("Reg term is:", 0.5*sum(sparse_Portfolio_Cplex_bigM.w[i]^2/thePortfolioData.γ[sparse_Portfolio_Cplex_bigM.indices[i]] for i=1:size(sparse_Portfolio_Cplex_bigM.indices,1)))
    open("results/Table3ResultsOutput.csv","a") do fp
    println(fp, "Instance, cardinality, regularizer, method, solvetime, nodecount")
    println(fp, "beasley"*string(portfolioNumber), ",", thePortfolioData.k, ",",
                     "100/sqrt(n)", ",", "CPLEX Big-M", ",",  sparse_Portfolio_Cplex_bigM.solveTime, ",", sparse_Portfolio_Cplex_bigM.nodeCount)
    end

    # # Test against MISOCP for correctness
    sparse_Portfolio_Cplex=cplex_raw_MISOCP(thePortfolioData, 300.0, 1e-8, 1)
    @show @test(abs.(sum(sparse_Portfolio_Cplex.w)-1.0)<1e-4)
    println("Variance is:", sparse_Portfolio_Cplex.portfolioVariance)
    println("Return is:", sparse_Portfolio_Cplex.expectedReturn)
    println("Reg term is:", 0.5*sum(sparse_Portfolio_Cplex.w[i]^2/thePortfolioData.γ[sparse_Portfolio_Cplex.indices[i]] for i=1:size(sparse_Portfolio_Cplex.indices,1)))
    open("results/Table3ResultsOutput.csv","a") do fp
    println(fp, "Instance, cardinality, regularizer, method, solvetime, nodecount")
    println(fp, "beasley"*string(portfolioNumber), ",", thePortfolioData.k, ",",
                    "100/sqrt(n)", ",", "CPLEX MISOCP", ",",  sparse_Portfolio_Cplex.solveTime, ",", sparse_Portfolio_Cplex.nodeCount)
    end
end


text:
function Base.show(io::IO,::MIME"text/plain",net::FastNet)
    n=countnodes(net)
    k=countlinks(net)
    print(io,"Network of $n nodes and $k links" )
end

function Base.show(io::IO,::MIME"text/html",net::FastNet)
    n=countnodes(net)
    k=countlinks(net)
    print(io,"Network of $n nodes and $k links" )
end

function Base.show(io::IO,::MIME"text/plain",sim::FastSim)
    print(io,"Simulation run, currently at time $(sim.t)" )
end

function Base.show(io::IO,::MIME"text/html",sim::FastSim)
    print(io,"Simulation run, currently at time $(sim.t)" )
end

function Base.show(io::IO,::MIME"text/plain",lt::LinkType)
    printlinktype(io,lt)
end

function Base.show(io::IO,::MIME"text/html",lt::LinkType)
    printlinktype(io,lt)
end