text:
# Excursions in formal power series (fps)
MaxSeriesSize=10
MaxRange = 50
MatrixSize=15
TT=Int64

(n1, n2, n3) = int(rand(3)*MaxSeriesSize)

X = FormalPowerSeries{TT}(int(rand(n1)*MaxRange))
Y = FormalPowerSeries{TT}(int(rand(n2)*MaxRange))
Z = FormalPowerSeries{TT}(int(rand(n3)*MaxRange))

c = int(rand()*MatrixSize) #Size of matrix representation to generate

nzeros = int(rand()*MaxSeriesSize)
@assert X == trim(X)
XX = deepcopy(X)
for i=1:nzeros
    idx = int(rand()*MaxRange)
    if !haskey(XX.c, idx)
        XX.c[idx] = convert(TT, 0)
    end
end
@assert trim(XX) == X

#Test addition, p.15, (1.3-4)
@assert X+X == 2X
@assert X+Y == Y+X
@assert MatrixForm(X+Y,c) == MatrixForm(X,c)+MatrixForm(Y,c)

#Test subtraction, p.15, (1.3-4)
@assert X-X == 0X
@assert X-Y == -(Y-X)
@assert MatrixForm(X-Y,c) == MatrixForm(X,c)-MatrixForm(Y,c)

#Test multiplication, p.15, (1.3-5)
@assert X*Y == Y*X
@assert MatrixForm(X*X,c) == MatrixForm(X,c)*MatrixForm(X,c)
@assert MatrixForm(X*Y,c) == MatrixForm(X,c)*MatrixForm(Y,c)
@assert MatrixForm(Y*X,c) == MatrixForm(Y,c)*MatrixForm(X,c)
@assert MatrixForm(Y*Y,c) == MatrixForm(Y,c)*MatrixForm(Y,c)

@assert X.*Y == Y.*X

