text:
# Wrapper classes have a pyo field that holds a PyObject
abstract type Wrapper end

# Helper function for use with map(pyo, list-of-Wrappers)
pyo(obj) = obj.pyo

#
# Wrap Python objects in a thin julia wrapper so we can use the type system
#
struct Model <: Wrapper
    pyo::PyObject

    function Model(f)
        return new(py"JuliaModel($f)")
    end
end

struct Parameter <: Wrapper
    pyo::PyObject

    function Parameter(name::AbstractString, default_value::Any=nothing)
        return new(rhodium.Parameter(name, default_value=default_value))
    end
end

Parameter(name::Symbol, default_value::Any=nothing) = Parameter(String(name), default_value)

Parameter(pair::Pair{Symbol, Any}) = Parameter(pair.first, pair.second)

struct Response <: Wrapper
    pyo::PyObject

    function Response(name::AbstractString, kind::Symbol)
        resp = if kind==:MAXIMIZE
            rhodium.Response.MAXIMIZE
        elseif kind==:MINIMIZE
            rhodium.Response.MINIMIZE
        elseif kind==:INFO
            rhodium.Response.INFO
        else
            error("The kind argument must be either :MAXIMIZE, :MINIMIZE or :INFO")
        end

        return new(rhodium.Response(name, resp))
    end

end

abstract type Lever  <: Wrapper end

struct IntegerLever <: Lever
    pyo::PyObject

    function IntegerLever(name::AbstractString, min_value::Int, max_value::Int; length::Int=1)
        return new(rhodium.IntegerLever(name, min_value, max_value, length=length))
    end
end

struct RealLever <: Lever
    pyo::PyObject

    function RealLever(name::AbstractString, min_value::Float64, max_value::Float64; length::Int=1)
        return new(rhodium.RealLever(name, min_value, max_value, length=length))
    end
end

struct CategoricalLever <: Lever
    pyo::PyObject

    function CategoricalLever(name::AbstractString, categories)
        return new(rhodium.CategoricalLever(name, categories))
    end
end

struct PermutationLever <: Lever
    pyo::PyObject

    function PermutationLever(name::AbstractString, options)
        return new(rhodium.PermutationLever(name, options))
    end
end

struct SubsetLever <: Lever
    pyo::PyObject

    function SubsetLever(name::AbstractString, options, size)
        return new(rhodium.SubsetLever(name, options, size))
    end
end

struct Constraint <: Wrapper
    pyo::PyObject

    function Constraint(con::AbstractString)
        return new(rhodium.Constraint(con))
    end

end

Constraint(con::Symbol) = Constraint(String(con))

struct Brush <: Wrapper
    pyo::PyObject

    function Brush(def::AbstractString)
        return new(rhodium.Brush(def))
    end
end

# In rhodium, a dataset is a subclass of list that holds only dicts.
struct DataSet <: Wrapper
    pyo::PyObject   # a python DataSet

    function DataSet(pyo::PyObject)
        return new(pyo)
    end

    # A string argument is interpreted in the python func as a file to load
    function DataSet(data::Union{AbstractString, AbstractArray}=nothing)
        new(rhodium.DataSet(data))
    end
end

Base.length(ds::DataSet) = length(ds.pyo)

function Base.iterate(ds::DataSet, state=1)
    if state>length(ds)
        return nothing
    else
        return get(ds.pyo, state-1), state + 1
    end
end

Base.getindex(ds::DataSet, i::Int) = get(ds.pyo, i-1)

function Base.findmax(ds::DataSet, key::Symbol)
    return pycall(ds.pyo.find_max, PyDict, String(key))
end

function Base.findmin(ds::DataSet, key::Symbol)
    return pycall(ds.pyo.find_min, PyDict, String(key))
end

# TODO This was a method on Base.find previously, maybe it should extend
# some other Base method?
function find(ds::DataSet, expr; inverse=false)
    return pycall(ds.pyo.find, PyObject, expr, inverse=inverse)
end

# Create a NamedTuple type expression from the contents of the given dict,
# returning the type or, if evaluate == false, the type expression.
function make_NT_type(dict::Union{Dict, PyDict})
    names = [Symbol(i) for i in keys(dict)]
    types = [typeof(i) for i in values(dict)]
    return NamedTuple{tuple(names...), Tuple{types...}}
end

function named_tuple(d::Union{Dict, PyDict})
    T = make_NT_type(d)
    return T(tuple(values(d)...))
end

function named_tuples(ds::DataSet)
    T = make_NT_type(ds[1])
    output = [T(tuple(values(dict)...)) for dict in ds]
end

#
# To some, "type piracy". To others, useful conversion methods. ;~)
#
DataFrames.DataFrame(ds::DataSet) = DataFrames.DataFrame(named_tuples(ds))

function DataStructures.Dict(nt::NamedTuple)
    return DataStructures.Dict{Symbol, Any}(collect(k => v for (k, v) in zip(keys(nt), values(nt))))
end

function DataStructures.Dict(pydict::PyDict)
    return DataStructures.Dict(Symbol(k) => v for (k, v) in pydict)
end

function pandas_dataframe(ds::DataSet; include=nothing, exclude=nothing)
    # TBD: Convert directly to np.recarray? That's the end result
    # on the python side, so we'd avoid building 2 intermediate DFs.
    df = DataFrame(ds)

    # However... numpy's drop_names is failing, complaining about
    # 'data type not understood',). We process the include/exclude
    # args here before passing the data to python.
    if include != nothing
        if ! (include isa AbstractArray)
            include = [include]
        end

        colnames = [Symbol(name) for name in include]
        df = df[colnames]
    end

    if exclude != nothing
        if ! (exclude isa AbstractArray)
            exclude = [exclude]
        end

        colnames = collect(setdiff(Set(names(df)), Set(map(Symbol, exclude))))
        df = df[colnames]
    end

    dict = Dict(k => df[k] for k in names(df))
    pandasDF = pd.DataFrame(dict)
    return pandasDF
end

"""
    set_parameters!(m::Model, parameters::Vector{Parameter})
    set_parameters!{T<:Union{Symbol,Pair{Symbol,Any}}}(m::Model, parameters::Vector{T})

Set model parameters using one of these forms:

  set_parameters!(m, [Parameter("a"), Parameter("b")...])

  # create parameters with the given names. Default values are `nothing`
  set_parameters!(m, [:a, :b, "c", ...])

  # create parameters with default values
  set_parameters!(m, [:name => 1, :name2 => 10.6])

"""
function set_parameters!(m::Model, parameters::Vector{Parameter})
    m.pyo.parameters = map(pyo, parameters)
    return nothing
end

set_parameters!(m::Model, v::Vector{T}) where {T<:Union{Symbol,Pair{Symbol,Any}}} = set_parameters!(m, map(Parameter, v))

function set_responses!(m::Model, responses::Vector{Response})
    m.pyo.responses = map(pyo, responses)
    return nothing
end

function set_responses!(m::Model, responses::Vector{Pair{Symbol,Symbol}})
    m.pyo.responses = map(responses) do i
        resp = if i.second==:MAXIMIZE
            rhodium.Response.MAXIMIZE
        elseif i.second==:MINIMIZE
            rhodium.Response.MINIMIZE
        elseif i.second==:INFO
            rhodium.Response.INFO
        else
            error("The kind argument must be either :MAXIMIZE, :MINIMIZE or :INFO")
        end
        return rhodium.Response(String(i.first), resp)
    end
    nothing
end

function set_levers!(m::Model, levers::Vector{T}) where T <: Lever
    m.pyo.levers = map(pyo, levers)
    return nothing
end

function set_constraints!(m::Model, constraints::Vector{Constraint})
    m.pyo.constraints = map(pyo, constraints)
    return nothing
end

set_constraints!(m::Model, v::Vector) = set_constraints!(m, map(Constraint, v))

function set_uncertainties!(m::Model, uncertainties::Vector{Pair{Symbol,T}} where T)
    m.pyo.uncertainties = map(uncertainties) do i
        if i.second isa Uniform{Float64}
            rhodium.UniformUncertainty(string(i.first), i.second.a, i.second.b)
        elseif i.second isa DiscreteUniform
            rhodium.IntegerUncertainty(string(i.first), i.second.a, i.second.b)
        else
            error("Distribution type $(typeof(i.second)) is not currently supported by Rhodium")
        end
    end
    return nothing
end

function sample_lhs(m::Model, nsamples::Int)
    # returns a rhodium DataSet (a subclass of list), which holds (python) OrderedDicts
    py_output = pycall(rhodium.sample_lhs, PyObject, m.pyo, nsamples)
    return DataSet(py_output)
end

function optimize(m::Model, algorithm, trials)
    py_output = pycall(rhodium.optimize, PyObject, m.pyo, algorithm, trials)
    return DataSet(py_output)
end

function evaluate(m::Model, policy::Dict)
    return pycall(rhodium.evaluate, PyDict, m.pyo, policy)
end

# If passed a named tuple, return result as one, too
evaluate(m::Model, policy::NamedTuple) = named_tuple(evaluate(m, Dict(policy)))

function evaluate(m::Model, policies::Vector)
    py_output = pycall(rhodium.evaluate, PyObject, m.pyo, policies)
    return DataSet(py_output)
end

function apply(ds::DataSet, expr; update=true)
    res = pycall(ds.pyo.apply, PyVector, expr, update)
    return collect(res)
end

# function apply(m::Model, results::Vector{T} where T<:NamedTuple, expr; update=false)
#     x = [Dict(k => v for (k,v) in zip(keys(result), values(result))) for result in results])
#     res = apply(whatever, expr, update=update)
#     return res
# end

function _add_brush(kwargs, brush)
    if brush != nothing
        return [kwargs...; (:brush, map(pyo, brush))]
    else
        return kwargs
    end
end

function scatter2d(m::Model, ds::DataSet; brush=nothing, kwargs...)
    kwargs2 = _add_brush(kwargs, brush)
    return rhodium.scatter2d(m.pyo, ds.pyo; kwargs2...)
end

function scatter3d(m::Model, ds::DataSet; brush=nothing, kwargs...)
    kwargs2 = _add_brush(kwargs, brush)
    return rhodium.scatter3d(m.pyo, ds.pyo; kwargs2...)
end

function pairs(m::Model, ds::DataSet; brush=nothing, kwargs...)
    kwargs2 = _add_brush(kwargs, brush)
    return rhodium.pairs(m.pyo, ds.pyo; kwargs2...)
end

function parallel_coordinates(m::Model, ds::DataSet; brush=nothing, kwargs...)
    kwargs2 = _add_brush(kwargs, brush)
    return rhodium.parallel_coordinates(m.pyo, ds.pyo; kwargs2...)
end

function use_seaborn(style="darkgrid")
    seaborn.set()
    seaborn.set_style(style)
end


text:
#   Unit "hilbert" of the FourierAnalysis Package for julia language
#
#   MIT License
#   Copyright (c) 2019, Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This unit implements the Hilbert transform and analytic signal
#   estimations for signals of arbitrary length, using FFTW.

