text:
module Options

using ..IonBase: dot_ion, ion_toml, init_dot_ion, templates
using Configurations
using Comonicon.Tools: cmd_error
using Pkg.Types: VersionSpec

function default_install_dir()
    haskey(ENV, "JULIA_INSTALL_PATH") && return ENV["JULIA_INSTALL_PATH"]
    install_dir = dot_ion("packages")
    if !ispath(install_dir)
        mkpath(install_dir)
    end
    return install_dir
end

@option mutable struct Julia
    active::Union{String, Nothing} = nothing
    stable::Union{VersionNumber, Nothing} = nothing
    nightly::Union{String, Nothing} = nothing
    install_dir::String = default_install_dir()
    # NOTE: we store nightly here too
    versions::Dict{VersionNumber, String} = Dict{VersionNumber, String}()
end

@option struct Plugins
    template::Vector{String} = String[]
    registry::Vector{String} = String[]
end

@option mutable struct Ion
    username::String = ""
    julia::Julia = Julia()
    plugins::Plugins = Plugins()
    default_template::String = "basic"
    templates::Dict{String, String} = Dict{String, String}(
        "basic" => templates("basic.toml"),
        "package" => templates("package.toml"),
        "academic" => templates("academic.toml"),
        "comonicon" => templates("comonicon.toml"),
        "comonicon-sysimg" => templates("comonicon-sysimg.toml"),
    )
end

function Configurations.option_convert(::Type{Julia}, ::Type{Dict{VersionNumber, String}}, x::Dict{String})
    d = Dict{VersionNumber, String}()
    for (k, v) in x
        d[VersionNumber(k)] = string(v)
    end
    return d
end

function load_plugins(m::Module, type::Symbol, ion::Ion=read())
    ex = Expr(:using)
    for each in getfield(ion.plugins, type)
        push!(ex.args, Expr(:(.), Symbol(each)))
    end
    Base.eval(m, ex)
    return
end

function read()::Ion
    if !ispath(ion_toml())
        return Ion()
    else
        return from_toml(Ion, ion_toml())
    end
end

# let's just pirate this
Base.String(x::VersionNumber) = string(x)
Base.convert(::Type{VersionNumber}, x::String) = VersionNumber(x)

toml_convert(x) = x

function toml_convert(x::VersionNumber)
    string(x)
end

function dump(option::Ion)
    init_dot_ion() # make sure .ion is avialable
    toml(toml_convert, ion_toml(), option)
    return
end

function withion(f)
    ion = read()
    ret = f(ion)
    dump(ion)
    return ret
end

function find_julia_bin(version::String="stable", option::Ion=read())
    find_julia_bin(version, option.julia)
end

function find_version(spec::String, versions)
    version_spec = VersionSpec(spec)
    curr = nothing
    for each in versions
        if each in version_spec && (curr === nothing || curr <= each)
            curr = each
        end
    end

    return curr
end

function find_julia_bin(version_string::String, option::Julia)
    if version_string == "stable"
        option.stable === nothing && return
        version_string = string(option.stable)
    elseif version_string == "latest" || version_string == "nightly"
        return option.nightly
    end

    match = find_version(version_string, keys(option.versions))
    match === nothing && return
    return option.versions[match]
end

function active_julia_bin(ion::Ion=read())::String
    # 1. use user specified julia binary
    if haskey(ENV, "JULIA_EXECUTABLE_PATH")
        return ENV["JULIA_EXECUTABLE_PATH"]
    end

    if !isnothing(ion.julia.active)
        return ion.julia.active
    end

    cmd_error(
        "cannot detect julia binary, " *
        "please install julia using: ion install julia [--version=stable] " *
        "or specify via environment variable JULIA_EXECUTABLE_PATH"
    )
end

function active_julia_depots(ion::Ion=read())::Vector{String}
    julia_bin = active_julia_bin(ion)
    return withenv("JULIA_PROJECT"=>nothing, "JULIA_LOAD_PATH"=>nothing, "JULIA_DEPOT_PATH"=>nothing) do
        Base.eval(Meta.parse(readchomp(`$julia_bin -E 'using Pkg; Pkg.Types.depots()'`)))
    end
end

function active_julia_version(ion::Ion=read())
    julia_bin = active_julia_bin(ion)
    return Base.eval(Meta.parse(readchomp(`$julia_bin -E 'VERSION'`)))
end

end # Options

"""
config Ion manager.

# Args

- `item`: the configuration item you want to set or view.
"""
@cast function config(item::String)
    ion = Options.read()
    @smatch split(item, "=") begin
        ["julia.install_dir"] => println(ion.julia.install_dir)
        ["julia.install_dir", value] => (ion.julia.install_dir = value;)
        ["username"] => println(ion.username)
        ["username", value] => (ion.username = value;)
        _ => cmd_error("invalid configuration item")
    end
    Options.dump(ion)
    return
end


text:
type InterpolatedHazardRateCurve{DC <: DayCount, P <: Interpolation, B <: BusinessCalendar} <: InterpolatedDefaultProbabilityCurve{P}
  settlementDays::Int
  referenceDate::Date
  dc::DC
  interp::P
  cal::B
  dates::Vector{Date}
  times::Vector{Float64}
  data::Vector{Float64}
end

function InterpolatedHazardRateCurve(dates::Vector{Date},
                                    hazardRates::Vector{Float64},
                                    dc::DayCount,
                                    interpolator::Interpolation)
  ihc = InterpolatedHazardRateCurve(0, dates[1], dc, interpolator, NullCalendar(), dates, zeros(length(dates)), hazardRates)

  initialize!(ihc)

  return ihc
end

function initialize!(ihc::InterpolatedHazardRateCurve)
  length(ihc.dates) == length(ihc.data) || error("dates / data mismatch")

  ihc.times[1] = 0.0
  for i = 2:length(ihc.dates)
    ihc.times[i] = year_fraction(ihc.dc, ihc.dates[1], ihc.dates[i])
  end

  # initialize interpolator
  Math.initialize!(ihc.interp, ihc.times, ihc.data)

  Math.update!(ihc.interp)

  return ihc
end

function hazard_rate_impl(ihc::InterpolatedHazardRateCurve, t::Float64)
  if t <= ihc.times[end]
    return Math.value(ihc.interp, t)
  end

  # flat hazard rate extrap
  return ihc.data[end]
end

function survival_probability_impl(ihc::InterpolatedHazardRateCurve, t::Float64)
  if t == 0.0
    return 1.0
  end

  if t <= ihc.times[end]
    integral = get_primitive(ihc.interp, t, true)
  else
    # flat hazard rate extrapolation
    integral = get_primitive(ihc.interp, ihc.times[end], true) + ihc.data[end] * (t - ihc.times[end])
  end

  return exp(-integral)
end

ext: jl


text:
n = 4

srcMap = BlockMap(4*n, comm)
desMap = BlockMap(n*numProc(comm), collect((1:n) .+ n*(pid%4)), comm)

function basicMPITest(impor)
    if isa(impor, Import)
        data = impor.importData
    else
        data = impor.exportData
    end
    @test srcMap == data.source
    @test desMap == data.target
    @test 0 == data.numSameIDs
    @test isa(data.distributor, Distributor{UInt64, UInt16, UInt32})
    @test [] == data.permuteToLIDs
    @test [] == data.permuteFromLIDs
    #TODO test remoteLIDs, exportLIDs, exportPIDs
    @test true == data.isLocallyComplete
end

# basic import
basicMPITest(Import(srcMap, desMap))
basicMPITest(Import(srcMap, desMap, nothing))

# basic export
basicMPITest(Export(srcMap, desMap))
basicMPITest(Export(srcMap, desMap, nothing))

max_stars_repo_path: test/MPIimport-export Tests.jl

max_issues_repo_path: test/MPIimport-export Tests.jl

max_forks_repo_path: test/MPIimport-export Tests.jl


text:
# This is the main module that encompasses everything.

# ---------------------------------------------------
module Deuron
# ---------------------------------------------------
using JSON

include("./core_interfaces.jl")

include("./core.jl")
include("./config.jl")
include("./file_io/logs.jl")

include("./sim/sim.jl")

core = Core()

# A paused simulation means we can resume when we return to the App
# A stopped simulation means the simulation starts fresh and ignores
#   any previous state. Quiting while a simulation is running means
#   we have stopped the simulation.
paused = false
stopped = false
reset = false
running = false
suspended = false
crashed = false

function boot()
    println("--> Welcome to Deuron simulation <--")
    msg = string("Previous run state: ", Config.exit_state(Deuron.core.config))
    println(msg)
    print_help()

    Logs.log_to_app(Deuron.core.logs, "", true)
    Logs.log_to_app(Deuron.core.logs, "---------------------------------------------", true)
    Logs.log_to_app(Deuron.core.logs, "******************* Ready *******************", true)
    Logs.log_to_app(Deuron.core.logs, "---------------------------------------------", true)
    Logs.log_to_app(Deuron.core.logs, msg, true)

    while true
        # println("input: ", input)

        msg = get_input()

        if msg == "q"  # quit
            if running && !paused
                println("?? Simulation is running! Are you sure you want to stop and quit ??")
                println("[y/n]?")
                print(">")
                msg = get_input()
                if msg == "y"
                    Logs.log_to_app(Deuron.core.logs, "Stopping...", true)
                    global stopped = true
                end
            elseif paused
                Logs.log_to_app(Deuron.core.logs, "Suspending...", true)
                global stopped = true
                global suspended = true
            else
                Logs.log_to_app(Deuron.core.logs, "Quit", true)
                println("Goodbye.")
                Base.exit(0)
            end
        elseif msg == "r"   # run/start
            if running || paused
                println("** Sim is paused or already running **")
            else
                global paused = false
                global running = true
                @async Sim.entry()
            end
        elseif msg == "p"   # pause
            if paused
                println("** Already paused **")
            else
                Logs.log_to_app(Deuron.core.logs, "Pausing...", true)
                global paused = true
            end
        elseif msg == "u"  # resume
            if !paused
                println("** Sim isn't paused **")
            else
                Logs.log_to_app(Deuron.core.logs, "Resuming...", true)
                global paused = false
            end
        elseif msg == "a"   # status
            if paused
                println("-- Simulation is paused --")
            elseif running
                println("-- Simulation is running --")
            elseif !running
                println("-- Simulation is idle --")
            end
        elseif msg == "s"   # reset
            global paused = true
            global reset = true
        end

        print_help()
    end
end

function get_input()
    input = readavailable(stdin)
    # Strip newline off
    String(input[1:length(input) - 1])
end

function print_help()
    # "r" either resumes any previous sim that was in progress or starts a new one.
    # The config file identifies which.
    if !stopped
        println("--------------------------------------------")
        println("Commands:")
        println("  q: quit")
        println("  r: run simulation")
        println("  p: pause simulation")
        println("  u: resume simulation")
        println("  s: reset simulation")
        println("  a: status of simulation")
        println("  h: this help menu")
        println("--------------------------------------------")
        print(">")
    end
end

function resume()
    global paused = false
end

function reset_complete()
    global reset = false
end

function stop()
    global running = false
end

function set_crashed()
    global crashed = true
end

end # End module -------------------------------------


text:
@testset "Satellite" begin
    @testset "Auxiliarly functions" begin
        doa = Spherical(1.0, 0.0, π / 2) # Zenith
        distance_from_earth_center = EARTH_RADIUS + 40000m
        distance = @inferred GNSSSimulator.calc_sat_user_distance(
            distance_from_earth_center,
            doa
        )

        @test distance == 40000m

        @test @inferred(GNSSSimulator.calc_carrier_phase(0m, 20Hz)) == 0

        doppler = @inferred GNSSSimulator.calc_doppler(
            distance_from_earth_center,
            doa,
            14_000.0m / 3.6s,
            1.157542e6Hz
        )
        @test doppler ≈ 0Hz rtol = 1
    end

    @testset "Create satellite with integers" begin
        system = GPSL1()
        sat = ConstantDopplerSatellite(
            system,
            1,
            carrier_doppler = 1000Hz,
            carrier_phase = 1,
            code_phase = 100,
            cn0 = 45dBHz
        )

        @test @inferred(GNSSSimulator.get_carrier_phase(sat)) == 1
        @test @inferred(GNSSSimulator.get_code_phase(sat)) == 100
        @test @inferred(GNSSSimulator.get_carrier_doppler(sat)) == 1000Hz
    end

    @testset "Satellite signal" begin
        system = GPSL1()
        sat = ConstantDopplerSatellite(
            system,
            1,
            carrier_doppler = 1000.0Hz,
            carrier_phase = π / 2,
            code_phase = 100.0,
            cn0 = 45dBHz
        )

        @test @inferred(GNSSSimulator.get_carrier_phase(sat)) == π / 2
        @test @inferred(GNSSSimulator.get_code_phase(sat)) == 100

        rng = MersenneTwister(1234)
        signal = @inferred GNSSSimulator.gen_signal!(
            sat,
            2.5e6Hz,
            10.0Hz,
            1/Hz,
            2500,
            rng
        )

        code_doppler = 1000.0 * 1023e3 / 1.57542e9
        reference_signal = get_code.(
            system,
            (0:2499) .* (1023e3 .+ code_doppler) ./ 2.5e6 .+ 100,
            1
        ) .* cis.(2π .* (0:2499) .* (1000.0 + 10.0) ./ 2.5e6 .+ π / 2) .* sqrt(10^(45 / 10))
        @test sqrt(sum(abs2.(signal .- reference_signal)) / 2500) < 1e-2 * sqrt(10^(45 / 10))

        next_sat = @inferred GNSSSimulator.propagate(
            sat,
            2500,
            100.0Hz,
            2.5e6Hz,
            1.0/Hz,
            Random.GLOBAL_RNG
        )

        @test @inferred(get_carrier_doppler(next_sat)) == 1000Hz
        @test @inferred(get_code_doppler(next_sat)) ≈ 1000Hz / 1540
        @test @inferred(get_carrier_phase(next_sat)) ≈ mod2pi(
            π / 2 + 2π * 1100.0Hz * 2500 / 2.5e6Hz + π
        ) - π
        @test @inferred(get_code_phase(next_sat)) ≈ mod(
            100.0 + (1023e3Hz + 1000.0Hz / 1540) * 2500 / 2.5e6Hz,
            1023
        )
        @test @inferred(get_existence(next_sat)) == true
        @test @inferred(get_carrier_to_noise_density_ratio(next_sat)) == 45dBHz
        @test @inferred(get_prn(next_sat)) == 1
        @test @inferred(get_gnss_system(next_sat)) == system

    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ e53e9d5a-6091-11eb-3582-398334c98d77