#The reciprocal series has associated matrix that is the matrix inverse
#of the original series
#Force reciprocal to exist
X.c[0] = 1
discrepancy = (norm(inv(float(MatrixForm(X,c)))[1, :]'[:, 1] - tovector(reciprocal(X, c),[0:c-1])))
if discrepancy > c*sqrt(eps())
    error(@sprintf("Error %f exceeds tolerance %f", discrepancy, c*sqrt(eps())))
end
#@assert norm(inv(float(MatrixForm(X,c)))[1, :]'[:, 1] - tovector(reciprocal(X, c),c)) < c*sqrt(eps())

#Test differentiation
XX = derivative(X)
for (k, v) in XX.c
    k==0 ? continue : true
    @assert X.c[k+1] == v/(k+1)
end

#Test product rule [H, Sec.1.4, p.19]
@assert derivative(X*Y) == derivative(X)*Y + X*derivative(Y)

#Test right distributive law of composition [H, Sec.1.6, p.38]
@assert compose(X,Z)*compose(Y,Z) == compose(X*Y, Z)

#Test chain rule [H, Sec.1.6, p.40]
@assert derivative(compose(X,Y)) == compose(derivative(X),Y)*derivative(Y)


text:
# Implements the method of PPT relaxations of MX operations

using Convex
using SCS

export pptRelax, pptRelaxCopies;

""" `(problem, F, p_succ) = pptRelax(rho,n,k,delta, verbose, eps, max_iters)`

Implements the PPT relaxation computing an upper bound on the fidelity achievable using measure and exchange operations for a given input state and fixed succes probability.

Inputs:
- *rho* quantum state to be distilled
- *n* number of qubits on one side, assuming dimensions nA=nB =2^n in the input.
- *k* desired output dimension of the maximally entangled state
- *delta* desired success probability
- *verbose* flag indicating whether the SCS solver should be verbose (default false)
- *eps* desired precision (default 1e-4) in SCS solver
- *max_iters* maximum number of iterations in SCS solver (default 2e4)

Outputs:
- *problem* SDP object from Convex.jl allowing further processing.
- *F* Fidelity achieved
- *p_succ* success probability attained in SDP (should equal *delta*)
"""

function pptRelax(rho::AbstractMatrix, n::Int, k::Int, delta::Real; verbose::Bool = false, eps::Number = 1e-4, max_iters = 2e4)
	# Presume the input dimensions are equally divided and we're using qubits
	nA = 2^n;
	nB = nA;

	return pptRelax(rho, nA, nB, k, delta, verbose = verbose, eps = eps, max_iters = max_iters)
end

""" `(problem, F, p_succ) = pptRelax(rho,nA, nB, k,delta, verbose, eps, max_iters)`

Implements the PPT relaxation computing an upper bound on the fidelity achievable using measure and exchange operations for a given input state and fixed succes probability.

Inputs:
- *rho* quantum state to be distilled on A and B
- *nA* dimension of the A system
- *nB* dimension of the B system
- *k* desired output dimension of the maximally entangled state
- *delta* desired success probability
- *verbose* flag indicating whether the SCS solver should be verbose (default true)
- *eps* desired precision (default 1e-4) in SCS solver
- *max_iters* maximum number of iterations in SCS solver (default 2e4)

Outputs:
- *problem* SDP object from Convex.jl allowing further processing.
- *F* Fidelity achieved
- *p_succ* success probability attained in SDP (should equal *delta*)
"""

function pptRelax(rho::AbstractArray, nA::Int, nB::Int, k::Int, delta::Number; verbose::Bool = true, eps::Number = 1e-4, max_iters::Number = 2e4)

	# Check whether rho is a quantum state
	@assert isQuantumState(rho) "The input state must be a valid quantum state."

	# Check whether dimensions match
	(d, db) = size(rho);
	@assert d == nA * nB "Input state doesn't match given dimensions." , d , " != ", nA*nB

	# Check whether delta is a valid probability
	@assert delta <= 1 "Success probability must be less than one."
	@assert 0 <= delta "Success probability must be positive."

	# define the identity matrix on the whole space
	id = eye(d);

	# define the variables - if the input is real it is sufficient to optimize over reals
	if real(rho) == rho
		D = Semidefinite(d);
		E = Semidefinite(d);
	else
		# complex variable support advertises HermitianSemidefinite
		# but it's buggy
		D = ComplexVariable(d,d);
		E = ComplexVariable(d,d);
	end

	# define the objective
	tv = Variable(1);
	problem = maximize(tv);
	problem.constraints += tv == trace(nA * nB * D * rho.');

	problem.constraints += D in :SDP;
	problem.constraints += E in :SDP;
	problem.constraints += (id/d - (D+E)) in :SDP;

	# Constraints relating to the PPT Condition
	EPT = ptranspose(E);
	DPT = ptranspose(D);
	problem.constraints += (DPT + EPT/(k+1)) in :SDP;
	problem.constraints += (- DPT + EPT/(k-1)) in :SDP;
	problem.constraints += (id/d - (DPT+EPT)) in :SDP;

	# Constraint coming from the success probability
	problem.constraints += nA * nB * trace(rho.'*(D+E)) == delta;

	solve!(problem, SCSSolver(verbose = verbose, eps = eps, max_iters = max_iters),verbose=verbose);
	p_succ = nA * nB * trace(rho.'*(D.value + E.value));
	F = problem.optval / p_succ;

	return (problem, F, p_succ)
end

""" `(problem, F, p) = pptRelaxCopies(rho, n, nA, nB, k, delta, verbose, eps, max_iters)`

Calls pptRelax for an input state of the form *rho&*^(tensor *n*), reordered appropriately.

"""

function pptRelaxCopies(rho::AbstractMatrix,n::Int,nA::Int,nB::Int,k::Int,delta::Real; verbose::Bool = false, eps::Number = 1e-4, max_iters::Number = 2e4)

	@assert n > 1 "Need at least two copies."

	# copy rho n times
	rhoNew = copies(rho, n);

  	# sort all the entries, such that all the systems on A are first and all
	# the systems on B are second
	rhoSorted = sortAB(rhoNew, nA, n)

	return pptRelax(rhoSorted,nA^n,nB^n, k, delta, verbose=verbose, eps=eps,max_iters=max_iters);
end

ext: jl


text:
# ******************************************************************************************
# Notices:
#
# Copyright © 2022 United States Government as represented by the Administrator of the
# National Aeronautics and Space Administration.  All Rights Reserved.
#
# Disclaimers
#
# No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND,
# EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY
# THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY
# WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION,
# IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER,
# CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS,
# RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM
# USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND
# LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND
# DISTRIBUTES IT "AS IS."
#
# Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED
# STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.
# IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES,
# EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON,
# OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND
# HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL
# AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY
# SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
# ******************************************************************************************

"""
Delete least significant keys if length exceeds `k`. Facilitates top-k priority queueing.
"""
function DataStructures.enqueue!(pq::PriorityQueue, key, value, k::Int64)
    pq[key] = value
    while length(pq) > k
        delete!(pq, first(keys(pq)))
    end
end

function replay!(mdp::Interface.AbstractASTMDP, result::LocalResult)
    as = getproperty(result, fieldnames(typeof(result))[1])
    Interface.reset!(mdp)
    for a in as
        Interface.act!(mdp, a)
    end
    return Interface.isevent(mdp.sim)
end

function replay!(mdp::Interface.AbstractASTMDP, result::GlobalResult)
    Interface.reset!(mdp)
    while !Interface.terminated(mdp)
        a = result(Interface.observe(mdp))
        Interface.act!(mdp, a)
    end
    return Interface.isevent(mdp.sim)
end

struct RandomPolicy end

function replay!(mdp::Interface.AbstractASTMDP, ::RandomPolicy)
    Interface.reset!(mdp)
    while !Interface.terminated(mdp)
        a = rand(Interface.actions(mdp))
        Interface.act!(mdp, a)
    end
    return Interface.isevent(mdp.sim)
end

struct NullPolicy end

function replay!(mdp::Interface.AbstractASTMDP, ::NullPolicy)
    Interface.reset!(mdp)
    while !Interface.terminated(mdp)
        a = zero(rand(Interface.actions(mdp)))
        Interface.act!(mdp, a)
    end
    return Interface.isevent(mdp.sim)
end


text:
immutable GTypeInfo
  size :: UInt16
  base_init :: Ptr{Void}
  base_finalize :: Ptr{Void}

  class_init :: Ptr{Void}
  class_finalize :: Ptr{Void}
  class_data :: Ptr{Void}

  instance_size :: UInt16
  n_preallocs :: UInt16
  instance_init :: Ptr{Void}

  value_table :: Ptr{Void}
end

immutable GInterfaceInfo
  interface_init :: Ptr{Void}
  interface_finalize :: Ptr{Void}
  interface_data :: Ptr{Void}
end
immutable GObjectClass
  typ :: Gtk.GLib.GType
  construct_props :: Ptr{Void}
  constructor :: Ptr{Void}
  set_property :: Ptr{Void}
  get_property :: Ptr{Void}
  dispose :: Ptr{Void}
  finalize :: Ptr{Void}
  dispatch_properties_changed :: Ptr{Void}
  notify :: Ptr{Void}
  constructed :: Ptr{Void}
  flags :: Cssize_t
  p1 :: Ptr{Void}
  p2 :: Ptr{Void}
  p3 :: Ptr{Void}
  p4 :: Ptr{Void}
  p5 :: Ptr{Void}
  p6 :: Ptr{Void}
end
immutable GObjectImpl
  g_type_instance :: Ptr{Void}
  ref_count :: Cuint
  qdata :: Ptr{Void}
end
immutable GTypeInterface
  g_type :: Gtk.GLib.GType
  g_instance_type :: Gtk.GLib.GType
end
function JuliaTypeToGType(typ :: Type)
  if typ <: Int8
    return g_type_from_name(:gchar)
  elseif typ <: UInt8
    return g_type_from_name(:guchar)
  elseif typ <: Bool
    return g_type_from_name(:gboolean)
  elseif typ <: Cint
    return g_type_from_name(:gint)
  elseif typ <: Cuint
    return g_type_from_name(:guint)
  elseif typ <: Clong
    return g_type_from_name(:glong)
  elseif typ <: Culong
    return g_type_from_name(:gulong)
  elseif typ <: Int64
    return g_type_from_name(:gint64)
  elseif typ <: UInt64
    return g_type_from_name(:guint64)
  elseif typ <: Float32
    return g_type_from_name(:gfloat)
  elseif typ <: Float64
    return g_type_from_name(:gdouble)
  elseif typ <: AbstractString
    return g_type_from_name(:gchararray)
  else
    return g_type_from_name(:invalid)
  end
end
function g_type_register_static(parent_type :: Gtk.GLib.GType, typename :: ASCIIString, type_info :: GTypeInfo, flags :: Cint)
  ccall((:g_type_register_static, Gtk.libgobject), Gtk.GLib.GType, (Gtk.GLib.GType, Cstring, Ref{GTypeInfo}, Cint),
        parent_type, typename, type_info, flags)
end
function g_type_add_interface_static(instance_type :: Gtk.GLib.GType, interface_type :: Gtk.GLib.GType, info :: GInterfaceInfo)
  ccall((:g_type_add_interface_static, Gtk.libgobject), Void, (Gtk.GLib.GType, Gtk.GLib.GType, Ref{GInterfaceInfo}),
        instance_type, interface_type, info)
end
function g_object_new(typ :: Gtk.GLib.GType)
  hnd = ccall((:g_object_new, Gtk.libgobject), Ptr{GObject}, (Gtk.GLib.GType, Ptr{Void}), typ, C_NULL)
  Gtk.GLib.wrap_gobject(hnd)
end

max_stars_repo_name: barcharcraz/eve-julia-gui

max_issues_repo_name: barcharcraz/eve-julia-gui

max_forks_repo_name: barcharcraz/eve-julia-gui


text:
using MetabolicSim
using Test

#generate some parameters
N,M = 2,2
Rm,ρ = [1.0, 1.0],[0.4, 0.2]
u = [1.0 1.0 ; 1.0 0.0]
l = [0.1 0.0 ; 0.2 0.2]
l_sum = sum(l, dims=2)[:]
#make parameter object
p_1 = MetabolicSim.Parameters(N,M,u,Rm,l,l_sum,ρ)
p_2 = MetabolicSim.make_parameters(N,M,u,Rm,l,ρ)

#function to test if parameter structs are equal
function isequal(a::MetabolicSim.Parameters, b::MetabolicSim.Parameters)
    itt = fieldnames(MetabolicSim.Parameters)
    for i = itt
        if getfield(a, i) != getfield(b, i)
            return false
        end
    end
    return true
end

@testset begin
    #test Parameters and make_parameters give the same objects
    @test isequal(p_1,p_2)

    #test make_parameters assertion error for leakage
    l = ones(M,M)
   @test_throws AssertionError MetabolicSim.make_parameters(N,M,u,Rm,l,ρ)
end


text:
"""
An instance of the uniform-matroid problem. It is sometimes called the m-set
problem.

Any set of `m` values or fewer is a solution. More formally, the independent
sets are the sets containing at most `m` elements.

It can be formalised as follows:

``\\max \\sum_i \\mathrm{rewards}_i x_i``
``\\mathrm{s.t.} \\sum_i x_i \\leq m, \\quad x \\in \\{0, 1\\}^d``
"""
struct UniformMatroidInstance{T <: Real, O <: CombinatorialObjective} <:
       CombinatorialInstance
    rewards::Vector{T} # TODO: rename as rewards for more consistency.
    m::Int
    objective::O

    function UniformMatroidInstance(
        rewards::Vector{T},
        m::Int,
        objective::O=Maximise(),
    ) where {T <: Real, O <: CombinatorialObjective}
        # Error checking.
        if m < 0
            error("m is less than zero: there is no solution.")
        end

        if m == 0
            error("m is zero: the only solution is to take no items.")
        end

        # Return a new instance.
        return new{T, O}(rewards, m, objective)
    end
end

dimension(i::UniformMatroidInstance) = length(i.rewards)

function copy(
    i::UniformMatroidInstance{T, O};
    rewards::Vector{T}=i.rewards,
    m::Int=i.m,
    objective::CombinatorialObjective=i.objective,
) where {T <: Real, O <: CombinatorialObjective}
    return UniformMatroidInstance(rewards, m, objective)
end

fastest_exact(::UniformMatroidInstance) = GreedyAlgorithm()

# Solution.

struct UniformMatroidSolution{T <: Real, O <: CombinatorialObjective} <:
       CombinatorialSolution
    instance::UniformMatroidInstance{T, O}
    variables::Vector{Int} # Indices to the chosen items.
end

function make_solution(
    i::UniformMatroidInstance{T, O},
    item::Dict{Int, Float64},
) where {T <: Real, O <: CombinatorialObjective}
    items_vector = Int[]
    for (k, v) in item
        if v >= 0.5
            push!(items_vector, k)
        end
    end

    return UniformMatroidSolution(i, items_vector)
end

# Budgeted solution.

struct MinBudgetedUniformMatroidSolution{T <: Real, U <: Real} <:
       MultipleMinBudgetedSolution
    instance::MinimumBudget{UniformMatroidInstance{T, Maximise}, U}
    variables::Vector{Int} # Indices to the chosen items for the min_budget.
    solutions::Dict{Int, Vector{Int}}
end

function MinBudgetedUniformMatroidSolution(
    instance::MinimumBudget{UniformMatroidInstance{T, O}, U},
    items::Vector{Int},
) where {T, O <: CombinatorialObjective, U}
    return MinBudgetedUniformMatroidSolution(
        instance,
        items,
        Dict{Int, Vector{Int}}(),
    )
end

function _check_has_solution_for_each_budget(s::MinBudgetedUniformMatroidSolution{T, U}) where {T, U}
    if length(s.solutions) == 0
        error("The provided solution object does not have one solution " *
              "per value of the budget.")
    end
    return
end

function value(s::MinBudgetedUniformMatroidSolution{T, U}) where {T, U}
    if -1 in s.variables || length(s.variables) == 0
        return -Inf
    end

    return sum(s.instance.instance.rewards[i] for i in s.variables)
end

function items(
    s::MinBudgetedUniformMatroidSolution{T, U},
    budget::Int,
) where {T, U}
    _check_has_solution_for_each_budget(s)
    return s.solutions[budget]
end

function items_all_budgets(
    s::MinBudgetedUniformMatroidSolution{T, U},
    max_budget::Int,
) where {T, U}
    _check_has_solution_for_each_budget(s)
    sol = Dict{Int, Vector{Int}}()
    m = s.instance.instance.m
    for budget in 0:max_budget
        sol[budget] = s.solutions[budget]
    end
    return sol
end

function value(
    s::MinBudgetedUniformMatroidSolution{T, U},
    budget::Int,
) where {T, U}
    _check_has_solution_for_each_budget(s)

    its = items(s, budget)
    if -1 in its || length(its) == 0
        return -Inf
    end
    return sum(s.instance.instance.rewards[i] for i in its)
end

function make_solution(
    i::MinimumBudget{UniformMatroidInstance{T, O}, U},
    item::Dict{Int, Float64},
) where {T, O, U}
    items_vector = Int[]
    for (k, v) in item
        if v >= 0.5
            push!(items_vector, k)
        end
    end

    return MinBudgetedUniformMatroidSolution(i, items_vector)
end


text:
begin
	f500kHz = 500.0kHz
	w100mus = 100.0μs
	wmus100mus=w100mus/μs # window in μs
	tlaser500kHz = uconvert(μs, 1.0/f500kHz)
	tlus2 = tlaser500kHz/μs
end

md" - Parameters:
 - laser frequency= $f500kHz
 - DAQ window= $w100mus
 - Laser pulses every $tlaser500kHz
"

md"# BOLD062\_A2\_500kHz\_730nm\_100us"

md"## Read waveform"

fn8n5k1m = 1  # test file

begin
	fls8n5k1m = glob("../labdata/BOLD062_A2_500kHz_730nm_100us/C1*.csv")
	n8n5k1m = length(fls8n5k1m)

	ffl8n5k1m = 200.0e+6
	wvfm8n5k1m = read_waveform(fls8n5k1m,fn8n5k1m; wmus=wmus100mus, wsym=true)
	first(wvfm8n5k1m,3)
end

md"
- Files in directory = $n8n5k1m
- sampling rate= $(round(sr8n5k1m, sigdigits=2)) GHz
- sampling period = $(round(sp8n5k1m, sigdigits=2) )ns
"

begin
	sr8n5k1m = sampling_rate(wvfm8n5k1m,wmus100mus)/GHz
	sp8n5k1m = sampling_period(wvfm8n5k1m,wmus100mus)/ns
end

md"
Unfiltered wvfm:
- mean = $(round(mean8n5k1m, sigdigits=4))
- std = $(round(std8n5k1m, sigdigits=4))
"
begin
	mean8n5k1m, std8n5k1m, thrp8n5k1m, thrn8n5k1m = thrx(wvfm8n5k1m; nsigma=3.0, waveform="raw");
	pfw8n5k1m = plot_filtered_waveform(wvfm8n5k1m; window=false, waveform="raw", thrp=thrp8n5k1m, thrn=thrn8n5k1m, scatter=false)
	pfft8n5k1m =plot_fft(wvfm8n5k1m,wmus100mus)
	plot(pfw8n5k1m, pfft8n5k1m)
end

md"### Filter signal
- Filter frequency (in Hz) = $ffl8n5k1m

Filtered wvfm:
- mean = $fmean8n5k1m
- std = $fstd8n5k1m
"
begin
    filter_signal_lp!(wvfm8n5k1m, wmus100mus; filtertype="Butterworth", flhz=ffl8n5k1m, n=4)
    first(wvfm8n5k1m,3)
end

begin
    fmean8n5k1m, fstd8n5k1m, fthrp8n5k1m, fthr8n5k1m = thrx(wvfm8n5k1m; nsigma=2.0, waveform="filtered")
    plot_filtered_waveform(wvfm8n5k1m; window=false, thrp=fthrp8n5k1m, thrn=fthr8n5k1m)
end

md"### Select peaks"

sprom8n5k1m = 2.0  #prominence cut

md"
- Peak Prominence cut = $sprom8n5k1m
"

begin
	r8n5k1m = select_filtered_peaks(wvfm8n5k1m, fthrp8n5k1m; promsel = sprom8n5k1m, wsel=0.0);
	if r8n5k1m != Nothing
		wvflt8n5k1m, spk8n5k1m = r8n5k1m
		plot_filtered_peaks(wvflt8n5k1m, spk8n5k1m)
	end
end

md"### Optimizing prominence cut"

begin
	P8n5k1m, W8n5k1m = proms(fls8n5k1m, 1, n8n5k1m;
	                 wmus=wmus100mus, wsym=true, flhz=ffl8n5k1m, promsel=0.5)
	hP8n5k1m = lfi.LaserLab.hist1d(P8n5k1m, 20, 0.0, 10.0)
	hW8n5k1m = lfi.LaserLab.hist1d(W8n5k1m, 20, 5.0, 15.0)
end

lfi.LaserLab.hist1d(hP8n5k1m, "Prominence (promsel=0.5)")

lfi.LaserLab.hist1d(hW8n5k1m, "Width")

md"### Selection"

begin
    spks8n5k1m = select_peaks(fls8n5k1m, 1, n8n5k1m; wmus=wmus100mus, wsym=true, flhz=ffl8n5k1m, promsel=sprom8n5k1m)
    t8n5k1m = reduce(vcat,[subtract_laser(spks8n5k1m[i].xs, tlus2) for i in 1:length(spks8n5k1m)])
    hp8n5k1m = lfi.LaserLab.hist1d(t8n5k1m, 20, 0.0, 2.0)
end
lfi.LaserLab.hist1d(hp8n5k1m, "times")

md"### Fits"

begin
    cofe8n5k1m, stder8n5k1m, tft8n5k1m = fit_peaks(hp8n5k1m; pa0=[100.0, 0.5], i0=1)
    pp8n5k1m = plot_fit(hp8n5k1m, cofe8n5k1m, tft8n5k1m; savefig=true,
	              fn="b62A2_500kHz_730nm_100mus_fit.png")
	plot(pp8n5k1m)
end

max_stars_repo_name: proscrite/LaserLab

max_issues_repo_name: proscrite/LaserLab

max_forks_repo_name: proscrite/LaserLab


text:
type FreeListEntry
  array::AFArray
  isResult::Bool
end

FreeListEntry(arr) = FreeListEntry(arr, false)

type FreeList <: AFImpl
  list::Vector{Vector{FreeListEntry}}
end

FreeList() = FreeList(Vector{Vector{FreeListEntry}}())

current(fl::FreeList) =
  length(fl.list) > 0 ?
      Nullable{Vector{FreeListEntry}}(last(fl.list)) :
      Nullable{Vector{FreeListEntry}}()

function previous(fl::FreeList)
    len = length(fl.list)
    len > 1 ?
        Nullable{Vector{FreeListEntry}}(fl.list[len - 1]) :
        Nullable{Vector{FreeListEntry}}()
end

newScope!(fl::FreeList) = push!(fl.list, Vector{FreeListEntry}())

function register!(fl::FreeList, arr::AFArray)
    curr = current(fl)
    if !isnull(curr)
        push!(get(curr), FreeListEntry(arr))
        return true
    end
    false
end

raiseNoScope() = error("There is no active scope.")

function markResult!(fl::FreeList, arr::AFArray)
  currN = current(fl)
  if !isnull(currN)
    curr = get(currN)
    idx = findfirst(x -> x.array == arr, curr)
    idx == 0 && error("Array isn't registered in the current scope.")
    curr[idx].isResult = true
    # prev = previous(fl)
    # if !isnull(prev)
    # 	push!(get(prev), FreeListEntry(arr))
    # end
    arr
  else
    raiseNoScope()
  end
end

function endScope!(fl::FreeList)
  curr = current(fl)
  if !isnull(curr)
    curr = get(curr)
    for entry in filter(x -> !x.isResult, curr)
      release!(entry.array)
    end
    pop!(fl.list)
  else
    raiseNoScope()
  end
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
# BED Record
# ==========

mutable struct Record
    # data and filled range
    data::Vector{UInt8}
    filled::UnitRange{Int}
    # number of columns
    ncols::Int
    # indexes
    chrom::UnitRange{Int}
    chromstart::UnitRange{Int}
    chromend::UnitRange{Int}
    name::UnitRange{Int}
    score::UnitRange{Int}
    strand::Int
    thickstart::UnitRange{Int}
    thickend::UnitRange{Int}
    itemrgb::UnitRange{Int}
    blockcount::UnitRange{Int}
    blocksizes::Vector{UnitRange{Int}}
    blockstarts::Vector{UnitRange{Int}}
end

"""
    BED.Record()

Create an unfilled BED record.
"""
function Record()
    return Record(
        UInt8[], 1:0, 0,
        # chrom-score
        1:0, 1:0, 1:0, 1:0, 1:0,
        # strand-itemrgb
        0, 1:0, 1:0, 1:0,
        # blockcount-blockstarts
        1:0, UnitRange{Int}[], UnitRange{Int}[])
end

"""
    BED.Record(data::Vector{UInt8})

Create a BED record object from `data`.

This function verifies and indexes fields for accessors.
Note that the ownership of `data` is transferred to a new record object.
"""
function Record(data::Vector{UInt8})
    return convert(Record, data)
end

function Base.convert(::Type{Record}, data::Vector{UInt8})
    record = Record(
        data, 1:0, 0,
        # chrom-score
        1:0, 1:0, 1:0, 1:0, 1:0,
        # strand-itemrgb
        0, 1:0, 1:0, 1:0,
        # blockcount-blockstarts
        1:0, UnitRange{Int}[], UnitRange{Int}[])
    index!(record)
    return record
end

"""
    BED.Record(str::AbstractString)

Create a BED record object from `str`.

This function verifies and indexes fields for accessors.
"""
function Record(str::AbstractString)
    return convert(Record, str)
end

function Base.convert(::Type{Record}, str::AbstractString)
    return convert(Record, Vector{UInt8}(str))
end

function initialize!(record::Record)
    record.filled = 1:0
    record.ncols = 0
    record.chrom = 1:0
    record.chromstart = 1:0
    record.chromend = 1:0
    record.name = 1:0
    record.score = 1:0
    record.strand = 0
    record.thickstart = 1:0
    record.thickend = 1:0
    record.itemrgb = 1:0
    record.blockcount = 1:0
    empty!(record.blocksizes)
    empty!(record.blockstarts)
    return record
end

function GenomicFeatures.Interval(record::Record)
    name = BioCore.seqname(record)
    lpos = BioCore.leftposition(record)
    rpos = BioCore.rightposition(record)
    strd = hasstrand(record) ? GenomicFeatures.strand(record) : GenomicFeatures.STRAND_BOTH
    return GenomicFeatures.Interval(name, lpos, rpos, strd, record)
end

function Base.convert(::Type{GenomicFeatures.Interval}, record::Record)
    return GenomicFeatures.Interval(record)
end

function Base.convert(::Type{GenomicFeatures.Interval{Record}}, record::Record)
    return convert(GenomicFeatures.Interval, record)
end

function isfilled(record::Record)
    return !isempty(record.filled)
end

function Base.:(==)(record1::Record, record2::Record)
    if isfilled(record1) == isfilled(record2) == true
        r1 = record1.filled
        r2 = record2.filled
        return length(r1) == length(r2) && memcmp(pointer(record1.data, first(r1)), pointer(record2.data, first(r2)), length(r1)) == 0
    else
        return isfilled(record1) == isfilled(record2) == false
    end
end

function Base.copy(record::Record)
    return Record(
        record.data[record.filled],
        record.filled,
        record.ncols,
        record.chrom,
        record.chromstart,
        record.chromend,
        record.name,
        record.score,
        record.strand,
        record.thickstart,
        record.thickend,
        record.itemrgb,
        record.blockcount,
        copy(record.blocksizes),
        copy(record.blockstarts))
end

function Base.write(io::IO, record::Record)
    return unsafe_write(io, pointer(record.data, first(record.filled)), length(record.filled))
end

function Base.print(io::IO, record::Record)
    write(io, record)
    return nothing
end

function Base.show(io::IO, record::Record)
    print(io, summary(record), ':')
    if isfilled(record)
        println(io)
        println(io, "   chromosome: ", chrom(record))
        println(io, "        start: ", chromstart(record))
        print(io, "          end: ", chromend(record))
        if hasname(record)
            println(io)
            print(io, "         name: ", name(record))
        end
        if hasscore(record)
            println(io)
            print(io, "        score: ", score(record))
        end
        if hasstrand(record)
            println(io)
            print(io, "       strand: ", strand(record))
        end
        if hasthickstart(record)
            println(io)
            print(io, "  thick start: ", thickstart(record))
        end
        if hasthickend(record)
            println(io)
            print(io, "    thick end: ", thickend(record))
        end
        if hasitemrgb(record)
            println(io)
            print(io, "     item RGB: ", itemrgb(record))
        end
        if hasblockcount(record)
            bcount = blockcount(record)
            println(io)
            print(io, "  block count: ", bcount)
            for (i, (bsize, bstart)) in enumerate(zip(blocksizes(record), blockstarts(record)))
                println(io)
                print(io, "      [$i]: size=$(bsize), start=$(bstart)")
            end
        end
    else
        print(io, " <not filled>")
    end
end

# Accessor functions
# ------------------

"""
    chrom(record::Record)::String

Get the chromosome name of `record`.
"""
function chrom(record::Record)::String
    checkfilled(record)
    return String(record.data[record.chrom])
end

function haschrom(record::Record)
    return isfilled(record)
end

function BioCore.seqname(record::Record)
    return chrom(record)
end

function BioCore.hasseqname(record::Record)
    return haschrom(record)
end

"""
    chromstart(record::Record)::Int

Get the starting position of `record`.

Note that the first base is numbered 1.
"""
function chromstart(record::Record)::Int
    checkfilled(record)
    return unsafe_parse_decimal(Int, record.data, record.chromstart) + 1
end

function haschromstart(record::Record)
    return isfilled(record)
end

function BioCore.leftposition(record::Record)
    return chromstart(record)
end

function BioCore.hasleftposition(record::Record)
    return haschromstart(record)
end

"""
    chromend(record::Record)::Int

Get the end position of `record`.
"""
function chromend(record::Record)::Int
    checkfilled(record)
    return unsafe_parse_decimal(Int, record.data, record.chromend)
end

function haschromend(record::Record)
    return isfilled(record)
end

function BioCore.rightposition(record::Record)
    return chromend(record)
end

function BioCore.hasrightposition(record::Record)
    return haschromend(record)
end

"""
    name(record::Record)::String

Get the name of `record`.
"""
function name(record::Record)::String
    checkfilled(record)
    if !hasname(record)
        missingerror(:name)
    end
    return String(record.data[record.name])
end

function hasname(record::Record)
    return record.ncols ≥ 4
end

"""
    score(record::Record)::Int

Get the score between 0 and 1000.
"""
function score(record::Record)::Int
    checkfilled(record)
    if !hasscore(record)
        missingerror(:score)
    end
    return unsafe_parse_decimal(Int, record.data, record.score)
end

function hasscore(record::Record)
    return record.ncols ≥ 5
end

"""
    strand(record::Record)::GenomicFeatures.Strand

Get the strand of `record`.
"""
function strand(record::Record)::GenomicFeatures.Strand
    checkfilled(record)
    if !hasstrand(record)
        missingerror(:strand)
    end
    return convert(GenomicFeatures.Strand, Char(record.data[record.strand]))
end

function hasstrand(record::Record)
    return record.ncols ≥ 6
end

function GenomicFeatures.strand(record::Record)
    return strand(record)
end

"""
    thickstart(record::Record)::Int

Get the starting position at which `record` is drawn thickly.

Note that the first base is numbered 1.
"""
function thickstart(record::Record)::Int
    checkfilled(record)
    if !hasthickstart(record)
        missingerror(:thickstart)
    end
    return unsafe_parse_decimal(Int, record.data, record.thickstart) + 1
end

function hasthickstart(record::Record)
    return record.ncols ≥ 7
end

"""
    thickend(record::Record)::Int

Get the end position at which `record` is drawn thickly.
"""
function thickend(record::Record)::Int
    checkfilled(record)
    if !hasthickend(record)
        missingerror(:thickend)
    end
    return unsafe_parse_decimal(Int, record.data, record.thickend)
end

function hasthickend(record::Record)
    return record.ncols ≥ 8
end

"""
    itemrgb(record::Record)::ColorTypes.RGB

Get the RGB value of `record`.

The return type is defined in [ColorTypes.jl](https://github.com/JuliaGraphics/ColorTypes.jl).
"""
function itemrgb(record::Record)::ColorTypes.RGB
    checkfilled(record)
    if !hasitemrgb(record)
        missingerror(:itemrgb)
    end
    return parse_rgbcolor(record.data, record.itemrgb)
end

function hasitemrgb(record::Record)
    return record.ncols ≥ 9
end

function parse_rgbcolor(data::Vector{UInt8}, range::UnitRange{Int})
    function searchcomma(s)
        i = s
        while i ≤ last(range)
            if data[i] == UInt8(',')
                break
            end
            i += 1
        end
        return i
    end
    lo = first(range)
    hi = searchcomma(lo)
    r = unsafe_parse_byte(data, lo:hi - 1)
    if hi > last(range)
        # single value
        g = b = r
    else
        # triplet
        lo = hi + 1
        hi = searchcomma(lo)
        g = unsafe_parse_byte(data, lo:hi - 1)
        lo = hi + 1
        hi = searchcomma(lo)
        b = unsafe_parse_byte(data, lo:hi - 1)
    end
    return ColorTypes.RGB(reinterpret(N0f8, r), reinterpret(N0f8, g), reinterpret(N0f8, b))
end

function unsafe_parse_byte(data::Vector{UInt8}, range::UnitRange{Int})
    val::UInt8 = 0x00
    for i in range
        val = val * 0x0a + (data[i] - UInt8('0'))
    end
    return val
end

"""
    blockcount(record::Record)::Int

Get the number of blocks (exons) in `record`.
"""
function blockcount(record::Record)::Int
    checkfilled(record)
    if !hasblockcount(record)
        missingerror(:blockcount)
    end
    return unsafe_parse_decimal(Int, record.data, record.blockcount)
end

function hasblockcount(record::Record)
    return record.ncols ≥ 10
end

"""
    blocksizes(record::Record)::Vector{Int}

Get the block (exon) sizes of `record`.
"""
function blocksizes(record::Record)::Vector{Int}
    checkfilled(record)
    if !hasblocksizes(record)
        missingerror(:blocksizes)
    end
    return [unsafe_parse_decimal(Int, record.data, r) for r in record.blocksizes]
end

function hasblocksizes(record::Record)
    return record.ncols ≥ 11
end

"""
    blockstarts(record::Record)::Vector{Int}

Get the block (exon) starts of `record`.

Note that the first base is numbered 1.
"""
function blockstarts(record::Record)::Vector{Int}
    checkfilled(record)
    if !hasblockstarts(record)
        missingerror(:blockstarts)
    end
    return [unsafe_parse_decimal(Int, record.data, r) + 1 for r in record.blockstarts]
end

function hasblockstarts(record::Record)
    return record.ncols ≥ 12
end

function checkfilled(record::Record)
    if !isfilled(record)
        throw(ArgumentError("unfilled BED record"))
    end
end

# r"[-+]?[0-9]+" must match `data[range]`.
function unsafe_parse_decimal(::Type{T}, data::Vector{UInt8}, range::UnitRange{Int}) where T <: Signed
    lo = first(range)
    if data[lo] == UInt8('-')
        sign = T(-1)
        lo += 1
    elseif data[lo] == UInt8('+')
        sign = T(+1)
        lo += 1
    else
        sign = T(+1)
    end
    x = zero(T)
    @inbounds for i in lo:last(range)
        x = Base.Checked.checked_mul(x, 10 % T)
        x = Base.Checked.checked_add(x, (data[i] - UInt8('0')) % T)
    end
    return sign * x
end

function memcmp(p1::Ptr, p2::Ptr, n::Integer)
    return ccall(:memcmp, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), p1, p2, n)
end


text:
function filter(obj::MeasurementFilter, arg0::ObservedMeasurement, arg1::SpacecraftState)
    return jcall(obj, "filter", void, (ObservedMeasurement, SpacecraftState), arg0, arg1)
end


text:
# few utils that are useful

export meshgrid, CrossSection, ExtractSubvolume, SubtractHorizontalMean
export ParseColumns_CSV_File, AboveSurface, BelowSurface, VoteMap
export InterpolateDataOnSurface, InterpolateDataFields2D
export RotateTranslateScale!
export DrapeOnTopo, LithostaticPressure!

using NearestNeighbors

"""
    meshgrid(vx,vy,vz)

Computes an (x,y,z)-grid from the vectors (vx,vy,vz).
For more information, see the MATLAB documentation.
"""
function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T},
                     vz::AbstractVector{T}) where {T}
    m, n, o = length(vy), length(vx), length(vz)
    vx = reshape(vx, 1, n, 1)
    vy = reshape(vy, m, 1, 1)
    vz = reshape(vz, 1, 1, o)
    om = ones(Int, m)
    on = ones(Int, n)
    oo = ones(Int, o)
    (vx[om, :, oo], vy[:, on, oo], vz[om, on, :])
end

"""
    CrossSection(Volume::GeoData; dims=(100,100), Interpolate=false, Depth_level=nothing; Lat_level=nothing; Lon_level=nothing; Start=nothing, End=nothing )

Creates a cross-section through a volumetric (3D) `GeoData` object.

- Cross-sections can be horizontal (map view at a given depth), if `Depth_level` is specified
- They can also be vertical, either by specifying `Lon_level` or `Lat_level` (for a fixed lon/lat), or by defining both `Start=(lon,lat)` & `End=(lon,lat)` points.
- `Interpolate` indicates whether we want to simply extract the data from the 3D volume (default) or whether we want to linearly interpolate it on a new grid, which has dimensions as specified in `dims`

# Example:
```julia-repl
julia> Lon,Lat,Depth   =   LonLatDepthGrid(10:20,30:40,(-300:25:0)km);
julia> Data            =   Depth*2;                # some data
julia> Vx,Vy,Vz        =   ustrip(Data*3),ustrip(Data*4),ustrip(Data*5);
julia> Data_set3D      =   GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon, Velocity=(Vx,Vy,Vz)));
julia> Data_cross      =   CrossSection(Data_set3D, Depth_level=-100km)
GeoData
  size  : (11, 11, 1)
  lon   ϵ [ 10.0 : 20.0]
  lat   ϵ [ 30.0 : 40.0]
  depth ϵ [ -100.0 km : -100.0 km]
  fields: (:Depthdata, :LonData, :Velocity)
```

"""
function CrossSection(V::GeoData; dims=(100,100), Interpolate=false, Depth_level=nothing, Lat_level=nothing, Lon_level=nothing, Start=nothing, End=nothing )

    CheckIsVolume(V);       # Check if it is a volume

    if !isnothing(Depth_level)    # Horizontal slice
        CheckBounds(V.depth, Depth_level)
        if Interpolate
            Lon,Lat,Depth = LonLatDepthGrid(    LinRange(minimum(V.lon.val), maximum(V.lon.val), dims[1]),
                                                LinRange(minimum(V.lat.val), maximum(V.lat.val), dims[2]),
                                                Depth_level)
        else
            ind_z   =   argmin(abs.(V.depth.val[1,1,:] .- Depth_level))
            iDepth  =   ind_z:ind_z;
            iLon    =   1:size(V.lon.val,1);
            iLat    =   1:size(V.lat.val,2);
        end
    end

    if !isnothing(Lat_level)   # vertical slice @ given latitude
        CheckBounds(V.lat, Lat_level)
        if Interpolate
            Lon,Lat,Depth = LonLatDepthGrid(    LinRange(minimum(V.lon.val), maximum(V.lon.val), dims[1]),
                                                Lat_level,
                                                LinRange(minimum(V.depth.val), maximum(V.depth.val), dims[2]))
        else
            ind_l   =   argmin(abs.(V.lat.val[1,:,1] .- Lat_level))
            iDepth  =   1:size(V.depth.val,3)
            iLon    =   1:size(V.lon.val,1);
            iLat    =   ind_l:ind_l
        end
    end

    if !isnothing(Lon_level)   # vertical slice @ given longitude
        CheckBounds(V.lon, Lon_level)
        if Interpolate
            Lon,Lat,Depth = LonLatDepthGrid(    Lon_level,
                                                LinRange(minimum(  V.lat.val), maximum(  V.lat.val), dims[1]),
                                                LinRange(minimum(V.depth.val), maximum(V.depth.val), dims[2]))
        else
            ind_l   =   argmin(abs.(V.lon.val[:,1,1] .- Lon_level))
            iDepth  =   1:size(V.depth.val,3)
            iLat    =   1:size(V.lat.val,2);
            iLon    =   ind_l:ind_l
        end
    end

    # diagonal profile defined by start and end lon/lat points
    if !isnothing(Start)
        if isnothing(End)
            error("Also define End coordinates if you indicate starting lon/lat value")
        end
        Interpolate = true; # we must interpolate in this case

        Lon_dum,Lat_p,Depth_p = LonLatDepthGrid(    Start[1],
                                                LinRange(Start[2], End[2], dims[1]),
                                                LinRange(minimum(V.depth.val), maximum(V.depth.val), dims[2]))

        Lon_p,Lat_dum,Depth = LonLatDepthGrid(    LinRange(Start[1], End[1], dims[1]),
                                                Start[2],
                                                LinRange(minimum(V.depth.val), maximum(V.depth.val), dims[2]))

        Lon             =   zeros(dims[1],dims[2],1)
        Lat             =   zeros(dims[1],dims[2],1)
        Depth           =   zeros(dims[1],dims[2],1)*Depth_p[1]

        # We need 3D matrixes for the paraview writing routine to know we are in 3D
        Lon[:,:,1]      =   Lon_p[:,1,:]
        Lat[:,:,1]      =   Lat_p[1,:,:]
        Depth[:,:,1]    =   Depth_p[1,:,:]

    end

    if Interpolate
        # Interpolate data on profile
        DataProfile = InterpolateDataFields(V, Lon, Lat, Depth);
    else
        # extract data (no interpolation)
        DataProfile = ExtractDataSets(V, iLon, iLat, iDepth);
    end

    return DataProfile

end

"""
    ExtractSubvolume(V::GeoData; Interpolate=false, Lon_level=nothing, Lat_level=nothing, Depth_level=nothing, dims=(50,50,50))

Extract or "cuts-out" a piece of a 2D or 3D GeoData set, defined by `Lon`, `Lat` and `Depth` coordinates.

This is useful if you are only interested in a part of a much bigger larger data set.

- `Lon_level`,`Lat_level` and `Depth_level` should be tuples that indicate `(minimum_value, maximum_value)` along the respective direction. If not specified we use the full range.
- By default, `Interpolate=false` and we find the closest indices within the data set (so your new data set will not go exactly from minimum to maximum).
- Alternatively, if `Interpolate=true` we interpolate the data onto a new grid that has dimensions `dims`. This can be useful to compare data sets that are originally given in different resolutions.

# 3D Example with no interpolation:
```julia-repl
julia> Lon,Lat,Depth   =   LonLatDepthGrid(10:20,30:40,(-300:25:0)km);
julia> Data            =   Depth*2;                # some data
julia> Vx,Vy,Vz        =   ustrip(Data*3),ustrip(Data*4),ustrip(Data*5);
julia> Data_set3D      =   GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon, Velocity=(Vx,Vy,Vz)))
GeoData
  size  : (11, 11, 13)
  lon   ϵ [ 10.0 : 20.0]
  lat   ϵ [ 30.0 : 40.0]
  depth ϵ [ -300.0 km : 0.0 km]
  fields: (:Depthdata, :LonData, :Velocity)
julia> Data_extracted = ExtractSubvolume(Data_set3D,Lon_level=(10,12),Lat_level=(35,40))
GeoData
  size  : (3, 6, 13)
  lon   ϵ [ 10.0 : 12.0]
  lat   ϵ [ 35.0 : 40.0]
  depth ϵ [ -300.0 km : 0.0 km]
  fields: (:Depthdata, :LonData, :Velocity)
```
By default it extracts the data points closest to the area defined by Lon_level/Lat_level/Depth_level.

# 3D Example with interpolation:
Alternatively, you can also interpolate the data onto a new grid:
```julia
julia> Data_extracted = ExtractSubvolume(Data_set3D,Lon_level=(10,12),Lat_level=(35,40), Interpolate=true, dims=(50,51,52))
GeoData
  size  : (50, 51, 52)
  lon   ϵ [ 10.0 : 12.0]
  lat   ϵ [ 35.0 : 40.0]
  depth ϵ [ -300.0 km : 0.0 km]
  fields: (:Depthdata, :LonData, :Velocity)
```

"""
function ExtractSubvolume(V::GeoData; Interpolate=false, Lon_level=nothing, Lat_level=nothing, Depth_level=nothing, dims=(50,50,50))

    if isnothing(Lon_level)
        Lon_level = (minimum(V.lon.val), maximum(V.lon.val))
    end
    if isnothing(Lat_level)
        Lat_level = (minimum(V.lat.val), maximum(V.lat.val))
    end
    if isnothing(Depth_level)
        Depth_level = (minimum(V.depth.val), maximum(V.depth.val))
    end
    if Interpolate
        Lon,Lat,Depth   = LonLatDepthGrid(  LinRange(Lon_level[1],      Lon_level[2],   dims[1]),
                                            LinRange(Lat_level[1],      Lat_level[2],   dims[2]),
                                            LinRange(Depth_level[1],    Depth_level[2], dims[3]) );
        Data_extract    =   InterpolateDataFields(V, Lon, Lat, Depth)

    else
        # Don't interpolate
        i_s, i_e    =   argmin(abs.(V.lon.val[:,1,1] .- Lon_level[1])), argmin(abs.(V.lon.val[:,1,1] .- Lon_level[2]))
        iLon        =   i_s:i_e;

        i_s, i_e    =   argmin(abs.(V.lat.val[1,:,1] .- Lat_level[1])), argmin(abs.(V.lat.val[1,:,1] .- Lat_level[2]))
        iLat        =   i_s:i_e;

        i_s, i_e    =   argmin(abs.(V.depth.val[1,1,:] .- Depth_level[1])), argmin(abs.(V.depth.val[1,1,:] .- Depth_level[2]))
        step        =   1;
        if i_e<i_s
            step=-1
        end
        iDepth      =   i_s:step:i_e;
        Data_extract =  ExtractDataSets(V, iLon, iLat, iDepth);
    end

    return Data_extract
end

function CheckBounds(Data::GeoUnit, Data_Cross)

    min_Data, max_Data = minimum(Data.val), maximum(Data.val);
    if Data_Cross < min_Data || Data_Cross>max_Data
        error("Outside bounds [$min_Data : $max_Data]; $Data_Cross")
    end
end

function CheckIsVolume(Volume::GeoData)
    if any(size(Volume.lon).==1)
        error("It appears your GeoData structure is not a volume; can't compute a cross-section")
    end
end

"""
    InterpolateDataFields(V::GeoData, Lon, Lat, Depth)

Interpolates a data field `V` on a grid defined by `Lon,Lat,Depth`
"""
function InterpolateDataFields(V::GeoData, Lon, Lat, Depth)

    Lon_vec     =  V.lon.val[:,1,1];
    Lat_vec     =  V.lat.val[1,:,1];
    Depth_vec   =  V.depth.val[1,1,:];
    if Depth_vec[1]>Depth_vec[end]
        ReverseData = true
    else
        ReverseData = false
    end

    fields_new  = V.fields;
    field_names = keys(fields_new);
    for i = 1:length(V.fields)
        if typeof(V.fields[i]) <: Tuple
            # vector or anything that contains more than 1 field
            data_tuple = fields_new[i]      # we have a tuple (likely a vector field), so we have to loop
            data_array = zeros(size(Lon,1),size(Lon,2),size(Lon,3),length(data_tuple));     # create a 3D array that holds the 2D interpolated values
            unit_array = zeros(size(data_array));

            for j=1:length(data_tuple)
                if ReverseData
                    ndim        =   length(size(data_tuple[j]))
                    interpol    =   LinearInterpolation((Lon_vec, Lat_vec, reverse(Depth_vec)), reverse(ustrip.(data_tuple[j]), dims=ndim) ,extrapolation_bc = Flat());      # create interpolation object
                else
                    interpol    =   LinearInterpolation((Lon_vec, Lat_vec, Depth_vec), ustrip.(data_tuple[j]),extrapolation_bc = Flat());      # create interpolation object
                end
                data_array[:,:,:,j] =   interpol.(Lon, Lat, Depth);
            end
            data_new    = tuple([data_array[:,:,:,c] for c in 1:size(data_array,4)]...)     # transform 3D matrix to tuple

        else
            # scalar field
            if ReverseData
                ndim        =   length(size(V.fields[i]))
                interpol    =   LinearInterpolation((Lon_vec, Lat_vec, reverse(Depth_vec)), reverse(V.fields[i], dims=ndim), extrapolation_bc = Flat(),);            # create interpolation object
            else
                interpol    =   LinearInterpolation((Lon_vec, Lat_vec, Depth_vec), V.fields[i], extrapolation_bc = Flat());            # create interpolation object
            end
            data_new    =   interpol.(Lon, Lat, Depth);                                                 # interpolate data field
        end

        # replace the one
        new_field   =   NamedTuple{(field_names[i],)}((data_new,))                          # Create a tuple with same name
        fields_new  =   merge(fields_new, new_field);                                       # replace the field in fields_new

    end

    # Create a GeoData struct with the newly interpolated fields
    Data_profile = GeoData(Lon, Lat, Depth, fields_new);

    return Data_profile
end

"""
    InterpolateDataFields(V::UTMData, EW, NS, Depth)

Interpolates a data field `V` on a grid defined by `UTM,Depth`
"""
function InterpolateDataFields(V::UTMData, EW, NS, Depth)

    EW_vec      =  V.EW.val[:,1,1];
    NS_vec      =  V.NS.val[1,:,1];
    Depth_vec   =  V.depth.val[1,1,:];
    if Depth_vec[1]>Depth_vec[end]
        ReverseData = true
    else
        ReverseData = false
    end

    fields_new  = V.fields;
    field_names = keys(fields_new);
    for i = 1:length(V.fields)
        if typeof(V.fields[i]) <: Tuple
            # vector or anything that contains more than 1 field
            data_tuple = fields_new[i]      # we have a tuple (likely a vector field), so we have to loop
            data_array = zeros(size(EW,1),size(EW,2),size(EW,3),length(data_tuple));     # create a 3D array that holds the 2D interpolated values
            unit_array = zeros(size(data_array));

            for j=1:length(data_tuple)
                if ReverseData
                    ndim        =   length(size(data_tuple[j]))
                    interpol    =   LinearInterpolation((EW_vec, NS_vec, reverse(Depth_vec)), reverse(ustrip.(data_tuple[j]), dims=ndim) ,extrapolation_bc = Flat());      # create interpolation object
                else
                    interpol    =   LinearInterpolation((EW_vec, NS_vec, Depth_vec), ustrip.(data_tuple[j]),extrapolation_bc = Flat());      # create interpolation object
                end
                data_array[:,:,:,j] =   interpol.(EW, NS, Depth);
            end
            data_new    = tuple([data_array[:,:,:,c] for c in 1:size(data_array,4)]...)     # transform 3D matrix to tuple

        else
            # scalar field
            if ReverseData
                ndim        =   length(size(V.fields[i]))
                interpol    =   LinearInterpolation((EW_vec, NS_vec, reverse(Depth_vec)), reverse(V.fields[i], dims=ndim), extrapolation_bc = Flat(),);            # create interpolation object
            else
                interpol    =   LinearInterpolation((EW_vec, NS_vec, Depth_vec), V.fields[i], extrapolation_bc = Flat());            # create interpolation object
            end
            data_new    =   interpol.(EW, NS, Depth);                                                 # interpolate data field
        end

        # replace the one
        new_field   =   NamedTuple{(field_names[i],)}((data_new,))                          # Create a tuple with same name
        fields_new  =   merge(fields_new, new_field);                                       # replace the field in fields_new

    end

    # Create a GeoData struct with the newly interpolated fields
    Data_profile = UTMData(EW, NS, Depth, fields_new);

    return Data_profile
end

"""
    InterpolateDataFields2D(V::GeoData, Lon, Lat)

Interpolates a data field `V` on a 2D grid defined by `Lon,Lat`. Typically used for horizontal surfaces
"""
function InterpolateDataFields2D(V::GeoData, Lon, Lat)

    Lon_vec     =  V.lon.val[:,1,1];
    Lat_vec     =  V.lat.val[1,:,1];

    fields_new  = V.fields;
    field_names = keys(fields_new);
    for i = 1:length(V.fields)
        if typeof(V.fields[i]) <: Tuple
            # vector or anything that contains more than 1 field
            data_tuple = fields_new[i]      # we have a tuple (likely a vector field), so we have to loop
            data_array = zeros(size(Lon,1),size(Lon,2),size(Lon,3),length(data_tuple));     # create a 3D array that holds the 2D interpolated values
            unit_array = zeros(size(data_array));

            for j=1:length(data_tuple)
                if length(size(data_tuple[j]))==3
                    interpol    =   LinearInterpolation((Lon_vec, Lat_vec), ustrip.(data_tuple[j][:,:,1]),extrapolation_bc = Flat());      # create interpolation object
                else
                    interpol    =   LinearInterpolation((Lon_vec, Lat_vec), ustrip.(data_tuple[j]),extrapolation_bc = Flat());      # create interpolation object
                end
                data_array[:,:,1,j] =   interpol.(Lon, Lat);
            end
            if length(size(data_tuple[1]))==3
                data_new    = tuple([data_array[:,:,:,c] for c in 1:size(data_array,4)]...)     # transform 3D matrix to tuple
            else
                data_new    = tuple([data_array[:,:,1,c] for c in 1:size(data_array,4)]...)     # transform 3D matrix to tuple
            end
        else
            # scalar field
            if length(size(V.fields[i]))==3
                interpol    =   LinearInterpolation((Lon_vec, Lat_vec), V.fields[i][:,:,1], extrapolation_bc = Flat());            # create interpolation object
            else
                interpol    =   LinearInterpolation((Lon_vec, Lat_vec), V.fields[i], extrapolation_bc = Flat());            # create interpolation object
            end

            data_new    =   interpol.(Lon, Lat);                                                 # interpolate data field
        end

        # replace the one
        new_field   =   NamedTuple{(field_names[i],)}((data_new,))                          # Create a tuple with same name
        fields_new  =   merge(fields_new, new_field);                                       # replace the field in fields_new

    end

    # Interpolate z-coordinate as well
    if length(size(V.lon))==3
        interpol    =   LinearInterpolation((Lon_vec, Lat_vec), V.depth.val[:,:,1], extrapolation_bc = Flat());            # create interpolation object
    else
        interpol    =   LinearInterpolation((Lon_vec, Lat_vec), V.depth.val, extrapolation_bc = Flat());            # create interpolation object
    end
    depth_new =  interpol.(Lon, Lat);

    # Create a GeoData struct with the newly interpolated fields
    # Data_profile = GeoData(Lon, Lat, Depth*0, fields_new);

    return depth_new, fields_new
end

"""
    InterpolateDataFields2D(V::UTMData, EW, NS)

Interpolates a data field `V` on a 2D grid defined by `UTM`. Typically used for horizontal surfaces
"""
function InterpolateDataFields2D(V::UTMData, EW, NS)

    EW_vec      =  V.EW.val[:,1,1];
    NS_vec      =  V.NS.val[1,:,1];

    fields_new  = V.fields;
    field_names = keys(fields_new);
    for i = 1:length(V.fields)
        if typeof(V.fields[i]) <: Tuple
            # vector or anything that contains more than 1 field
            data_tuple = fields_new[i]      # we have a tuple (likely a vector field), so we have to loop
            data_array = zeros(size(EW,1),size(EW,2),size(EW,3),length(data_tuple));     # create a 3D array that holds the 2D interpolated values
            unit_array = zeros(size(data_array));

            for j=1:length(data_tuple)
                interpol    =   LinearInterpolation((EW_vec, NS_vec), ustrip.(data_tuple[j]),extrapolation_bc = Flat());      # create interpolation object
                data_array[:,:,1,j] =   interpol.(EW, NS);
            end
            data_new    = tuple([data_array[:,:,1,c] for c in 1:size(data_array,4)]...)     # transform 3D matrix to tuple

        else
            # scalar field
            if length(size(V.fields[i]))==3
                interpol    =   LinearInterpolation((EW_vec, NS_vec), V.fields[i][:,:,1], extrapolation_bc = Flat());            # create interpolation object
            else
                interpol    =   LinearInterpolation((EW_vec, NS_vec), V.fields[i], extrapolation_bc = Flat());            # create interpolation object
            end

            data_new    =   interpol.(EW, NS);                                                 # interpolate data field
        end

        # replace the one
        new_field   =   NamedTuple{(field_names[i],)}((data_new,))                          # Create a tuple with same name
        fields_new  =   merge(fields_new, new_field);                                       # replace the field in fields_new

    end

    # Interpolate z-coordinate as well
    if length(size(V.depth))==3
        interpol    =   LinearInterpolation((EW_vec, NS_vec), V.depth.val[:,:,1], extrapolation_bc = Flat());            # create interpolation object
    else
        interpol    =   LinearInterpolation((EW_vec, NS_vec), V.depth.val, extrapolation_bc = Flat());            # create interpolation object
    end
    depth_new =  interpol.(EW, NS);

    # Create a UTMData struct with the newly interpolated fields
    # Data_profile = UTMData(EW, NS, Depth*0, fields_new);

    return depth_new, fields_new
end

"""
    Surf_interp = InterpolateDataOnSurface(V::ParaviewData, Surf::ParaviewData)

Interpolates a 3D data set `V` on a surface defined by `Surf`. nex
# Example
```julia
julia> Data
ParaviewData
  size  : (33, 33, 33)
  x     ϵ [ -3.0 : 3.0]
  y     ϵ [ -2.0 : 2.0]
  z     ϵ [ -2.0 : 0.0]
  fields: (:phase, :density, :visc_total, :visc_creep, :velocity, :pressure, :temperature, :dev_stress, :strain_rate, :j2_dev_stress, :j2_strain_rate, :plast_strain, :plast_dissip, :tot_displ, :yield, :moment_res, :cont_res)
julia> surf
ParaviewData
  size  : (96, 96, 1)
  x     ϵ [ -2.9671875 : 3.2671875]
  y     ϵ [ -1.9791666666666667 : 1.9791666666666667]
  z     ϵ [ -1.5353766679763794 : -0.69925457239151]
  fields: (:Depth,)
julia> Surf_interp = InterpolateDataOnSurface(Data, surf)
  ParaviewData
    size  : (96, 96, 1)
    x     ϵ [ -2.9671875 : 3.2671875]
    y     ϵ [ -1.9791666666666667 : 1.9791666666666667]
    z     ϵ [ -1.5353766679763794 : -0.69925457239151]
    fields: (:phase, :density, :visc_total, :visc_creep, :velocity, :pressure, :temperature, :dev_stress, :strain_rate, :j2_dev_stress, :j2_strain_rate, :plast_strain, :plast_dissip, :tot_displ, :yield, :moment_res, :cont_res)
```
"""
function InterpolateDataOnSurface(V::ParaviewData, Surf::ParaviewData)

    # Create GeoData structure:
    V_geo               =   GeoData(V.x.val, V.y.val, V.z.val, V.fields)
    V_geo.depth.val     =   ustrip(V_geo.depth.val);

    Surf_geo            =   GeoData(Surf.x.val, Surf.y.val, Surf.z.val, Surf.fields)
    Surf_geo.depth.val  =   ustrip(Surf_geo.depth.val);

    Surf_interp_geo     =   InterpolateDataOnSurface(V_geo, Surf_geo)
    Surf_interp         =   ParaviewData(Surf_interp_geo.lon.val, Surf_interp_geo.lat.val, ustrip.(Surf_interp_geo.depth.val), Surf_interp_geo.fields)

    return Surf_interp

end

"""
    Surf_interp = InterpolateDataOnSurface(V::GeoData, Surf::GeoData)

Interpolates a 3D data set `V` on a surface defined by `Surf`
"""
function InterpolateDataOnSurface(V::GeoData, Surf::GeoData)

    Surf_interp = InterpolateDataFields(V, Surf.lon.val, Surf.lat.val, Surf.depth.val)

    return Surf_interp
end

# Extracts a sub-data set using indices
function ExtractDataSets(V::GeoData, iLon, iLat, iDepth)

    Lon     =   zeros(typeof(V.lon.val[1]), length(iLon),length(iLat),length(iDepth));
    Lat     =   zeros(typeof(V.lat.val[1]), length(iLon),length(iLat),length(iDepth));
    Depth   =   zeros(typeof(V.depth.val[1]), length(iLon),length(iLat),length(iDepth));

    iLo                 =   1:length(iLon);
    iLa                 =   1:length(iLat);
    iDe                 =   1:length(iDepth)
    Lon[iLo,iLa,iDe]    =     V.lon.val[iLon, iLat, iDepth];
    Lat[iLo,iLa,iDe]    =     V.lat.val[iLon, iLat, iDepth];
    Depth[iLo,iLa,iDe]  =   V.depth.val[iLon, iLat, iDepth];

    fields_new  = V.fields;
    field_names = keys(fields_new);
    for i = 1:length(V.fields)
        if typeof(V.fields[i]) <: Tuple
            # vector or anything that contains more than 1 field
            data_tuple = fields_new[i]      # we have a tuple (likely a vector field), so we have to loop
            data_array = zeros(typeof(data_tuple[1][1]),length(iLon),length(iLat),length(iDepth),length(data_tuple));     # create a 3D array that holds the 2D interpolated values
            unit_array = zeros(size(data_array));

            for j=1:length(data_tuple)
                data_field           =   data_tuple[j];
                data_array[:,:,:,j]  =   data_field[iLon, iLat, iDepth];
            end
            data_new    = tuple([data_array[:,:,:,c] for c in 1:size(data_array,4)]...)       # transform 4D matrix to tuple

        else
            # scalar field
            data_new                =   zeros(typeof(V.fields[i][1]), length(iLon),length(iLat),length(iDepth));
            data_new[iLo,iLa,iDe]   =   V.fields[i][iLon, iLat, iDepth]                                 # interpolate data field
        end

        # replace the one
        new_field   =   NamedTuple{(field_names[i],)}((data_new,))                          # Create a tuple with same name
        fields_new  =   merge(fields_new, new_field);                                       # replace the field in fields_new

    end

    # Create a GeoData struct with the newly interpolated fields
    Data_profile = GeoData(Lon, Lat, Depth, fields_new);

end

"""
    V_sub = SubtractHorizontalMean(V::AbstractArray{T, 3}; Percentage=false)

Subtracts the horizontal average of the 3D data array V.

If `Percentage=true`, the result is given as percentage; otherwise absolute values are returned

"""
function SubtractHorizontalMean(V::AbstractArray{T, 3}; Percentage=false) where T

    nx        = size(V,1);
    ny        = size(V,2);
    NumLayers = size(V,3); # get the number of depth levels

    if Percentage
        V_sub     = zeros(size(V));                 # no units
    else
        V_sub     = zeros(typeof(V[1]), size(V));
    end

    for iLayer = 1:NumLayers
        average             =   mean(filter(!isnan, vec(V[:,:,iLayer])));

        if Percentage
            V_sub[:,:,iLayer]   =   ustrip(V[:,:,iLayer]) .- ustrip(average);
            V_sub[:,:,iLayer]   =   V_sub[:,:,iLayer]./ustrip(average)*100.0;     # the result is normalized
        else
            V_sub[:,:,iLayer]   =   V[:,:,iLayer] .- average;
        end
    end

    return V_sub
end

"""
    V_sub = SubtractHorizontalMean(V::AbstractArray{T, 2}; Percentage=false)

Subtracts the horizontal average of the 2D data array V.

If `Percentage=true`, the result is given as percentage; otherwise absolute values are returned

"""
function SubtractHorizontalMean(V::AbstractArray{T, 2}; Percentage=false) where T

    nx        = size(V,1);
    NumLayers = size(V,2); # get the number of depth levels

    if Percentage
        V_sub     = zeros(size(V));                 # no units
    else
        V_sub     = zeros(typeof(V[1]), size(V));
    end

    for iLayer = 1:NumLayers
        average             =   mean(filter(!isnan, vec(V[:,iLayer])));

        if Percentage
            V_sub[:,iLayer]   =   ustrip(V[:,iLayer]) .- ustrip(average);
            V_sub[:,iLayer]   =   V_sub[:,iLayer]./ustrip(average)*100.0;     # the result is normalized
        else
            V_sub[:,iLayer]   =   V[:,iLayer] .- average;
        end
    end

    return V_sub
end

"""
    ParseColumns_CSV_File(data_file, num_columns)

This parses numbers from CSV file that is read in with `CSV.File`.
That is useful in case the CSV files has tables that contain both strings (e.g., station names) and numbers (lat/lon/height) and you are only intested in the numbers

# Example
This example assumes that the data starts at line 18, that the colums are separated by spaces, and that it contains at most 4 columns with data:
```julia-repl
julia> using CSV
julia> data_file        =   CSV.File("FileName.txt",datarow=18,header=false,delim=' ')
julia> data = ParseColumns_CSV_File(data_file, 4)
```

"""
function ParseColumns_CSV_File(data_file, num_columns)
    data                =   zeros(size(data_file,1), num_columns);
    for (row_num,row) in enumerate(data_file)
        num         =   0;
        for i=1:length(row)
            if typeof(row[i])==Float64
                num          +=  1;
                data[row_num,num] = row[i]
            else

                try parse(Float64,row[i])
                    num          +=  1;
                    data[row_num,num] = parse(Float64,row[i])
                catch
                end
            end

        end
    end
    return data
end

"""
    AboveSurface(Data::GeoData, DataSurface::GeoData; above=true)

Returns a boolean array of size(Data.Lon), which is true for points that are above the surface DataSurface (or for points below if above=false).

This can be used, for example, to mask points above/below the Moho in a volumetric dataset or in a profile.

# Example
First we create a 3D data set and a 2D surface:
```julia
julia> Lon,Lat,Depth   =   LonLatDepthGrid(10:20,30:40,(-300:25:0)km);
julia> Data            =   Depth*2;
julia> Data_set3D      =   GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon))
GeoData
  size  : (11, 11, 13)
  lon   ϵ [ 10.0 : 20.0]
  lat   ϵ [ 30.0 : 40.0]
  depth ϵ [ -300.0 km : 0.0 km]
  fields: (:Depthdata, :LonData)
julia> Lon,Lat,Depth   =   LonLatDepthGrid(10:20,30:40,-40km);
julia> Data_Moho       =   GeoData(Lon,Lat,Depth+Lon*km, (MohoDepth=Depth,))
  GeoData
    size  : (11, 11, 1)
    lon   ϵ [ 10.0 : 20.0]
    lat   ϵ [ 30.0 : 40.0]
    depth ϵ [ -30.0 km : -20.0 km]
    fields: (:MohoDepth,)
```
Next, we intersect the surface with the data set:
```julia
julia> Above       =   AboveSurface(Data_set3D, Data_Moho);
```
Now, `Above` is a boolean array that is true for points above the surface and false for points below and at the surface.

"""
function AboveSurface(Data::GeoData, DataSurface::GeoData; above=true)

    if size(DataSurface.lon)[3]!=1
        error("It seems that DataSurface is not a surface")
    end

    # Create interpolation object for surface
    Lon_vec     =  DataSurface.lon.val[:,1,1];
    Lat_vec     =  DataSurface.lat.val[1,:,1];
    interpol    =  LinearInterpolation((Lon_vec, Lat_vec), ustrip.(DataSurface.depth.val[:,:,1]));            # create interpolation object

    DepthSurface = interpol.(Data.lon.val,Data.lat.val);
    DepthSurface = DepthSurface*unit(DataSurface.depth.val[1])

    if above
        Above       =   Data.depth.val .> DepthSurface;
    else
        Above       =   Data.depth.val .< DepthSurface;
    end

    return Above
end

"""
    Below = BelowSurface(Data::GeoData, DataSurface::GeoData)

Determines if points within the 3D `Data` structure are below the GeoData surface `DataSurface`
"""
function BelowSurface(Data::GeoData, DataSurface::GeoData)
    return AboveSurface(Data::GeoData, DataSurface::GeoData; above=false)
end

"""
    Above = AboveSurface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData; above=true)

Determines if points within the 3D `Data_Cart` structure are above the Cartesian surface `DataSurface_Cart`
"""
function AboveSurface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData; above=true)

    Data            =   GeoData(ustrip.(Data_Cart.x.val),       ustrip.(Data_Cart.y.val),        ustrip.(Data_Cart.z.val), Data_Cart.fields)
    DataSurface     =   GeoData(ustrip.(DataSurface_Cart.x.val),ustrip.(DataSurface_Cart.y.val), ustrip.(DataSurface_Cart.z.val), DataSurface_Cart.fields )

    return Above    =   AboveSurface(Data, DataSurface; above=above)
end

"""
    Above = AboveSurface(Data_Cart::CartData, DataSurface_Cart::CartData; above=true)

Determines if points within the 3D `Data_Cart` structure are above the Cartesian surface `DataSurface_Cart`
"""
function AboveSurface(Data_Cart::CartData, DataSurface_Cart::CartData; above=true)

    Data            =   GeoData(ustrip.(Data_Cart.x.val),       ustrip.(Data_Cart.y.val),        ustrip.(Data_Cart.z.val), Data_Cart.fields)
    DataSurface     =   GeoData(ustrip.(DataSurface_Cart.x.val),ustrip.(DataSurface_Cart.y.val), ustrip.(DataSurface_Cart.z.val), DataSurface_Cart.fields )

    return Above    =   AboveSurface(Data, DataSurface; above=above)
end

"""
    Below = BelowSurface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData)

Determines if points within the 3D Data_Cart structure are below the Cartesian surface DataSurface_Cart
"""
function BelowSurface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData)
    return AboveSurface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData; above=false)