"""
```
   (1)
   function analyticsignal( X  :: Union{Vector{T}, Matrix{T}},
                            wl :: Int     = size(X, 1);
                     nonlinear :: Bool    =  false,
                     planner   :: Planner =  getplanner,
                     ⏩       :: Bool    =  true) where T<:Real

   (2)
   function analyticsignal( 𝐗      :: Vector{Matrix{T}},
                            wl     :: Int;
                        nonlinear  ::  Bool      =  false,
                        planner    ::  Planner   =  getplanner,
                        ⏩        ::  Bool      =  true) where T<:Real
```

(1)

Compute the analytic signal(AS) of vector `X` or of
all column vectors of matrix `X` via the FFT and iFFT procedure,
as explained in Marple(1999).
If `wl`=size(X, 1) (default), use the standard method passing to the FFT
and iFFT all samples in `X` altogether, whereas if `wl`<size(X, 1) a Welch-like
method is used (see below).

Return the analytic signal `𝑌`, a complex vector if `X` is a vector
or a complex matrix holding in its columns the analytic signal of the
columns of `X` if `X` is a matrix. `𝑌` has the same number of samples (rows) as `X`.
Contrarely to what is done in the [DSP](https://github.com/JuliaDSP/DSP.jl)
package, the DC level of the signal is removed, therefore,
if the input signal features a non-zero DC level,
the real part of the AS will be equal to the input signal with the
DC level removed. The imaginary part of `𝑌` is the Hilbert transform
of such no-DC `X`.

The Welch-like AS allows an efficient estimation of the AS for vectors and
matrices of any length, that is, even if they are very large; it proceeds
computing the AS on 50% sliding overlapping windows and forming the AS
by retaining the central half of each window. The number of points effectively
used to obtain the final estimation is ``wl÷2`` (integer division)
if `wl` is even, ``wl÷2+1`` otherwise. Using the central half of the
AS computed at each window eliminates edge effects in the range ``wl÷2:end-wl÷2``.
Before running the procedure, ``wl÷2`` zeros are padded at the beginning
and at the end of `X` and trimmed on exiting the function,
in order to reduce edge effects at the beginning
and end of `X`. This procedure does not eliminate these edge effects
completely, thus the first and last ``wl÷2`` samples of the AS estimation
should be discarded. In practice, one requires the AS of a larger data
segment and trims at least ``wl÷2`` samples at the beginning and end of
the estimation. This is done automatically by the [`TFanalyticsignal`](@ref)
function.

!!! note "Nota Bene"
    In order to avoid FFT computation of very long epochs,
    if `wl` > 2^14, then `wl` is set to 2^10.

    See [window length in FFTW](@ref) for setting efficiently argument `wl`.

    The input signal should be previously band-pass or high-pass filtered
    so as not to contain frequency components below the first discrete
    Fourier frequency obtained with windows of `wl` samples, that is,
    below sr/wl Hz, where sr is the sampling rate.

**Optional Keyword Arguments**:

`nonlinear`, if true, the analytic signal is normalized so that its amplitude
is ``1.0`` at all points. This allow non-linear univariate and bivariate estimations
(see [timefrequencyuni.jl](@ref) and [timefrequencybi.jl](@ref)).

`planner` is an instance of the [`Planner`](@ref) object, holding the forward
and backward FFTW plans used to compute the FFTs and the iFFTs.
By default the planner is computed, but it can be passed as an
argumet here if it is pre-computed. This is interesting if the
`analyticsignal` function is to be invoked repeatedly.

if `⏩` is true, the method is run in multi-threaded mode across the series in `X`
if the number of series is at least twice the number of threads Julia
is instructed to use. See [Threads](@ref).

(2)

Compute the analytic signal for all ``k``
multivariate data matrices given as a vector of matrices `𝐗`.
Return a vector of matrices hodling the corresponding analytic signals
as in method (1). The FFT and iFFT plans are computed only once.
The ``k`` matrices in `𝐗` may have different number of columns (i.e., different
number of series) and different number of rows (samples).
However, the number of rows must be larger than `wl` for all of them.
if `⏩` is true, this method run in multi-threaded mode across the
matrices in `𝐗` if the number of matrices is at least twice
the number of threads Julia is instructed to use, otherwise it
tries to run each analytic signal estimation in multi-threaded mode
as per method (1). See [Threads](@ref).

This function is called by the following functions operating
on time-frequency reprsentations: [`TFanalyticsignal`](@ref),
[`TFamplitude`](@ref), [`TFphase`](@ref), [`meanAmplitude`](@ref),
[`concentration`](@ref), [`meanDirection`](@ref), [`comodulation`](@ref),
[`coherence`](@ref).

**References**
Marple S.L. (1999)
Computing the Discrete-Time Analytic Signal via FFT.
IEEE Transactions on Signal Processing 47(9), 2600-3.

**Examples**:
```
using FourierAnalysis, FFTW, LinearAlgebra, Statistics, Plots, DSP
t=128; lab=["x", "real(y)", "imag(y)"]

# Analytic signal of one vector
x=sinusoidal(10, 2, 128, t, π/2; DC=10) # cosine
y=analyticsignal(x)
# note that real(y) is x without the DC level, i.e., x=real(y)+DC
plot([x, real(y), imag(y)]; labels=lab)

# make a check
s=sinusoidal(10, 2, 128, t, 0) # sine
norm(s-imag(y)) # should be zero

# Analytic Signal by DSP.jl
y2=hilbert(x)
norm(s-imag(y2)) # should be zero
# DSP.jl does not remove the DC level
# thus x=real(y2) in this case
plot([x, real(y2), imag(y2)]; labels=lab)

# Analytic signal of multiple vectors
x=hcat(x, sinusoidal(10, 3, 128, t, π/2; DC=10))
y=analyticsignal(x)

# welch-like analytic signal of one vector
# (note edge effects)
x=sinusoidal(10, 2, 128, t*4, π/2; DC=0)
y=analyticsignal(x, t)
plot([x, real(y), imag(y)]; labels=lab)

# Welch-like analytic signal of multiple vectors
x=hcat(x, sinusoidal(10, 3, 128, t*4, π/2; DC=0))
y=analyticsignal(x, t)
```
"""
function analyticsignal( X  :: Union{Vector{T}, Matrix{T}},
                         wl :: Int     = size(X, 1);
                  nonlinear :: Bool    =  false,
                  planner   :: Planner =  getplanner,
                  ⏩       :: Bool    =  true) where T<:Real

   # get all needed paramters (see below the _paramHT! function)
   X_, 𝚙, i𝚙, e, tₐₗₗ, n, wl½, two_wl⁻¹, cT, f, g, ζ = _paramHT!(X, wl, planner)

   𝑌=zeros(cT, tₐₗₗ, n) # Preallocate memory for the Analytic Signal of X

   # Add to `𝑌` the analytic signal of a vector `X_` or of the column vectors of
   # matrix `X_`, computed on the epoch(s) of duration `t` samples starting at
   # sample `at`. It does not return anything.
   # ARGUMENTS:
   # `n` (Int) is the number of columns in `X_`; it must be 1 if `X_` is a vector.
   # `wl½` = wl÷2
   # `two_wl⁻¹` = 2/wl
   # `e` (Int) the number of sliding windows
   # `𝚙` is the forward FFTW plan performing the FFT.
   # `i𝚙` is the backward FFTW plan performing the iFFT (Hilber transform).
   # `ζ`, =zeros(cT, wl-wl½_), a zero-vector append to the iFFT vectors, since the second half of the FFT is not computed
   # `f`, the lower limit (in samples) of the central half of each Hilbert transform to cumulate with respect to the FFT window
   # `g`, the upper limit (in samples) of the central half of each Hilbert transform to cumulate with respect to the FFT window
   function _analyticsignal!(at, 𝚗)
      y = 𝚙*X_[at:at+wl-1, 𝚗]        # FFT
      @inbounds begin                # transform
         y[1]=0.
         for i=2:wl½ y[i] *= two_wl⁻¹ end
      end
      e==1 ? 𝑌[:, 𝚗]+=i𝚙*vcat(y, ζ) : 𝑌[at+f-1:at+g-1, 𝚗]+=(i𝚙*vcat(y, ζ))[f:g]  # iFFT, cumulate
   end

   # add the Analytic Signal
   as!(𝚎, 𝚗)=_analyticsignal!((𝚎*wl½)+1, 𝚗)
   _thread(⏩, n) ? (for 𝚎=0:e-1 @threads for 𝚗=1:n as!(𝚎, 𝚗) end end) : (for 𝚗=1:n, 𝚎=0:e-1 as!(𝚎, 𝚗) end)

   # return a vector if `X` is a vector, a matrix if `X` is a matrix.
   # for the high-resolution version eliminate the first and last t½ samples that have been previously padded.
   𝑌_ = n==1 ? (e>1 ? 𝑌[:][wl½+1:end-wl½] : 𝑌[:]) : (e>1 ? 𝑌[wl½+1:end-wl½, 1:n] : 𝑌)
   if nonlinear @inbounds for i in eachindex(𝑌_) 𝑌_[i]/=abs(𝑌_[i]) end end
   return 𝑌_
end

function analyticsignal( 𝐗      :: Vector{Matrix{T}},
                         wl     :: Int;
                     nonlinear  ::  Bool      =  false,
                     planner    ::  Planner   =  getplanner,
                     ⏩        ::  Bool      =  true) where T<:Real

   cT, k=typeof(Complex(T(0))), length(𝐗)

   planner ≠ getplanner ? plan=planner : plan=Planner(plan_estimate, -1.0, wl, T, true)

   𝒀=Vector{Matrix{cT}}(undef, k)

   # function to add the analytic signal
   as(i, ⏩)=analyticsignal(𝐗[i], wl; nonlinear=nonlinear, planner=plan, ⏩=⏩)
   _thread(⏩, k) ? (@threads for i=1:k 𝒀[i]=as(i, false) end) : (for i=1:k 𝒀[i]=as(i, ⏩) end)

   return 𝒀
end

# Internal function: prepare all parameters that are needed for
# Analytic Signal and High-Resolution Analytic Signal estimations
#
# ARGUMENTS:
# `X` is the input data vector or matrix (n time-series)
# `inrange` is the range (UnitRange{Int}, e.g., 128:1280-1) of the vectors
#      or of rows of `X` on which the analysis is to be performed
# `wl` (Int) is the epoch length for FFT
# `plan` a FFTW.rFFTWPlan forward plan for computing the FFT
# `plan` a FFTW.cFFTWPlan backward plan for computing the FFT
# `planflag` flags for the computation of the plan by FFYW. See doc of `spectra`
# `plantime` (Int) is the maximum time allowed (in seconds) to optimize the FFTW plan
# `⏩` if true multi-threaded computations are requested
#
# RETURN
# `X_`: if a standard anaytic signal is requested (t==size(X, 1)) then X_=X
#       else X_ is the data in `X` with `t`-1 zeros prepended and appended
# `𝚙` return the forward FFTW plan. It is computed if argument `plan`=`computeplan`
# `i𝚙` return the backward FFTW plan. It is computed if argument `iplan`=`computeiplan`
# `e` (Int) the number of sliding windows for the validated `inrange`, always with 1-sample step
# `n` (Int) 1 if `X` is a Vector, the number of columns of `X` otherwise
# `from` the first sample in the validated `inrange`
# `wl½` `t`÷2 computed by the right aritmetic bit-wise operation `wl`>>1
# `wl½_`, the length of FFT vectors:  t½+1 if DC=true, t½ is DC=false
# `two_wl⁻¹` = 2/wl
# `cT` the complex type corresponding to type `T` in function definition (e.g., ComplexF64 if T=Float64)
# `thrn`: true if ⏩=true and # of threads is >1 and is >= 2*`n`, false otherwise
# `navg`, the number of analytic signal estimations for each sample. This is
#       different from `e` because only the central half of each Hilbert transform is used for each window
# `f`, the lower limit (in samples) of the central half of each Hilbert transform to cumulate with respect to the FFT window
# `g`, the upper limit (in samples) of the central half of each Hilbert transform to cumulate with respect to the FFT window
# `ζ`, =zeros(cT, wl-wl½_), a zero-vector append to the iFFT vectors, since the second half of the FFT is not computed
# NB if wl > 2^14 then t is set to 2^10. This affects the default behavior of all AS functions
function _paramHT!( X       :: AbstractArray{T},
                    wl      :: Int,
                    planner :: Planner) where T<:Real

    size(X, 1) < wl && @error 📌*", the number of samples in input matrix is smaller than the desired window length."
    wl > 2^14 ? wl = 2^10 : nothing
    n, wl½, two_wl⁻¹ = size(X, 2), (wl>>1), T(2/wl)
    if size(X, 1)>wl
        X_=[zeros(T, wl½, n); X; zeros(T, wl½, n)] else X_=X
        isodd(wl) && @error 📌*", for Welch-like Analytic Signal estimation `wl` must be even."
    end
    tₐₗₗ = size(X_, 1) # for high-resolution anaytic signal, this is not size(X, 1)
    cT = typeof(Complex(T(0)))
    if planner ≠ getplanner
       𝚙=planner.p
       i𝚙=planner.ip
    else
       𝚙=plan_rfft(zeros(T, wl), flags=plan_estimate, timelimit=-1.0)
       i𝚙=plan_bfft(zeros(cT, wl), flags=plan_estimate, timelimit=-1.0)
    end
    f=wl½÷2+1 # lower limit of central region to copy
    g=f+wl½-1 # upper limit of central region to copy
    e = (tₐₗₗ-wl)÷wl½+1 # number of 50% overlapping epochs
    ζ=zeros(cT, wl-wl½-1)

    return X_, 𝚙, i𝚙, e, tₐₗₗ, n, wl½, two_wl⁻¹, cT, f, g, ζ
end

ext: jl


text:
# ghg_read.jl
#
#   Load a GHG file
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 1.7.0
# 18.11.2014
# Last Edit: 02.12.2021

"# ghg_read(source::String,verbose::Bool,filetype::String)