begin
	# Load required packages
	using Pkg, HTTP, CSV, DataFrames, Random, Distances,
        UMAP, Clustering, Gadfly, DelimitedFiles,
        MultivariateStats, Cairo, Fontconfig, StatsBase

	# Set default plot size
    set_default_plot_size(50cm, 20cm)

	# Set working directory
	cd("/Users/martintreppner/Desktop")
end;

# ╔═╡ f1a368be-6091-11eb-0da3-77b12261da8c
begin
	# Install and Load single-cell deep Boltzmann machine (scDBM)
	# Pkg.add(PackageSpec(url="https://github.com/MTreppner/scDBM.jl", rev="master"))
	using scDBM
end;

# ╔═╡ c5044126-60b8-11eb-1d18-33c34995a8d5
begin
	using PyCall

	n_genes = 2000
	# Load python packages
	random = pyimport("random")
	os = pyimport("os")
	np = pyimport("numpy")
	pd = pyimport("pandas")
	scvi = pyimport("scvi")
	scvi_dataset = pyimport("scvi.dataset")
	scvi_models = pyimport("scvi.models")
	scvi_inference = pyimport("scvi.inference")

	dir = "notebook_test/segerstolpe_hvg.csv";

	countmatrix_scvi = scvi_dataset.CsvDataset(dir, save_path = "", new_n_genes = n_genes);
	# Set hyperparameters
	n_latent=14
	n_hidden=64
	n_layers=2
	n_epochs=200
	lr_scvi=0.01
	use_batches = false
	use_cuda = true

	# Train the model and output model likelihood every epoch
	vae = scvi_models.VAE(countmatrix_scvi.nb_genes,
		n_batch=countmatrix_scvi.n_batches * use_batches,
		n_latent=n_latent,
		n_hidden=n_hidden,
		n_layers=n_layers,
		reconstruction_loss="nb",
		dispersion = "gene"
	);

	trainer = scvi_inference.UnsupervisedTrainer(vae,
		countmatrix_scvi,
		train_size=0.7,
		use_cuda=use_cuda,
		frequency=1,
		n_epochs_kl_warmup=75,
		batch_size=36
	);

	trainer.train(n_epochs=n_epochs,lr=lr_scvi)

	full = trainer.create_posterior(trainer.model, countmatrix_scvi);

	batch_size = 32
	n_cells = 16

	# Generate Data from posterior
	gen_data_scvi = full.generate(batch_size=batch_size,n_samples=10, sample_prior=false, n_cells=n_cells);
	gen_data_scvi_orig = Array{Float64,2}(gen_data_scvi[1][:,:,2]);

####### scVI Zero data

	n_genes = 2000
	# Load python packages
	random = pyimport("random")
	os = pyimport("os")
	np = pyimport("numpy")
	pd = pyimport("pandas")
	scvi = pyimport("scvi")
	scvi_dataset = pyimport("scvi.dataset")
	scvi_models = pyimport("scvi.models")
	scvi_inference = pyimport("scvi.inference")

	dir = "notebook_test/segerstolpe_hvg_zeros.csv";

	countmatrix_scvi = scvi_dataset.CsvDataset(dir, save_path = "", new_n_genes = n_genes);

	# Set hyperparameters
	n_latent=14
	n_hidden=64
	n_layers=2
	n_epochs=200
	lr_scvi=0.01
	use_batches = false
	use_cuda = true

	# Train the model and output model likelihood every epoch
	vae = scvi_models.VAE(countmatrix_scvi.nb_genes,
		n_batch=countmatrix_scvi.n_batches * use_batches,
		n_latent=n_latent,
		n_hidden=n_hidden,
		n_layers=n_layers,
		reconstruction_loss="nb",
		dispersion = "gene"
	);

	trainer = scvi_inference.UnsupervisedTrainer(vae,
		countmatrix_scvi,
		train_size=0.7,
		use_cuda=use_cuda,
		frequency=1,
		n_epochs_kl_warmup=75,
		batch_size=36
	);

	trainer.train(n_epochs=n_epochs,lr=lr_scvi)

	full = trainer.create_posterior(trainer.model, countmatrix_scvi);

	batch_size = 32
	n_cells = 16

	# Generate Data from posterior
	gen_data_scvi = full.generate(batch_size=batch_size,n_samples=10, sample_prior=false, n_cells=n_cells);
	gen_data_scvi_zeros = Array{Float64,2}(gen_data_scvi[1][:,:,2]);

end;

# ╔═╡ a6314f18-60a5-11eb-0307-450764cbf160
md"
#### Supplementary Material:

##### Synthetic Single-Cell RNA-Sequencing Data from Small Pilot Studies using Deep Generative Models
"

# ╔═╡ 81a8c5a8-6092-11eb-3414-69625042ef17
md"
###### scDBM on original data
"