end

"""
    Below = BelowSurface(Data_Cart::CartData, DataSurface_Cart::CartData)

Determines if points within the 3D Data_Cart structure are below the Cartesian surface DataSurface_Cart
"""
function BelowSurface(Data_Cart::CartData, DataSurface_Cart::CartData)
    return AboveSurface(Data_Cart::CartData, DataSurface_Cart::CartData; above=false)
end

"""
    VoteMap(DataSets::Vector{GeoData}, criteria::Vector{String}, dims=(50,50,50))

Creates a Vote map which shows consistent features in different 2D/3D tomographic datasets.

The way it works is:
- Find a common region between the different GeoData sets (overlapping lon/lat/depth regions)
- Interpolate the fields of all DataSets to common coordinates
- Filter data points in one model (e.g., areas with a velocity anomaly > 2 percent). Set everything that satisfies this criteria to 1 and everything else to 0.
- Sum the results of the different datasets

If a feature is consistent between different datasets, it will have larger values.

# Example
We assume that we have 2 seismic velocity datasets `Data_Zhao_Pwave` and `DataKoulakov_Alps`:
```julia
julia> Data_Zhao_Pwave
GeoData
  size  : (121, 94, 101)
  lon   ϵ [ 0.0 : 18.0]
  lat   ϵ [ 38.0 : 51.95]
  depth ϵ [ -1001.0 km : -1.0 km]
  fields: (:dVp_Percentage,)
julia> DataKoulakov_Alps
  GeoData
    size  : (108, 81, 35)
    lon   ϵ [ 4.0 : 20.049999999999997]
    lat   ϵ [ 37.035928143712574 : 49.01197604790419]
    depth ϵ [ -700.0 km : -10.0 km]
    fields: (:dVp_percentage, :dVs_percentage)
```
You can create a VoteMap which combines the two data sets with:
```julia
julia> Data_VoteMap = VoteMap([Data_Zhao_Pwave,DataKoulakov_Alps],["dVp_Percentage>2.5","dVp_percentage>3.0"])
GeoData
  size  : (50, 50, 50)
  lon   ϵ [ 4.0 : 18.0]
  lat   ϵ [ 38.0 : 49.01197604790419]
  depth ϵ [ -700.0 km : -10.0 km]
  fields: (:VoteMap,)
```

You can also create a VoteMap of a single dataset:
```julia
julia> Data_VoteMap = VoteMap(Data_Zhao_Pwave,"dVp_Percentage>2.5", dims=(50,51,52))
GeoData
  size  : (50, 51, 52)
  lon   ϵ [ 0.0 : 18.0]
  lat   ϵ [ 38.0 : 51.95]
  depth ϵ [ -1001.0 km : -1.0 km]
  fields: (:VoteMap,)
```

"""
function VoteMap(DataSets::Vector{GeoData}, criteria::Vector{String}; dims=(50,50,50))

    numDataSets = length(DataSets)

    if length(criteria) != numDataSets
        error("Need the same number of criteria as the number of data sets")
    end

    # Determine the overlapping lon/lat/depth regions of all datasets
    lon_limits  = [minimum(DataSets[1].lon.val);        maximum(DataSets[1].lon.val)];
    lat_limits  = [minimum(DataSets[1].lat.val);        maximum(DataSets[1].lat.val)];
    z_limits    = [minimum(DataSets[1].depth.val);      maximum(DataSets[1].depth.val)];
    for i=1:numDataSets
        lon_limits[1]   =   maximum([lon_limits[1]  minimum(DataSets[i].lon.val)]);
        lon_limits[2]   =   minimum([lon_limits[2]  maximum(DataSets[i].lon.val)]);

        lat_limits[1]   =   maximum([lat_limits[1]  minimum(DataSets[i].lat.val)]);
        lat_limits[2]   =   minimum([lat_limits[2]  maximum(DataSets[i].lat.val)]);

        z_limits[1]     =   maximum([z_limits[1]    minimum(DataSets[i].depth.val)]);
        z_limits[2]     =   minimum([z_limits[2]    maximum(DataSets[i].depth.val)]);
    end

    # Loop over all datasets, and interpolate the data set to the new (usually smaller) domain
    VoteMap             =   zeros(Int64,dims)
    for i=1:numDataSets
        VoteMap_Local   =   zeros(Int64,dims)

        # Interpolate data set to smaller domain
        DataSet         =   ExtractSubvolume(DataSets[i]; Interpolate=true, Lon_level=lon_limits, Lat_level=lat_limits, Depth_level=z_limits, dims=dims);

        # Extract the criteria to evaluate
        expr            =   Meta.parse(criteria[i]);     # the expression, such as Vs>1.0

        # Extract data field
        if !haskey(DataSet.fields,expr.args[2])
            error("The GeoData set does not have the field: $(expr.args[2])")
        end

        Array3D         =   ustrip.(DataSet.fields[expr.args[2]]);                  # strip units, just in case

        # Modify the value, to be Array3D
        expr_mod        =   Expr(:call, expr.args[1], :($Array3D), expr.args[3]);      # modify the original expression to use Array3D as variable name

        # The expression should have a ".", such as Array .> 1.0. If not, it will not apply this in a pointwise manner
        #   Here, we add this dot if it is not there yet
        if cmp(String(expr_mod.args[1])[1],Char('.'))==1
            expr_mod.args[1] = Symbol(".",expr_mod.args[1]);
        end

        ind                 = eval(expr_mod);    # evaluate the modified expression
        VoteMap_Local[ind] .= 1;                 # assign vote-map

        VoteMap = VoteMap + VoteMap_Local;       # Sum
    end

    DataSet     =   ExtractSubvolume(DataSets[1], Interpolate=true, Lon_level=lon_limits, Lat_level=lat_limits, Depth_level=z_limits, dims=dims);

    # Construct GeoData set that holds the VoteMap (makes it easier to write paraview files)
    VoteData    =   GeoData(DataSet.lon.val,DataSet.lat.val,DataSet.depth.val, (VoteMap=VoteMap,));

    return VoteData