`time,data = ghg_read(source)` Load a single GHG file\n
* **source**::String = Single GHG file or a directory of GHG files

---

#### Keywords:\n
* verbose::Bool = Display information as the function runs, TRUE is default
* filetype::String = Data file type to load, \"primary\" is default and loads the primary high frequency file. \"biomet\" would load the BIOMET file if it is present.
* errorlog::String = Log each loading error, \"\" is default
* select::Vector{Symbol} (optional) = List of column names to be loaded, names not included will not be loaded. [] is default\n\n"
function ghg_read(source::String;verbose::Bool=false,filetype::String="primary",errorlog::String="",select::Vector{Symbol}=Symbol[])
	####################
	##  SELECT Check  ##
	####################
	if !isempty(select)
		if filetype == "primary"
			select = [setdiff([:Date,:Time]);select] # Date and Time must be included
		elseif filetype == "biomet"
			select = [setdiff([:DATE,:TIME]);select] # Date and Time must be included
		end
	end

	##################
	##  Initialize  ##
	##################
	epoch = DateTime(1970,1,1,1) # Licor epoch
	D = DataFrame()
	t = DateTime[]
	dfmt = Dates.DateFormat("yyyy-mm-dd HH:MM:SS:sss")
	rootnamereg = r"\d{4}\-\d{2}\-\d{2}T\d{6}_AIU\-\d{4}" # Regular expression for the basic file name

	#################
	##  Load Data  ##
	#################
	# Load Header
	list = []
	try
		list = zipList(source) # List the files in the GHG file
	catch
		@warn("Failed to read " * source)
		if !isempty(errorlog)
			try
				fid3 = open(errorlog,"a+")
				write(fid3,filetype * ", failed to read: " * source * "\r\n")
				close(fid3)
			catch
				@warn("Writing problem \"failed to read\"")
			end
		end
		return Array{DateTime}(undef,0), DataFrame(), ""
	end
	iData = 0
	hf = false # High frequency?
	for i=1:1:length(list)
		if length(list[i]) >= 31
			if occursin(rootnamereg,list[i][1:26])
				if filetype == "primary" && length(list[i]) == 31
					iData = i
					hf = true
				elseif filetype == "biomet" && list[i][end-10:end] == "biomet.data"
					iData = i
					hf = false
				end
			end
		end
	end

	#############
	##  Unzip  ##
	#############
	verbose ? println("\t   Unzipping") : nothing
	temp_dir = tempdir() # Temporary directory for unzipped files
	if iData > 0
		if Sys.iswindows()
			if filetype == "biomet"
				subfile = splitext(basename(source))[1] * "-biomet.data"
			else
				subfile = splitext(basename(source))[1] * ".data"
			end

			#read(`$exe7z x $source -y -o$temp_dir $subfile`); # Temp
			zipExtractFile(source,temp_dir,subfile)
		elseif Sys.islinux()
			if filetype == "primary"
				ext = ".data"
			else
				ext = "-biomet.data"
			end
			try
				run(pipeline(`unzip -p $[source] $[splitext(splitdir(source)[2])[1]]$[ext]`,joinpath(temp_dir,splitext(splitdir(source)[2])[1] * ext)))
			catch
				@warn("Failed to unzip " * source)
				if !isempty(errorlog)
					try
						fid3 = open(errorlog,"a+")
						write(fid3,filetype * ", failed to unzip: " * source * "\r\n")
						close(fid3)
					catch
						@warn("Writing problem \"failed to unzip\"")
					end
				end
				ghg_cleanup(list,temp_dir)
				return Array{DateTime}(undef,0), DataFrame(), ""
			end
		end
	else
		@warn("No Biomet file found: " * source)
		ghg_cleanup(list,temp_dir)
		return Array{DateTime}(undef,0), DataFrame(), ""
	end

	#################
	##  Load Data  ##
	#################
	verbose ? println("\t   Loading " * joinpath(temp_dir,list[iData])) : nothing
	t = Array{DateTime}(undef,0)
	D = DataFrame()
	cols = Array{Any}(undef,0,0)
	fid = open(joinpath(temp_dir,list[iData]),"r")
	if iData != 0
		# Prepare header
		if filetype == "primary"
			header_line = 9
		else
			header_line = 7
		end
		for j=1:1:header_line - 2
			readline(fid)
		end
		if eof(fid)
			# Erroneous Biomet file, it has no column names or data
			@warn("No data: " * joinpath(temp_dir,list[iData]))
			if !isempty(errorlog)
				try
					fid3 = open(errorlog,"a+")
					write(fid3,filetype * ", no data: " * source * "\r\n")
					close(fid3)
				catch
					@warn("Writing problem \"no data\"")
				end
			end
			close(fid)
			ghg_cleanup(list,temp_dir)
			return Array{DateTime}(undef,0), DataFrame(), ""
		end
		l_columns = readline(fid)
		cols = split(l_columns,'\t') # Split column line into column names
		new_names = Array{String}(undef,length(cols))
		temp_name = ""
		for i=1:1:length(cols)
			temp_name = replace(cols[i]," " => "_")
			temp_name = replace(replace(temp_name,"(" => ""),")" => "")
			temp_name = replace(replace(temp_name,"%" => ""),"^" => "")
			temp_name = replace(temp_name,"_-" => "")
			temp_name = replace(temp_name,"/" => "_per_")
			new_names[i] = temp_name
		end

		# Prepare column types
		col_types = fill!(Array{DataType}(undef,length(cols)),Float64)
		if hf
			stringCols = ["DATAH","Date","Time"]
		else
			stringCols = ["DATAH","DATE","TIME"]
		end
		intCols = ["Seconds","Nanoseconds","Diagnostic Value","Diagnostic Value 2","CHK"]
		col_types[cols .∈ (stringCols,)] .= String
		col_types[cols .∈ (intCols,)] .= Int

		# Load Data
		if isempty(select)
			D = CSV.read(joinpath(temp_dir,list[iData]),DataFrame,types = col_types,header = new_names,delim = '\t',skipto = header_line) # # ZipFile.jl temporary fix
		else
			D = CSV.read(joinpath(temp_dir,list[iData]),DataFrame,types = col_types,header = new_names,delim = '\t',skipto = header_line,select = select)
		end

		# Convert Diagnostic_Value Column
		if "Diagnostic_Value" ∈ names(D)
			D.Diagnostic_Value = LI7200Diagnostic.(D.Diagnostic_Value)
		end

		# Convert Time
		t = Array{DateTime}(undef,size(D,1)) # Preallocate
		if hf
			for j=1:1:length(t)
				# Normal Data
				t[j] = DateTime.(D.Date[j] * " " * D.Time[j],dfmt)
			end
		else
			for j=1:1:length(t)
				# Biomet
				t[j] = DateTime.(D.DATE[j] * " " * D.TIME[j],dfmt)
			end
		end
	end
	close(fid)
	#if Sys.iswindows()
	#	# Garbage collection ensures the temporary file is closed so that it can be deleted.
	#	# Windows does not seem to close it in time whereas Linux does
	#	GC.gc()
	#end

	##########################################################
	##  Delete Temporary Files (until ZipFile.jl is fixed)  ##
	##########################################################
	ghg_cleanup(list,temp_dir)

	##############
	##  Output  ##
	##############
	return t,D,cols
end

function ghg_cleanup(list::Array{String,1},temp_dir::String)
	junk = ""
	for i = 1:1:length(list)
		junk = joinpath(temp_dir,list[i])
		if isfile(junk)
			rm(junk)
		end
	end
end

max_stars_repo_name: JeremyRueffer/ClimateDataIO

max_issues_repo_name: JeremyRueffer/ClimateDataIO

max_forks_repo_name: JeremyRueffer/ClimateDataIO


text:
B_data = [  -2.5         0.5         0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   5.26957  -2.45701   -0.812567   0.0        0.0         0.0         0.0         0.0        0.0         0.0         0.0        0.0          0.0         0.0         0.0          0.0        0.0        0.0          0.0          0.0
1.0        -3.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   5.26957  -2.45701   -0.812567   0.0        0.0         0.0         0.0         0.0        0.0         0.0         0.0        0.0          0.0         0.0         0.0          0.0        0.0        0.0          0.0          0.0
0.0         0.0        -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0       0.0        0.0        0.0        0.0         0.0         0.0         0.0        0.0         0.728864   -0.283275   0.0544112    0.0         0.0         0.0          0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0  -2.5  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0       2.0        0.0        0.0        0.0         0.0         0.0         0.0        0.0         1.05443    -0.549114  -0.00531726   0.0         0.0         0.0          0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0       0.0        2.0        0.0        0.0         0.0         0.0         0.0        0.0         0.564472   -0.211292   0.14682      0.0         0.0         0.0          0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        2.0        0.0         0.0         0.0         0.0        0.0         0.0         0.0        0.0          0.728864   -0.283275    0.0544112    0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0  -2.5  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0       0.0        0.0        0.0        2.0         0.0         0.0         0.0        0.0         0.0         0.0        0.0          1.05443    -0.549114   -0.00531726   0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        0.0        0.0         2.0         0.0         0.0        0.0         0.0         0.0        0.0          0.564472   -0.211292    0.14682      0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         2.0         0.0        0.0         0.0         0.0        0.0          0.0         0.0         0.0          0.728864  -0.283275   0.0544112    0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0       0.0        0.0        0.0        0.0         0.0         0.0         2.0        0.0         0.0         0.0        0.0          0.0         0.0         0.0          1.05443   -0.549114  -0.00531726   0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         0.0         0.0        2.0         0.0         0.0        0.0          0.0         0.0         0.0          0.564472  -0.211292   0.14682      0.0          0.0
0.0         0.0         0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.91546  -1.1331     0.217645   0.0        0.0         0.0         0.0         0.0        0.0         0.5         0.0        0.0          0.0         0.0         0.0          0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   4.21772  -2.19645   -0.021269   0.0        0.0         0.0         0.0         0.0        0.0         0.0         0.5        0.0          0.0         0.0         0.0          0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   2.25789  -0.845166   0.587279   0.0        0.0         0.0         0.0         0.0        0.0         0.0         0.0        0.5          0.0         0.0         0.0          0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        2.91546   -1.1331      0.217645    0.0         0.0        0.0         0.0         0.0        0.0          0.5         0.0         0.0          0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5  -0.0   0.0   0.0   0.0   0.0       0.0        0.0        4.21772   -2.19645    -0.021269    0.0         0.0        0.0         0.0         0.0        0.0          0.0         0.5         0.0          0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0       0.0        0.0        2.25789   -0.845166    0.587279    0.0         0.0        0.0         0.0         0.0        0.0          0.0         0.0         0.5          0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         2.91546    -1.1331     0.217645    0.0         0.0        0.0          0.0         0.0         0.0          0.5        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5  -0.0   0.0       0.0        0.0        0.0        0.0         0.0         4.21772    -2.19645   -0.021269    0.0         0.0        0.0          0.0         0.0         0.0          0.0        0.5        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0       0.0        0.0        0.0        0.0         0.0         2.25789    -0.845166   0.587279    0.0         0.0        0.0          0.0         0.0         0.0          0.0        0.0        0.5          0.0          0.0
0.0         0.0         1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -3.50473   0.0        0.0        3.96463   -1.84856    -0.611345    0.0         0.0        0.0         1.45773    -0.56655    0.108822     0.0         0.0         0.0          0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0      -3.50473   -1.55879    8.07172   -3.76354    -1.24466     0.0         0.0        0.0         2.10886    -1.09823   -0.0106345    0.0         0.0         0.0          0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0       1.55879   -3.50473   -0.142449   0.0664188   0.0219656   0.0         0.0        0.0         1.12894    -0.422583   0.29364      0.0         0.0         0.0          0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0       -6.51418    0.0         0.0        11.8939     -5.54569   -1.83404     0.0         0.0        0.0          1.45773    -0.56655     0.108822     0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0       0.0        0.0        0.0       -6.51418    -4.67637    24.2152    -11.2906    -3.73397     0.0         0.0        0.0          2.10886    -1.09823    -0.0106345    0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        0.0        4.67637    -6.51418    -0.427348    0.199257   0.0658969   0.0         0.0        0.0          1.12894    -0.422583    0.29364      0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0        -2.56427     0.0        0.0         0.0         0.0        0.0          0.0         0.0         0.0          1.66769   -0.66445    0.0764458    0.0821063    0.402476
0.0         0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0       0.0        0.0        0.0        0.0         0.0         0.0        -2.56427   -0.584547    0.0         0.0        0.0          0.0         0.0         0.0          2.53634   -1.29754   -0.0765513    0.167163     0.819413
0.0         0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         0.0         0.584547  -2.56427     0.0         0.0        0.0          0.0         0.0         0.0          1.1214    -0.419066   0.294803    -0.00295008  -0.014461
1.63747     0.34258     0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   3.64539  -1.47344    0.105089   0.0        0.0         0.0         0.0         0.0        0.0        -5.25709     0.0        0.0          0.0         0.0         0.0          0.0        0.0        0.0          0.0          0.0
3.33378     0.697469    0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   5.70382  -2.88937   -0.250424   0.0        0.0         0.0         0.0         0.0        0.0         0.0        -5.25709   -2.33819      0.0         0.0         0.0          0.0        0.0        0.0          0.0          0.0
-0.0588344  -0.0123089   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   2.23166  -0.832938   0.591324   0.0        0.0         0.0         0.0         0.0        0.0         0.0         2.33819   -5.25709      0.0         0.0         0.0          0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        2.91546   -1.1331      0.217645    0.0         0.0        0.0        17.8409     -8.31853   -2.75105     -9.77127     0.0         0.0          0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        4.21772   -2.19645    -0.021269    0.0         0.0        0.0        36.3227    -16.9359    -5.60095      0.0        -9.77127    -7.01456      0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0       0.0        0.0        2.25789   -0.845166    0.587279    0.0         0.0        0.0        -0.641022    0.298885   0.0988454    0.0         7.01456    -9.77127      0.0        0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         2.91546    -1.1331     0.217645    0.0         0.0        0.0          2.23011    -1.03982    -0.343882    -3.84641    0.0        0.0          0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         4.21772    -2.19645   -0.021269    0.0         0.0        0.0          4.54034    -2.11699    -0.700119     0.0       -3.84641   -0.87682      0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0       0.0        0.0        0.0        0.0         0.0         2.25789    -0.845166   0.587279    0.0         0.0        0.0         -0.0801277   0.0373606   0.0123557    0.0        0.87682   -3.84641      0.0          0.0
0.0         0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         0.0         0.0        0.0         0.0         0.0        0.0          0.0         0.0         0.0          1.31739   -0.614252  -0.203142    -2.5          2.0
0.0         0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         0.0         0.0        0.0         0.0         0.0        0.0          0.0         0.0         0.0          2.63479   -1.2285    -0.406283     1.0         -2.0]