# ╔═╡ f6e0abde-6091-11eb-2a1c-010ab8a51452
begin
	global countmatrix = CSV.read("notebook_test/segerstolpe_hvg.csv")
	global countmatrix = Array{Float64,2}(Array{Float64,2}(countmatrix[:,2:end])')

	# Sample pilot dataset
	pilot_data = countmatrix;

	# Train- Test-split
	Random.seed!(101);
	data, datatest = splitdata(pilot_data, 0.3);
	datadict = DataDict("Training data" => data, "Test data" => datatest);

	# Set hyperparameters
	epochs = 60;
	init_disp = (ones(size(data,2)) .* 0.5);
	lr = [1.0e-11*ones(6000);0.00005*ones(1000)];
	regularization = 2.5;

	# Train scDBM
	Random.seed!(59);
	global dbm = fitdbm(data,
		epochs = 10,
		learningrate = 1.0e-11,
		batchsizepretraining = 16,
		  pretraining = [
				# Negative-Binomial RBM
				TrainLayer(nhidden = 64,
				learningrates = lr,
				epochs = epochs,
				rbmtype = NegativeBinomialBernoulliRBM,
				inversedispersion = init_disp,
				fisherscoring = 1,
				lambda = regularization),
				# Bernoulli RBM
				TrainLayer(nhidden = 16,
				learningrate = 1.0e-4,
				epochs = 40)]
		)

	# Generate synthetic cells
	Random.seed!(559);
	number_gensamples = 2090#size(pilot_data,1);
	synthetic_cells = initparticles(dbm, number_gensamples);
	gibbssamplenegativebinomial!(pilot_data,
		synthetic_cells,
		dbm,
		30
	);
end;

# ╔═╡ a02127d2-6092-11eb-26fc-e749e38caaca
begin
	using Distributions
	global countmatrix_impute = CSV.read("notebook_test/segerstolpe_hvg.csv")
	global countmatrix_impute = Array{Float64,2}(Array{Float64,2}(countmatrix_impute[:,2:end])')

	Random.seed!(112);
	tmp = rand.(Uniform(0,1), size(countmatrix_impute,1), size(countmatrix_impute,2))
	tmp1 = tmp .> 0.2
	tmp_zero_proportion = mapslices(sum, tmp1, dims=1) ./ size(countmatrix_impute,1)

	countmatrix_impute = countmatrix_impute .* tmp1

	tmp = DataFrame(countmatrix_impute')
	CSV.write("notebook_test/segerstolpe_hvg_zeros.csv", tmp, writeheader = true)

	# Sample pilot dataset
	pilot_data_zeros = countmatrix_impute;

	# Train- Test-split
	Random.seed!(101);
	data_zeros, datatest_zeros = splitdata(pilot_data_zeros, 0.3);
	datadict_zeros = DataDict("Training data" => data_zeros,
		"Test data" => datatest_zeros);

	# Set hyperparameters
	epochs_zeros = 60;
	init_disp_zeros = (ones(size(data_zeros,2)) .* 0.5);
	lr_zeros = [1.0e-11*ones(6000);0.00005*ones(1000)];
	regularization_zeros = 2.5;

	# Train scDBM
	Random.seed!(59);
	global dbm_zeros = fitdbm(data_zeros,
		epochs = 10,
		learningrate = 1.0e-11,
		batchsizepretraining = 16,
		  pretraining = [
				# Negative-Binomial RBM
				TrainLayer(nhidden = 64,
				learningrates = lr_zeros,
				epochs = epochs_zeros,
				rbmtype = NegativeBinomialBernoulliRBM,
				inversedispersion = init_disp_zeros,
				fisherscoring = 1,
				lambda = regularization_zeros),
				# Bernoulli RBM
				TrainLayer(nhidden = 16,
				learningrate = 1.0e-4,
				epochs = 40)]
		)

	# Generate synthetic cells
	Random.seed!(559);
	number_gensamples_zeros = 2090
	synthetic_cells_zeros = initparticles(dbm_zeros, number_gensamples);
	gibbssamplenegativebinomial!(pilot_data_zeros,
		synthetic_cells_zeros,
		dbm_zeros,
		30
	);
end;

# ╔═╡ 9509ec9e-6092-11eb-2157-a540547604cb
md"
###### scDBM on original data with added zeros
"

# ╔═╡ e40e07ce-6094-11eb-1137-b5acb0fec21a
md"
###### scVI
"

# ╔═╡ e5f44554-60a2-11eb-33ab-e199bec0c726
begin
	cd("/Users/martintreppner/Desktop")

	# scDBM
	scdbm_out = DataFrame(synthetic_cells[1])
	CSV.write("notebook_test/imputation_scdbm.csv", scdbm_out, writeheader = true)

	scdbm_zeros_out = DataFrame(synthetic_cells_zeros[1])
	CSV.write("notebook_test/imputation_scdbm_zeros.csv",
		scdbm_zeros_out,
		writeheader = true)

	# scVI
	scvi_out = DataFrame(gen_data_scvi_orig)
	CSV.write("notebook_test/imputation_scvi.csv",
		scvi_out,
		writeheader = true)

	scvi_zeros_out = DataFrame(gen_data_scvi_zeros)
	CSV.write("notebook_test/imputation_scvi_zeros.csv",
		scvi_zeros_out,
		writeheader = true)
end;

# ╔═╡ Cell order:
# ╟─a6314f18-60a5-11eb-0307-450764cbf160
# ╠═e53e9d5a-6091-11eb-3582-398334c98d77
# ╠═f1a368be-6091-11eb-0da3-77b12261da8c
# ╟─81a8c5a8-6092-11eb-3414-69625042ef17
# ╠═f6e0abde-6091-11eb-2a1c-010ab8a51452
# ╟─9509ec9e-6092-11eb-2157-a540547604cb
# ╠═a02127d2-6092-11eb-26fc-e749e38caaca
# ╟─e40e07ce-6094-11eb-1137-b5acb0fec21a
# ╠═c5044126-60b8-11eb-1d18-33c34995a8d5
# ╠═e5f44554-60a2-11eb-33ab-e199bec0c726

max_stars_repo_name: MTreppner/scDBM-paper

max_issues_repo_name: MTreppner/scDBM-paper

max_forks_repo_name: MTreppner/scDBM-paper


text:
Xc1, yc1 = gen_classif(classes=["M", "F"])
Xc2, yc2 = gen_classif(classes=["A", "B", "C"])
Xr, yr   = gen_reg()

@testset "DummyReg" begin
    m, f = simple_test_reg(DummyRegressor(), Xr, yr)
    fp = fitted_params(m, f)
    @test keys(fp) == (:constant, :n_outputs)
    infos = info_dict(m)
    @test infos[:input_scitype] == MLJBase.Table(MLJBase.Continuous)
    @test infos[:target_scitype] == AbstractVector{MLJBase.Continuous}
end

@testset "DummyClf" begin
    m, f = simple_test_classif_prob(DummyClassifier(), Xc2, yc2; nodummy=true)
    fp = fitted_params(m, f)
    @test keys(fp) == (:classes, :n_classes, :n_outputs)
    infos = info_dict(m)
    @test infos[:input_scitype] == MLJBase.Table(MLJBase.Continuous)
    @test infos[:target_scitype] == AbstractVector{<:MLJBase.Finite}
    @test !isempty(infos[:docstring])
end

@testset "GaussianNBClf" begin
    m, f = simple_test_classif_prob(GaussianNBClassifier(), Xc2, yc2)
    fp =  fitted_params(m, f)
    @test keys(fp) == (:class_prior, :class_count, :theta, :sigma, :epsilon)
    infos = info_dict(m)
    @test infos[:input_scitype] == MLJBase.Table(MLJBase.Continuous)
    @test infos[:target_scitype] == AbstractVector{<:MLJBase.Finite}
    @test !isempty(infos[:docstring])
end

@testset "KNNReg" begin
    m, f = simple_test_reg(KNeighborsRegressor(), Xr, yr)
    fp = fitted_params(m, f)
    @test keys(fp) == (:effective_metric, :effective_metric_params)
    infos = info_dict(m)
    @test infos[:input_scitype] == MLJBase.Table(MLJBase.Continuous)
    @test infos[:target_scitype] == AbstractVector{MLJBase.Continuous}
    @test !isempty(infos[:docstring])
end

@testset "KNNClf" begin
    m, f = simple_test_classif_prob(KNeighborsClassifier(), Xc2, yc2)
    fp = fitted_params(m, f)
    @test keys(fp) == (:classes, :effective_metric, :effective_metric_params, :outputs_2d)
    infos = info_dict(m)
    @test infos[:input_scitype] == MLJBase.Table(MLJBase.Continuous)
    @test infos[:target_scitype] == AbstractVector{<:MLJBase.Finite}
    @test !isempty(infos[:docstring])
end

X3 = 5 * abs.(MLJBase.matrix(Xc2))
X3 .= ceil.(X3)
Xc3 = MLJBase.table(Int.(X3))

@testset "BernNBClf" begin
    m, f = simple_test_classif_prob(BernoulliNBClassifier(binarize=2.0), Xc3, yc2; nodummy=true)
    fp =  fitted_params(m, f)
    @test keys(fp) == (:class_log_prior, :feature_log_prob, :class_count, :feature_count)
    infos = info_dict(m)
    @test infos[:input_scitype] == MLJBase.Table(MLJBase.Count)
    @test infos[:target_scitype] == AbstractVector{<:MLJBase.Finite}
    @test !isempty(infos[:docstring])
end

@testset "MultiNBClf" begin
    m, f = simple_test_classif_prob(MultinomialNBClassifier(), Xc3, yc2; nodummy=true)
    fp =  fitted_params(m, f)
    @test keys(fp) == (:class_log_prior, :intercept, :feature_log_prob, :coef, :class_count, :feature_count)
    infos = info_dict(m)
    @test infos[:input_scitype] == MLJBase.Table(MLJBase.Count)
    @test infos[:target_scitype] == AbstractVector{<:MLJBase.Finite}
    @test !isempty(infos[:docstring])
end

@testset "ComplNBClf" begin
    m, f = simple_test_classif_prob(ComplementNBClassifier(), Xc3, yc2; nodummy=true)
    fp =  fitted_params(m, f)
    @test keys(fp) == (:class_log_prior, :feature_log_prob, :class_count, :feature_count, :feature_all)
    infos = info_dict(m)
    @test infos[:input_scitype] == MLJBase.Table(MLJBase.Count)
    @test infos[:target_scitype] == AbstractVector{<:MLJBase.Finite}
    @test !isempty(infos[:docstring])
end


text:
include("dg_tree/dg.jl")
include("dg_curved/dg.jl")
include("dg_unstructured_quad/dg.jl")
include("dg_p4est/dg.jl")
include("dg_common.jl")
include("fdsbp_tree/fdsbp_2d.jl")


text:
@testset "PlutoUI" begin

    df = DataFrame(a=1:10, b=rand(10), c="hello".* string.(1:10), d=Bool.((1:10) .% 2), e=Date("2020-08-15") .+ Day.(1:10), f="world!" .* string.(1:10))

    empty_file_picker = Dict{Any, Any}("name" => "", "data" => UInt8[], "type" => "")
    @test_throws ErrorException read_table(empty_file_picker)

    filenames = ["test.csv", "test.xlsx", "test_array.json", "test.arrow", "test.zip", "test3.zip"]
    # "test_obj.json" not checked here - strange behaviour with Julia 1.6 nightly. To be checked after Julia 1.6 release.

    for filename ∈ filenames
        @info "testing import of file $filename from PlutoUI.FilePicker"
        file_picker = Dict{Any, Any}("name" => filename,
            "data" => read(joinpath(testpath, filename)),
            "type" => "")
        df_recovered = DataFrame(read_table(file_picker); copycols=false)
        df_recovered[!, :e] = Date.(df_recovered.e) # some impots do not convert date columns automatically

        @test df_recovered == df
    end

end


text:
__precompile__()
module Gallium
    using ASTInterpreter
    using Base.Meta
    using DWARF
    using ObjFileBase
    using ELF
    using MachO
    using COFF
    using AbstractTrees
    import ASTInterpreter: @enter

    # Debugger User Interface
    export breakpoint, @enter, @breakpoint, @conditional

    include("remote.jl")
    include("registers.jl")
    include("x86_64/registers.jl")
    include("x86_32/registers.jl")
    include("win64seh.jl")
    include("modules.jl")
    include("unwind.jl")
    include("Hooking/Hooking.jl")
    include("ptrace.jl")

    using .Registers
    using .Registers: ip, get_dwarf
    using .Hooking

    type JuliaStackFrame
        oh
        file
        ip
        sstart
        line::Int
        linfo::LambdaInfo
        variables::Dict
        env::Environment
    end

    type CStackFrame
        ip::Ptr{Void}
        file::AbstractString
        line::Int
        declfile::AbstractString
        declline::Int
        stacktop::Bool
    end

    # Fake "Interpreter" that is just a native stack
    immutable NativeStack
        stack
        RCs
        modules
        session
    end
    NativeStack(stack::Vector) = NativeStack(stack,Any[],active_modules,LocalSession())

    function ASTInterpreter.language_specific_prompt(state, stack::JuliaStackFrame)
        state.julia_prompt
    end
    function ASTInterpreter.language_specific_prompt(state, stack::NativeStack)
        ASTInterpreter.language_specific_prompt(state, stack.stack[end])
    end

    ASTInterpreter.done!(stack::NativeStack) = nothing

    function ASTInterpreter.print_frame(_, io, num, x::JuliaStackFrame)
        print(io, "[$num] ")
        linfo = x.linfo
        ASTInterpreter.print_linfo_desc(io, linfo)
        println(io)
        ASTInterpreter.print_locals(io, linfo, x.env, (io,name)->begin
            if haskey(x.variables, name)
                if x.variables[name] == :available
                    println(io, "<undefined>")
                else
                    println(io, "<not available here>")
                end
            else
                println(io, "<optimized out>")
            end
        end)
    end
    ASTInterpreter.print_frame(_, io, num, x::NativeStack) =
        ASTInterpreter.print_frame(_, io, num, x.stack[end])

    function ASTInterpreter.print_status(state, x::JuliaStackFrame; kwargs...)
        if x.line < 0
            println("Got a negative line number. Bug?")
        elseif (!isa(x.file,AbstractString) || isempty(x.file)) || x.line == 0
            println("<No file found. Did DWARF parsing fail?>")
        else
            dfile = Base.find_source_file(string(x.file))
            ASTInterpreter.print_sourcecode(x.linfo,
                ASTInterpreter.readfileorhist(
                    dfile == nothing ? x.file : dfile), x.line)
        end
    end

    function ASTInterpreter.print_status(state, x::NativeStack; kwargs...)
        ASTInterpreter.print_status(state, x.stack[end]; kwargs...)
    end

    function ASTInterpreter.get_env_for_eval(x::JuliaStackFrame)
        copy(x.env)
    end
    ASTInterpreter.get_env_for_eval(x::NativeStack) =
        ASTInterpreter.get_env_for_eval(x.stack[end])
    ASTInterpreter.get_linfo(x::JuliaStackFrame) = x.linfo
    ASTInterpreter.get_linfo(x::NativeStack) =
        ASTInterpreter.get_linfo(x.stack[end])

    const GalliumFrame = Union{NativeStack, JuliaStackFrame, CStackFrame}
    using DWARF: CallFrameInfo
    function ASTInterpreter.execute_command(state, x::GalliumFrame, ::Val{:cfi}, command)
        modules = state.top_interp.modules
        mod, base, ip = modbaseip_for_stack(state, x)
        modrel = UInt(ip - base)
        loc, fde = Unwinder.find_fde(mod, modrel)
        cie = realize_cie(fde)
        target_delta = modrel - loc
        out = IOContext(STDOUT, :reg_map =>
            isa(Gallium.getarch(state.top_interp.session),Gallium.X86_64.X86_64Arch) ?
            Gallium.X86_64.dwarf_numbering : Gallium.X86_32.dwarf_numbering)
        drs = CallFrameInfo.RegStates()
        CallFrameInfo.dump_program(out, cie, ptrT = ObjFileBase.intptr(ObjFileBase.handle(mod)),
            target = UInt(target_delta), rs = drs); println(out)
        CallFrameInfo.dump_program(out, fde,
            cie = cie, target = UInt(target_delta), rs = drs)
        return false
    end

    function obtain_linetable(state, stack)
        mod, base, ip = modbaseip_for_stack(state, stack)
        dbgs = debugsections(dhandle(mod))
        lip = compute_ip(dhandle(mod), base, ip)
        cu = DWARF.searchcuforip(dbgs, lip)
        line_offset = get(DWARF.extract_attribute(cu, DWARF.DW_AT_stmt_list))
        seek(dbgs.debug_line, UInt(line_offset.value))
        DWARF.LineTableSupport.LineTable(handle(dbgs).io), lip
    end

    function ASTInterpreter.execute_command(state, stack::GalliumFrame, ::Val{:lip}, command)
        mod, base, ip = modbaseip_for_stack(state, stack)
        dbgs = debugsections(dhandle(mod))
        lip = compute_ip(dhandle(mod), base, ip)
        show(STDOUT, lip); println(STDOUT)
        return false
    end

    function ASTInterpreter.execute_command(state, stack::GalliumFrame, ::Union{Val{:linetabprog}, Val{:linetab}}, command)
        linetab = obtain_linetable(state, stack)[1]
        (command == Val{:linetabprog} ? DWARF.LineTableSupport.dump_program :
            DWARF.LineTableSupport.dump_table)(STDOUT, linetab)
        return false
    end

    function find_module(session, modules, frame::CStackFrame)
        ip = frame.ip - (frame.stacktop ? 0 : 1)
        base, mod = find_module(session, modules, ip)
    end

    function modbaseip_for_stack(state, stack)
        modules = state.top_interp.modules
        frame = isa(stack, NativeStack) ? stack.stack[end] : stack
        session = isa(state.top_interp, NativeStack) ? state.top_interp.session : LocalSession()
        ip = frame.ip - (frame.stacktop ? 0 : 1)
        base, mod = find_module(session, modules, ip)
        mod, base, ip
    end

    function ASTInterpreter.execute_command(state, stack::GalliumFrame, ::Val{:handle}, command)
        Base.LineEdit.transition(state.s, :abort)
        mod, _, __ = modbaseip_for_stack(state, stack)
        eval(Main,:(h = $(handle(mod))))
        return false
    end

    function compute_ip(handle, base, theip)
        if isa(handle, ELF.ELFHandle)
            phs = ELF.ProgramHeaders(handle)
            idx = findfirst(p->p.p_offset==0&&p.p_type==ELF.PT_LOAD, phs)
            UInt(theip + (phs[idx].p_vaddr-base))
        elseif isa(handle, COFF.COFFHandle)
            # Map from loaded address to file address
            UInt(theip) + (UInt(COFF.readoptheader(handle).windows.ImageBase) - UInt(base))
        elseif ObjFileBase.isrelocatable(handle)
            UInt(theip)
        end
    end

    function compute_symbol_value(handle, base, symbol)
        if isa(handle, COFF.COFFHandle)
            value = COFF.symbolvalue(symbol, COFF.Sections(ObjFileBase.handle(handle)))
            # Map from file address to loaded address
            return UInt(value) + (Int(base) - Int(COFF.readoptheader(h).windows.ImageBase))
        else
            value = ObjFileBase.symbolvalue(symbol, ObjFileBase.Sections(ObjFileBase.handle(handle)))
            return base + value
        end
    end

    function ASTInterpreter.execute_command(state, stack::GalliumFrame,
            cmd::Union{Val{:cu},Val{:sp}}, command)
        mod, base, ip = modbaseip_for_stack(state, stack)
        dbgs = debugsections(dhandle(mod))
        lip = compute_ip(dhandle(mod), base, ip)
        unit = DWARF.searchcuforip(dbgs, UInt(lip))
        (cmd == Val{:sp}()) && (unit = DWARF.searchspforip(unit, UInt(lip)))
        AbstractTrees.print_tree(show, IOContext(STDOUT,:strtab=>StrTab(dbgs.debug_str)), unit)
        return false
    end

    function ASTInterpreter.execute_command(state, stack::GalliumFrame,
            cmd::Val{:seh}, command)
        mod, base, ip = modbaseip_for_stack(state, stack)
        entry = Unwinder.find_seh_entry(mod, ip-base)
        i = 1
        while i <= length(entry.opcodes)
            i += print_op(STDOUT, entry.opcodes, i)
        end
        return false
    end

    # Use this hook to expose extra functionality
    function ASTInterpreter.execute_command(x::JuliaStackFrame, command)
        lip = UInt(x.ip)-x.sstart-1
        if isrelocatable(handle(x.oh))
            lip = UInt(x.ip)-1
        end
        if command == "ip"
            println(x.ip)
            return
        elseif command == "sstart"
            println(x.sstart)
            return
        elseif startswith(command, "bp")
            subcmds = split(command,' ')[2:end]
            if subcmds[1] == "list"
                list_breakpoints()
            elseif subcmds[1] == "disable"
                bp = breakpoints[parse(Int,subcmds[2])]
                disable(bp)
                println(bp)
            elseif subcmds[1] == "enable"
                bp = breakpoints[parse(Int,subcmds[2])]
                enable(bp)
                println(bp)
            end
        elseif startswith(command, "b")
            nothing
        end
    end

    function ASTInterpreter.execute_command(x::NativeStack, command)
        ASTInterpreter.execute_command(x.stack[end], command)
    end

    export breakpoint

    # Move this somewhere better
    function ObjFileBase.getSectionLoadAddress(LOI::Dict, sec)
        return LOI[Symbol(sectionname(sec))]
    end

    function search_linetab(linetab, ip)
        local last_entry
        first = true
        for entry in linetab
            if !first
                if entry.address > reinterpret(UInt64,ip)
                    return last_entry
                end
            end
            first = false
            last_entry = entry
        end
        last_entry
    end

    function rec_backtrace(RC)
        ips = Array(UInt64, 0)
        rec_backtrace(RC->(push!(ips,ip(RC)); true), RC)
        ips
    end

    function rec_backtrace_hook(RC)
        ips = Array(UInt64, 0)
        rec_backtrace_hook(RC->(push!(ips,ip(RC)); true), RC)
        ips
    end

    global active_modules = LazyJITModules()
    global allow_bad_unwind = true
    function rec_backtrace(callback, RC, session = LocalSession(), modules = active_modules, ip_only = false, cfi_cache = nothing; stacktop = true)
        callback(RC) || return
        while true
            (ok, RC) = try
                Unwinder.unwind_step(session, modules, RC, cfi_cache;
                    stacktop = stacktop, ip_only = ip_only, allow_frame_based=allow_bad_unwind)
            catch e # e.g. invalid memory access, invalid unwind info etc.
                if !allow_bad_unwind::Bool
                    println("Was at ip: ", ip(RC))
                    rethrow(e)
                end
                break
            end
            stacktop = false
            ok || break
            callback(RC) || break
        end
    end

    function step_first!(RC)
        set_ip!(RC,unsafe_load(convert(Ptr{UInt},get_dwarf(RC,:rsp)[])))
        set_sp!(RC,get_dwarf(RC,:rsp)[]+sizeof(Ptr{Void}))
    end

    function rec_backtrace_hook(callback, RC, session = LocalSession(), modules = active_modules, ip_only = false)
        callback(RC) || return
        step_first!(RC)
        rec_backtrace(callback, RC, session, modules, ip_only; stacktop = false)
    end

    # Validate that an address is a valid location in the julia heap
    function heap_validate(ptr)
        typeptr = Ptr{Ptr{Void}}(UInt64(ptr)-sizeof(Ptr{Void}))
        Hooking.mem_validate(typeptr,sizeof(Ptr)) || return false
        T = UInt(unsafe_load(typeptr))&(~UInt(0x3))
        typetypeptr = Ptr{Ptr{Void}}(T-sizeof(Ptr))
        Hooking.mem_validate(typetypeptr,sizeof(Ptr)) || return false
        UInt(unsafe_load(typetypeptr))&(~UInt(0x3)) == UInt(pointer_from_objref(DataType))
    end

    function get_ipinfo(session::LocalSession, theip)
        Base.StackTraces.lookup(theip)
    end

    function get_ipinfo(session, theip)
        sf = StackFrame(:unknown, "", 0, Nullable{LambdaInfo}(),
            true, false, theip)
        [sf]
    end

    function getreg(RC, fbreg, reg)
        if reg == DWARF.DW_OP_fbreg
            (fbreg == -1) && error("fbreg requested but not found")
            Gallium.get_dwarf(RC, fbreg)
        else
            Gallium.get_dwarf(RC, reg)
        end
    end

    function iterate_variables(RC, found_cb, not_found_cb, dbgs, lip)
        cu = DWARF.searchcuforip(dbgs, lip)
        sp = DWARF.searchspforip(cu, lip)

        culow = DWARF.extract_attribute(cu,DWARF.DW_AT_low_pc)

        fbreg = DWARF.extract_attribute(sp, DWARF.DW_AT_frame_base)
        # Array is for DWARF 2 support.
        fbreg = isnull(fbreg) ? -1 :
            (isa(get(fbreg).value, Array) ? get(fbreg).value[1] : get(fbreg).value.expr[1]) - DWARF.DW_OP_reg0

        local getreg
        getreg = (reg)->Gallium.getreg(RC, fbreg, reg)
        getword(addr) = unsafe_load(reinterpret(Ptr{UInt64}, addr))
        addr_func(x) = x

        for vardie in (filter(children(sp)) do child
                    tag = DWARF.tag(child)
                    tag == DWARF.DW_TAG_formal_parameter ||
                    tag == DWARF.DW_TAG_variable
                end)
            name = DWARF.extract_attribute(vardie,DWARF.DW_AT_name)
            loc = DWARF.extract_attribute(vardie,DWARF.DW_AT_location)
            (isnull(name) || isnull(loc)) && continue
            name = Symbol(bytestring(get(name).value,StrTab(dbgs.debug_str)))
            loc = get(loc)
            if loc.spec.form == DWARF.DW_FORM_exprloc
                sm = DWARF.Expressions.StateMachine{typeof(loc.value).parameters[1]}()
                val = DWARF.Expressions.evaluate_simple_location(
                    sm, loc.value.expr, getreg, getword, addr_func, :NativeEndian)
                found_cb(dbgs, vardie, getreg, name, val)
            elseif loc.spec.form == DWARF.DW_FORM_sec_offset
                T = UInt64
                seek(dbgs.debug_loc, loc.value.offset)
                list = read(dbgs.debug_loc, DWARF.LocationList{T})
                sm = DWARF.Expressions.StateMachine{T}()
                found = false
                for entry in list.entries
                    if entry.first <= lip - UInt(culow.value) <= entry.last
                        val = DWARF.Expressions.evaluate_simple_location(
                            sm, entry.data, getreg, getword, addr_func, :NativeEndian)
                        found = true
                        found_cb(dbgs, vardie, getreg, name, val)
                        break
                    end
                end
                found || not_found_cb(dbgs, vardie, name)
            else
                not_found_cb(dbgs, vardie, name)
            end
        end
    end

    function make_base_stackframe(session, modules, theip; firstframe = false)
        sstart, h = try
            find_module(session, modules, theip)
        catch # Unwind got it wrong, but still include at least one stack frame
            return Base.StackFrame(Symbol("???"),Symbol("???"),0,
                Nullable{LambdaInfo}(),true,false,theip)
        end
        dh = dhandle(h)
        lip = compute_ip(dh, sstart, theip-(firstframe?0:1))
        dbgs = debugsections(dh)
        file, line, linetab = linetable_entry(dbgs, lip)
        name = symbolicate(session, modules, theip-(firstframe?0:1))
        Base.StackFrame(Symbol(name),Symbol(file),line,Nullable{LambdaInfo}(),
            true, false, theip)
    end

    function linetable_entry(dbgs, lip)
        cu = DWARF.searchcuforip(dbgs, lip)
        sp = DWARF.searchspforip(cu, lip)
        # Process Compilation Unit to get line table
        line_offset = DWARF.extract_attribute(cu, DWARF.DW_AT_stmt_list)
        line_offset = isnull(line_offset) ? 0 : convert(UInt, get(line_offset).value)

        seek(dbgs.debug_line, line_offset)
        linetab = DWARF.LineTableSupport.LineTable(handle(dbgs).io)
        entry = search_linetab(linetab, lip)
        line = entry.line
        fileentry = linetab.header.file_names[entry.file]
        if fileentry.dir_idx == 0
            found_file = Base.find_source_file(fileentry.name)
            file = found_file != nothing ? found_file : fileentry.name
        else
            file = joinpath(linetab.header.include_directories[fileentry.dir_idx],
                fileentry.name)
        end
        file, line, linetab
    end

    function frameinfo(RC, session, modules; rich_c = false, firstframe = false)
        theip = reinterpret(Ptr{Void},UInt(ip(RC)))
        ipinfo = get_ipinfo(session, theip)[end]
        fromC = ipinfo.from_c
        file = ""
        line = 0
        local declfile="", declline=0
        if fromC
            sstart, h = try
                find_module(session, modules, theip)
            catch err # Unwind got it wrong, but still include at least one stack frame
                allow_bad_unwind::Bool || rethrow(err)
                return false, Nullable(CStackFrame(theip, file, line, declfile, declline, firstframe))
            end
            if rich_c
                dh = dhandle(h)
                lip = compute_ip(dh, sstart, theip-(firstframe?0:1))
                local dbgs=nothing, cu=nothing, sp=nothing

                try
                    dbgs = debugsections(dh)
                    cu = DWARF.searchcuforip(dbgs, lip)
                    sp = DWARF.searchspforip(cu, lip)
                catch
                end
                if sp !== nothing
                    file, line, linetab = linetable_entry(dbgs, lip)

                    declfile = DWARF.extract_attribute(sp, DWARF.DW_AT_decl_file)
                    declfile = isnull(declfile) ? "" : linetab.header.file_names[convert(UInt,get(declfile))].name
                    declline = DWARF.extract_attribute(sp, DWARF.DW_AT_decl_line)
                    declline = isnull(declline) ? 0 : convert(UInt, get(declline))
                end
            end
            return (true, Nullable(CStackFrame(theip, file, line, declfile, declline, firstframe)))
        else
            sstart, h = find_module(session, modules, theip)
            isnull(ipinfo.linfo) && return (false, Nullable{JuliaStackFrame}())
            tlinfo = get(ipinfo.linfo)
            env = ASTInterpreter.prepare_locals(tlinfo)
            copy!(env.sparams, tlinfo.sparam_vals)
            slottypes = tlinfo.slottypes
            tlinfo = tlinfo.def.lambda_template
            variables = Dict()
            dbgs = debugsections(dhandle(h))
            try
                lip = UInt(theip)-sstart-1
                if isrelocatable(handle(h))
                    lip = UInt(theip)-1
                end

                file, line, _ = linetable_entry(dbgs, lip)

                # Process Subprogram to extract local variables
                strtab = ObjFileBase.load_strtab(dbgs.debug_str)
                vartypes = Dict{Symbol,Type}()
                if slottypes === nothing
                    for name in tlinfo.slotnames
                        vartypes[name] = Any
                    end
                else
                    for (name, ty) in zip(tlinfo.slotnames, slottypes)
                        vartypes[name] = ty
                    end
                end
                variables = Dict()
                function found_cb(dbgs, vardie, getreg, name, val)
                    haskey(vartypes, name) || return
                    variables[name] = :found
                    if isa(val, DWARF.Expressions.MemoryLocation)
                        if isbits(vartypes[name])
                            ptr = reinterpret(Ptr{Void}, val.i)
                            if ptr != C_NULL && heap_validate(ptr)
                                val = unsafe_load(reinterpret(Ptr{vartypes[name]},
                                    ptr))
                            else
                                val = Nullable{vartypes[name]}()
                            end
                        else
                            ptr = reinterpret(Ptr{Void}, val.i)
                            if ptr == C_NULL
                                val = Nullable{Ptr{vartypes[name]}}()
                            elseif heap_validate(ptr)
                                val = unsafe_pointer_to_objref(ptr)
                            # This is a heuristic. Should update to check
                            # whether the variable is declared as jl_value_t
                            elseif Hooking.mem_validate(ptr, sizeof(Ptr{Void}))
                                ptr2 = unsafe_load(Ptr{Ptr{Void}}(ptr))
                                if ptr2 == C_NULL
                                    val = Nullable{Ptr{vartypes[name]}}()
                                elseif heap_validate(ptr2)
                                    val = unsafe_pointer_to_objref(ptr2)
                                end
                            end
                        end
                        variables[name] = :available
                    elseif isa(val, DWARF.Expressions.RegisterLocation)
                        # The value will generally be in the low bits of the
                        # register. This should give the appropriate value
                        val = getreg(val.i)[]
                        if !isbits(vartypes[name])
                            if val == 0
                                val = Nullable{vartypes[name]}()
                            elseif heap_validate(val)
                                val = unsafe_pointer_to_objref(Ptr{Void}(val))
                            else
                                return
                            end
                        else
                            val = reinterpret(vartypes[name],[val])[]
                        end
                        variables[name] = :available
                    end
                    varidx = findfirst(tlinfo.slotnames, name)
                    if varidx != 0
                        env.locals[varidx] = Nullable{Any}(val)
                    end
                end
                iterate_variables(RC,found_cb,(dbgs, vardie, name)->variables[name] = :found,dbgs,lip)
            catch err
                if !isa(err, ErrorException) || err.msg != "Not found"
                    @show err
                    Base.show_backtrace(STDOUT, catch_backtrace())
                end
            end
            # process SP.variables here
            #h = readmeta(IOBuffer(data))
            return (true, Nullable(JuliaStackFrame(h, file, UInt(ip(RC)), sstart, line, get(ipinfo.linfo), variables, env)))
        end
    end

    function stackwalk(RC, session = LocalSession(), modules = active_modules;
            fromhook = false, rich_c = false, ip_only = false, collectRCs=false)
        stack = Any[]
        RCs = Any[]
        firstframe = true
        (fromhook ? rec_backtrace_hook : rec_backtrace)(RC, session, modules, ip_only) do RC
            keep_walking, frame = frameinfo(RC, session, modules, rich_c = rich_c, firstframe = firstframe)
            !isnull(frame) && (push!(stack, get(frame));
                collectRCs && push!(RCs, RC))
            firstframe = false
            return keep_walking
        end
        (reverse!(stack), reverse!(RCs))
    end

    function matches_condition(interp, condition)
        condition == nothing && return true
        if isa(condition, Expr)
            ok, res = ASTInterpreter.eval_in_interp(interp, condition)
            !ok && println("Conditional breakpoint errored. Breaking.")
            return !ok || res
        else
            error("Unexpected condition kind")
        end
    end

    function process_lowlevel_conditionals(loc, RC)
        !haskey(bps_at_location, loc) && return true
        stop = true
        for bp in bps_at_location[loc]
            for cond in bp.conditions
                isa(cond, Function) && cond(loc, RC) && return true
                # Will get evaluated later
                isa(cond, Expr) && return true
                stop = false
            end
        end
        stop
    end

    macro conditional(bp, condition)
        esc(:(let bp = $bp
            push!(bp.conditions, $(Expr(:quote, condition)))
            bp
        end))
    end

    function conditional(f, bp)
        push!(bp.conditions, f)
        bp
    end

    """
        Suspends other tasks that may want to access the terminal.
        In an out of process debugger, this would be the place
        where we take control of the tty for the debugger. Since
        we're not out-of-process, we instead just prevent those
        tasks from being notified of any reads by modifying internal
        datastructures.
    """
    function suspend_other_tasks()
        STDINwaitq = STDIN.readnotify.waitq
        GlobalWaitq = copy(Base.Workqueue)
        STDIN.readnotify.waitq = Any[]
        empty!(Base.Workqueue)
        STDINwaitq, GlobalWaitq
    end
    restore_other_tasks(state) = (STDIN.readnotify.waitq = state[1]; append!(Base.Workqueue, state[2]))

    function breakpoint_hit(hook, RC)
        if !process_lowlevel_conditionals(Location(LocalSession(), hook.addr), RC)
            return
        end
        stack = stackwalk(RC; fromhook = true)[1]
        stacktop = pop!(stack)
        linfo = stacktop.linfo
        fT = linfo.def.sig.parameters[1]
        def_linfo = linfo.def.lambda_template
        argnames = [Symbol("#target_self#");def_linfo.slotnames[2:def_linfo.nargs]]
        spectypes = [fT;linfo.specTypes.parameters[2:end]...]
        bps = bps_at_location[Location(LocalSession(),hook.addr)]
        target_line = minimum(map(bps) do bp
            idx = findfirst(s->isa(s, FileLineSource), bp.sources)
            idx != 0 ? bp.sources[idx].line : def_linfo.def.line
        end)
        conditions = reduce(vcat,map(bp->bp.conditions, bps))
        thunk = Expr(:->,Expr(:tuple,map(x->Expr(:(::),x[1],x[2]),zip(argnames,spectypes))...),Expr(:block,
            :(linfo = $(quot(linfo))),
            :((loctree, code) = ASTInterpreter.reparse_meth(linfo)),
            :(__env = ASTInterpreter.prepare_locals(linfo.def.lambda_template)),
            :(copy!(__env.sparams, linfo.sparam_vals)),
            [ :(__env.locals[$i] = Nullable{Any}($(argnames[i]))) for i = 1:length(argnames) ]...,
            :(interp = ASTInterpreter.enter(linfo,__env,
                $(collect(filter(x->!isa(x,CStackFrame),stack)));
                    loctree = loctree, code = code)),
            (target_line != linfo.def.line ?
                :(ASTInterpreter.advance_to_line(interp, $target_line)) :
                :(nothing)),
            :(tty_state = Gallium.suspend_other_tasks()),
            :((isempty($conditions) ||
                any(c->Gallium.matches_condition(interp,c),$conditions)) &&
                ASTInterpreter.RunDebugREPL(interp)),
            :(Gallium.restore_other_tasks(tty_state)),
            :(ASTInterpreter.finish!(interp)),
            :(return interp.retval::$(linfo.rettype))))
        f = eval(thunk)
        t = Tuple{spectypes...}
        faddr = Hooking.get_function_addr(f, t)
        Hooking.Deopt(faddr)
    end
    abstract LocationSource
    immutable Location
        vm
        addr::UInt64
    end
    type Breakpoint
        active_locations::Vector{Location}
        inactive_locations::Vector{Location}
        sources::Vector{LocationSource}
        disable_new::Bool
        conditions::Vector{Any}
    end
    Breakpoint(locations::Vector{Location}) = Breakpoint(locations, Location[], LocationSource[], false, Any[])
    Breakpoint() = Breakpoint(Location[], Location[], LocationSource[], false, Any[])

    function print_location(io::IO, vm::LocalSession, loc)
        ipinfo = Base.StackTraces.lookup(loc.addr+1)[end]
        if ipinfo.from_c
            println(io, "At address ", loc.addr)
        else
            linfo = get(ipinfo.linfo)
            ASTInterpreter.print_linfo_desc(io, linfo, true)
            println(io)
        end
    end

    function print_locations(io::IO, locations, prefix = " - ")
        for loc in locations
            print(io,prefix)
            print_location(io, loc.vm, loc)
        end
    end

    function Base.show(io::IO, b::Breakpoint)
        if isempty(b.active_locations) && isempty(b.inactive_locations) &&
            isempty(b.sources)
            println(io, "Empty Breakpoint")
            return
        end
        println(io, "Locations (+: active, -: inactive, *: source):")
        !isempty(b.active_locations) && print_locations(io, b.active_locations, " + ")
        !isempty(b.inactive_locations) && print_locations(io, b.inactive_locations, " - ")
        for source in b.sources
            print(io," * ")
            println(io,source)
        end
    end

    const bps_at_location = Dict{Location, Set{Breakpoint}}()
    # session => (addr => loc)
    disable(s::LocalSession, loc) = unhook(Ptr{Void}(loc.addr))
    disable(loc::Location) = disable(loc.vm, loc)
    function disable(bp::Breakpoint, loc::Location)
        pop!(bps_at_location[loc],bp)
        if isempty(bps_at_location[loc])
            disable(loc)
            delete!(bps_at_location,loc)
        end
    end
    function disable(b::Breakpoint)
        locs = copy(b.active_locations)
        empty!(b.active_locations)
        for loc in locs
            disable(b, loc)
            push!(b.inactive_locations, loc)
        end
        b.disable_new = true
    end
    remove(b::Breakpoint) = (disable(b); deleteat!(breakpoints, findfirst(breakpoints, b)); nothing)

    enable(s::LocalSession, loc) = hook(breakpoint_hit, Ptr{Void}(loc.addr))
    enable(loc::Location) = enable(loc.vm, loc)
    function enable(bp::Breakpoint, loc::Location)
        if !haskey(bps_at_location, loc)
            enable(loc)
            bps_at_location[loc] = Set{Breakpoint}()
        end
        push!(bps_at_location[loc], bp)
    end
    function enable(b::Breakpoint)
        locs = copy(b.inactive_locations)
        empty!(b.inactive_locations)
        for loc in locs
            enable(b, loc)
            push!(b.active_locations, loc)
        end
        b.disable_new = false
    end

    function _breakpoint_spec(spec::LambdaInfo, bp)
        llvmf = ccall(:jl_get_llvmf, Ptr{Void}, (Any, Bool, Bool), spec.specTypes, false, true)
        @assert llvmf != C_NULL
        fptr = ccall(:jl_get_llvm_fptr, UInt64, (Ptr{Void},), llvmf)
        @assert fptr != 0
        loc = Location(LocalSession(), fptr)
        add_location(bp, loc)
    end

    function _breakpoint_method(meth::Method, bp::Breakpoint, predicate = linfo->true)
        cache = meth.sig.parameters[1].name.mt.cache
        cache != nothing || return
        Base.visit(cache) do spec
            spec.def == meth || return
            predicate(spec) || return
            _breakpoint_spec(spec, bp)
        end
    end

    type SpecSource <: LocationSource
        bp::Breakpoint
        meth::Method
        predicate
        function SpecSource(bp::Breakpoint, meth::Method, predicate)
            !haskey(TracedMethods, meth) && (TracedMethods[meth] = Set{SpecSource}())
            ccall(:jl_trace_method, Void, (Any,), meth)
            this = new(bp, meth, predicate)
            push!(TracedMethods[meth], this)
            finalizer(this,function (this)
                pop!(TracedMethods[this.meth], this)
                if isempty(TracedMethods[this.meth])
                    ccall(:jl_untrace_method, Void, (Any,), this.meth)
                    delete!(TracedMethods, this.meth)
                end
            end)
            this
        end
    end
    function fire(s::SpecSource, linfo::LambdaInfo)
        s.predicate(linfo) || return
        _breakpoint_spec(linfo, s.bp)
    end

    const TracedMethods = Dict{Method, Set{SpecSource}}()
    function Base.show(io::IO, source::SpecSource)
        print(io,"Any matching specialization of ")
        ASTInterpreter.print_linfo_desc(io, source.meth.lambda_template, true)
    end

    function rebreak_tracer(x::Ptr{Void})
        linfo = unsafe_pointer_to_objref(x)::LambdaInfo
        !haskey(TracedMethods, linfo.def) && return nothing
        for s in TracedMethods[linfo.def]
            fire(s, linfo)
        end
        nothing
    end

    function add_meth_to_bp!(bp::Breakpoint, meth::Union{Method, TypeMapEntry}, predicate = linfo->true)
        isa(meth, TypeMapEntry) && (meth = meth.func)
        _breakpoint_method(meth, bp, predicate)
        push!(bp.sources, SpecSource(bp, meth, predicate))
        bp
    end

    function breakpoint(meth::Union{Method, TypeMapEntry})
        bp = add_meth_to_bp!(Breakpoint(), meth)
        push!(breakpoints, bp)
        bp
     end

    const breakpoints = Vector{Breakpoint}()

    function list_breakpoints()
        for (i, bp) in enumerate(breakpoints)
            println("[$i] $bp")
        end
    end

    function breakpoint(addr::Ptr{Void})
        hook(breakpoint_hit, addr)
    end

    function add_location(bp, loc)
        if bp.disable_new
            push!(bp.inactive_locations, loc)
        else
            push!(bp.active_locations, loc)
            enable(bp, loc)
        end
    end

    function _breakpoint_concrete(bp, t)
        addr = try
            Hooking.get_function_addr(t)
        catch err
            error("no method found for the specified argument types")
        end
        add_location(bp, Location(LocalSession(),addr))
    end

    function breakpoint(func, args::Union{Tuple,DataType})
        argtt = Base.to_tuple_type(args)
        t = Tuple{typeof(func), argtt.parameters...}
        bp = Breakpoint()
        if Base.isleaftype(t)
            _breakpoint_concrete(bp, t)
        else
            spec_predicate(linfo) = linfo.specTypes <: t
            meth_predicate(meth) = t <: meth.lambda_template.specTypes || meth.lambda_template.specTypes <: t
            for meth in methods(func, argtt)
                add_meth_to_bp!(bp, meth, spec_predicate)
            end
            push!(bp.sources, MethSource(bp, typeof(func), meth_predicate, spec_predicate))
        end
        push!(breakpoints, bp)
        bp
    end

    include("breakfile.jl")

    function method_tracer(x::Ptr{Void})
        ccall(:jl_trace_linfo, Void, (Ptr{Void},), x)
        nothing
    end

    function __init__()
        ccall(:jl_register_linfo_tracer, Void, (Ptr{Void},), cfunction(rebreak_tracer,Void,(Ptr{Void},)))
        ccall(:jl_register_method_tracer, Void, (Ptr{Void},), cfunction(method_tracer,Void,(Ptr{Void},)))
        ccall(:jl_register_newmeth_tracer, Void, (Ptr{Void},), cfunction(newmeth_tracer, Void, (Ptr{Void},)))
        update_shlibs!(LocalSession(), active_modules)
    end

    function breakpoint(f)
        bp = Breakpoint()
        for meth in methods(f)
            add_meth_to_bp!(bp, meth)
        end
        unshift!(bp.sources, MethSource(bp, typeof(f)))
        push!(breakpoints, bp)
        bp
    end

    macro breakpoint(ex0)
        Base.gen_call_with_extracted_types(:(Gallium.breakpoint),ex0)
    end

    # For now this is a very simple implementation. A better implementation
    # would trap and reuse logic. That will become important once we actually
    # support optimized code to avoid cloberring registers. For now do the dead
    # simple, stupid thing.
    function breakpoint()
        RC = Hooking.getcontext()
        # -1 to skip breakpoint (getcontext is inlined)
        stack, RCs = stackwalk(RC; fromhook = false)
        ASTInterpreter.RunDebugREPL(NativeStack(collect(filter(x->isa(x,JuliaStackFrame),stack[1:end-1]))))
    end

    const bp_on_error_conditions = Any[]
    function breakpoint_on_error_hit(thehook, RC)
        unhook(thehook)
        err = unsafe_pointer_to_objref(Ptr{Void}(get_dwarf(RC, :rdi)[]))
        stack = stackwalk(RC; fromhook = true)[1]
        stack = NativeStack(filter(x->isa(x,JuliaStackFrame),stack))
        if !isempty(bp_on_error_conditions) &&
            !any(c->Gallium.matches_condition(stack,c),bp_on_error_conditions)
            hook(thehook)
            rethrow(err)
        end
        ips = [x.ip-1 for x in stack]
        Base.with_output_color(:red, STDERR) do io
            print(io, "ERROR: ")
            Base.showerror(io, err, reverse(ips); backtrace=false)
            println(io)
        end
        println(STDOUT)
        ASTInterpreter.RunDebugREPL()
        # Get a somewhat sensible backtrace when returning
        try; throw(); catch; end
        hook(thehook)
        rethrow(err)
    end

    # Compiling these function has an error thrown/caught in type inference.
    # Precompile them here, to make sure we make it throught
    precompile(breakpoint_on_error_hit,(Hooking.Hook,X86_64.BasicRegs))
    precompile(Hooking.callback,(Ptr{Void},))

    function breakpoint_on_error(enable = true)
        addr = cglobal(:jl_throw)
        if enable
            hook(breakpoint_on_error_hit, addr)
        else
            unhook(addr)
        end
    end

    include("precompile.jl")
end


text:
import LinearAlgebra: I

"""
    get_data(params::Dict{String, Any}; eliminate_identical_gates = true)

Given the user input `params` dictionary, this function returns a dictionary of processed data which contains all the
necessary information to formulate the optimization model for the circuit design problem.
"""
function get_data(params::Dict{String, Any}; eliminate_identical_gates = true)

    # Number of qubits
    if "num_qubits" in keys(params)
        if params["num_qubits"] < 2
            Memento.error(_LOGGER, "Minimum of 2 qubits is necessary")
        end
        num_qubits = params["num_qubits"]
    else
        Memento.error(_LOGGER, "Number of qubits has to be specified by the user")
    end

    # Depth
    if "maximum_depth" in keys(params)
        if params["maximum_depth"] < 2
            Memento.error(_LOGGER, "Minimum depth of 2 is necessary")
        end
        maximum_depth = params["maximum_depth"]
    else
        Memento.error(_LOGGER, "Maximum depth of the decomposition has to be specified by the user")
    end

    # Elementary gates
    if !("elementary_gates" in keys(params)) || isempty(params["elementary_gates"])
        Memento.error(_LOGGER, "Specify at least two unique unitary elementary gates")
    end

    # Input Circuit
    if "input_circuit" in keys(params)
        input_circuit = params["input_circuit"]
    else
        # default value
        input_circuit = []
    end

    input_circuit_dict = Dict{String,Any}()

    if length(input_circuit) > 0 && (length(input_circuit) <= params["maximum_depth"])

        input_circuit_dict = QCO.get_input_circuit_dict(input_circuit, params)

    else
        (length(input_circuit) > 0) && (Memento.warn(_LOGGER, "Neglecting the input circuit as it's depth is greater than the allowable depth"))
    end

    # Decomposition type
    if "decomposition_type" in keys(params)
        decomposition_type = params["decomposition_type"]
    else
        decomposition_type = "exact"
    end

    # Objective function
    if "objective" in keys(params)
        objective = params["objective"]
    else
        objective = "minimize_depth"
    end

    elementary_gates = unique(params["elementary_gates"])

    if length(elementary_gates) < length(params["elementary_gates"])
        Memento.warn(_LOGGER, "Eliminating non-unique gates in the input elementary gates")
    end

    gates_dict, are_elementary_gates_real = QCO.get_elementary_gates_dictionary(params, elementary_gates)

    target_real, is_target_real = QCO.get_target_gate(params, are_elementary_gates_real)

    gates_dict_unique, M_real_unique, identity_idx, cnot_idx = QCO.eliminate_nonunique_gates(gates_dict, are_elementary_gates_real, eliminate_identical_gates = eliminate_identical_gates)

    # Initial gate
    if "initial_gate" in keys(params)
        if params["initial_gate"] == "Identity"
            if are_elementary_gates_real && is_target_real
                initial_gate = real(QCO.IGate(num_qubits))
            else
                initial_gate = QCO.complex_to_real_gate(QCO.IGate(num_qubits))
            end
        else
            Memento.error(_LOGGER, "Currently, only \"Identity\" is supported as an initial gate")
            # Add code here to support non-identity as an initial gate.
        end
    else
        if are_elementary_gates_real && is_target_real
            initial_gate = real(QCO.IGate(num_qubits))
        else
            initial_gate = QCO.complex_to_real_gate(QCO.IGate(num_qubits))
        end
    end

    data = Dict{String, Any}("num_qubits" => num_qubits,
                             "maximum_depth" => maximum_depth,
                             "gates_dict" => gates_dict_unique,
                             "gates_real" => M_real_unique,
                             "initial_gate" => initial_gate,
                             "identity_idx" => identity_idx,
                             "cnot_idx" => cnot_idx,
                             "elementary_gates" => elementary_gates,
                             "target_gate" => target_real,
                             "are_gates_real" => (are_elementary_gates_real && is_target_real),
                             "objective" => objective,
                             "decomposition_type" => decomposition_type
                             )

    if data["are_gates_real"]
        Memento.info(_LOGGER, "Detected all-real elementary and target gates")
    end

    # Rotation and Universal gate angle discretizations
    data = QCO._populate_data_angle_discretization!(data, params)

    # Determinant test for input gates
    QCO._determinant_test_for_infeasibility(data)

    # Input circuit
    if length(keys(input_circuit_dict)) > 0
        data["input_circuit"] = input_circuit_dict
    end

    # CNOT lower/upper bound
    data = QCO._get_cnot_bounds!(data, params)

    return data
end

"""
    eliminate_nonunique_gates(gates_dict::Dict{String, Any})

"""
function eliminate_nonunique_gates(gates_dict::Dict{String, Any}, are_elementary_gates_real::Bool; eliminate_identical_gates = false)

    num_gates = length(keys(gates_dict))

    # This identifies if all the elementary gates have only real entries and returns compact matrices
    if are_elementary_gates_real
        M_real = zeros(size(gates_dict["1"]["matrix"])[1], size(gates_dict["1"]["matrix"])[2], num_gates)
    else
        M_real = zeros(2*size(gates_dict["1"]["matrix"])[1], 2*size(gates_dict["1"]["matrix"])[2], num_gates)
    end

    for i=1:num_gates
        if are_elementary_gates_real
            M_real[:,:,i] = real(gates_dict["$i"]["matrix"])
        else
            M_real[:,:,i] = QCO.complex_to_real_gate(gates_dict["$i"]["matrix"])
        end
    end

    M_real_unique = M_real
    M_real_idx = collect(1:size(M_real)[3])

    if eliminate_identical_gates
        M_real_unique, M_real_idx = QCO.unique_matrices(M_real)
    end

    gates_dict_unique = Dict{String, Any}()

    if size(M_real_unique)[3] < size(M_real)[3]

        Memento.info(_LOGGER, "Detected $(size(M_real)[3]-size(M_real_unique)[3]) non-unique gates (after discretization)")

        for i = 1:length(M_real_idx)
            gates_dict_unique["$i"] = gates_dict["$(M_real_idx[i])"]
        end

    else
        gates_dict_unique = gates_dict
    end

    identity_idx = QCO._get_identity_idx(M_real_unique)

    for i_id = 1:length(identity_idx)
        if !("Identity" in gates_dict_unique["$(identity_idx[i_id])"]["type"])
            push!(gates_dict_unique["$(identity_idx[i_id])"]["type"], "Identity")
        end
    end

    cnot_idx = QCO._get_cnot_idx(gates_dict_unique)

    return gates_dict_unique, M_real_unique, identity_idx, cnot_idx
end

function _populate_data_angle_discretization!(data::Dict{String, Any}, params::Dict{String, Any})

    one_angle_gates, three_angle_gates = QCO._get_angle_gates_idx(data["elementary_gates"])

    if !isempty(union(one_angle_gates, three_angle_gates))

        if length(findall(x -> (occursin("discretization", x)), collect(keys(params)))) == 0
            Memento.error(_LOGGER, "Enter valid discretization angles in the imput params")
        end

        data["discretization"] = Dict{String, Any}()

        if !isempty(one_angle_gates)
            for i in one_angle_gates
                gate_type = QCO._parse_gate_string(data["elementary_gates"][i], type = true)

                data["discretization"][gate_type] = Float64.(params["$(gate_type)_discretization"])
            end
        end

        if !isempty(three_angle_gates)
            for i in three_angle_gates
                gate_type = QCO._parse_gate_string(data["elementary_gates"][i], type = true)

                data["discretization"]["$(gate_type)_θ"] = Float64.(params["$(gate_type)_θ_discretization"])
                data["discretization"]["$(gate_type)_ϕ"] = Float64.(params["$(gate_type)_ϕ_discretization"])
                data["discretization"]["$(gate_type)_λ"] = Float64.(params["$(gate_type)_λ_discretization"])

            end
        end

    end

    return data
end

"""
    get_target_gate(params::Dict{String, Any}, are_elementary_gates_real::Bool)

Given the user input `params` dictionary and a boolean if all the input elementary gates are real,
this function returns the corresponding real version of the target gate.
"""
function get_target_gate(params::Dict{String, Any}, are_elementary_gates_real::Bool)

    if !("target_gate" in keys(params)) || isempty(params["target_gate"])
        Memento.error(_LOGGER, "Target gate not found in the input data")
    end

    if (size(params["target_gate"])[1] != size(params["target_gate"])[2]) || (size(params["target_gate"])[1] != 2^params["num_qubits"])
        Memento.error(_LOGGER, "Dimensions of target gate do not match the input num_qubits")
    end

    # This identifies if all the target gate has only real entries and returns a compact matrix
    is_target_real = QCO.is_gate_real(params["target_gate"])

    if are_elementary_gates_real
        if !is_target_real
            Memento.error(_LOGGER, "Infeasible decomposition: all elementary gates have zero imaginary parts")
        else
            return real(params["target_gate"]), is_target_real
        end
    else
        return QCO.complex_to_real_gate(params["target_gate"]), is_target_real
    end

end

function get_elementary_gates_dictionary(params::Dict{String, Any}, elementary_gates::Array{String,1})

    num_qubits = params["num_qubits"]

    kron_gates_idx  = QCO._get_kron_gates_idx(elementary_gates)

    one_angle_gates, three_angle_gates = QCO._get_angle_gates_idx(elementary_gates)

    one_angle_gates_dict   = Dict{String,Any}()
    three_angle_gates_dict = Dict{String,Any}()

    if !isempty(one_angle_gates)
        one_angle_gates_dict   = QCO.get_all_one_angle_gates(params, elementary_gates, one_angle_gates)
    end
    if !isempty(three_angle_gates)
        three_angle_gates_dict = QCO.get_all_three_angle_gates(params, elementary_gates, three_angle_gates)
    end

    all_angle_gates_dict = merge(one_angle_gates_dict, three_angle_gates_dict)

    gates_dict = Dict{String, Any}()

    counter = 1

    for i=1:length(elementary_gates)

        if i in union(one_angle_gates, three_angle_gates)

            M_elementary_dict = all_angle_gates_dict[elementary_gates[i]]

            for k in keys(M_elementary_dict) # Angle
                for l in keys(M_elementary_dict[k]["$(num_qubits)qubit_rep"]) # qubits (which will now be 1)

                    gates_dict["$counter"] = Dict{String, Any}("type" => [elementary_gates[i]],
                                                                "angle" => Any,
                                                                "qubit_loc" => l,
                                                                "matrix" => M_elementary_dict[k]["$(num_qubits)qubit_rep"][l])

                    if i in one_angle_gates
                        gates_dict["$counter"]["angle"] = M_elementary_dict[k]["angle"]

                    elseif i in three_angle_gates
                        gates_dict["$counter"]["angle"] = Dict{String, Any}("θ" => M_elementary_dict[k]["θ"],
                                                                            "ϕ" => M_elementary_dict[k]["ϕ"],
                                                                            "λ" => M_elementary_dict[k]["λ"],)
                    end

                    counter += 1
                end
            end

        # Elementary gates which contain Kronecker symbols
        elseif i in kron_gates_idx
            M = QCO.get_full_sized_kron_symbol_gate(elementary_gates[i], num_qubits)

            gates_dict["$counter"] = Dict{String, Any}("type" => [elementary_gates[i]],
                                                       "matrix" => M)
            counter += 1
        else

            M = QCO.get_full_sized_gate(elementary_gates[i], num_qubits)

            gates_dict["$counter"] = Dict{String, Any}("type" => [elementary_gates[i]],
                                                       "matrix" => M)
            counter += 1
        end

    end

    are_elementary_gates_real = true

    for i in keys(gates_dict)
        if !(QCO.is_gate_real(gates_dict[i]["matrix"]))
            are_elementary_gates_real = false
            continue
        end
    end

    return gates_dict, are_elementary_gates_real
end

function get_all_one_angle_gates(params::Dict{String, Any}, elementary_gates::Array{String,1}, gates_idx::Vector{Int64})

    gates_complex = Dict{String, Any}()

    if length(gates_idx) >= 1

        for i=1:length(gates_idx)

            input_gate = elementary_gates[gates_idx[i]]

            gate_name = QCO._parse_gate_string(input_gate, type=true)

            if isempty(params[string(gate_name,"_discretization")])
                Memento.error(_LOGGER, "Empty discretization angles for $(input_gate) gate. Input a valid angle")
            end

            angle_disc = Float64.(params[string(gate_name,"_discretization")])

            gates_complex[input_gate] = Dict{String, Any}()
            gates_complex[input_gate] = QCO.get_discretized_one_angle_gates(input_gate, gates_complex[input_gate], angle_disc, params["num_qubits"])

        end

    end

    return gates_complex
end

function get_discretized_one_angle_gates(gate_type::String, M1::Dict{String, Any}, discretization::Array{Float64,1}, num_qubits::Int64)

    if length(discretization) >= 1

        for i=1:length(discretization)
            angles = discretization[i]
            M1["angle_$i"] = Dict{String, Any}("angle" => angles,
                                             "$(num_qubits)qubit_rep" => Dict{String, Any}() )

            qubit_loc = QCO._parse_gate_string(gate_type, qubits=true)
            if length(qubit_loc) == 1
                qubit_loc_str = string(qubit_loc[1])
            elseif length(qubit_loc) == 2
                qubit_loc_str = string(qubit_loc[1], qubit_separator, qubit_loc[2])
            end

            M1["angle_$i"]["$(num_qubits)qubit_rep"]["qubit_$(qubit_loc_str)"] = QCO.get_full_sized_gate(gate_type, num_qubits, angle = angles)
        end
    end

    return M1
end

function get_all_three_angle_gates(params::Dict{String, Any}, elementary_gates::Array{String,1}, gates_idx::Vector{Int64})

    gates_complex = Dict{String, Any}()

    if length(gates_idx) >= 1

        for i=1:length(gates_idx)

            input_gate = elementary_gates[gates_idx[i]]

            gate_name = QCO._parse_gate_string(input_gate, type = true)

            gates_complex[input_gate] = Dict{String, Any}()

            for angle in ["θ", "ϕ", "λ"]
                if isempty(params["$(gate_name)_$(angle)_discretization"])
                    Memento.error(_LOGGER, "Empty $(angle) discretization angle for $input_gate gate. Input atleast one valid angle")
                end
            end

            θ_disc = Vector{Float64}(params["$(gate_name)_θ_discretization"])
            ϕ_disc = Vector{Float64}(params["$(gate_name)_ϕ_discretization"])
            λ_disc = Vector{Float64}(params["$(gate_name)_λ_discretization"])

            gates_complex[input_gate] = QCO.get_discretized_three_angle_gates(input_gate, gates_complex[input_gate], θ_disc, ϕ_disc, λ_disc, params["num_qubits"])

        end

    end

    return gates_complex
end

function get_discretized_three_angle_gates(gate_type::String, M3::Dict{String, Any}, θ_discretization::Array{Float64,1}, ϕ_discretization::Array{Float64,1}, λ_discretization::Array{Float64,1}, num_qubits::Int64)

    counter = 1

    for i=1:length(θ_discretization)
        for j=1:length(ϕ_discretization)
            for k=1:length(λ_discretization)
                angles = [θ_discretization[i], ϕ_discretization[j], λ_discretization[k]]

                M3["angle_$(counter)"] = Dict{String, Any}("θ" => angles[1],
                                                           "ϕ" => angles[2],
                                                           "λ" => angles[3],
                                                           "$(num_qubits)qubit_rep" => Dict{String, Any}()
                                                          )
                qubit_loc = QCO._parse_gate_string(gate_type, qubits=true)
                if length(qubit_loc) == 1
                    qubit_loc_str = string(qubit_loc[1])
                elseif length(qubit_loc) == 2
                    qubit_loc_str = string(qubit_loc[1], qubit_separator, qubit_loc[2])
                end

                M3["angle_$(counter)"]["$(num_qubits)qubit_rep"]["qubit_$(qubit_loc_str)"] = QCO.get_full_sized_gate(gate_type, num_qubits, angle = angles)

                counter += 1
            end
        end
    end

    return M3
end

"""
    get_full_sized_gate(input::String, num_qubits::Int64; angle = nothing)

Given an input string representing the gate and number of qubits of the circuit, this function returns a full-sized
gate with respect to the input number of qubits. For example, if `num_qubits = 3` and the input gate in `H_3`
(Hadamard on third qubit), then this function returns `IGate ⨷ IGate ⨷ HGate`, where IGate and HGate are single qubit Identity and Hadamard gates, respectively.
Note that `angle` vector is an optional input which is necessary when the input gate is parametrized by Euler angles.
"""
function get_full_sized_gate(input::String, num_qubits::Int64; angle = nothing)

    if num_qubits > 10
        Memento.error(_LOGGER, "Greater than 10 qubits is currently not supported")
    end

    if input == "Identity"
        return QCO.IGate(num_qubits)
    end

    gate_type, qubit_loc = QCO._parse_gate_string(input, type = true, qubits = true)

    if !(gate_type in union(QCO.ONE_QUBIT_GATES, QCO.TWO_QUBIT_GATES))
        Memento.error(_LOGGER, "Specified $input gate does not exist in the predefined set of gates")
    end

    QCO._catch_input_gate_errors(gate_type, qubit_loc, num_qubits, input)

    #----------------------;
    #   One qubit gates    ;
    #----------------------;
    if length(qubit_loc) == 1

        if gate_type in QCO.ONE_QUBIT_GATES_CONSTANTS

            return QCO.kron_single_qubit_gate(num_qubits, getfield(QCO, Symbol(gate_type, "Gate"))(), "q$(qubit_loc[1])")

        elseif gate_type in QCO.ONE_QUBIT_GATES_ANGLE_PARAMETERS

            if (angle != nothing) && (length(angle) > 0)

                if length(angle) == 1
                    return QCO.kron_single_qubit_gate(num_qubits, getfield(QCO, Symbol(gate_type, "Gate"))(angle), "q$(qubit_loc[1])")
                elseif length(angle) == 3
                    return QCO.kron_single_qubit_gate(num_qubits, getfield(QCO, Symbol(gate_type, "Gate"))(angle[1], angle[2], angle[3]), "q$(qubit_loc[1])")
                end

            else
                Memento.error(_LOGGER, "Enter a valid angle parameter for the input $input gate")
            end

        end

    #----------------------;
    #   Two qubit gates    ;
    #----------------------;
    elseif length(qubit_loc) == 2

        if gate_type in QCO.TWO_QUBIT_GATES_CONSTANTS

            if (qubit_loc[1] < qubit_loc[2]) || (gate_type in QCO.TWO_QUBIT_GATES_CONSTANTS_SYMMETRIC)
                return QCO.kron_two_qubit_gate(num_qubits, getfield(QCO, Symbol(gate_type, "Gate"))(), "q$(qubit_loc[1])", "q$(qubit_loc[2])")
            else
                return QCO.kron_two_qubit_gate(num_qubits, getfield(QCO, Symbol(gate_type, "RevGate"))(), "q$(qubit_loc[1])", "q$(qubit_loc[2])")
            end

        elseif gate_type in QCO.TWO_QUBIT_GATES_ANGLE_PARAMETERS

            if (angle != nothing) && (length(angle) > 0)

                if length(angle) == 1
                    if (qubit_loc[1] < qubit_loc[2])
                        return QCO.kron_two_qubit_gate(num_qubits, getfield(QCO, Symbol(gate_type, "Gate"))(angle), "q$(qubit_loc[1])", "q$(qubit_loc[2])")
                    else
                        return QCO.kron_two_qubit_gate(num_qubits, getfield(QCO, Symbol(gate_type, "RevGate"))(angle), "q$(qubit_loc[1])", "q$(qubit_loc[2])")
                    end
                elseif length(angle) == 3
                    if (qubit_loc[1] < qubit_loc[2])
                        return QCO.kron_two_qubit_gate(num_qubits, getfield(QCO, Symbol(gate_type, "Gate"))(angle[1], angle[2], angle[3]), "q$(qubit_loc[1])", "q$(qubit_loc[2])")
                    else
                        return QCO.kron_two_qubit_gate(num_qubits, getfield(QCO, Symbol(gate_type, "RevGate"))(angle[1], angle[2], angle[3]), "q$(qubit_loc[1])", "q$(qubit_loc[2])")
                    end
                end

            else
                Memento.error(_LOGGER, "Enter a valid angle parameter for the input $input gate")
            end

        end
    end

end

"""
    get_full_sized_kron_symbol_gate(input::String, num_qubits::Int64)

Given an input string with kronecker symbols representing the gate and number of qubits of the circuit, this function returns a full-sized
gate with respect to the input number of qubits. For example, if `num_qubits = 3` and the input gate in `I_1xT_1xH_3`,
then this function returns `IGate ⨷ TGate ⨷ HGate`, where IGate, TGate and HGate are single-qubit Identity, T and Hadamard gates, respectively.
Note that this function currently does not support an input gate parametrized with Euler angles.
"""
function get_full_sized_kron_symbol_gate(input::String, num_qubits::Int64)

    kron_gates = QCO._parse_gates_with_kron_symbol(input)

    M = 1

    for i = 1:length(kron_gates)

        gate_type, qubit_loc = QCO._parse_gate_string(kron_gates[i], type = true, qubits = true)

        if !(gate_type in union(QCO.ONE_QUBIT_GATES_CONSTANTS, QCO.TWO_QUBIT_GATES_CONSTANTS))
            Memento.error(_LOGGER, "Specified $input gate is not supported in conjunction with the Kronecker product operation")
        end

        QCO._catch_input_gate_errors(gate_type, qubit_loc, num_qubits, input)

        # One qubit gates
        if (length(qubit_loc) == 1) && (gate_type in QCO.ONE_QUBIT_GATES_CONSTANTS)
            if (gate_type == "I") || (gate_type == "Identity")
                M = kron(M, getfield(QCO, Symbol(gate_type, "Gate"))(1))
            else
                M = kron(M, getfield(QCO, Symbol(gate_type, "Gate"))())
            end

        # Two qubit gates
        elseif (length(qubit_loc) == 2) && (gate_type in QCO.TWO_QUBIT_GATES_CONSTANTS)
            if (qubit_loc[1] < qubit_loc[2]) || (gate_type in QCO.TWO_QUBIT_GATES_CONSTANTS_SYMMETRIC)
                M = kron(M, getfield(QCO, Symbol(gate_type, "Gate"))())
            else
                M = kron(M, getfield(QCO, Symbol(gate_type, "RevGate"))())
            end
        end
    end

    if size(M)[1] == 2^(num_qubits)
        return M
    else
        Memento.error(_LOGGER, "Dimensions mismatch in evaluation of elementary gates with Kronecker symbol")
    end

end

"""
    get_input_circuit_dict(input_circuit::Vector{Tuple{Int64,String}}, params::Dict{String,Any})

Given the user input circuit which serves as a warm-start to the optimization model, and user input params dictionary,
this function outputs the post-processed dictionary of the input circuit which is used by the optimization model.
"""
function get_input_circuit_dict(input_circuit::Vector{Tuple{Int64,String}}, params::Dict{String,Any})

    input_circuit_dict = Dict{String, Any}()

    status = true
    gate_type = []

    for i = 1:length(input_circuit)
        if !(input_circuit[i][2] in params["elementary_gates"])
            status = false
            gate_type = input_circuit[i][2]
        end
    end

    if status
        for i = 1:length(input_circuit)

            if i == input_circuit[i][1]
                input_circuit_dict["$i"] = Dict{String, Any}("depth" => input_circuit[i][1],
                                                             "gate" => input_circuit[i][2])
                # Later: add support for universal and rotation gates here
            else
                input_circuit_dict = Dict{String, Any}()
                Memento.warn(_LOGGER, "Neglecting the input circuit as multiple gates cannot be input at the same depth")
                break
            end

        end
    else
        Memento.warn(_LOGGER, "Neglecting the input circuit as gate $gate_type is not in input elementary gates")
    end

    return input_circuit_dict
end

"""
    _catch_input_gate_errors(gate_type::String, qubit_loc::Vector{Int64}, num_qubits::Int64, input_gate::String)

Given an input gate string, number of qubits of the circuit and the qubit locations for the input gate,
this function catches and throws any errors, should the input gate type and qubits are invalid.
"""
function _catch_input_gate_errors(gate_type::String, qubit_loc::Vector{Int64}, num_qubits::Int64, input_gate::String)

    if (gate_type in QCO.TWO_QUBIT_GATES) && (length(qubit_loc) != 2)
        Memento.error(_LOGGER, "Specify two qubits for $input_gate, which is a 2-qubit gate")
    elseif (gate_type in QCO.ONE_QUBIT_GATES) && (length(qubit_loc) != 1)
        Memento.error(_LOGGER, "Specify only one qubit for $input_gate, which is a 1-qubit gate")
    end

    if isempty(qubit_loc)
        Memento.error(_LOGGER, "Specify a valid qubit location(s) for the input $input_gate gate")

    elseif !issubset(qubit_loc, 1:num_qubits)
        Memento.error(_LOGGER, "Specified qubit(s) for the $input_gate gate ∉ [1,...,$num_qubits]")

    elseif (length(qubit_loc) == 2) && (isapprox(qubit_loc[1], qubit_loc[2]))
        Memento.error(_LOGGER, "Specified $input_gate gate cannot have identical control and target qubits")

    elseif length(qubit_loc) > 2
        Memento.error(_LOGGER, "Only 1- and 2-qubit elementary gates are currently supported")
    end

end

function _get_R_gates_idx(elementary_gates::Array{String,1})

    return findall(x -> (startswith(x, "RX") || startswith(x, "RY") || startswith(x, "RZ")) && !(occursin(kron_symbol, x)), elementary_gates)
end

function _get_U3_gates_idx(elementary_gates::Array{String,1})

    return findall(x -> (startswith(x, "U3")) && !(occursin(kron_symbol, x)), elementary_gates)
end

function _get_CR_gates_idx(elementary_gates::Array{String,1})

    return findall(x -> (startswith(x, "CRX") || startswith(x, "CRY") || startswith(x, "CRZ")) && !(occursin(kron_symbol, x)), elementary_gates)
end

function _get_CU3_gates_idx(elementary_gates::Array{String, 1})

    return findall(x -> (startswith(x, "CU3")) && !(occursin(kron_symbol, x)), elementary_gates)
end

function _get_kron_gates_idx(elementary_gates::Array{String, 1})

    return findall(x -> occursin(kron_symbol, x), elementary_gates)
end

function _get_Phase_gates_idx(elementary_gates::Array{String, 1})

    return findall(x -> (startswith(x, "Phase")) && !(occursin(kron_symbol, x)), elementary_gates)
end

function _get_identity_idx(M::Array{Float64,3})

    identity_idx = Int64[]

    for i=1:size(M)[3]
        if isapprox(M[:,:,i], Matrix(I, size(M)[1], size(M)[2]), atol=1E-5)
            push!(identity_idx, i)
        end
    end

    return identity_idx
end

function _get_cnot_idx(gates_dict::Dict{String, Any})

    cnot_idx = Int64[]

    # Counts for both (CNot_i_j and CNot_j_i) or (CX_i_j and CX_j_i)
    for i in keys(gates_dict)
        # Input gates with kron symbols
        if occursin(kron_symbol, gates_dict[i]["type"][1])
            gate_types = QCO._parse_gates_with_kron_symbol(gates_dict[i]["type"][1])

            for j = 1:length(gate_types)
                gate_type = QCO._parse_gate_string(gate_types[j], type = true)
                if gate_type in ["CNot", "CX"]
                    push!(cnot_idx, parse(Int64, i))
                end
            end
        else

        # Single input gate per depth
            if "Identity" in gates_dict[i]["type"]
                continue
            end

            gate_type = QCO._parse_gate_string(gates_dict[i]["type"][1], type = true)
            if gate_type in ["CNot", "CX"]
                push!(cnot_idx, parse(Int64, i))
            end
        end
    end

    return cnot_idx
end

function _get_angle_gates_idx(elementary_gates::Array{String,1})

    # Update this list as new gates with angle parameters are added in gates.jl
    R_gates_idx     = QCO._get_R_gates_idx(elementary_gates)
    Phase_gates_idx = QCO._get_Phase_gates_idx(elementary_gates)
    CR_gates_idx    = QCO._get_CR_gates_idx(elementary_gates)
    U3_gates_idx    = QCO._get_U3_gates_idx(elementary_gates)
    CU3_gates_idx   = QCO._get_CU3_gates_idx(elementary_gates)

    one_angle_gates   = union(R_gates_idx, Phase_gates_idx, CR_gates_idx)
    three_angle_gates = union(U3_gates_idx, CU3_gates_idx)

    return one_angle_gates, three_angle_gates
end

function _get_cnot_bounds!(data::Dict{String, Any}, params::Dict{String, Any})

    cnot_lb = 0
    cnot_ub = data["maximum_depth"]

    if "set_cnot_lower_bound" in keys(params)
        cnot_lb = params["set_cnot_lower_bound"]
    end

    if "set_cnot_upper_bound" in keys(params)
        cnot_ub = params["set_cnot_upper_bound"]
    end

    if cnot_lb < cnot_ub
        if cnot_lb > 0
            data["cnot_lower_bound"] = params["set_cnot_lower_bound"]
        end
        if cnot_ub < data["maximum_depth"]
            data["cnot_upper_bound"] = params["set_cnot_upper_bound"]
        end
    elseif isapprox(cnot_lb, cnot_ub, atol=1E-6)
        data["cnot_lower_bound"] = cnot_lb
        data["cnot_upper_bound"] = cnot_ub
    else
        Memento.warn(_LOGGER, "Invalid CNot-gate lower/upper bound")
    end

    return data
end


text:
"""
The ElecChain supertype corresponds to the chain of
    electronic components involved in the DAQ setup
    that appear before the trigger.
"""
abstract type ElecChain end

"""
Generic electronics chain consisting of a preamplifier
    and an FADC module
"""
mutable struct GenericElecChain <: ElecChain
    preamp::PreAmp

    fadc::FADC
end

"""
    GenericElecChain(sim_conf)

PropDict -> GenericElecChain

Construct electronics components based on simulation
    configuration <sim_conf> and create a GenericChain instance
    based on these components.
"""
function GenericElecChain(sim_conf::PropDict)
    preamp = PreAmp(sim_conf)
    fadc = FADC(sim_conf)

    GenericElecChain(preamp, fadc)
end

"""
    ElecChain(sim_conf)

PropDict -> <ElecChain>

Construct an ElecChain supertype struct based on given simulation configuration.
Type of returned instance depends on settings in <sim_conf>
Currently only one type of ElecChain available (GenericElecChain),
    rendering this function temporarily redundant.
"""
function ElecChain(sim_conf::PropDict)
    GenericElecChain(sim_conf)
end

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

"""
    simulate(wf, elec_chain, ::NoiseFromSim)

RDWaveform, ElecChain -> RDWaveform

Simulate effects of the electronics chain on the waveform
    modelling noise based on each component.

"""
function simulate(wf::RDWaveform, elec_chain::ElecChain, ::NoiseFromSim)
    # simply simulate elec chain the way it is
    # noise simulation of components is included in the simulation of components
    simulate(wf, elec_chain)
end

"""
    simulate(wf, elec_chain, noise_model)

RDWaveform, GenericElecChain, NoiseFromData -> RDWaveform

Simulate effects of the electronics chain on the waveform
    modelling noise based on baselines extracted from data
    (note: also takes preamp offset into account)
"""
function simulate(wf::RDWaveform, elec_chain::GenericElecChain, noise_model::NoiseFromData)
    ## simulate elec chain

    # check if noise and offset are double counted
    if (elec_chain.preamp.noise_σ != 0u"keV") || (elec_chain.preamp.offset != 0u"keV")
        @info "WARNING!\n
        You are simulating noise and offset based on data baselines,
        but your PreAmp settings are:\n
        -- noise sigma = $(elec_chain.preamp.noise_σ)\n
        -- offset = $(elec_chain.preamp.offset)"
        @info "Proceeding with the simulation, but there might be double counting"
    end

    # simulate elec chain
    wf_elec = simulate(wf, elec_chain)

    ## simulate noise and offset via slapping data baselines on top of the resulting waveform
    simulate_noise(wf_elec, noise_model)
end

"""
    simulate(wf, elec_chain)

RDWaveform, ElecChain -> RDWaveform

Simulate effects of the electronics chain components on the waveform.
"""
function simulate(wf::RDWaveform, elec_chain::ElecChain)
    # elec chain components
    component_names = fieldnames(typeof(elec_chain))

    # start with originally given waveform
    wf_elec = wf
    # simulate effect of each component
    for name in component_names
        component = getfield(elec_chain, name)
        wf_elec = simulate(wf_elec, component)
    end

    wf_elec
end


text:
#examples of using jld22sac functions to convert JLD2 files to sac format.
using SeisIO, SeisConvert, Sockets, JLD2
# using SeisIO, Sockets, JLD2
# include("/Users/xiaotaoyang/SOFT/mySeisJL/code/SeisConvert.jl/src/jld22sac.jl")
# include("/Users/xiaotaoyang/SOFT/mySeisJL/code/SeisConvert.jl/src/sac2jld2.jl")
# include("/Users/xiaotaoyang/SOFT/mySeisJL/code/SeisConvert.jl/src/corr2seis.jl")
# include("/Users/xiaotaoyang/SOFT/mySeisJL/code/SeisConvert.jl/src/utilities.jl")
# Saves individual sac file (in) to JLD2 file.
# The output file name will be the same as the input, with new extension *.jld2.

# 1. prepare a SeisData jld2 file first.
println("preparing JLD2 file using test dataset")
sacdirlist = ls("sac/Raw")
timestamplist=["2010.001T00.00.00","2010.002T00.00.00"]
sac2jld2(sacdirlist,timestamplist,"testSeisData.jld2")

# 2. Convert it back to SAC using jld22sac
println(" Convert it back to SAC using jld22sac")
jld22sac("testSeisData.jld2","SeisData_SAC_all_TG",informat="TD",outformat="TG",verbose=true)

# 3. Convert it back to SAC using jld22sac, output as SRC format/structure
println(" Convert it back to SAC using jld22sac, output as SRC format/structure")
jld22sac("testSeisData.jld2","SeisData_SAC_all_SRC",informat="TD",outformat="SRC",verbose=true)

# 4. Convert it back to SAC using jld22sac, subset timestamplist
println(" Convert it back to SAC using jld22sac, subset timestamplist")
subsetcondition = Dict("timestamplist" => ["2010.001T00.00.00"])
jld22sac("testSeisData.jld2","SeisData_SAC_subset_TG",subset=subsetcondition,informat="TD",outformat="TG",verbose=true)

"""
convert CorrData to SAC
"""
jld22sac("SeisXcorr_testout.jld2","CorrData_SAC_all_SRC",informat="TD",outformat="SRC",verbose=true)