end

# Make this work for single data sets as well
function VoteMap(DataSets::GeoData, criteria::String; dims=(50,50,50))
    VoteMap([DataSets], [criteria]; dims=dims)
end

"""
    RotateTranslateScale!(Data::ParaviewData; Rotate=0, Translate=(0,0,0), Scale=(1.0,1.0,1.0))

Does an in-place rotation, translation and scaling of the Cartesian dataset `Data`.

# Parameters
Note that we apply the transformations in exactly this order:
-   `Scale`:        scaling applied to the `x,y,z` coordinates of the data set
-   `Rotate`:       rotation around the `x/y` axis (around the center of the box)
-   `Translate`:    translation

# Example
```julia
julia> X,Y,Z   =   XYZGrid(10:20,30:40,-50:-10);
julia> Data_C  =   ParaviewData(X,Y,Z,(Depth=Z,))
ParaviewData
  size  : (11, 11, 41)
  x     ϵ [ 10.0 : 20.0]
  y     ϵ [ 30.0 : 40.0]
  z     ϵ [ -50.0 : -10.0]
  fields: (:Depth,)
julia> RotateTranslateScale!(Data_C, Rotate=30);
julia> Data_C
ParaviewData
  size  : (11, 11, 41)
  x     ϵ [ 8.169872981077807 : 21.83012701892219]
  y     ϵ [ 28.16987298107781 : 41.83012701892219]
  z     ϵ [ -50.0 : -10.0]
  fields: (:Depth,)
```
"""
function RotateTranslateScale!(Data::ParaviewData; Rotate=0, Translate=(0,0,0), Scale=(1.0,1.0,1.0))

    X,Y,Z       = Data.x.val,   Data.y.val,     Data.z.val;         # Extract coordinates
    Xr,Yr,Zr    = X,Y,Z;                                            # Rotated coordinates

    # 1) Scaling
    if length(Scale)==1
        Scale = [Scale Scale Scale];
    end
    Xr .*= Scale[1];
    Yr .*= Scale[2];
    Zr .*= Scale[3];

    # 2) 2D rotation around X/Y axis, around center of box
    Xm,Ym = mean(X), mean(Y);
    R = [cosd(Rotate[1]) -sind(Rotate[1]); sind(Rotate[1]) cosd(Rotate[1])]; # 2D rotation matrix

    for i in eachindex(X)
        Rot_XY = R*[X[i]-Xm; Y[i]-Ym];
        Xr[i]  = Rot_XY[1] + Xm;
        Yr[i]  = Rot_XY[2] + Ym;
    end

    # 3) Add translation
    Xr .+= Translate[1];
    Yr .+= Translate[2];
    Zr .+= Translate[3];

    # Modify original structure
    Data.x.val = Xr;
    Data.y.val = Yr;
    Data.z.val = Zr;