text:
## Motion Planning Benchmark
## =========================

# Construct a path from the start region to the target region
# which does not pass trhoug nay of the obstacles

typealias Point Lifted{Vector{Float64}}
typealias Vec Lifted{Vector{Float64}}
typealias Mat Lifted{Matrix{Float64}}
typealias PointPair Lifted{Matrix{Float64}}
typealias Edge Lifted{Matrix{Float64}}  #Parametric form of line is p1 + (p2 - p1)
typealias Scalar Lifted{Float64}

points_to_vec(p1::Point, p2::Point) = p1 - p2
points_to_vec(edge::PointPair) = points_to_vec(edge[:,1], edge[:,2])

# Parametric form of line is p1 + (p2 - p1)
points_to_parametric(p1::Point,p2::Point) = liftedarray([p1 points_to_vec(p2,p1)])
points_to_parametric(edge::PointPair) = points_to_parametric(edge[:,1], edge[:,2])
parametric_to_point(p::Edge, s::Scalar) = p[:,1] + s * p[:,2]

# Is the point in the box?
function ispointinpoly(point, box)
  r = (point[1] >= box[1,1]) &
    (point[1] <= box[2,1]) &
    (point[2] >= box[1,2]) &
    (point[2] <= box[2,2])
  r
end

# Linear Check
function avoid_obstacles(points, obs)
  observations = true
  for i = 1:size(points,2),ob in obs
    observations &= !point_in_poly(points[:,i],ob)
#     @show i, observations
  end
  observations
end

# Where if anywhere, along p does it intersect segment
function intersect_segments(ppos::Point, pdir::Vec, qpos::Point, qdir::Vec)
  w = ppos - qpos
  u = pdir
  v = qdir
  (v[2] * w[1] - v[1] * w[2]) / (v[1] * u[2] - v[2] * u[1])
end

function short_lines(points, length::Float64)
  true
end

function isvalid_path(points, obstacles, origin, dest)
  point_in_poly(points[:,1],origin) &
  ispointinpoly(points[:,end], dest) &
  avoid_obstacles(points,obstacles) &
  short_lines(points, .1)
end

function linearplan(points)
  origin = [0.0 0.0
           0.2 0.2]
  dest = [.9 .9
          1.0 1.0]

  box1 = [.3  .1
          .65 .7]
  box2 = [.05 .7
          .25 .99]

  obstacles = Matrix[box1,box2]
  good_path = isvalid_path(points,obstacles,origin,dest)
  points,good_path
end

## Nonlinear
## =========
function isvalid_pathnl(points, obstacles, origin, dest)
  ispointinpoly(points[:,1],origin) &
  ispointinpoly(points[:,end], dest) &
  avoid_obstacles_nl(points,obstacles) &
  short_lines(points, .1)
end

# Nonlinear Check
function avoid_obstacles_nl(points, obs)
  observations = true
  for i = 1:size(points,2)-1, ob in obs
    p = points[:,i]
    pdir = points[:,i+1] - points[:,i]

    obpos = ob[:,1]
    obdir = ob[:,2]
    s = intersect_segments(p, pdir, obpos, obdir)
    observations &= (s < 0) | (s > 1)
  end
  observations
end

function nonlinearplan(points)
  origin = [0.0 0.0
           0.2 0.2]
  dest = [.9 .9
          1.0 1.0]
  obstacles = Array[[8.01 3.01; 1.02 9],
                    [0.5 3.08; 2.02 9.04],
                    [0.0 9.99; 3 5.04]]
  obs = map(points_to_parametric, obstacles)
  good_path = isvalid_pathnl(points,obs,origin,dest)
  points,good_path
end

points = mvuniform(0.0, 10.0, 2,6)
points, condition = nonlinearplan(points)
@show condition.smt
# s = cond_sample_tlmh(points,condition,1)

# # ## Draw
# obstacles = Array[[8.01 3.01; 1.02 9],
#                   [0.5 3.08; 2.02 9.04],
#                   [0.0 9.99; 3 5.04]]
# a = make_compose_lines(obstacles)
# points = s[1]
# b = [Compose.line([pair(points[:,i]),pair(points[:,i+1])]) for i = 1:(size(points,2)-1)]
# b
# draw_lines(a,b)


text:
#generate synthetic psfs

function synth_psf_img(psf::Psf{T,N}, imgradii::NTuple{N,Int}; shift=(zeros(N)...,), bead_max=one(T), bias=zero(T)) where {T,N}
    imginds = roiranges(imgradii)
    imgdat = fill(T(bias), (2 .* imgradii .+ 1))
    img = OffsetArray(imgdat, imginds)
    synth_psf_img!(img, psf; shift=shift, bead_max=bead_max)
end

mayberound(t::Val{T}, v) where {T<:Integer} = round(T, v)
mayberound(t, v) = v

#adds a psf onto existing input img
function synth_psf_img!(img::OffsetArray{TA,N}, psf::Psf{T,N}; shift=(zeros(N)...,), bead_max=one(T)) where {TA,T,N}
    peak0 = psf[ntuple(i->zero(T), N)...]
    for I in CartesianIndices(img)
        img[I] += mayberound(Val(TA), bead_max*psf[(Tuple(I).+shift)...]/peak0)
    end
    return img
end

rng_subtr(r::UnitRange, v::Int) = (first(r)-v):(last(r)-v)
rng_subtr(r, v::Int) = r-v

#because subtracting a tuple of ints from a tuple of ranges doesn't seem to work...
tuple_subtr(t1::Tuple, t2::Tuple) = (rng_subtr(first(t1), first(t2)), tuple_subtr(Base.tail(t1), Base.tail(t2))...)
tuple_subtr(t1::Tuple{}, t2::Tuple{}) = ()

#noise will be uniformly distributed, +- noise_max
function fake_bead_sample(img_size::NTuple{3,Int}, bead_sigma, bead_shift, bead_ctrs::AbstractVector, bead_max, noise_max, bias; img_eltype=UInt16)
    fullimg = zeros(img_eltype, img_size...)
    psf0 = Psf(bead_sigma)
    roi_radii = 2 .* ceil.(Int, bead_sigma)
    fake_bead_sample!(fullimg, psf0, roi_radii, bead_shift, bead_ctrs, bead_max, noise_max, bias)
end

function fake_bead_sample!(fullimg::AbstractArray, psf0::Psf, roi_radii, bead_shift, bead_ctrs, bead_max, noise_max, bias)
    for i in eachindex(fullimg)
        fullimg[i] = round(eltype(fullimg), bias + (rand()-0.5)*2*noise_max)
    end
    rngs0 = ExtractPSF.roiranges(roi_radii) #indices at which to index psf if its roi is contained in image
    for c in bead_ctrs
        crngs = ExtractPSF.roiranges(c, roi_radii) #parent indices we'd like to use
        valid_rngs = intersect.(crngs, axes(fullimg)) #valid indices in parent array
        vshifted = tuple_subtr(valid_rngs, c)
        psf_rngs = intersect.(rngs0, vshifted) #evaluation indices for psf
        imgv = view(fullimg, valid_rngs...)
        oa = OffsetArray(imgv, psf_rngs)
        ExtractPSF.synth_psf_img!(oa, psf0; shift=bead_shift, bead_max=bead_max)
    end
    return fullimg
end

max_stars_repo_name: HolyLab/ExtractPSF

max_forks_repo_name: HolyLab/ExtractPSF


text:
using Theon
using Test

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


text:
# Step 2 Run DICE
using MimiDICE2010

m = MimiDICE2010.get_model()
run(m)

# Step 3: Altering Parameters

using Mimi

# update the forcings of equilibrium CO2
update_param!(m, :fco22x, 3.000)
run(m)

# update the time horizon
const ts = 10 # timestep 10 years
const years = collect(1995:ts:2505)
nyears = length(years)
set_dimension!(m, :time, years)

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance


text:
using TIFF
using FileIO
using Printf
include("julia2ims.jl")

"""
20201216: Package Spinning disk tif stack to imaris
"""

src_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_again/20201215/sCMOS"
out_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_again/20201215"
t_len = 28 #1-28
z_len = 20
pos_len = 40

buffer = zeros(N0f16, 1900, 1300, z_len*t_len*pos_len)
img_num = 866

Threads.@threads for index in 0:24
    println(index)
    filename = "HE7-1NLS-uw-PWM_1_MMStack_$index.ome.tif"
    raw = TIFF.load("$src_dir/$filename")
    id = index*img_num + 1
    buffer[:, :, id:id+img_num-1] = raw.data
    GC.gc()
end

#Threads.@threads for index in 25:54
#    filename = "HE7-1NLS-uw-PWM_1_MMStack_$index.ome.tif"
#    raw = load(File(format"TIFF", "$src_dir/$filename"))
#    id = index*img_num + 1
#    buffer[:, :, id:id+img_num-1] = raw.data
#    GC.gc()
#end

t0 = 25*866
remain = t_len*z_len*pos_len - t0 # 750
raw = load(File(format"TIFF", "$src_dir/HE7-1NLS-uw-PWM_1_MMStack_25.ome.tif"))
buffer[:, :, t0+1:end] = raw[:, :, 1:remain]

buffer = reshape(buffer, (1900, 1300, z_len, t_len*pos_len))
Threads.@threads for pos in 1:pos_len
    imsname = @sprintf("%s/pos%02d-1.tiff", out_dir, pos)
    #save2ims(reinterpret(UInt16, buffer[:, :, :, :, pos]),
    #         imsname, dx=0.108, compression=1)
    save(imsname, reshape(buffer[:, :, :, pos:40:end],
        (1900, 1300, z_len*t_len)))
end


text:
include("solve_method.jl")
include("centrality_corr.jl")
include("regularization.jl")
include("solvers/K2LDL.jl")
include("solvers/K2_5LDL.jl")

function compute_α_dual(v, dir_v)
    n = length(v)
    T = eltype(v)
    if n == 0
        return one(T)
    end
    α = one(T)
    @inbounds @simd for i=1:n
        if dir_v[i] < zero(T)
            α_new = -v[i] * T(0.999) / dir_v[i]
            if α_new < α
                α = α_new
            end
        end
    end
    return α
end