end

"""
    Topo = DrapeOnTopo(Topo::GeoData, Data::GeoData)

This drapes fields of a data set `Data` on the topography `Topo`

"""
function DrapeOnTopo(Topo::GeoData, Data::GeoData)

    Lon,Lat,Depth    =   LonLatDepthGrid( Topo.lon.val[:,1,1], Topo.lat.val[1,:,1],Topo.depth.val[1,1,:]);

    # use nearest neighbour to interpolate data
    coord       =   [vec(Data.lon.val)'; vec(Data.lat.val)'];
    kdtree      =   KDTree(coord; leafsize = 10);
    points      =   [vec(Lon)';vec(Lat)'];
    idx,dist    =   nn(kdtree, points);

    idx_out     = findall(  (Lon .<  minimum(Data.lon.val)) .| (Lon .>  maximum(Data.lon.val)) .|
                            (Lat .<  minimum(Data.lat.val)) .| (Lat .>  maximum(Data.lat.val)) )

    fields_new  = Topo.fields;
    field_names = keys(Data.fields);

    for i = 1:length(Data.fields)

        if typeof(Data.fields[i]) <: Tuple

            # vector or anything that contains more than 1 field
            data_tuple = Data.fields[i]      # we have a tuple (likely a vector field), so we have to loop

            data_array = zeros(typeof(data_tuple[1][1]),size(Topo.lon.val,1),size(Topo.lon.val,2),size(Topo.lon.val,3),length(Data.fields[i]));     # create a 3D array that holds the 2D interpolated values
            unit_array = zeros(size(data_array));

            for j=1:length(data_tuple)
                data_field           =   data_tuple[j];
                tmp                  =   data_array[:,:,:,1];
                tmp                  =   data_field[idx]
                data_array[:,:,:,j]  =   tmp
            end

            data_new    = tuple([data_array[:,:,:,c] for c in 1:size(data_array,4)]...)       # transform 4D matrix to tuple

            # remove points outside domain
            for j=1:length(data_tuple)
                tmp           =   data_new[j];
                tmp[idx_out] .= NaN
                data_array[:,:,:,j]  =   tmp
            end
            data_new    = tuple([data_array[:,:,:,c] for c in 1:size(data_array,4)]...)       # transform 4D matrix to tuple

        else

            # scalar field
            data_new        =   zeros(typeof(Data.fields[i][1]), size(Topo.lon.val,1),size(Topo.lon.val,2),size(Topo.lon.val,3));
            data_new        =   Data.fields[i][idx]                                 # interpolate data field

        end

        # replace the one
        new_field   =   NamedTuple{(field_names[i],)}((data_new,))                          # Create a tuple with same name
        fields_new  =   merge(fields_new, new_field);                                       # replace the field in fields_new

    end

    Topo_new        =   GeoData(Topo.lon.val,Topo.lat.val,Topo.depth.val, fields_new)

    return Topo_new

end

"""
    DrapeOnTopo(Topo::CartData, Data::CartData)

Drapes Cartesian Data on topography
"""
function DrapeOnTopo(Topo::CartData, Data::CartData)
    Topo_lonlat = GeoData(ustrip.(Topo.x.val),ustrip.(Topo.y.val), ustrip.(Topo.z.val), Topo.fields )
    Data_lonlat = GeoData(ustrip.(Data.x.val),ustrip.(Data.y.val), ustrip.(Data.z.val), Data.fields )

    Topo_new_lonlat = DrapeOnTopo(Topo_lonlat, Data_lonlat)

    Topo_new = CartData(Topo_new_lonlat.lon.val, Topo_new_lonlat.lat.val, Topo_new_lonlat.depth.val, Topo_new_lonlat.fields)

    return Topo_new
end

"""
    LithostaticPressure!(Plithos::Array, Density::Array, dz::Number; g=9.81)

Computes lithostatic pressure from a 3D density array, assuming constant soacing `dz` in vertical direction. Optionally, the gravitational acceleration `g` can be specified.

"""
function LithostaticPressure!(Plithos::Array{T,N}, Density::Array{T,N}, dz::Number; g=9.81) where {T,N}

    Plithos[:] = Density*dz*g;

    selectdim(Plithos,N,size(Plithos)[N]) .= 0      # set upper row to zero

    Plithos[:] = reverse!(cumsum(reverse!(Plithos),dims=N))

    return nothing