function compute_α_primal(v, dir_v, lvar, uvar)
    n = length(v)
    T = eltype(v)
    α_l, α_u = one(T), one(T)
    @inbounds @simd for i=1:n
        if dir_v[i] > zero(T)
            α_u_new = (uvar[i] - v[i]) * T(0.999) / dir_v[i]
            if α_u_new < α_u
                α_u = α_u_new
            end
        elseif dir_v[i] < zero(T)
            α_l_new = (lvar[i] - v[i]) * T(0.999) / dir_v[i]
            if α_l_new < α_l
                α_l = α_l_new
            end
        end
    end
    return min(α_l, α_u)
end

@inline function compute_αs(x, s_l, s_u, lvar, uvar, Δxy, Δs_l, Δs_u, nvar)
    α_pri = @views compute_α_primal(x, Δxy[1:nvar], lvar, uvar)
    α_dual_l = compute_α_dual(s_l, Δs_l)
    α_dual_u = compute_α_dual(s_u, Δs_u)
    return α_pri, min(α_dual_l, α_dual_u)
end

@inline function compute_μ(x_m_lvar, uvar_m_x, s_l, s_u, nb_low, nb_upp)
    return (dot(s_l, x_m_lvar) + dot(s_u, uvar_m_x)) / (nb_low + nb_upp)
end

function update_pt!(x, y, s_l, s_u, α_pri, α_dual, Δxy, Δs_l, Δs_u, ncon, nvar)
    x .= @views x .+ α_pri .* Δxy[1:nvar]
    y .= @views y .+ α_dual .* Δxy[nvar+1: ncon+nvar]
    s_l .= s_l .+ α_dual .* Δs_l
    s_u .= s_u .+ α_dual .* Δs_u
end

# "security" if x is too close from lvar or uvar
function boundary_safety!(x_m_lvar, uvar_m_x, nlow, nupp, T)
    if 0 in x_m_lvar
        @inbounds @simd for i=1:nlow
            if x_m_lvar[i] == 0
                x_m_lvar[i] = eps(T)^2
            end
        end
    end
    if 0 in uvar_m_x
        @inbounds @simd for i=1:nupp
            if uvar_m_x[i] == 0
                uvar_m_x[i] = eps(T)^2
            end
        end
    end
end