end


text:
using Schistoxpkg
using Random
using JLD
using Distributions
#-using Plots

# include parameters and add a file name
include("parameters.jl")
filename = "equ.jld"

# make parameters structure to hold all the parameters for the simulation
pars = Parameters(N, time_step, N_communities, community_probs, community_contact_rate,
        density_dependent_fecundity, average_worm_lifespan,
        max_age, initial_worms, initial_miracidia, initial_miracidia_days, init_env_cercariae,
        worm_stages, contact_rate, max_fec_contact_rate_product, max_fecundity, age_contact_rates,
        ages_for_contacts, contact_rate_by_age_array, mda_adherence, mda_access,  female_factor, male_factor, miracidia_maturity,
        birth_rate, human_cercariae_prop, predis_aggregation, cercariae_survival, miracidia_survival,
        death_prob_by_age, ages_for_death, r, vaccine_effectiveness, drug_effectiveness,
        spec_ages, ages_per_index, record_frequency, use_kato_katz, kato_katz_par, heavy_burden_threshold,
        rate_acquired_immunity, M0, human_larvae_maturity_time, egg_sample_size, egg_production_distribution)
pars = make_age_contact_rate_array(pars, scenario, [4,9,20,100],[0.001,0.6,0.4,0.08]);

# create the larvae variables along with the human structure
humans, miracidia, cercariae = create_population_specified_ages(pars)

# update the ages and death ages, so that there aren't any people with the death age lower than actual age
humans = generate_ages_and_deaths(20000, humans, pars)
humans = update_contact_rate(humans,  pars)

#
number_years = 50
num_time_steps = trunc(Int, 365*number_years / time_step)
num_repeats = 1

# create mda information

mda_info = []

vaccine_info = []
#
@time humans, miracidia, cercariae, record =
        update_env_no_births_deaths_human_larvae(num_time_steps, humans,  miracidia, cercariae, pars, mda_info, vaccine_info)
# save population
save_population_to_file(filename, humans,  miracidia, cercariae, pars)
eggs = (p->p.eggs).(humans)
mean(eggs)
# plot school age children prevalence and heavy burden
pp = (p->p.sac_burden[1]).(record)
pp2 = (p->p.sac_burden[3]).(record)
times = (p->p.time).(record);

plot(times,pp)
plot(times,pp2)

#
# # do simulations with mda for school age children every year for 20 years
#
#         num_repeats = 5 #number of simulations to run
#         number_years = 20
#         drug_efficacy = 0.863 #Toor et al. JID paper in SI: drug efficacy 86.3% for S. mansoni and 94% for S. haematobium
#         num_time_steps = 365*number_years / time_step
#
#         mda_info = create_mda(0, .85, 0.5, 1, number_years, 1, [0,1], [0,1], [0,1], pars.drug_effectiveness)
#
#
#
#         vaccine_info = []
#         # times, prev, sac_prev, high_burden, high_burden_sac, adult_prev, high_adult_burden =
#         # run_repeated_sims_no_births_deaths(filename, num_time_steps, mda_info, vaccine_info, num_repeats)
#         times, prev, sac_prev, high_burden, high_burden_sac, adult_prev, high_adult_burden =
#                 #run_repeated_sims_no_births_deaths_human_larvae(filename, num_time_steps, mda_info, vaccine_info, num_repeats)
# run_repeated_sims_no_population_change(filename, num_time_steps, mda_info, vaccine_info, num_repeats)
#         clf()
#         plt.plot(times, mean.(sac_prev))
#         plt.plot(times, mean.(high_burden_sac))
#         gcf()


text:
"""
MST_PRIM
--------
    compute a minimum spanning tree with Prim's algorithm.
    T = mst_prim_matrix(A) computes a minimum spanning tree T using Prim's algorithm
    for the spanning tree of a graph with non-negative edge weights.

    T = mst_prim_matrix(A,false) produces an MST for just the component at A containing
    vertex 1.  T = mst_prim_matrix(A,0,u) produces the MST for the component
    containing vertex u.

    (ti,tj,tv,nverts) = mst_prim(A) returns the edges from the matrix and does not
    convert to a sparse matrix structure.  This saves a bit of work and is
    required when there are 0 edge weights.

Functions
---------
- (ti,tj,tv,nverts) = mst_prim(A::MatrixNetwork,full::Bool,u::Int64)
- (ti,tj,tv,nverts) = mst_prim(A::MatrixNetwork)
- (ti,tj,tv,nverts) = mst_prim{T}(A::SparseMatrixCSC{T,Int64},full::Bool,u::Int64)
- (ti,tj,tv,nverts) = mst_prim{T}(A::SparseMatrixCSC{T,Int64})
- T = mst_prim_matrix(A::MatrixNetwork,full::Bool,u::Int64) #output modifier
- T = mst_prim_matrix(A::MatrixNetwork) #output modifier
- T = mst_prim_matrix{T}(A::SparseMatrixCSC{T,Int64},full::Bool,u::Int64) #output modifier
- T = mst_prim_matrix{T}(A::SparseMatrixCSC{T,Int64}) #output modifier

Example
-------
~~~
A = load_matrix_network("airports")
A = -A #convert to travel time
A = max.(A,A')
A = sparse(A)
(ti,tj,tv,nverts) = mst_prim(A)
~~~
"""
function mst_prim(A::MatrixNetwork,full::Bool,u::Int64)

    (rp,ci,ai) = (A.rp,A.ci,A.vals)
    if any(ai.<0)
        error("prim''s algorithm cannot handle negative edge weights.")
    end
    nverts = length(rp) - 1
    d = Inf*ones(Float64,nverts)
    T = zeros(Int64,nverts)
    L = zeros(Int64,nverts)
    pred = zeros(Int64,length(rp)-1)

    # enter the main dijkstra loop
    for iter = 1:nverts
        if iter == 1
            root = u
        else
            root = mod(u+iter-1,nverts) + 1
            if L[root] > 0
                continue
            end
        end
        n = 1
        T[n] = root
        L[root] = n # oops, n is now the size of the heap
        d[root] = 0
        while n > 0
            v = T[1]
            L[v] = -1
            ntop = T[n]
            T[1] = ntop
            n = n-1
            if n > 0
                L[ntop] = 1
            end         # pop the head off the heap
            k = 1
            kt = ntop     # move element T[1] down the heap
            while true
                i = 2*k
                if i > n
                    break
                end       # end of heap
                if i == n
                    it = T[i]     # only one child, so skip
                else              # pick the smallest child
                    lc = T[i]
                    rc = T[i+1]
                    it = lc
                    if d[rc] < d[lc]
                        i = i+1
                        it = rc
                    end # right child is smaller
                end
                if d[kt] < d[it]
                    break     # at correct place, so end
                else
                    T[k] = it
                    L[it] = k
                    T[i] = kt
                    L[kt] = i
                    k=i # swap
                end
            end       # end heap down

            # for each vertex adjacent to v, relax it
            for ei = rp[v]:rp[v+1]-1      # ei is the edge index
                w = ci[ei]
                ew = ai[ei]          # w is the target, ew is the edge weight
                if L[w] < 0
                    continue
                end      # make sure we don't visit w twice
                # relax edge (v,w,ew)
                if d[w] > ew
                    d[w] = ew
                    pred[w] = v
                    # check if w is in the heap
                    k = L[w]
                    onlyup = false
                    if k == 0
                        # element not in heap, only move the element up the heap
                        n = n + 1
                        T[n] = w
                        L[w] = n
                        k = n
                        kt = w
                        onlyup = true
                    else
                        kt = T[k]
                    end
                    # update the heap, move the element down in the heap
                    while !onlyup
                        i = 2*k
                        if i > n
                            break
                        end          # end of heap
                        if i == n
                            it = T[i]    # only one child, so skip
                        else             # pick the smallest child
                            lc = T[i]
                            rc = T[i+1]
                            it = lc
                            if d[rc] < d[lc]
                                i = i+1
                                it = rc
                                end # right child is smaller
                            end
                            if d[kt] < d[it]
                                break    # at correct place, so end
                            else
                                T[k] = it
                                L[it] = k
                                T[i] = kt
                                L[kt] = i
                                k = i # swap
                            end
                        end
                        # move the element up the heap
                        j = k
                        tj = T[j]
                        while j > 1                      # j==1 => element at top of heap
                            j2 = convert(Int64,floor(j/2))
                            tj2 = T[j2]                  # parent element
                            if d[tj2] < d[tj]
                                break                    # parent is smaller, so done
                            else                         # parent is larger, so swap
                            T[j2] = tj
                            L[tj] = j2
                            T[j] = tj2
                            L[tj2] = j
                            j = j2
                        end
                    end
                end
            end
        end
        if !full
            break
        end
    end
    nmstedges = 0
    for i = 1:nverts
        if pred[i] > 0
            nmstedges = nmstedges + 1
        end
    end
    ti = zeros(Int64,nmstedges)
    tj = copy(ti)
    tv = zeros(Int64,nmstedges)
    k = 1
    for i = 1:nverts
        if pred[i] > 0
            j = pred[i]
            ti[k] = i
            tj[k] = j
            for rpi = rp[i]:rp[i+1]-1
                if ci[rpi] == j
                    tv[k] = ai[rpi]
                    break
                end
            end
            k = k + 1
        end
    end

    return (ti,tj,tv,nverts)

end

### Support more input:
mst_prim(A::MatrixNetwork) = mst_prim(A,false,1)
mst_prim(A::MatrixNetwork,full::Bool) = mst_prim(A,full,1)

### Support sparse matrices:

mst_prim(A::SparseMatrixCSC{T,Int64},full::Bool,u::Int64) where {T} = mst_prim(MatrixNetwork(A),full,u)
mst_prim(A::SparseMatrixCSC{T,Int64},full::Bool) where {T} = mst_prim(MatrixNetwork(A),full,1)
mst_prim(A::SparseMatrixCSC{T,Int64}) where {T} = mst_prim(MatrixNetwork(A),false,1)

### Output modifiers:
function mst_prim_matrix(A::MatrixNetwork,full::Bool,u::Int64)
    (ti,tj,tv,nverts) = mst_prim(A,full,u)
    T = sparse(ti,tj,tv,nverts,nverts)
    T = T + T'
    return T
end
mst_prim_matrix(A::MatrixNetwork) = mst_prim_matrix(A,false,1)
mst_prim_matrix(A::MatrixNetwork,full::Bool) = mst_prim_matrix(A,full,1)

mst_prim_matrix(A::SparseMatrixCSC{T,Int64},full::Bool,u::Int64) where {T} = mst_prim_matrix(MatrixNetwork(A),full,u)
mst_prim_matrix(A::SparseMatrixCSC{T,Int64},full::Bool) where {T} = mst_prim_matrix(MatrixNetwork(A),full,1)
mst_prim_matrix(A::SparseMatrixCSC{T,Int64}) where {T} = mst_prim_matrix(MatrixNetwork(A),false,1)


text:
function f(x,y)
	x+y
end

mf1 = methods(f, (Any, Any))[end]
ast1 = Base.uncompressed_ast(mf1.func.code)
@show ast1

h(x) = x+1

function f(x::Int64, y)
	a = h(x+y)
    b = x*y
    # Insert a line in between just to see how lines numbered below...
    a/b
end

mf2 = methods(f, (Int64, Any))[end]
ast2 = Base.uncompressed_ast(mf2.func.code)
@show ast2

body = ast2.args[3]

# Now when we have the methods body Expr we can for example instrument it, here using Match...
#using Match
#function instrument(e::Expr)
#	@match e begin
#		Expr(:block, , _)     => name
#		Expr(exprtype, _...)     => error("Cannot instrument $e ($exprtype)")
#	end
#end
#@show instrument(body)


text:
@doc raw"""
    alternating_gradient_descent(M, F, gradF, x)

perform an alternating gradient descent

# Input

* `M` – the product manifold ``\mathcal M = \mathcal M_1 × \mathcal M_2 × ⋯ ×\mathcal M_n``
* `F` – the objective function (cost) defined on `M`.
* `gradF` – a gradient, that can be of two cases
  * is a single function returning a `ProductRepr` or
  * is a vector functions each returning a component part of the whole gradient
* `x` – an initial value ``x ∈ \mathcal M``

# Optional
* `evaluation` – ([`AllocatingEvaluation`](@ref)) specify whether the gradient(s) works by
   allocation (default) form `gradF(M, x)` or [`MutatingEvaluation`](@ref) in place, i.e.
   is of the form `gradF!(M, X, x)` (elementwise).
* `evaluation_order` – (`:Linear`) – whether
  to use a randomly permuted sequence (`:FixedRandom`), a per
  cycle permuted sequence (`:Random`) or the default `:Linear` one.
* `inner_iterations`– (`5`) how many gradient steps to take in a component before alternating to the next
* `stopping_criterion` ([`StopAfterIteration`](@ref)`(1000)`)– a [`StoppingCriterion`](@ref)
* `stepsize` ([`ArmijoLinesearch`](@ref)`()`) a [`Stepsize`](@ref)
* `order` - (`[1:n]`) the initial permutation, where `n` is the number of gradients in `gradF`.
* `retraction_method` – (`default_retraction_method(M)`) a `retraction(M, p, X)` to use.

# Output
* `x_opt` – the resulting (approximately critical) point of gradientDescent
OR
* `options` - the options returned by the solver (see `return_options`)

!!! note

    This Problem requires the `ProductManifold` from `Manifolds.jl`, so `Manifolds.jl` to be loaded.

!!! note

    The input of each of the (component) gradients is still the whole vector `x`,
    just that all other then the `i`th input component are assumed to be fixed and just
    the `i`th components gradient is computed / returned.

"""
function alternating_gradient_descent(
    M::ProductManifold, F, gradF::Union{TgF,AbstractVector{<:TgF}}, x; kwargs...
) where {TgF}
    x_res = copy(M, x)
    return alternating_gradient_descent!(M, F, gradF, x_res; kwargs...)
end
@doc raw"""
    alternating_gradient_descent!(M, F, gradF, x)

perform a alternating gradient descent in place of `x`.

# Input

* `M` a manifold ``\mathcal M``
* `F` – the objective functioN (cost)
* `gradF` – a gradient function, that either returns a vector of the subgradients
  or is a vector of gradients
* `x` – an initial value ``x ∈ \mathcal M``

for all optional parameters, see [`alternating_gradient_descent`](@ref).
"""
function alternating_gradient_descent!(
    M::ProductManifold,
    F,
    gradF::Union{TgF,AbstractVector{<:TgF}},
    x;
    evaluation::AbstractEvaluationType=AllocatingEvaluation(),
    inner_iterations::Int=5,
    stopping_criterion::StoppingCriterion=StopAfterIteration(100) |
                                          StopWhenGradientNormLess(1e-9),
    stepsize::Stepsize=ArmijoLinesearch(M),
    order_type::Symbol=:Linear,
    order=collect(1:(gradF isa Function ? length(gradF(M, x)) : length(gradF))),
    retraction_method::AbstractRetractionMethod=default_retraction_method(M),
    return_options=false,
    kwargs...,
) where {TgF}
    p = AlternatingGradientProblem(M, F, gradF; evaluation=evaluation)
    o = AlternatingGradientDescentOptions(
        M,
        x;
        inner_iterations=inner_iterations,
        stopping_criterion=stopping_criterion,
        stepsize=stepsize,
        order_type=order_type,
        order=order,
        retraction_method=retraction_method,
    )
    o = decorate_options(o; kwargs...)
    resultO = solve(p, o)
    if return_options
        return resultO
    else
        return get_solver_result(resultO)
    end
end
function initialize_solver!(
    ::AlternatingGradientProblem, o::AlternatingGradientDescentOptions
)
    o.k = 1
    o.i = 1
    (o.order_type == :FixedRandom || o.order_type == :Random) && (shuffle!(o.order))
    return o
end
function step_solver!(
    p::AlternatingGradientProblem, o::AlternatingGradientDescentOptions, iter
)
    s, o.gradient = o.direction(p, o, iter)
    j = o.order[o.k]
    retract!(p.M[j], o.x[p.M, j], o.x[p.M, j], -s * o.gradient[p.M, j])
    o.i += 1
    if o.i > o.inner_iterations
        o.k = ((o.k) % length(o.order)) + 1
        (o.order_type == :Random) && (shuffle!(o.order))
        o.i = 1
    end
    return o
end
get_solver_result(o::AlternatingGradientDescentOptions) = o.x


text:
module CubeSatSim

using StaticArrats, DifferentialEquations

include("config.jl")
include("simulation.jl")

export setup

end # module


text:
# initialise.jl
# Define functions to prepare global constants and plotting styles

#=
CONSTANTS

Atomic units:
-------------------------
kb = 1
hbar = 1
4*pi*e0 = 1
a_0 = 1
e = 1
-------------------------
distance = bohr (a_0)
energy = hartree (27.211385 eV)
time = hbar/E_h (2.418884326505(16)×10−17 s)
dipole moment = ea_0 (8.47835326(19)×10−30 C·m = 2.541746 Debye)
polarisability = e2 a_02/E_h (1.65E-41 C^2 m^2 /J)
temperature = E_h/kB (3.1577464E5 K)
electric field = E_h/ea_0 (5.14220652(11)×10^11 V/m)
=#

"""
    init!(num,Edis, Jdis, J, Temp, radius)

Initialise global parameters to define the system

N - (int)  Number of sites in model.
Edisorder - (eV)  Energetic disorder, Gaussian form.
Jdisorder - (eV)  Transfer integral disorder.
J0 - (eV) Transfer Integral.
T - (K) Temperature
radius - (Bohr)

Edits const global:
N::Int
Edisorder - (Hartree)
Jdisorder - (Hartree)
J0 - (Hartree)
T - (E_h/kB)
B - (1/E_h)
r - (Bohr)
"""
function init!(num=50, Edis=0.0, Jdis=0.0, J=0.5, Temp=300, radius=1)
    const global N=num
    const global Edisorder=Edis
    const global Jdisorder=Jdis
    const global J0=J/27.211385
    const global modelJ = J0*cos(π/180.0).^2
    const global T=Temp/3.1577464E5
    const global B=1/T
    const global r=radius
    # This effectively reduces down to the 'alpha' parameter in the Frohlich polaron Hamiltonian
    # Contains the exponential decay of the dipole response (exp(-TimeStep*DecayRate))
    const global dipolestrength=0.2
    const global field_ext = 5E-8#0.6/5.1422065E11

    return N, Edisorder, Jdisorder, modelJ, B, dipolestrength, r
end

"""
    randH(SiteEnergy, Edisorder, Jdisorder, modelJ, N)

Originally borrowed from 'Sturm': https://github.com/jarvist/Teclo/blob/master/Sturm.jl

Determine on-site and nearest neighbour energies for TightBinding
Hamiltonian

SiteEnergy - self-energy of sites (Hartree)
Edisorder  - scale of disorder of on-site energies (Hartree)
Jdisorder  - scale of disorder of transfer integral (Hartree)
modelJ     - (function) for Transfer Integral for nearest neighbours (Hartree)
N::Int     - Number of sites in the system

returns S,E

S - Diagonal elements of TightBinding Hamiltonian (Hartree)
E - Off-diagonal elements of TightBinding Hamiltonian (Hartree)
"""
function randH(SiteEnergy, Edisorder, Jdisorder, modelJ, N)

    S=SiteEnergy + Edisorder*randn(N)
    E=modelJ + Jdisorder*randn(N-1)
    return (S,E)
end

"""
    prepare_model()

Convenience function to put together everything for initiating a 1D model.

Returns S,E,H,psi,dipoles

S       - Diagonal elements of TightBinding Hamiltonian (Hartree)
E       - Off-Diagonal elements of TightBinding Hamiltonian (Hartree)
H       - TightBinding Hamiltonian (Hartree)
psi     - wavefunction
dipoles - initial dipole matrix of zeros (e*a_0)
"""
function prepare_model()
    S,E = randH(0.05, Edisorder, Jdisorder, modelJ, N)
    H = diagm(E,-1)+diagm(S)+diagm(E,1) #build full matrix from diagonal elements
    psi = eigvecs(H)[:,1] # 1=gnd state

    #    plot_H(H)
    #    TODO: UPDATE THIS FOR PLOTS.jl

    ## Construct initial dipoles structure
    dipoles=zeros(N)

    return S,E,H,psi,dipoles
end

"""
    plot_model(S,psi,density,dipoles;title="",verbose::Bool=false)

Wrapper function to pretty-print and plot time evolved state. .
Currently  displays energy, wavefunction, electron density and dipole

S - On-site energies (Hartrees)
psi - Wavefunction
density - Electron Density
dipoles - Site Dipoles

Directly outputs Plots.jl plots.
"""
function plot_model(S,psi,density,dipoles;title="",verbose::Bool=false)
    if verbose
        println("Site energies: ",S)
        println("Psi: ",psi)
        println("Electron density: ",density," Sum: ",sum(density))
        println("Dipoles: ",dipoles," Sum: ",sum(dipoles))
    end

    # Using Plots interface

    plot(S/norm(S),label="Site Energies",color=:red, width=2, yaxis = (" ", (-1,1)), size=(1024,768)) # makes it larger, but also slower!
    plot!(density,label="Electon Density",color=:orange, width=4, fill=true, fill=(0,0.3,:orange))
    plot!(dipoles/norm(dipoles),label="Dipoles",color=:blue, width=2)

    plt = twinx()
    psioffset=-0.0 # Shift wavefunctions down
    plot!(plt,real(psi)+psioffset,label="Re[Psi]", legend = :bottomright, color=:green , width=2, fill=true, fill=(0,0.3,:green),yaxis=(" ",(-0.5,1.5)))
    plot!(plt,imag(psi)+psioffset,label="Im[Psi]",color=:pink, width=2, fill=true, fill=(0,0.3,:pink))

    xaxis!("Site")
    title!(title)

    gui()

end