function update_IterData!(itd, pt, fd, id, safety)

    T = eltype(itd.x_m_lvar)
    itd.x_m_lvar .= @views pt.x[id.ilow] .- fd.lvar[id.ilow]
    itd.uvar_m_x .= @views fd.uvar[id.iupp] .- pt.x[id.iupp]
    safety && boundary_safety!(itd.x_m_lvar, itd.uvar_m_x, id.nlow, id.nupp, T)

    itd.μ = compute_μ(itd.x_m_lvar, itd.uvar_m_x, pt.s_l, pt.s_u, id.nlow, id.nupp)
    itd.Qx = mul!(itd.Qx, Symmetric(fd.Q, :U), pt.x)
    itd.xTQx_2 =  dot(pt.x, itd.Qx) / 2
    itd.ATy = mul!(itd.ATy, fd.AT, pt.y)
    itd.Ax = mul!(itd.Ax, fd.AT', pt.x)
    itd.cTx = dot(fd.c, pt.x)
    itd.pri_obj = itd.xTQx_2 + itd.cTx + fd.c0
    itd.dual_obj = @views dot(fd.b, pt.y) - itd.xTQx_2 + dot(pt.s_l, fd.lvar[id.ilow]) -
                        dot(pt.s_u, fd.uvar[id.iupp]) + fd.c0
    itd.pdd = abs(itd.pri_obj - itd.dual_obj ) / (one(T) + abs(itd.pri_obj))
end

function update_data!(pt :: Point{T}, α_pri :: T, α_dual :: T, itd :: IterData{T}, pad :: PreallocatedData{T},
                      res :: Residuals{T}, fd :: QM_FloatData{T}, id :: QM_IntData) where {T<:Real}

    # (x, y, s_l, s_u) += α * Δ
    update_pt!(pt.x, pt.y, pt.s_l, pt.s_u, α_pri, α_dual, itd.Δxy, itd.Δs_l, itd.Δs_u, id.ncon, id.nvar)
    update_IterData!(itd, pt, fd, id, true)

    #update Residuals
    res.n_Δx = @views α_pri * norm(itd.Δxy[1:id.nvar])
    res.rb .= itd.Ax .- fd.b
    res.rc .= itd.ATy .- itd.Qx .- fd.c
    res.rc[id.ilow] .+= pt.s_l
    res.rc[id.iupp] .-= pt.s_u
    # update stopping criterion values:
#         rcNorm, rbNorm = norm(rc), norm(rb)
#         xNorm = norm(x)
#         yNorm = norm(y)
#         optimal = pdd < ϵ_pdd && rbNorm < ϵ_rb * max(1, bNorm + ANorm * xNorm) &&
#                     rcNorm < ϵ_rc * max(1, cNorm + QNorm * xNorm + ANorm * yNorm)
    res.rcNorm, res.rbNorm = norm(res.rc, Inf), norm(res.rb, Inf)
end

function iter!(pt :: Point{T}, itd :: IterData{T}, fd :: Abstract_QM_FloatData{T}, id :: QM_IntData, res :: Residuals{T},
               sc :: StopCrit{Tc}, dda :: DescentDirectionAllocs{T}, pad :: PreallocatedData{T}, ϵ :: Tolerances{T},
               cnts :: Counters, T0 :: DataType, display :: Bool) where {T<:Real, Tc<:Real}

    @inbounds while cnts.k < sc.max_iter && !sc.optimal && !sc.tired # && !small_μ && !small_μ

        out = update_pad!(pad, dda, pt, itd, fd, id, res, cnts, T0) # update data for the solver! function used
        out == 1 && break

        # Solve system to find a direction of descent
        out = update_dd!(dda, pt, itd, fd, id, res, pad, cnts, T0)
        out == 1 && break

        α_pri, α_dual = compute_αs(pt.x, pt.s_l, pt.s_u, fd.lvar, fd.uvar, itd.Δxy, itd.Δs_l, itd.Δs_u, id.nvar)

        if cnts.kc > 0   # centrality corrections
            α_pri, α_dual = multi_centrality_corr!(dda, pad, pt, α_pri, α_dual, itd, fd, id, cnts, res, T0)
            ## TODO replace by centrality_corr.jl, deal with α
        end

        update_data!(pt, α_pri, α_dual, itd, pad, res, fd, id) # update point, residuals, objectives...

        sc.optimal = itd.pdd < ϵ.pdd && res.rbNorm < ϵ.tol_rb && res.rcNorm < ϵ.tol_rc
        sc.small_Δx, sc.small_μ = res.n_Δx < ϵ.Δx, itd.μ < ϵ.μ

        cnts.k += 1
        if T == Float32
            cnts.km += 1
        elseif T == Float64
            cnts.km += 4
        else
            cnts.km += 16
        end

        sc.Δt = time() - sc.start_time
        sc.tired = sc.Δt > sc.max_time

        if display == true
            @info log_row(Any[cnts.k, itd.pri_obj, itd.pdd, res.rbNorm, res.rcNorm, res.n_Δx, α_pri, α_dual, itd.μ])
        end
    end
end


text:
struct niw_hyperparams <: distribution_hyper_params
    κ::Float64
    m::AbstractArray{Float64}
    ν::Float64
    ψ::AbstractArray{Float64}
end

mutable struct niw_sufficient_statistics <: sufficient_statistics
    N::Float64
    points_sum::AbstractArray{Float64,1}
    S::AbstractArray{Float64,2}
end

function calc_posterior(prior:: niw_hyperparams, suff_statistics::niw_sufficient_statistics)
    if suff_statistics.N == 0
        return prior
    end
    κ = prior.κ + suff_statistics.N
    ν = prior.ν + suff_statistics.N
    m = (prior.m.*prior.κ + suff_statistics.points_sum) / κ
    ψ = (prior.ν * prior.ψ + prior.κ*prior.m*prior.m' -κ*m*m'+ suff_statistics.S) / ν
    ψ = Matrix(Hermitian(ψ))
    if  isposdef(ψ) == false
        println(ψ)
        println(m)
        println(κ*m*m')
        println(suff_statistics)
        println(prior)
    end
    return niw_hyperparams(κ,m,ν,ψ)
end

#function calc_posterior(prior:: niw_hyperparams, suff_statistics::niw_sufficient_statistics)
#    κ = prior.κ + suff_statistics.N
#    m = (κ * prior.m + suff_statistics.points_sum) / κ
#    ν = prior.ν + suff_statistics.N
#    ψ = prior.ψ + suff_statistics.S + prior.κ*prior.m*prior.m'-κ*m*m'
#    return niw_hyperparams(κ,m,ν,ψ)
#end

function sample_distribution(hyperparams::niw_hyperparams)
    Σ = rand(Distributions.InverseWishart(hyperparams.ν, hyperparams.ν* hyperparams.ψ))
    μ = rand(Distributions.MvNormal(hyperparams.m, Σ/hyperparams.κ))
    return mv_gaussian(μ,Σ,inv(Σ),logdet(Σ))
end

function create_sufficient_statistics(hyper::niw_hyperparams,posterior::niw_hyperparams,points::AbstractArray{Float64,2}, pts_to_group = 0)
    if size(points,2) == 0
        return niw_sufficient_statistics(size(points,2),zeros(length(hyper.m)),zeros(length(hyper.m),length(hyper.m)))
    end
    pts = Array(points)
    points_sum = @view sum(pts, dims = 2)[:]
    S = pts * pts'
    #println(size(points))
    #println(points_sum)
    #println(S)
    return niw_sufficient_statistics(size(points,2),points_sum,S)
end

function log_marginal_likelihood(hyper::niw_hyperparams, posterior_hyper::niw_hyperparams, suff_stats::niw_sufficient_statistics)
    D = length(suff_stats.points_sum)
    logpi = log(pi)
    return -suff_stats.N*D/2*logpi +
        log_multivariate_gamma(posterior_hyper.ν/2, D)-
        log_multivariate_gamma(hyper.ν/2, D) +
         (hyper.ν/2)*logdet(hyper.ψ*hyper.ν)-
         (posterior_hyper.ν/2)*logdet(posterior_hyper.ψ*posterior_hyper.ν) +
         (D/2)*(log(hyper.κ))-(D/2)*log(posterior_hyper.κ)
end

function aggregate_suff_stats(suff_l::niw_sufficient_statistics, suff_r::niw_sufficient_statistics)
    return niw_sufficient_statistics(suff_l.N+suff_r.N, suff_l.points_sum + suff_r.points_sum, suff_l.S+suff_r.S)
end


text:
# 3x3 matrix A
a = [5 -3 11; 5 6 22; 8 19 0]
println("Matrix A = ", a)

# 3x3 matrix B
b = [12 3 -4; 5 16 9; 8 -3 4]
println("Matrix B = ", b)

# Task 1: addition and subtraction
println("\n============ Task 1: matrix addition and subtraction ============\n")
# expected results
c_add_exp = [17 0 7; 10 22 31; 16 16 4]
c_sub_exp = [-7 -6 15; 0 -10 13; 0 22 -4]

c_add = a + b
print("A+B = ", c_add)
c_sub = a - b
print("A-B = ", c_sub)

println("Does addition work as expected? ", c_add == c_add_exp)
println("Does subtracttion work as expected? ", c_sub == c_sub_exp)

# Task 2: matrix multiplications
println("\n============ Task 2: matrix multiplication ============\n")
c_mult_ab = a * b
c_mult_ba = b * a

println("A * B = ", c_mult_ab)
println("B * A = ", c_mult_ba)
println("==> A*B != B*A because rows from the first matrix will be multiplied with the columns from the seconds matrix. If those rows and columns are in different order, different values will result.")

# Task 3: matrix division
println("\n============ Task 3: matrix division ============\n")
c_diff_1 = a / b
println("A/B = ", c_diff_1)
c_diff_2 = a \ b
println("A\\B = ", c_diff_2)
println("==> A/B != A\\B, because the '\\' operator divides inverse which means that A\\B == B/A")

# Task 4: Single matrix operations
println("\n============ Task 4: single matrix operations ============\n")
a_incr = a + 1
println("A+1 = ", a_incr)
a_decr = a - 1
println("A-1 = ", a_decr)
a_double = a * 2
println("A*2 = ", a_double)
a_half = a / 2
println("A/2 = ", a_half)

# Task 5: matrix and vector multiplication
println("\n============ Task 5: matrix and vector multiplication ============\n")
# 3x4 matrix D
d = [5 -3 11 5; 6 22 8 19; 1 5 8 11]
println("Matrix D = ", d)
# 4x vector V
v = [2, 6, 9, 12]
println("Vector V = ", v)
result = d * v
println("D*V = ", result)
println("==> The result from the multiplication of a 3x4 matrix with a 4-dimensional vector is a 3-dimensional vector!")

max_stars_repo_name: DominicSchiller/osmi-module-datascience

max_issues_repo_name: DominicSchiller/osmi-module-datascience

max_forks_repo_name: DominicSchiller/osmi-module-datascience


text:
module EnsembleKalmanProcesses

using Distributions, Statistics, LinearAlgebra, DocStringExtensions

# auxiliary modules
include("ParameterDistributions.jl")
include("DataContainers.jl")
include("Observations.jl")
include("Localizers.jl")

# algorithmic updates
include("EnsembleKalmanProcess.jl")

end # module


text:
#
# function expand_star_syntax(s::String)
#     s = strip(s)
#     if !occursin("*", s)
#         return s
#     end
#     # Handle e.g `pawecutdg*`
#     if typeof(Meta.parse(string(s[1]))) == Symbol && s[end - 1] == '*'
#         return s
#     end
#     flag, rest = split(s)
#     spl = strip_split(rest, "*")
#
#     # Handle "*2" case i.e. return "*2"
#     if length(spl) == 1
#         @assert s[1] != '*' "error in expanding start syntax"
#         return string(s[1])
#     end
#
#     val = spl[2]
#     out = "$flag $val"
#     len = Meta.parse(string(spl[1]))
#     for i = 1:len-1
#         out *= " $val"
#     end
#     return out
# end
#
# @inline function eval_abinit_operator(s)
#     spl = split(s)
#     out = "$(spl[1])"
#     for p in spl[2:end]
#         out *= " $(eval(Meta.parse(string(p))))"
#     end
#     return out
# end
#
# "series like 10 ecut+ 5 are not supported"
# function read_abi_datasets(filename::String, T=Float64)
#     lines = readlines(filename) .|>
#         strip                   .|>
#         lowercase               .|>
#         filter_comment
#     lines = filter(x -> length(x) > 1, lines) .|>
#         expand_star_syntax                    .|>
#         eval_abinit_operator
#
#     lines = split(join(lines," "))
#
#     datasets = Vector{Dict{Symbol, Any}}()
#     dataset = Dict{Symbol,Any}()
#     flag = Symbol("start")
#     for (l, line) in enumerate(lines)
#         if occursin("?", line) || occursin(":", line)
#             continue
#         end
#         #convert to the correct abinit units
#         if line in ABI_UNIT_NAMES
#             dataset[flag] = convert_2abi(dataset[flag], line)
#             continue
#         end
#         if typeof(Meta.parse(line)) == Symbol
#             j = tryparse(Int, string(line[end]))
#             if !isnull(j)
#                 dtset = get(j)+1
#                 while dtset > length(datasets)
#                     push!(datasets, Dict{Symbol, Any}())
#                 end
#
#                 dataset = datasets[dtset]
#                 flag = Symbol(line[1:end-1])
#             else
#                 dtset = 1
#                 while dtset > length(datasets)
#                     push!(datasets, Dict{Symbol, Any}())
#                 end
#                 dataset = datasets[dtset]
#                 flag = Symbol(line)
#             end
#         else
#             if flag == :ndtset
#                 continue
#             end
#             if haskey(dataset,flag)
#                 push!(dataset[flag], Meta.parse(line))
#             else
#                 dataset[flag] = [Meta.parse(line)]
#             end
#         end
#     end
#     #make flags that are length one just be their value and that they are in the abivars
#     for (d,data) in enumerate(datasets)
#         for (flag, val) in data
#             flag_type = abi_flag_type(flag)
#             if flag_type == Nothing
#                 error("Flag $flag in dataset $d not found in abinit variables.")
#             elseif flag_type != typeof(val)
#                 try
#                     val = convert.(flag_type, val)
#                 catch
#                     error("Flag type of flag $flag not correct, found: $(typeof(val)) expected:$flag_type.")
#                 end
#             elseif flag == :spgroup || flag == :nobj
#                 error("Abinit spgroup builder is not supported. Structure must be given explicitly!")
#             end
#             if length(val) == 1
#                 data[flag] = val[1]
#             end
#         end
#     end
#
#     return datasets
# end
#
# #Everything is in angstrom, bohr lengths begone foul beasts!
# function extract_structures!(abi_datasets...; structure_name = "NoName")
#     #possible different structures for different calculations
#     out_structures = Vector{Vector{AbinitDataBlock}}()
#     out_structures = Vector{Structure}()
#     natom = 1
#     znucl = 1
#     ntypat =1
#     typat = [1]
#     for data in abi_datasets
#
#         acell = pop!(data, :acell, [1.0, 1.0, 1.0])
#         rprimd = zeros(3, 3)
#         if haskey(data, :angdeg)
#             if haskey(data, :rprim)
#                 error("rprim and angdeg cannot be used at the same time.")
#             end
#             angdeg = pop!(data, :angdeg)
#             if !all(x -> x > 0, angdeg)
#                 error("Angles must be > 0 but got $angdeg")
#             elseif sum(angdeg)
#                 error("The sum of angdeg must be lower that 360, angdeg $angdeg")
#             end
#             tol12 = 1e-12
#             rprim = zeros(3,3)
#             if abs(angdeg[1] -angdeg[2]) < tol12 && abs(angdeg[2] - angdeg[3]) < tol12 &&
#             abs(angdeg[1]-90.) + abs(angdeg[2]-90.) + abs(angdeg[3] -90.) > tol12
#
#             # Treat the case of equal angles (except all right angles):
#             # generates trigonal symmetry wrt third axis
#                 cosang = cos(pi * angdeg[1]/180.0)
#                 a2 = 2.0 / 3.0 * (1.0 - cosang)
#                 aa = sqrt(a2)
#                 cc = sqrt(1.0-a2)
#                 rprim[1,1] = aa
#                 rprim[1,2] = 0.0
#                 rprim[1,3] = cc
#                 rprim[2,1] = -0.5 * aa
#                 rprim[2,2] = sqrt(3.0) * 0.5 * aa
#                 rprim[2,3] = cc
#                 rprim[3,1] = -0.5 * aa
#                 rprim[3,2] = -sqrt(3.0) * 0.5 * aa
#                 rprim[3,3] = cc
#             else
#                 # Treat all the other cases
#                 rprim[1,1] = 1.0
#                 rprim[2,1] = cos(pi * angdeg[3] / 180.)
#                 rprim[2,2] = sin(pi * angdeg[3] / 180.)
#                 rprim[3,1] = cos(pi * angdeg[2] / 180.)
#                 rprim[3,2] = (cos(pi * angdeg[1] / 180.0) - rprim[2,1] * rprim[3,1]) / rprim[2,2]
#                 rprim[3,3] = sqrt(1.0 - rprim[3,1]^2 - rprim[3,2]^2)
#             end
#             for i = 1:3
#                 rprimd[i,:] *= acell[i]
#             end
#         else
#             rprimd = reshape(pop!(data, :rprim, eye(3)), (3,3))
#             for i=1:3
#                 rprimd[i,:] *= acell[i]
#             end
#         end
#         rprimd *= conversions[:bohr2ang]
#         rprimd = Mat3(rprimd)
#         ntypat = pop!(data, :ntypat, ntypat)
#         natom = pop!(data, :natom, natom)
#         znucl = pop!(data, :znucl, znucl)
#         #some kind of npsp not supported
#         typat = pop!(data, :typat, typat)
#         typat = typat[1:natom]
#
#         s_atoms = Atom{eltype(rprimd)}[]
#         if haskey(data, :xred)
#             xred = pop!(data, :xred)
#             atoms = reshape(xred, (div(length(xred), 3), 3))
#
#             for i=1:length(typat)
#                 typ = typat[i]
#                 atom = Array(atoms[i,:])
#                 z = znucl[typ]
#                 nz = length(filter( x-> x==z, znucl))
#                 _element = element(z)
#                 z_sym = nz == 1 ? _element.symbol : Symbol(String(_element.symbol) * "$typ")
#                 push!(s_atoms, Atom(z_sym, _element, rprimd' * Point3(atom)))
#             end
#
#
#         elseif haskey(data, :xcart)
#             xcart = pop!(data, :xcart)
#             atoms = reshape(xcart, (div(length(xcart), 3), 3))
#             for (typ, atom) in zip(typat, rows(atoms))
#                 z = znucl[typ]
#                 nz = length(filter( x-> x==z, znucl))
#                 _element = element(z)
#                 z_sym = nz == 1 ? _element.symbol : Symbol(String(_element.symbol) * "$typ")
#                 push!(s_atoms, Atom(z_sym, _element, conversions[:bohr2ang] * Point3(atom)))
#             end
#
#         elseif haskey(data, :xangst)
#             xangst = pop!(data, :xangst)
#             atoms = reshape(xangst, (div(length(xangst), 3), 3))
#             for (typ, atom) in zip(typat, rows(atoms))
#                 z = znucl[typ]
#                 nz = length(filter( x-> x==z, znucl))
#                 element = element(z)
#                 z_sym = nz == 1 ? _element.symbol : Symbol(String(_element.symbol) * "$typ")
#
#                 push!(s_atoms, Atom(z_sym, _element, Point3(atom)))
#             end
#
#         else
#             try
#                 rprimd = out_structures[1].cell
#                 atoms = out_structures[1].atoms
#             catch
#                 error("xred|xcart|xangs must be given in calculation.")
#             end
#         end
#         # up till here everything was in bohr, not anymore.
#         if !isempty(s_atoms)
#             push!(out_structures, Structure(rprimd, s_atoms))
#         end
#     end
#     return out_structures
# end
#
# """
#     read_abi_calculation(filename::String, T=Float64)
#
# Returns an ABINIT calculation. We assume that jdtset is on a seperate line.
# If # DATASET is supplied as first line, it will make sure that amount of datasets are read no matter what ndtset and jdtset are.
# `ndtset` and `jdtset` will be taken into account to decide which calculations will be marked as 'should run'.
# Either all structures of the calculation are the same or all are different, otherwise there will be an error.
# """
# function read_abi_calculation(filename::String, T=Float64; runcommand= "", pseudos=[""], structure_name = "NoName")
#     datasets = read_abi_datasets(filename, T)
#     structures = extract_structures!(datasets..., structure_name = structure_name)
#     calculations = Array{AbinitInput,1}()
#     jdtset = Int[]
#     ndtset = 0
#     @assert length(datasets) == length(structures) || length(structures) == 1 "All structures of the calculation have to be the same or all different."
#     for (i, data) in enumerate(datasets)
#         structure = length(structures) == 1 ? structures[1] : structures[i]
#         file, ext= splitext(splitdir(filename)[2])
#         if length(datasets) > 1
#             newfile = file * "$i" * ext
#         else
#             newfile = file * ext
#         end
#         if haskey(data, :jdtset)
#             jdtset = pop!(data, :jdtset)
#         end
#         run = i in jdtset
#         if isempty(data)
#             continue
#         end
#         push!(calculations, AbinitInput(newfile, structure, data, [AbinitDataBlock(:pseudos, :pseudos, pseudos)], runcommand, run))
#     end
#     return calculations
# end
#
# function write_abi_structure(f, structure)
#     write(f, "acell 1.0 1.0 1.0\n")
#     write(f, "xangst\n")
#     for at in structure.atoms
#         write(f, at.position_cart)
#         write(f, "\n")
#     end
#     write(f, "rprim\n")
#     write_cell(f, structure.cell)
#
#     unique = unique(structure.atoms)
#     write(f, "nat $(length(structure.atoms))\n")
#     write(f, "ntypat $(length(unique))\n")
#     write(f, "typat\n")
#     i = 1
#     for at in structure.atoms
#         write(f, "$(findfirst(x-> x.name== at.name, unique)) ")
#         if i == 3
#             i=0
#             write(f,"\n")
#         end
#         i+=1
#     end
#     write(f, "\n")
#     write(f, "znucl ")
#     for at in unique
#         write(f, "$(at.element.Z) ")
#     end
#     write(f, "\n")
# end
#
# #question: Does it matter that we might be writing redundant data such as spinat etc?
# """
#     write_abi_datasets(calculations::Array{AbinitInput,1}, directory)
#
# Takes all the calculations, sees which ones have the same structure and constructs calculation files for each seperate sturcture.
# The filename of the first calculation of a certain structure is used as file for all the datasets.
# """
# function write_abi_datasets(calculations::Vector{Calculation{Abinit}}, directory)
#     calculation_groups = Vector{Vector{Calculation{Abinit}}}([[calculations[end]]])
#     for calculation in reverse(calculations)[2:end]
#         for group in calculation_groups
#             if group[1].structure == calculation.structure
#                 push!(group, calculation)
#                 break
#             end
#         end
#     end
#
#     filenames    = String[]
#     run_commands = String[]
#     pseudos      = Vector{Vector{String}}()
#
#     for group in calculation_groups
#         run_indices = Int[]
#         for (i, _calculation) in enumerate(reverse(group))
#             if _calculation.run
#                 push!(run_indices, i)
#             end
#         end
#         if length(run_indices) == length(group)
#             jdtset = ""
#             ndtset = length(group)
#         else
#             jdtset = join(["$c" for c in run_indices], " ")
#             ndtset = length(run_indices)
#         end
#         file, ext = splitext(group[end].filename)
#         push!(filenames, file[1:end-1] * ext)
#         push!(run_commands, group[end].runcommand)
#         push!(pseudos, data(group[end], :pseudos))
#         open(directory * file[1:end-1] * ext, "w") do f
#             write(f, "# DATASETS $(length(group))\n")
#             write(f, "ndtset $ndtset\n")
#             if jdtset != ""
#                 write(f, "jdtset $jdtset\n")
#             end
#             write_abi_structure(f, group[1].structure)
#             write(f, "\n")
#             for (t, dt) in enumerate(reverse(group))
#                 write_flag(flag_data) = write_flag_line(f, flag_data[1], flag_data[2], "", t)
#                 write(f, "#=============== BEGIN DATASET $t ===============#\n")
#                 write_flag.(collect(dt.flags))
#                 write(f, "#================ END DATASET $t ================#\n")
#             end
#         end
#     end
#     return zip(filenames, pseudos, run_commands)
# end
#
# #very stupid
function read_abi_output(filename::String, T = Float64)
    if occursin("FATBANDS", filename)
        return read_abi_fatbands(filename, T)
    elseif occursin("EBANDS.agr", filename)
        return read_abi_ebands(filename, T)
    elseif occursin("_EIG", filename)
        return read_abi_eig(filename, T)
    else
        error("Please supply a file with FATBANDS, EBANDS.agr or _EIG in the filename.")
    end
end
#
# "Reads an abinit FATBANDS output file and returns the found Bands, with pdos values in the data field. K-points are not given (they aren't present in the output file)."
# function read_abi_fatbands(filename::String, T=Float64)
#     bands = Band[]
#     open(filename, "r") do f
#         while !eof(f)
#             line = readline(f)
#             if occursin("BAND number", line)
#                 extra   = Dict{Symbol,Any}(:pdos => T[])
#                 eigvals = Vector{T}()
#                 line    = readline(f)
#                 while line != "" && line != "&"
#                     eigval, pdos = Meta.parse.(T, split(line)[2:end])
#                     push!(eigvals, eigval)
#                     push!(extra[:pdos], pdos)
#                     line = readline(f)
#                 end
#                 push!(bands, Band{T}([T[0.0, 0.0, 0.0] for i=1:length(eigvals)], [T[0.0, 0.0, 0.0] for i=1:length(eigvals)], eigvals, extra))
#             end
#         end
#     end
#     return bands
# end
#
# "Reads an abinit EBANDS.agr output file and returns the found Bands. K-points only given in crystalline coordinates."
function read_abi_ebands(filename::String, T = Float64)
    bands = Band[]
    open(filename, "r") do f
        k_points_cryst = Vector{Vector{T}}()
        while !eof(f)
            line = readline(f)
            if occursin("List of k-points", line)
                line = readline(f)
                while line[1] != '@'
                    k_point = parse.(T,
                                     replace.(strip.(strip.(strip.(split(line)[4:end], '['),
                                                            ']'), ','), 'E', 'e'))#ohno the replace here!
                    push!(k_points_cryst, k_point)
                    line = readline(f)
                end

            elseif occursin("target", line)
                readline(f)
                eigvals = T[]
                line = readline(f)
                while line[1] != '&'
                    push!(eigvals, parse(T, split(line)[2]))
                    line = readline(f)
                end
                push!(bands,
                      Band([T[0.0, 0.0, 0.0] for i in 1:length(eigvals)], k_points_cryst,
                           eigvals))
            end
        end
    end
    return bands
end
#
# "Reads and abinit _EIG output file and returns the found Bands. K-points are only given in crystalline coordinates."
# function read_abi_eig(filename::String, T=Float64)
#     bands = Band[]
#     k_points_array = Vector{Vector{T}}()
#     eigvals_array  = Vector{Vector{T}}()
#     open(filename, "r") do f
#         while !eof(f)
#             line = readline(f)
#             if occursin("kpt#", line)
#                 push!(k_points_array, Meta.parse.(T, split(line)[7:9]))
#                 read_abi_eig_block!(f, k_points_array, eigvals_array, T)
#             end
#         end
#     end
#     for eigvals in eigvals_array
#         push!(bands, Band([T[0.0, 0.0, 0.0] for i=1:length(eigvals)], k_points_array, eigvals))
#     end
#     return bands
# end
#
# function read_abi_eig_block!(f, k_points_array, eigvals_array, T=Float64)
#     first_run = isempty(eigvals_array)
#     if eof(f)
#         return
#     end
#     i = 1
#     while !eof(f)
#         line = readline(f)
#         if occursin("Eigenvalues", line)
#             continue
#         elseif occursin("kpt#", line)
#             push!(k_points_array, Meta.parse.(T, split(line)[7:9]))
#             read_abi_eig_block!(f, k_points_array, eigvals_array, T)
#         else
#             for val in Meta.parse.(T, split(line))
#                 if first_run
#                     push!(eigvals_array, [val])
#                 else
#                     push!(eigvals_array[i], val)
#                 end
#                 i += 1
#             end
#         end
#     end
# end


text:
# Example of GOKU-net model on friction-less pendulum data created with Luxor

using LatentDiffEq
using FileIO
using Parameters: @with_kw
using ProgressMeter: Progress, next!
using Random
using Statistics
using MLDataUtils
using BSON: @save
using Flux.Data: DataLoader
using Flux
using OrdinaryDiffEq
using StochasticDiffEq
using ModelingToolkit
using DiffEqSensitivity
using Images
using Plots
import GR
using CUDA
CUDA.allowscalar(false)

include("pendulum.jl")
include("create_data.jl")

################################################################################
## Arguments for the train function
@with_kw mutable struct Args
    ## Global model
    model_type = GOKU()

    ## Latent Differential Equations
    diffeq = Pendulum()
    # diffeq = SPendulum()

    ## Training params
    η = 1e-3                        # learning rate
    decay = 0.001f0                 # decay applied to weights during optimisation
    batch_size = 64                 # minibatch size
    seq_len = 50                    # sequence length for training samples
    epochs = 1500                   # number of epochs for training
    seed = 333                      # random seed
    cuda = true                     # GPU usage (not working well yet)
    dt = 0.05                       # timestep for ode solve
    variational = true              # variational or deterministic training

    ## Annealing schedule
    start_β = 0f0                   # start value
    end_β = 1f0                     # end value
    n_cycle = 4                     # number of annealing cycles
    ratio = 0.9                     # proportion used to increase β (and 1-ratio used to fix β)

    ## Progressive observation training
    progressive_training = false    # progressive training usage
    prog_training_duration = 200    # number of epochs to reach the final seq_len
    start_seq_len = 10              # training sequence length at first step

    ## Visualization
    vis_len = 60                    # number of test frames to visualize after each epoch
    save_figure = true              # true: save visualization figure in save_path folder
                                    # false: display image instead of saving it
end

################################################################################
################################################################################
## Training done manualy

function train(; kws...)
    ## Load hyperparameters and GPU config
    args = Args(; kws...)
    @unpack_Args args

    if cuda && has_cuda_gpu()
        device = gpu
        @info "Training on GPU"
    else
        device = cpu
        @info "Training on CPU"
    end

    ############################################################################
    ## Prepare training data
    root_dir = @__DIR__
    data_path = "$root_dir/data/data.bson"

    if ~isfile(data_path)
        @info "Generating data"
        latent_data, u0s, ps, high_dim_data = generate_dataset(diffeq = diffeq)
        data = (latent_data, u0s, ps, high_dim_data)
        mkpath("$root_dir/data")
        @save data_path data
    end

    seed > 0 && Random.seed!(seed)

    data_loaded = load(data_path, :data)
    train_data = data_loaded[4]
    latent_data = data_loaded[1]
    params_data = data_loaded[3]

    # Stack time for each sample
    train_data = Flux.stack.(train_data, 3)

    # Stack all samples
    train_data = Flux.stack(train_data, 4) # 28x28x400x450
    h, w, full_seq_len, observations = size(train_data)
    latent_data = Flux.stack(latent_data, 3)
    params_data = Flux.stack(params_data, 3)

    # Vectorize frames
    train_data = reshape(train_data, :, full_seq_len, observations) # input_dim, time_size, samples
    train_data = Float32.(train_data)

    # Split into train and validation sets
    train_set, val_set = Array.(splitobs(train_data, 0.9))
    train_set_latent, val_set_latent = Array.(splitobs(latent_data, 0.9))
    train_set_params, val_set_params = Array.(splitobs(params_data, 0.9))

    # Prepare data loader
    loader_train = DataLoader(train_set, batchsize=batch_size, shuffle=true, partial=false)

    val_set = permutedims(val_set, [1,3,2])
    t_val = range(0.f0, step=dt, length = size(val_set, 3))
    input_dim = size(train_set, 1)

    ############################################################################
    # Create model
    encoder_layers, decoder_layers = default_layers(model_type, input_dim, diffeq, device = device)
    model = LatentDiffEqModel(model_type, encoder_layers, decoder_layers)

    # Get parameters
    ps = Flux.params(model)

    ############################################################################
    ## Define optimizer
    #opt = ADAM(η)
    # opt = AdaBelief(η)
    opt = ADAMW(η, (0.9,0.999), decay)

    ############################################################################
    ## Various definitions
    if progressive_training
        prog_seq_lengths = range(start_seq_len, seq_len, step=(seq_len-start_seq_len)/(prog_training_duration-1))
        prog_seq_lengths = Int.(round.(prog_seq_lengths))
    else
        prog_training_duration = 0
    end

    # KL annealing scheduling
    annealing_schedule = frange_cycle_linear(epochs, start_β, end_β, n_cycle, ratio)

    # Preparation for saving best models weights
    mkpath("$root_dir/output")
    best_val_loss = Inf32
    val_loss = 0f0

    # mkpath("$root_dir/output")
    # args = struct2dict(args)
    # @save "$root_dir/output/args.bson" args

    ## Visualization options
    if save_figure
        vis_dir = "$root_dir/output/visualization"
        mkpath(vis_dir)
        @info "Visualizations at $vis_dir"
        GR.inline("pdf")
    end

    ############################################################################
    ## Main train loop
    @info "Start Training of $(typeof(model_type))-net, total $epochs epochs"
    for epoch = 1:epochs

        # Set annealing factor
        β = annealing_schedule[epoch]

        # Set a sequence length for training samples
        seq_len = epoch ≤ prog_training_duration ? prog_seq_lengths[epoch] : seq_len

        # Model evaluation length
        t = range(0.f0, step=dt, length=seq_len)

        @info "Epoch $epoch .. (Sequence training length $seq_len)"
        progress = Progress(length(loader_train))

        for x in loader_train

            # Permute dimesions for having (pixels, batch_size, time)
            x = PermutedDimsArray(x, [1,3,2])

            # Use only random sequences of length seq_len for the current minibatch
            x = time_loader(x, full_seq_len, seq_len)

            # Run the model with the current parameters and compute the loss
            loss, back = Flux.pullback(ps) do
                loss_batch(model, x |> device, t, β, variational)
            end

            # Backpropagate and update
            grad = back(1f0)
            Flux.Optimise.update!(opt, ps, grad)

            # Use validation set to get loss and visualisation
            val_loss = loss_batch(model, val_set |> device, t_val, β, false)

            # Progress meter
            next!(progress; showvalues=[(:loss, loss), (:val_loss, val_loss)])
        end

        visualize_val_image(model, val_set |> device, val_set_latent, val_set_params, vis_len, dt, h, w, device, save_figure, epoch)

        if val_loss < best_val_loss
            best_val_loss = deepcopy(val_loss)
            weights = Flux.params(model)
            @save "$root_dir/output/best_model_weights.bson" weights
            @info "Model saved"
        end
    end
end

################################################################################
## Loss definition

function loss_batch(model, x, t, β, variational)

    # Make prediction
    X̂, μ, logσ² = model(x, t, variational)
    x̂, ẑ, l̂ = X̂

    # Compute reconstruction loss
    reconstruction_loss = sum(mean((x .- x̂).^2, dims=(2,3)))

    # Compute KL losses for parameters and initial values
    kl_loss = vector_kl(μ, logσ²)

    return reconstruction_loss + β * kl_loss
end

################################################################################
## Visualization function

function visualize_val_image(model, val_set, val_set_latent, val_set_params, vis_len, dt, h, w, device, save_figure, epoch)

    # randomly pick a sample from val_set and a random time interval of length vis_len
    j = rand(1:size(val_set, 2))
    idxs = rand_time(size(val_set, 3), vis_len)
    x = val_set[:, j:j, idxs] |> device
    true_latent = val_set_latent[:, idxs, j]
    true_params = val_set_params[j]

    # create time range for the model diffeq solving
    t_val = range(0.f0, step=dt, length=vis_len)

    # run model with current parameters on the picked sample
    X̂, μ, logσ² = model(x, t_val)
    x̂, ẑ, l̂ = X̂
    ẑ₀, θ̂ = l̂

    θ̂ = cpu(θ̂)
    ẑ = cpu(ẑ)
    x̂ = cpu(x̂)
    x = cpu(x)
    θ̂ = θ̂[1]

    # plot actual and inferred angles
    plt1 = plot(ẑ[1,1,:], legend=false, ylabel="inferred angle", box = :on, color=:indigo, yforeground_color_axis=:indigo, yforeground_color_text=:indigo, yguidefontcolor=:indigo, rightmargin = 2.0Plots.cm)
    xlabel!("time")
    plt1 = plot!(twinx(), true_latent[1,:], color=:darkorange1, box = :on, xticks=:none, legend=false, ylabel="true angle", yforeground_color_axis=:darkorange1, yforeground_color_text=:darkorange1, yguidefontcolor=:darkorange1)
    title!("Sample from validation set")

    # downsample
    x = @view x[:, :, 1:6:end]
    x̂ = @view x̂[:, :, 1:6:end]

    # build frames vectors
    to_image(x) = Gray{N0f8}.(reshape(x, h, w))
    frames_val = [to_image(xₜ) for xₜ in eachslice(x, dims = 3)]
    frames_pred = [to_image(x̂ₜ) for x̂ₜ in eachslice(x̂, dims = 3)]

    # plot a mosaic view of the frames
    plt2 = mosaicview(frames_val..., frames_pred..., nrow=2, rowmajor=true)
    plt2 = plot(plt2, leg = false, ticks = nothing, border = :none)
    annotate!((208, -21, ("True Pendulum Length = $(round(true_params, digits = 2))", 9, :gray, :right)))
    annotate!((208, -11, ("Inferred Pendulum Length = $(round(θ̂, digits = 2))", 9, :gray, :right)))
    plt = plot(plt1, plt2, layout = @layout([a; b]))
    save_figure ? savefig(plt, "output/visualization/fig_$epoch.pdf") : display(plt)
    return nothing
end

train()

max_issues_repo_name: gabrevaya/LatentDiffEq

max_forks_repo_name: gabrevaya/LatentDiffEq


text:
using Celeste

function main()
    f = FITSIO.FITS(ENV["FIELD_EXTENTS"])

    hdu = f[2]::FITSIO.TableHDU

    # read in the entire table.
    all_run = read(hdu, "run")::Vector{Int16}
    all_camcol = read(hdu, "camcol")::Vector{UInt8}
    all_field = read(hdu, "field")::Vector{Int16}

    close(f)

    for i in eachindex(all_run)
        println("RUN=$(all_run[i]) CAMCOL=$(all_camcol[i]) FIELD=$(all_field[i])")
    end
end

main()


text:
#=
This is the main function that does the google overlay.
Uses the various functions below to:
    1. Work out the Lon | Lat coordinates from the bounding box
    2. Use those coordinates and the zoom to calc the inputs for
       the Googe API - center and size
    3. Calls the Google API with calculated and user provided Inputs
    4. Saves that to a temp picture
    5. Overlays the canvas with the Google image
=#
function google_overlay(canvas::Compose.Context, key::String, zoom::Int;
                        scale::Int=1, maptype::String=roadmap)

    # Get MBR from the canvas units
    MBR = canvas_to_MBR(canvas::Compose.Context)

    # Coordinates from the MBR
    cxEast, cyNorth, cxWest, cySouth = MBR.right, MBR.top, MBR.left, MBR.bottom

    # Use coordinates to get Google API parameters
    mid_Lon, mid_Lat, p_width, p_height = bbox_coords_and_zoom_to_center_and_size(cxEast, cyNorth, cxWest, cySouth, zoom)

    # Build a dictionary of the google options for the API call
    options = Dict(
        :center => "$(mid_Lat),$(mid_Lon)",
        :zoom => "$zoom",
        :size => "$(p_width)x$(p_height)",
        :scale => "$scale",
        :maptype => maptype,
        :key => key)

    # Fetch the Google map image
    google_img = load_google_image(options)

    # Overlay the canvas onto the Google map image
    overlay(canvas, google_img)
end

#=
Takes the Compose.Context canvas bounding box and converts the units back to
longitude and lattitude to use in fetching the right image from Google.
=#
function canvas_to_MBR(canvas::Compose.Context)

    # From UnitBox
    left = canvas.units.x0
    top = canvas.units.y0
    width = canvas.units.width
    height = canvas.units.height
    right = left + width
    bottom = top + height

    # Convert back to longitude and lattitude
    left, top = webmercator_to_lonlat(left, top)
    right, bottom = webmercator_to_lonlat(right, bottom)

    # Create an MBR
    Shapefile.Rect(left, bottom, right, top)
end

function bbox_coords_and_zoom_to_center_and_size(cxEast, cyNorth, cxWest, cySouth, zoom)

    # Convert Lon and Lat to pixel numbers
    pxEast, pyNorth = lonlat_to_pixels(cxEast, cyNorth, zoom)
    pxWest, pySouth = lonlat_to_pixels(cxWest, cySouth, zoom)

    # Calculate image size in pixels
    width = pxEast - pxWest
    height = pySouth - pyNorth

    # Google has max size of 640x640 so give error if choice exceeds this
    (width > 640 || height > 640) &&
    throw("Bounds Error - Google sets 640 as the upper limit on the width and height. Select a lower zoom.")

    # Calculate centre point in pixels
    mid_px = (pxEast + pxWest) / 2
    mid_py = (pyNorth + pySouth) / 2

    # Convert center coordinates to Lon/Lat
    mid_cxLon, mid_cyLat = pixels_to_lonlat(mid_px, mid_py, zoom)

    # Return center x, center y, width and height - these are used in the google API call
    mid_cxLon, mid_cyLat, ceil(Int, width), ceil(Int, height)
end

# Fetch the image from the Google static ... API
function load_google_image(params::Dict)
    img = load(fetch_google_image("temp_google_api_img.png", params))
end
function fetch_google_image(save_path::AbstractString, params::Dict)
    download(google_staticmap_url(params), save_path)
end

# Build the url string to be used in the Google API
# Inputs are the options ...
function google_staticmap_url(params::Dict)
    url = "https://maps.googleapis.com/maps/api/staticmap?"
    for (key, value) in params
        url *= "$key=$value"
        url *= "&"
    end
    return url[1:end-1]
end

#=
Overlays the shape plot canvas onto the Google image.
The transparency in the color used in the shape plot is used.
=#
function overlay(shape_plot, google_img)

    # Image dimensions
    height, width = size(google_img)

    # Create Cairo surface, draw polygon to surface, reinterpret as an array of RGB fixed point numbers
    surface = CairoARGBSurface(zeros(UInt32, height, width))
    draw(PNG(surface), shape_plot)
    overlay = reinterpret(BGRA{Normed{UInt8, 8}}, surface.data')

    # Create combined image
    outimg = similar(google_img)
    for i=1:height, j=1:width

        # Transparency value at pixel i,j
        α = overlay[i,j].alpha

        # Convert pixel value to same element type (RGB of FixedPointNumbers) as the image
        overlay_pixel = convert(eltype(google_img), overlay[i,j])

        # Merge the two pixels - weighting determined by the transparency
        β = one(α) - α
        outimg[i,j] = β*google_img[i,j] + α*overlay_pixel
    end

    # Return the combined image
    return outimg
end


text:
using Test
using Julog
using PDDL, PDDL.Parser

include("strips/test.jl")
include("typing/test.jl")
include("axioms/test.jl")
include("adl/test.jl")
include("fluents/test.jl")
include("functions/test.jl")


text:
############# Distribution types vectorizations   ################
# single parameter distributions

for d in [:Bernoulli, :TDist, :Exponential, :Poisson]
	@eval begin
		function ($d)(ps::Array)
			ds = Array($d, size(ps))
			for i in 1:length(ds)
				ds[i] = ($d)(ps[i])
			end
			ds
		end
	end
end

# two parameter distributions

for d in [:Normal, :Uniform, :Weibull, :Gamma, :Cauchy, :LogNormal, :Binomial, :Beta, :Laplace]
	@eval begin
		function ($d)(p1::Array, p2::Array)
			ds = Array($d, size(p1))
			for i in 1:length(ds)
				ds[i] = ($d)(p1[i], p2[i])
			end
			ds
		end

		function ($d)(p1::Array, p2::Real)
			ds = Array($d, size(p1))
			for i in 1:length(ds)
				ds[i] = ($d)(p1[i], p2)
			end
			ds
		end

		function ($d)(p1::Real, p2::Array)
			ds = Array($d, size(p2))
			for i in 1:length(ds)
				ds[i] = ($d)(p1, p2[i])
			end
			ds
		end
	end
end

############# logpdf vectorization on the distribution argument   ################

function logpdf{T<:Distribution}(ds::Array{T}, x::AbstractArray)
	res = Array(Float64, size(ds))
	size(ds) == size(x) || error("x and distributions sizes do not match")
	for i in 1:length(x)
		res[i] = logpdf(ds[i], x[i])
	end
	res
end

function logcdf{T<:Distribution}(ds::Array{T}, x::AbstractArray)
	res = Array(Float64, size(ds))
	size(ds) == size(x) || error("x and distributions sizes do not match")
	for i in 1:length(x)
		res[i] = logcdf(ds[i], x[i])
	end
	res
end

function logccdf{T<:Distribution}(ds::Array{T}, x::AbstractArray)
	res = Array(Float64, size(ds))
	size(ds) == size(x) || error("x and distributions sizes do not match")
	for i in 1:length(x)
		res[i] = logccdf(ds[i], x[i])
	end
	res
end


text:
function Min()
    return Min(())
end

function Min(arg0::Min)
    return Min((Min,), arg0)
end

function aggregate(obj::Min, arg0::Min)
    return jcall(obj, "aggregate", void, (Min,), arg0)
end

function clear(obj::Min)
    return jcall(obj, "clear", void, ())
end

function copy(obj::Min)
    return jcall(obj, "copy", Min, ())
end

function evaluate(obj::Min, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function get_n(obj::Min)
    return jcall(obj, "getN", jlong, ())
end

function get_result(obj::Min)
    return jcall(obj, "getResult", jdouble, ())
end

function increment(obj::Min, arg0::jdouble)
    return jcall(obj, "increment", void, (jdouble,), arg0)
end


text:
module Gts
export gts_point_new, GtsObject, GtsPoint, gts_point_class, gts_destroy

const libGTS = "libgts"

abstract AbstractGtsObject

immutable GtsObject <: AbstractGtsObject
    gtsObjectClass::Ptr{Void}
    reserved::Ptr{Void}
    flags::UInt32
end

immutable GtsPoint <: AbstractGtsObject
    gtsObject::GtsObject
    x::Float64
    y::Float64
    z::Float64
end

function gts_point_class()
    return ccall((:gts_point_class, libGTS), Ptr{Void}, ())
end

function gts_point_new(x, y, z)
    pointclass = gts_point_class()
    gtspointer = ccall((:gts_point_new, libGTS), Ref{GtsPoint}, (Ptr{Void}, Cfloat, Cfloat, Cfloat), pointclass, x, y, z)
    return gtspointer
end

function gts_destroy{T<:AbstractGtsObject}(ob::Ref{T})
    return ccall((:gts_object_destroy, libGTS), Void, (Ref{T},), ob)
end

end # module