text:
# Copyright (c) 2021 California Institute of Technology (“Caltech”).
# U.S. Government sponsorship acknowledged.
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# - Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# - Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# - Neither the name of Caltech nor its operating division, the Jet Propulsion
#   Laboratory, nor the names of its contributors may be used to endorse or
#   promote products derived from this software without specific prior written
#   permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
# PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
# OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
# OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
# ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

##################### Radial cryowindow ################################

function calc_Rb̃llnn(ℓ, rlo, rhi, keff)
    I, E = quadgk(r->sphericalbesselj(ℓ, r*keff), rlo, rhi)
    return I * keff * √(2/π) / (rhi - rlo)
end

function calc_Rb̃llnn(ℓ, rlo, rhi, keff, jl_spl)
    xlo = rlo * keff
    xhi = rhi * keff
    I = √(2/π) / (rhi - rlo) * integrate(jl_spl, xlo, xhi)
    if !isfinite(I) # || !(I ≈ I0)
        I0 = calc_Rb̃llnn(ℓ, rlo, rhi, keff)
        @error "Rb̃llnn error" ℓ rlo rhi keff xlo xhi jl_spl(xlo) jl_spl(xhi) I I0 I≈I0
        error("Rb̃llnn error")
    end
    return I
end

function get_kN(cfb, L)
    kk, ll = CryoFaBs.idx2mode(cfb, 1:size(cfb,2))
    kN = kk[ll .== L]
    # There are ~2L+1 sets of keff that are equal, choose one. (If we don't know which one it is, then we don't care which one it is.)
    rlo = cfb.rbounds[1:end-1]
    kN = kN[1:length(rlo)]
    return kN
end

function get_kj(cfb, L; nmult=10)
    kN = get_kN(cfb, L)
    #kj = kN
    Δr = (cfb.rmax - cfb.rmin) / length(cfb.reff)
    kmax = 2 / Δr
    #kj = kN
    kj = range(0.0, π/2*maximum(kN), length=nmult*length(kN))
    #kj = range(minimum(kN), maximum(kN), length=4*length(kN))
    return kj
end

function radial_cryo_window(cfb::AngRadCryoFaB, iL::Integer, ℓ::Integer, k::AbstractVector{Float64})
    rlo = cfb.rbounds[1:end-1]
    rhi = cfb.rbounds[2:end]
    Rb̃ll = calc_Rb̃llnn.(ℓ, rlo, rhi, k')
    Rinv = blocks(cfb.RadTransInv)[iL]
    b̃ll = Rinv * Rb̃ll
    return b̃ll::Matrix{Float64}
end

function radial_cryo_window2(cfb::AngRadCryoFaB, iL::Integer, ℓ::Integer, k::AbstractVector{Float64})
    rlo = cfb.rbounds[1:end-1]
    rhi = cfb.rbounds[2:end]
    xmin = k[1]*rlo[1]
    xmax = k[end]*rhi[end]
    nsamps = (xmax - xmin) * 20 / π
    x = range(xmin, xmax, length=ceil(Int, nsamps))
    jl_spl = Spline1D(x, sphericalbesselj.(ℓ, x))
    Rb̃ll = calc_Rb̃llnn.(ℓ, rlo, rhi, k', jl_spl)
    Rinv = blocks(cfb.RadTransInv)[iL]
    b̃ll = Rinv * Rb̃ll
    return b̃ll::Matrix{Float64}
end

# vim: set sw=4 et sts=4 :


text:
"""
    T2mapSEcorr([io = stderr,] image::Array{T,4}; <keyword arguments>)
    T2mapSEcorr([io = stderr,] image::Array{T,4}, opts::T2mapOptions{T})

Uses nonnegative least squares (NNLS) to compute T2 distributions in the presence of stimulated echos by optimizing the refocusing pulse flip angle.
Records parameter maps and T2 distributions for further partitioning.

# Arguments
- `image`: 4D array with intensity data as `(row, column, slice, echo)`
- A series of optional keyword argument settings which will be used to construct a [`T2mapOptions`](@ref) struct internally, or a [`T2mapOptions`](@ref) struct directly

# Outputs
- `maps`: dictionary containing parameter maps with the following fields:
    - **Default Fields**
        - `"echotimes"`     Echo times of time signal (length `nTE` 1D array)
        - `"t2times"`       T2 times corresponding to T2-distributions (length `nT2` 1D array)
        - `"refangleset"`   Refocusing angles used during flip angle optimization (length `nRefAngles` 1D array by default; scalar if `SetFlipAngle` is used)
        - `"decaybasisset"` Decay basis sets corresponding to `"refangleset"` (`nTE x nT2 x nRefAngles` 3D array by default; `nTE x nT2` 2D array if `SetFlipAngle` is used)
        - `"gdn"`:          Map of general density = sum(T2distribution) (`MatrixSize` 3D array)
        - `"ggm"`:          Map of general geometric mean of T2-distribution (`MatrixSize` 3D array)
        - `"gva"`:          Map of general variance (`MatrixSize` 3D array)
        - `"fnr"`:          Map of fit to noise ratio = gdn / sqrt(sum(residuals.^2) / (nTE-1)) (`MatrixSize` 3D array)
        - `"snr"`:          Map of signal to noise ratio = maximum(signal) / std(residuals) (`MatrixSize` 3D array)
        - `"alpha"`:        Map of optimized refocusing pulse flip angle (`MatrixSize` 3D array)
    - **Optional Fields**
        - `"resnorm"`:      ``\\ell^2``-norm of NNLS fit residuals; see `SaveResidualNorm` option (`MatrixSize` 3D array)
        - `"decaycurve"`:   Signal decay curve resulting from NNLS fit; see `SaveDecayCurve` option (`MatrixSize x nTE` 4D array)
        - `"mu"`:           Regularization parameter used during from NNLS fit; see `SaveRegParam` option (`MatrixSize` 3D array)
        - `"chi2factor"`:   ``\\chi^2`` increase factor relative to unregularized NNLS fit; see `SaveRegParam` option (`MatrixSize` 3D array)
        - `"decaybasis"`:   Decay bases resulting from flip angle optimization; see `SaveNNLSBasis` option (`MatrixSize x nTE x nT2` 5D array, or `nTE x nT2` 2D array if `SetFlipAngle` is used)
- `distributions`: T2-distribution array with data as `(row, column, slice, T2 amplitude)` (`MatrixSize x nT2` 4D array)

# Examples
```julia-repl
julia> image = DECAES.mock_image(MatrixSize = (100,100,1), nTE = 32); # mock image with size 100x100x1x32

julia> maps, dist = T2mapSEcorr(image; TE = 10e-3, nT2 = 40, T2Range = (10e-3, 2.0), Reg = "lcurve", Silent = true); # compute the T2-maps and T2-distribution

julia> maps
Dict{String, Any} with 10 entries:
  "echotimes"     => [0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1  …  0.23, …
  "t2times"       => [0.01, 0.0114551, 0.013122, 0.0150315, 0.0172188, 0.0197244, 0.022594…
  "refangleset"   => [50.0, 54.1935, 58.3871, 62.5806, 66.7742, 70.9677, 75.1613, 79.3548,…
  "gdn"           => [185326.0 1.26487e5 … 1.85322e5 1.025e5; 1.48254e5 1.91871e5 … 1.5159…
  "fnr"           => [460.929 398.485 … 415.599 428.055; 316.05 340.206 … 372.7 410.487; ……
  "alpha"         => [165.348 164.657 … 163.031 166.735; 165.212 164.994 … 165.6 165.094; …
  "gva"           => [0.438218 0.330969 … 0.367929 0.350691; 0.287478 0.307593 … 0.3892 0.…
  "ggm"           => [0.0501681 0.0518976 … 0.0509003 0.0519268; 0.0535984 0.0528403 … 0.0…
  "snr"           => [357.763 317.439 … 322.451 340.88; 250.674 275.185 … 294.411 324.464;…
  "decaybasisset" => [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0…
```

See also:
* [`T2partSEcorr`](@ref)
* [`lsqnonneg`](@ref)
* [`lsqnonneg_chi2`](@ref)
* [`lsqnonneg_gcv`](@ref)
* [`lsqnonneg_lcurve`](@ref)
* [`EPGdecaycurve`](@ref)
"""
T2mapSEcorr(image::Array{T,4}; kwargs...) where {T} = T2mapSEcorr(stderr, image; kwargs...)
T2mapSEcorr(io::IO, image::Array{T,4}; kwargs...) where {T} = T2mapSEcorr(io, image, T2mapOptions(image; kwargs...))
T2mapSEcorr(image::Array{T,4}, opts::T2mapOptions{T}) where {T} = T2mapSEcorr(stderr, image, opts)

function T2mapSEcorr(io::IO, image::Array{T,4}, opts::T2mapOptions{T}) where {T}
    # =========================================================================
    # Initialize output data structures and thread-local buffers
    # =========================================================================
    @assert size(image) == (opts.MatrixSize..., opts.nTE)

    # Print settings to terminal
    !opts.Silent && printbody(io, _show_string(opts))

    # =========================================================================
    # Initialization
    # =========================================================================
    # Initialize output maps and distributions
    maps = init_output_t2maps(opts)
    distributions = init_output_t2distributions(opts)

    # Initialize reusable temporary buffers
    thread_buffers = [thread_buffer_maker(opts) for _ in 1:Threads.nthreads()]

    # =========================================================================
    # Process all pixels
    # =========================================================================
    LinearAlgebra.BLAS.set_num_threads(1) # Prevent BLAS from stealing julia threads

    # Run analysis in parallel
    indices = filter(I -> image[I,1] > opts.Threshold, CartesianIndices(opts.MatrixSize))
    if opts.Silent
        tforeach(indices) do I
            thread_buffer = thread_buffers[Threads.threadid()]
            voxelwise_T2_distribution!(thread_buffer, maps, distributions, image, opts, I)
        end
    else
        bigblocks = Iterators.partition(indices, 8 * Threads.nthreads() * default_blocksize()) .|> copy
        progmeter = DECAESProgress(io, length(bigblocks), "Computing T2-Distribution: "; dt = 5.0)
        for bigblock in bigblocks
            tforeach(bigblock) do I
                thread_buffer = thread_buffers[Threads.threadid()]
                voxelwise_T2_distribution!(thread_buffer, maps, distributions, image, opts, I)
            end
            next!(progmeter)
        end
    end

    LinearAlgebra.BLAS.set_num_threads(Threads.nthreads()) # Reset BLAS threads

    return maps, distributions
end

function init_output_t2distributions(opts::T2mapOptions{T}) where {T}
    distributions = fill(T(NaN), opts.MatrixSize..., opts.nT2)
    return distributions
end

function init_output_t2maps(opts::T2mapOptions{T}) where {T}
    maps = Dict{String,Any}()
    init_output_t2maps!(thread_buffer_maker(opts), maps, opts)
    return maps
end

function init_output_t2maps!(thread_buffer, maps, opts::T2mapOptions{T}) where {T}
    @unpack T2_times, flip_angles, refcon_angles, decay_basis = thread_buffer

    # Misc. processing parameters
    maps["echotimes"]      = convert(Array{T}, copy(opts.TE .* (1:opts.nTE)))
    maps["t2times"]        = convert(Array{T}, copy(T2_times))
    maps["refangleset"]    = opts.SetFlipAngle === nothing ? convert(Array{T}, copy(flip_angles)) : T(opts.SetFlipAngle)
    maps["refconangleset"] = opts.SetRefConAngle === nothing ? convert(Array{T}, copy(refcon_angles)) : T(opts.SetRefConAngle)
    # maps["decaybasisset"]  = convert(Array{T}, copy(decay_basis_set))

    # Default output maps
    maps["gdn"]   = fill(T(NaN), opts.MatrixSize...)
    maps["ggm"]   = fill(T(NaN), opts.MatrixSize...)
    maps["gva"]   = fill(T(NaN), opts.MatrixSize...)
    maps["fnr"]   = fill(T(NaN), opts.MatrixSize...)
    maps["snr"]   = fill(T(NaN), opts.MatrixSize...)
    maps["alpha"] = fill(T(NaN), opts.MatrixSize...)
    maps["beta"]  = fill(T(NaN), opts.MatrixSize...)

    # Optional output maps
    if opts.SaveResidualNorm
        maps["resnorm"] = fill(T(NaN), opts.MatrixSize...)
    end

    if opts.SaveDecayCurve
        maps["decaycurve"] = fill(T(NaN), opts.MatrixSize..., opts.nTE)
    end

    if opts.SaveRegParam
        maps["mu"] = fill(T(NaN), opts.MatrixSize...)
        maps["chi2factor"] = fill(T(NaN), opts.MatrixSize...)
    end

    if opts.SaveNNLSBasis
        if opts.SetFlipAngle === nothing || opts.SetRefConAngle === nothing
            maps["decaybasis"] = fill(T(NaN), opts.MatrixSize..., opts.nTE, opts.nT2) # unique decay basis set for each voxel
        else
            maps["decaybasis"] = convert(Array{T}, copy(decay_basis)) # single decay basis set used for all voxels
        end
    end

    return nothing
end

# =========================================================
# Main loop function
# =========================================================
function voxelwise_T2_distribution!(thread_buffer, maps, distributions, image, opts::T2mapOptions, I::CartesianIndex)
    # Skip low signal voxels
    @inbounds if !(image[I,1] > opts.Threshold)
        return nothing
    end

    # Extract decay curve from the voxel
    @inbounds @simd for j in 1:opts.nTE
        thread_buffer.decay_data[j] = image[I,j]
    end

    # Find optimum flip angle and compute EPG decay basis
    if opts.SetFlipAngle === nothing || opts.SetRefConAngle === nothing
        optimize_flip_angle!(thread_buffer.flip_angle_work, opts)
    end

    # Calculate T2 distribution and map parameters
    t2_distribution!(thread_buffer.T2_dist_work)

    # Save loop results to outputs
    save_results!(thread_buffer, maps, distributions, opts, I)

    return nothing
end

# =========================================================
# EPG decay basis set construction
# =========================================================
struct EPGBasisSetFunctor{T, ETL, opt_vars, Tθ <: EPGParameterization{T,ETL}, W <: AbstractEPGWorkspace{T, ETL}, F <: EPGFunctor{T, ETL, opt_vars}, J <: EPGJacobianFunctor{T, ETL, opt_vars}}
    θ::Tθ
    T2_times::Vector{T}
    epg_work::W
    epg_functor!::F
    epg_jac_functor!::J
end

function EPGBasisSetFunctor(o::T2mapOptions{T}, θ::EPGParameterization{T,ETL}, ::Val{opt_vars}) where {T, ETL, opt_vars}
    epg_work = EPGdecaycurve_work(θ)
    epg_functor! = EPGFunctor(θ, opt_vars)
    epg_jac_functor! = EPGJacobianFunctor(θ, opt_vars)
    EPGBasisSetFunctor(θ, t2_times(o), epg_work, epg_functor!, epg_jac_functor!)
end

#### EPG basis set

function epg_decay_basis!(f::EPGBasisSetFunctor{T,ETL}, decay_basis::AbstractMatrix{T}, θ::EPGParameterization{T,ETL}) where {T,ETL}
    epg_decay_basis!(decay_basis, f.epg_work, θ, f.T2_times)
end
epg_decay_basis!(f::EPGBasisSetFunctor{T,ETL}, decay_basis::AbstractMatrix{T}, x::SVector{D,T}, ::Val{opt_vars}) where {D,T,ETL,opt_vars} = epg_decay_basis!(f, decay_basis, restructure(f.θ, x, Val(opt_vars)))
epg_decay_basis!(f::EPGBasisSetFunctor{T,ETL,opt_vars}, decay_basis::AbstractMatrix{T}, x::SVector{D,T}) where {D,T,ETL,opt_vars} = epg_decay_basis!(f, decay_basis, restructure(f.θ, x, Val(opt_vars)))

function epg_decay_basis!(decay_basis::AbstractMatrix{T}, decay_curve_work::AbstractEPGWorkspace{T,ETL}, θ::EPGParameterization{T,ETL}, T2_times::AbstractVector) where {T,ETL}
    # Compute the NNLS basis over T2 space
    @inbounds for j in 1:length(T2_times)
        decay_curve = uview(decay_basis, :, j) # `UnsafeArrays.uview` is a bit faster than `Base.view`
        θj = restructure(θ, (; T2 = T2_times[j])) # remake options with T2 of basis `j`
        EPGdecaycurve!(decay_curve, decay_curve_work, θj)
    end
    return decay_basis
end

function epg_decay_basis(θ::EPGParameterization{T,ETL}, T2_times::AbstractVector) where {T,ETL}
    decay_basis = zeros(T, ETL, length(T2_times))
    decay_curve_work = EPGdecaycurve_work(θ)
    epg_decay_basis!(decay_basis, decay_curve_work, θ, T2_times)
end

#### Jacobian of EPG basis set

function ∇epg_decay_basis!(f::EPGBasisSetFunctor{T,ETL}, ∇decay_basis::AbstractArray{T,3}, decay_basis::AbstractMatrix{T}, θ::EPGParameterization{T,ETL}) where {T,ETL}
    ∇epg_decay_basis!(∇decay_basis, decay_basis, f.epg_jac_functor!, θ, f.T2_times)
end
∇epg_decay_basis!(f::EPGBasisSetFunctor{T,ETL}, ∇decay_basis::AbstractArray{T,3}, decay_basis::AbstractMatrix{T}, x::SVector{D,T}, ::Val{opt_vars}) where {D,T,ETL,opt_vars} = ∇epg_decay_basis!(f, ∇decay_basis, decay_basis, restructure(f.θ, x, Val(opt_vars)))
∇epg_decay_basis!(f::EPGBasisSetFunctor{T,ETL,opt_vars}, ∇decay_basis::AbstractArray{T,3}, decay_basis::AbstractMatrix{T}, x::SVector{D,T}) where {D,T,ETL,opt_vars} = ∇epg_decay_basis!(f, ∇decay_basis, decay_basis, restructure(f.θ, x, Val(opt_vars)))

function ∇epg_decay_basis!(∇decay_basis::AbstractArray{T,3}, decay_basis::AbstractMatrix{T}, decay_curve_jac!::EPGJacobianFunctor{T,ETL}, θ::EPGParameterization{T,ETL}, T2_times::AbstractVector) where {T,ETL}
    # Compute the NNLS basis over T2 space
    @inbounds for j in 1:length(T2_times)
        decay_curve = uview(decay_basis, :, j) # `UnsafeArrays.uview` is a bit faster than `Base.view`
        ∇decay_curve = uview(∇decay_basis, :, j, :)
        θj = restructure(θ, (; T2 = T2_times[j])) # remake options with T2 of basis `j`
        decay_curve_jac!(∇decay_curve, decay_curve, θj)
    end
    return ∇decay_basis
end

function ∇epg_decay_basis(θ::EPGParameterization{T,ETL}, T2_times::AbstractVector, Fs::NTuple{N,Symbol}) where {T,ETL,N}
    decay_basis = zeros(T, ETL, length(T2_times))
    ∇decay_basis = zeros(T, ETL, length(T2_times), N)
    decay_curve_jac! = EPGJacobianFunctor(θ, Fs)
    ∇epg_decay_basis!(∇decay_basis, decay_basis, decay_curve_jac!, θ, T2_times)
    return decay_basis, ∇decay_basis
end

# =========================================================
# Ensemble of EPG decay basis sets for discrete parameter search
# =========================================================
struct EPGBasisSetEnsemble{D, T, ETL, opt_vars, A1<:AbstractArray{T}, A2<:AbstractArray{T}, F <: EPGBasisSetFunctor{T, ETL, opt_vars}, P <: NNLSDiscreteSurrogateSearch{D, T}}
    decay_basis_set::A1
    ∇decay_basis_set::A2
    epg_basis_functor!::F
    nnls_search_prob::P
end

function EPGBasisSetEnsemble(o::T2mapOptions{T}, θ::EPGOptions{T,ETL}, ::Val{opt_vars}, decay_data::AbstractVector{T}) where {T, ETL, opt_vars}
    opt_ranges = ntuple(length(opt_vars)) do i
        opt_vars[i] === :α ?
            flip_angles(o) :
        opt_vars[i] === :β ?
            refcon_angles(o) :
            error("Optimization variable must be one of :α or :β")
    end
    decay_basis_set    = zeros(T, ETL, o.nT2, length.(opt_ranges)...)
    ∇decay_basis_set   = zeros(T, ETL, o.nT2, length(opt_vars), length.(opt_ranges)...)
    epg_basis_functor! = EPGBasisSetFunctor(o, θ, Val(opt_vars))
    nnls_search_prob   = NNLSDiscreteSurrogateSearch(decay_basis_set, ∇decay_basis_set, opt_ranges, decay_data; legacy = o.legacy)
    EPGBasisSetEnsemble(decay_basis_set, ∇decay_basis_set, epg_basis_functor!, nnls_search_prob)
end

function epg_decay_basis!(work::EPGBasisSetEnsemble{D,T,ETL,opt_vars}, θ::EPGParameterization{T,ETL}) where {D,T,ETL,opt_vars}
    @inbounds for I in CartesianIndices(work.nnls_search_prob.αs)
        x = work.nnls_search_prob.αs[I]
        θ = restructure(θ, x, Val(opt_vars))
        decay_basis = uview(work.decay_basis_set, :, :, I)
        epg_decay_basis!(work.epg_basis_functor!, decay_basis, θ)
    end
    return work
end
epg_decay_basis!(work::EPGBasisSetEnsemble{D,T,ETL}, x::SVector{D,T}, ::Val{opt_vars}) where {D,T,ETL,opt_vars} = epg_decay_basis!(work, restructure(work.epg_basis_functor!.θ, x, Val(opt_vars)))
epg_decay_basis!(work::EPGBasisSetEnsemble{D,T,ETL,opt_vars}, x::SVector{D,T}) where {D,T,ETL,opt_vars} = epg_decay_basis!(work, restructure(work.epg_basis_functor!.θ, x, Val(opt_vars)))

function ∇epg_decay_basis!(work::EPGBasisSetEnsemble{D,T,ETL,opt_vars}, θ::EPGParameterization{T,ETL}) where {D,T,ETL,opt_vars}
    @inbounds for I in CartesianIndices(work.nnls_search_prob.αs)
        x = work.nnls_search_prob.αs[I]
        θ = restructure(θ, x, Val(opt_vars))
        decay_basis = uview(work.decay_basis_set, :, :, I)
        ∇decay_basis = uview(work.∇decay_basis_set, :, :, :, I)
        ∇epg_decay_basis!(work.epg_basis_functor!, ∇decay_basis, decay_basis, θ)
    end
    return work
end
∇epg_decay_basis!(work::EPGBasisSetEnsemble{D,T,ETL}, x::SVector{D,T}, ::Val{opt_vars}) where {D,T,ETL,opt_vars} = ∇epg_decay_basis!(work, restructure(work.epg_basis_functor!.θ, x, Val(opt_vars)))
∇epg_decay_basis!(work::EPGBasisSetEnsemble{D,T,ETL,opt_vars}, x::SVector{D,T}) where {D,T,ETL,opt_vars} = ∇epg_decay_basis!(work, restructure(work.epg_basis_functor!.θ, x, Val(opt_vars)))

# =========================================================
# Flip angle optimization
# =========================================================
struct FlipAngleOptimizationWorkspace{T, ETL, A1<:AbstractMatrix{T}, A2<:AbstractVector{T}, B <: EPGBasisSetFunctor{T,ETL}, E1 <: Union{Nothing, <:EPGBasisSetEnsemble{1, T, ETL}}, E2 <: Union{Nothing, <:EPGBasisSetEnsemble{1, T, ETL}}, E3 <: Union{Nothing, <:EPGBasisSetEnsemble{2, T, ETL}}, S1 <: Union{Nothing, AbstractSurrogate{1, T}}, S2 <: Union{Nothing, AbstractSurrogate{1, T}}, S3 <: Union{Nothing, AbstractSurrogate{2, T}}}
    decay_basis::A1
    decay_data::A2
    α::Base.RefValue{T}
    β::Base.RefValue{T}
    αβ_basis_set::B
    α_basis_set_ensemble::E1
    β_basis_set_ensemble::E2
    αβ_basis_set_ensemble::E3
    α_surrogate::S1
    β_surrogate::S2
    αβ_surrogate::S3
end

function FlipAngleOptimizationWorkspace(o::T2mapOptions{T}, decay_basis::AbstractMatrix{T}, decay_data::AbstractVector{T}) where {T}
    α = Ref(o.SetFlipAngle === nothing ? T(NaN) : o.SetFlipAngle)
    β = Ref(o.SetRefConAngle === nothing ? T(NaN) : o.SetRefConAngle)
    θ = EPGOptions((; α = α[], TE = o.TE, T2 = T(NaN), T1 = o.T1, β = β[]), Val(o.nTE))
    αβ_basis_set = EPGBasisSetFunctor(o, θ, Val((:α, :β)))
    α_basis_set_ensemble = β_basis_set_ensemble = αβ_basis_set_ensemble = nothing
    α_surrogate = β_surrogate = αβ_surrogate = nothing

    if o.SetFlipAngle !== nothing && o.SetRefConAngle !== nothing
        # Compute basis for fixed `SetFlipAngle` and `SetRefConAngle`
        epg_decay_basis!(αβ_basis_set, decay_basis, SA{T}[α[], β[]])
    else
        # Compute basis for each angle
        if o.SetFlipAngle === nothing
            α_basis_set_ensemble = EPGBasisSetEnsemble(o, θ, Val((:α,)), decay_data)
            ∇epg_decay_basis!(α_basis_set_ensemble, θ)
            α_surrogate = o.legacy ?
                CubicSplineSurrogate(α_basis_set_ensemble.nnls_search_prob; legacy = true) :
                HermiteSplineSurrogate(α_basis_set_ensemble.nnls_search_prob)
        end
        if o.SetRefConAngle === nothing
            β_basis_set_ensemble = EPGBasisSetEnsemble(o, θ, Val((:β,)), decay_data)
            ∇epg_decay_basis!(β_basis_set_ensemble, θ)
            β_surrogate = CubicSplineSurrogate(β_basis_set_ensemble.nnls_search_prob) do I
                if o.SetFlipAngle === nothing
                    βI = β_basis_set_ensemble.nnls_search_prob.αs[I][1]
                    decay_basis = uview(β_basis_set_ensemble.decay_basis_set, :, :, I)
                    epg_decay_basis!(β_basis_set_ensemble.epg_basis_functor!, decay_basis, SA[α[], βI], Val((:α, :β)))
                end
            end
        end
        if o.SetFlipAngle === nothing && o.SetRefConAngle === nothing
            αβ_basis_set_ensemble = EPGBasisSetEnsemble(o, θ, Val((:α, :β)), decay_data)
            ∇epg_decay_basis!(αβ_basis_set_ensemble, θ)
            αβ_surrogate = HermiteSplineSurrogate(αβ_basis_set_ensemble.nnls_search_prob)
        end
    end

    FlipAngleOptimizationWorkspace(decay_basis, decay_data, α, β, αβ_basis_set, α_basis_set_ensemble, β_basis_set_ensemble, αβ_basis_set_ensemble, α_surrogate, β_surrogate, αβ_surrogate)
end

function optimize_flip_angle!(work::FlipAngleOptimizationWorkspace, o::T2mapOptions)

    if o.SetFlipAngle === nothing
        empty!(work.α_surrogate)
        nnls_searcher = DiscreteSurrogateSearcher(work.α_surrogate; mineval = o.nRefAnglesMin, maxeval = o.nRefAngles)
        α_opt, _ = bisection_search(work.α_surrogate, nnls_searcher; maxeval = o.nRefAngles)
        work.α[] = α_opt[1]
    end

    if o.SetRefConAngle === nothing
        empty!(work.β_surrogate)
        nnls_searcher = DiscreteSurrogateSearcher(work.β_surrogate; mineval = o.nRefAnglesMin, maxeval = o.nRefAngles)
        β_opt, _ = bisection_search(work.β_surrogate, nnls_searcher; maxeval = o.nRefAngles)
        work.β[] = β_opt[1]
    end

    # Compute basis using optimized flip angles
    if o.SetFlipAngle === nothing || o.SetRefConAngle === nothing
        epg_decay_basis!(work.αβ_basis_set, work.decay_basis, SA[work.α[], work.β[]])
    end

    return nothing
end

# =========================================================
# T2-distribution fitting
# =========================================================
abstract type RegularizationMethod end
struct NoRegularization <: RegularizationMethod end
struct ChiSquared{T} <: RegularizationMethod; Chi2Factor::T; legacy::Bool; end
struct GCV <: RegularizationMethod end
struct LCurve <: RegularizationMethod end

function regularization_method(o::T2mapOptions)
    o.Reg == "none"   ? NoRegularization() : # Fit T2 distribution using unregularized NNLS
    o.Reg == "chi2"   ? ChiSquared(o.Chi2Factor, o.legacy) : # Fit T2 distribution using chi2-based regularized NNLS
    o.Reg == "gcv"    ? GCV() : # Fit T2 distribution using GCV-based regularized NNLS
    o.Reg == "lcurve" ? LCurve() : # Fit T2 distribution using L-curve-based regularized NNLS
    error("Unrecognized regularization method: $(o.Reg)")
end

nnls_workspace(::NoRegularization, decay_basis::AbstractMatrix{T}, decay_data::AbstractVector{T}) where {T} = lsqnonneg_work(decay_basis, decay_data)
nnls_workspace(::ChiSquared, decay_basis::AbstractMatrix{T}, decay_data::AbstractVector{T}) where {T} = lsqnonneg_chi2_work(decay_basis, decay_data)
nnls_workspace(::GCV, decay_basis::AbstractMatrix{T}, decay_data::AbstractVector{T}) where {T} = lsqnonneg_gcv_work(decay_basis, decay_data)
nnls_workspace(::LCurve, decay_basis::AbstractMatrix{T}, decay_data::AbstractVector{T}) where {T} = lsqnonneg_lcurve_work(decay_basis, decay_data)

struct T2DistWorkspace{Reg,T,W}
    reg::Reg
    nnls_work::W
    μ::Base.RefValue{T}
    χ²fact::Base.RefValue{T}
end
function T2DistWorkspace(reg::RegularizationMethod, decay_basis::AbstractMatrix{T}, decay_data::AbstractVector{T}, μ::Base.RefValue{T} = Ref(T(NaN)), χ²fact::Base.RefValue{T} = Ref(T(NaN))) where {T}
    T2DistWorkspace(reg, nnls_workspace(reg, decay_basis, decay_data), μ, χ²fact)
end

function nnls!(t2work::T2DistWorkspace{NoRegularization, T}) where {T}
    x, t2work.μ[], t2work.χ²fact[] = lsqnonneg!(t2work.nnls_work), zero(T), one(T)
end

function nnls!(t2work::T2DistWorkspace{ChiSquared{T}, T}) where {T}
    x, t2work.μ[], t2work.χ²fact[] = lsqnonneg_chi2!(t2work.nnls_work, t2work.reg.Chi2Factor; legacy = t2work.reg.legacy)
end

function nnls!(t2work::T2DistWorkspace{GCV, T}) where {T}
    x, t2work.μ[], t2work.χ²fact[] = lsqnonneg_gcv!(t2work.nnls_work)
end

function nnls!(t2work::T2DistWorkspace{LCurve, T}) where {T}
    x, t2work.μ[], t2work.χ²fact[] = lsqnonneg_lcurve!(t2work.nnls_work)
end

t2_distribution!(t2work::T2DistWorkspace) = nnls!(t2work)
solution(t2work::T2DistWorkspace) = solution(t2work.nnls_work)

# =========================================================
# Save thread local results to output maps
# =========================================================
function save_results!(thread_buffer, maps, distributions, o::T2mapOptions, I::CartesianIndex)
    @unpack T2_dist_work, flip_angle_work, logT2_times, decay_data, decay_basis, decay_calc, residuals, gva_buf = thread_buffer
    T2_dist = solution(T2_dist_work)

    # Compute and save parameters of distribution
    @inbounds begin
        @unpack gdn, ggm, gva, fnr, snr, alpha, beta = maps
        mul!(decay_calc, decay_basis, T2_dist)
        residuals .= decay_calc .- decay_data
        gdn[I]   = sum(T2_dist)
        ggm[I]   = exp(dot(T2_dist, logT2_times) / sum(T2_dist))
        gva_buf .= (logT2_times .- log(ggm[I])).^2 .* T2_dist
        gva[I]   = exp(sum(gva_buf) / sum(T2_dist)) - 1
        fnr[I]   = sum(T2_dist) / sqrt(sum(abs2, residuals)/(o.nTE-1))
        snr[I]   = maximum(decay_data) / std(residuals)
        alpha[I] = flip_angle_work.α[]
        beta[I]  = flip_angle_work.β[]
    end

    # Save distribution
    @inbounds @simd for j in 1:o.nT2
        distributions[I,j] = T2_dist[j]
    end

    # Optionally save regularization parameters
    if o.SaveRegParam
        @unpack mu, chi2factor = maps
        @inbounds mu[I], chi2factor[I] = T2_dist_work.μ[], T2_dist_work.χ²fact[]
    end

    # Optionally save ℓ²-norm of residuals
    if o.SaveResidualNorm
        @unpack resnorm = maps
        @inbounds resnorm[I] = sqrt(sum(abs2, residuals))
    end

    # Optionally save signal decay curve from fit
    if o.SaveDecayCurve
        @unpack decaycurve = maps
        @inbounds @simd for j in 1:o.nTE
            decaycurve[I,j] = decay_calc[j]
        end
    end

    # Optionally save NNLS basis
    if o.SaveNNLSBasis && o.SetFlipAngle === nothing
        @unpack decaybasis = maps
        @inbounds @simd for J in CartesianIndices((o.nTE, o.nT2))
            decaybasis[I,J] = decay_basis[J]
        end
    end

    return nothing
end

# =========================================================
# Utility functions
# =========================================================
function thread_buffer_maker(o::T2mapOptions{T}) where {T}
    decay_basis = zeros(T, o.nTE, o.nT2)
    decay_data  = zeros(T, o.nTE)
    buffer = (
        T2_times         = logrange(o.T2Range..., o.nT2),
        logT2_times      = log.(logrange(o.T2Range..., o.nT2)),
        flip_angles      = flip_angles(o),
        refcon_angles    = refcon_angles(o),
        decay_basis      = decay_basis,
        decay_data       = decay_data,
        decay_calc       = zeros(T, o.nTE),
        residuals        = zeros(T, o.nTE),
        gva_buf          = zeros(T, o.nT2),
        decay_curve_work = EPGdecaycurve_work(T, o.nTE),
        flip_angle_work  = FlipAngleOptimizationWorkspace(o, decay_basis, decay_data),
        T2_dist_work     = T2DistWorkspace(regularization_method(o), decay_basis, decay_data),
    )
    return buffer
end

ext: jl


text:
"""
`AllenCahnProblem1D` - structure to solve the one-dimensional Allen-Cahn equation, that is, solves

uₜ = Auₓₓ + f(u)

on (a, b), where A = Δ and f(u) = -u³ + u.
"""
struct AllenCahnProblem1D{TLBC, TRBC, TA}
    a::Real # Left endpoint
    b::Real # Right endpoint
    ϵ::Real # Model parameter
    nₓ::Int # Total number of interior spatial mesh points
    Δx::Real # Spatial mesh spacing
    x::Vector # Spatial mesh points including boundary
    nₜ::Int # Number of timesteps
    t_max::Real # Final time for integration
    Δt::Real # Length of timestep
    left_bc::TLBC # Boundary condition on the left endpoint
    right_bc::TRBC # Boundary condition on the right endpoint
    A::TA # Matrix approximation of differential operator for the problem
    rhs::Vector # Right hand side of the end system
    f::Function # Nonlinear part of semilinear ODE system (u - u³ for this problem)
    u::Vector # Solution of the problem
end

"""
`BoundaryCondition` - abstract type to include boundary conditions
"""
abstract type BoundaryCondition end

"""
`NeumannBC` - structure to house a Neumann boundary condition, i.e., a boundary condition for which

uₓ(a) = vₗ  or  uₓ(b) = vᵣ
"""
struct NeumannBC <: BoundaryCondition
    boundary_condition::Real
end

"""
`PeriodicBC` - structure to house a Periodic boundary condition, i.e., a boundary condition for which

u(a) = u(b), uₓ(a) = uₓ(b)
"""
struct PeriodicBC <: BoundaryCondition
end

"""
`TimeSteppingMethod` - abstract type to include time stepping methods
"""
abstract type TimeSteppingMethod end

"""
`BackwardEulerMethod` - structure for Backward Euler method
"""
struct BackwardEulerMethod <: TimeSteppingMethod
end

"""
`CrankNicolsonMethod` - structure for Crank-Nicolson method
"""
struct CrankNicolsonMethod <: TimeSteppingMethod
end


text:
function BSPTree()
    return BSPTree(())
end

function BSPTree(arg0::Object)
    return BSPTree((Object,), arg0)
end

function BSPTree(arg0::SubHyperplane, arg1::BSPTree, arg2::BSPTree, arg3::Object)
    return BSPTree((SubHyperplane, BSPTree, BSPTree, Object), arg0, arg1, arg2, arg3)
end

function copy_self(obj::BSPTree)
    return jcall(obj, "copySelf", BSPTree, ())
end

function get_attribute(obj::BSPTree)
    return jcall(obj, "getAttribute", Object, ())
end

function get_cell(obj::BSPTree, arg0::Point, arg1::jdouble)
    return jcall(obj, "getCell", BSPTree, (Point, jdouble), arg0, arg1)
end

function get_close_cuts(obj::BSPTree, arg0::Point, arg1::jdouble)
    return jcall(obj, "getCloseCuts", List, (Point, jdouble), arg0, arg1)
end

function get_cut(obj::BSPTree)
    return jcall(obj, "getCut", SubHyperplane, ())
end

function get_minus(obj::BSPTree)
    return jcall(obj, "getMinus", BSPTree, ())
end

function get_parent(obj::BSPTree)
    return jcall(obj, "getParent", BSPTree, ())
end

function get_plus(obj::BSPTree)
    return jcall(obj, "getPlus", BSPTree, ())
end

function insert_cut(obj::BSPTree, arg0::Hyperplane)
    return jcall(obj, "insertCut", jboolean, (Hyperplane,), arg0)
end

function insert_in_tree(obj::BSPTree, arg0::BSPTree, arg1::jboolean, arg2::BSPTree_VanishingCutHandler)
    return jcall(obj, "insertInTree", void, (BSPTree, jboolean, BSPTree_VanishingCutHandler), arg0, arg1, arg2)
end

function merge(obj::BSPTree, arg0::BSPTree, arg1::BSPTree_LeafMerger)
    return jcall(obj, "merge", BSPTree, (BSPTree, BSPTree_LeafMerger), arg0, arg1)
end

function prune_around_convex_cell(obj::BSPTree, arg0::Object, arg1::Object, arg2::Object)
    return jcall(obj, "pruneAroundConvexCell", BSPTree, (Object, Object, Object), arg0, arg1, arg2)
end

function set_attribute(obj::BSPTree, arg0::Object)
    return jcall(obj, "setAttribute", void, (Object,), arg0)
end

function split(obj::BSPTree, arg0::SubHyperplane)
    return jcall(obj, "split", BSPTree, (SubHyperplane,), arg0)
end

function visit(obj::BSPTree, arg0::BSPTreeVisitor)
    return jcall(obj, "visit", void, (BSPTreeVisitor,), arg0)
end


text:
# Very simple initial sentinel structure for Graff elements in DFG.
# TODO: Need to flesh out further in next release.

export User, Robot, Session

abstract type AbstractCGNode
end

mutable struct User <: AbstractCGNode
    id::Symbol
    name::String
    description::String
    # labels::Vector{Symbol}
    data::Dict{Symbol, String}
end

mutable struct Robot <: AbstractCGNode
    id::Symbol
    userId::Symbol
    name::String
    description::String
    # labels::Vector{Symbol}
    data::Dict{Symbol, String}
end

mutable struct Session <: AbstractCGNode
    id::Symbol
    robotId::Symbol
    userId::Symbol
    name::String
    description::String
    # labels::Vector{Symbol}
    data::Dict{Symbol, String}
end


text:
export instantaneous_zenith_angle, daily_zenith_angle

# true anomaly, radians
function true_anomaly(MA::FT, e::FT) where {FT <: Real}
    TA = MA + (2*e - FT(1/4)*e^3)*sin(MA)
            + FT(5/4)*e^2*sin(2*MA)
            + FT(13/12)*e^3*sin(3*MA)
    return mod(TA, FT(2π))
end

# equation of time, radians
function equation_of_time(e::FT, MA::FT, γ::FT, ϖ::FT) where {FT <: Real}
    _Δt = -2*e*sin(MA) + tan(γ/2)^2*sin(2*(MA + ϖ));
    return mod(_Δt+FT(π), FT(2π)) - FT(π)
end

# calculate the distance, declination, and hour angle (at lon=0)
function distance_declination_hourangle(::Type{FT},
                                        date::DateTime,
                                        param_set::APS,
                                        eot_correction::Bool,
                                        milankovitch::Bool) where {FT <: Real}
    Ya::FT = year_anom(param_set)
    day_length::FT = Planet.day(param_set)
    d0::FT = orbit_semimaj(param_set)
    M0::FT = mean_anom_epoch(param_set)

    epoch_string = "2000-01-01T11:58:56.816"
    # epoch_string::String = epoch(param_set)
    date0 = DateTime(epoch_string,dateformat"y-m-dTHH:MM:SS.s")

    days_per_year = Ya / day_length;
    Δt_years = FT(datetime2julian(date) - datetime2julian(date0)) / days_per_year;

    # mean anomaly given mean anomaly at epoch
    MA = mod(FT(2π) * (Δt_years) + M0, FT(2π))

    # calculate orbital parameters or take values at J2000
    if milankovitch
        ϖ = FT(ϖ_spline(Δt_years));
        γ = FT(γ_spline(Δt_years));
        e = FT(e_spline(Δt_years));
    else
        ϖ = FT(lon_perihelion_epoch(param_set));
        γ = FT(obliq_epoch(param_set));
        e = FT(eccentricity_epoch(param_set));
    end

    # true anomaly, radians
    TA = true_anomaly(MA, e)

    # true longitude, radians
    TL = mod(TA + ϖ, FT(2π))

    # declination, radians
    δ = mod(asin(sin(γ) * sin(TL)), FT(2π))

    # earth-sun distance
    d = d0 * (1 - e^2) / (1 + e*cos(TA))

    # hour angle, zero at local solar noon, radians
    if eot_correction
        Δt = equation_of_time(e, MA, γ, ϖ) / FT(2π) * day_length # radians to seconds
    else
        Δt = FT(0)
    end
    time_in_day = FT(mod(datetime2julian(date), 1)) * day_length; # s in this day
    η_UTC = mod(FT(2π) * (time_in_day + Δt) / day_length, FT(2π));

    return d, δ, η_UTC
end

"""
    instantaneous_zenith_angle(date::DateTime,
                               longitude::FT,
                               latitude::FT,
                               param_set::APS;
                               eot_correction::Bool=true,
                               milankovitch::Bool=true) where {FT <: Real}

Returns the zenith angle and earth-sun distance
at a particular longitude and latitude on the given date (and time UTC)
given orbital parameters: obliquity, longitude of perihelion, and eccentricity
param_set is an AbstractParameterSet from CLIMAParameters.jl.

`eot_correction` is an optional Boolean keyword argument that defaults to true
when set to true the equation of time correction is turned on.
This switch functionality is implemented for easy comparisons with reanalyses.

`milankovitch` is an optional Boolean keyword argument that defaults to true
when set to true the orbital parameters are calculated for the given DateTime
when set to false the orbital parameters at the J2000 epoch from CLIMAParameters are used.
"""
function instantaneous_zenith_angle(date::DateTime,
                                    longitude::FT,
                                    latitude::FT,
                                    param_set::APS;
                                    eot_correction::Bool=true,
                                    milankovitch::Bool=true) where {FT <: Real}
    λ = deg2rad(longitude)
    ϕ = deg2rad(latitude)

    d, δ, η_UTC = distance_declination_hourangle(FT, date, param_set, eot_correction, milankovitch)

    # hour angle
    η = mod(η_UTC + λ, FT(2π))

    # zenith angle, radians
    θ = mod(acos(cos(ϕ)*cos(δ)*cos(η) + sin(ϕ)*sin(δ)), FT(2π))

    # solar azimuth angle, ζ = 0 when due E and increasing CCW
    # ζ = 3π/2 (due S) when η=0 at local solar noon
    ζ = mod(FT(3π/2)- atan(sin(η), cos(η)*sin(ϕ) - tan(δ)*cos(ϕ)), FT(2π))

    return θ, ζ, d
end

"""
    daily_zenith_angle(date::DateTime,
                       latitude::FT,
                       param_set::APS;
                       eot_correction::Bool=true,
                       milankovitch::Bool=true) where {FT <: Real}

Returns the effective zenith angle corresponding to the diurnally averaged insolation
and earth-sun distance at a particular latitude given the date.

`param_set` is an AbstractParameterSet from CLIMAParameters.jl.

`eot_correction` is an optional Boolean keyword argument that defaults to true
when set to true the equation of time correction is turned on.
This switch functionality is implemented for easy comparisons with reanalyses.

`milankovitch` is an optional Boolean keyword argument that defaults to true
when set to true the orbital parameters are calculated for the given DateTime,
when set to false the orbital parameters at the J2000 epoch from CLIMAParameters are used.
"""
function daily_zenith_angle(date::DateTime,
                            latitude::FT,
                            param_set::APS;
                            eot_correction::Bool=true,
                            milankovitch::Bool=true) where {FT <: Real}
    ϕ = deg2rad(latitude)

    d, δ, _ = distance_declination_hourangle(FT, date, param_set, eot_correction, milankovitch)

    # sunrise/sunset angle
    T = tan(ϕ) * tan(δ)
    if T >= FT(1)
        ηd = FT(π)
    elseif T <= FT(-1)
        ηd = FT(0)
    else
        ηd = acos(FT(-1)*T)
    end

    # effective zenith angle to get diurnally averaged insolation (i.e., averaging cosine of zenith angle)
    daily_θ = mod(acos(FT(1/π)*(ηd*sin(ϕ)*sin(δ) + cos(ϕ)*cos(δ)*sin(ηd))), FT(2π))

    return daily_θ, d
end


text:
#=

  Discrete tomography in Julia ConstraintSolver.jl

  Note: The origin of the problem is from ECLiPSe,
  but this model has been transformed in this way
     MiniZinc -> SICStus Prolog -> ECLiPSe -> B-Prolog -> Picat -> ConstraintSolver.jl
  Here is my own take at the problem.

  Problem from http://eclipse-clp.org/examples/tomo.ecl.txt
  """
  This is a little "tomography" problem, taken from an old issue
  of Scientific American.

  A matrix which contains zeroes and ones gets "x-rayed" vertically and
  horizontally, giving the total number of ones in each row and column.
  The problem is to reconstruct the contents of the matrix from this
  information. Sample run:

  ?- go.
     0 0 7 1 6 3 4 5 2 7 0 0
  0
  0
  8      * * * * * * * *
  2      *             *
  6      *   * * * *   *
  4      *   *     *   *
  5      *   *   * *   *
  3      *   *         *
  7      *   * * * * * *
  0
  0

  Eclipse solution by Joachim Schimpf, IC-Parc
  """

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

#
# print_matrix_str(x,zstring)
#
# Print each row and convert each element to the corresponding position
# in zstring. "z" is for zero-based, i.e. first char corresponds to 0, etc
#
function print_matrix_str(x,zstring="")
    for row in eachrow(x)
        if zstring == ""
            println(row)
        else
            println(join(row.|>c->zstring[c+1]))
        end
    end
    println()
end

function discrete_tomography(problem,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            # "branch_strategy" => :IMPS, # default
                                                            "branch_strategy" => :ABS, # Activity Based Search
                                                            "activity.decay" => 0.999, # default 0.999
                                                            "activity.max_probes" => 100, # default, 10
                                                            "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    row_sums = problem[:r]
    col_sums = problem[:c]

    rows = length(row_sums)
    cols = length(col_sums)

    @variable(model, x[1:rows,1:cols], Bin)

    # Check rows and cols
    for i in 1:rows
        @constraint(model,row_sums[i] == sum(x[i,:]))
    end
    for j in 1:cols
        @constraint(model,col_sums[j] == sum(x[:,j]))
    end

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                # println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                print_matrix_str(x_val,".X")

            end
        end
    else
        println("status:$status")
    end

    return status
end

discrete_tomography_problems = Dict(

#
# The three first problems are from the ECLiPSe model:
#
# The above stated problem
:1  => Dict(
       :r => [0,0,8,2,6,4,5,3,7,0,0], # row sums
       :c => [0,0,7,1,6,3,4,5,2,7,0,0] # column sums
       ),

:2 => Dict(
       :r => [10,4,8,5,6],
       :c => [5,3,4,0,5,0,5,2,2,0,1,5,1]
),

# This give three slightly different solutions.
:3 => Dict(
        :r => [11,5,4],
        :c => [3,2,3,1,1,1,1,2,3,2,1]
),

# This is my own problem.
:4 => Dict(
        :r => [0,2,2,2,2,2,8,8,4,4,4,4,4,0],
        :c => [0,0,0,12,12,2,2,2,2,7,7,0,0,0]
),

)

for p in sort(collect(keys(discrete_tomography_problems)))
    println("problem $p")
    @time discrete_tomography(discrete_tomography_problems[p])
end

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
module TestConversion
using CustomTest
using Base.Test
using ControlSystems

f11 = tf(ss([-1 0;1 1],[1;0],[1 1],0))
f12 = tf([1,0],[1,0,-1])

#Direct term was previously missing
f21 = tf(ss([-1 0;1 1],[1;0],[1 1],1))
f22 = tf([1,1,-1],[1,0,-1])

#Test ss2tf
@test_approx_eq_eps f11 f12 1e-15
@test_approx_eq_eps f21 f22 1e-15

#Test tf2ss
#ss is not unique so go back and forth
@test_approx_eq_eps f11 tf(ss(f12)) 1e-15
@test_approx_eq_eps f21 tf(ss(f22)) 1e-15

end


text:
function find_files(folder, ext=".dat", delim=';', types = ["CV", "C&D", "EIS"];
        exclude_with=ext, select_with="", rename=true)

    dirs = dirs_in_folder(folder, true)
    data = Dict{String,Vector{DataFile}}()
    for dir in dirs
        files = readdir(dir, join=true)
        files = files[isfile.(files)]

        for file in files
            for type in types
                if occursin(type, file) && occursin(select_with, file) && !exclude(file, exclude_with)
                    if haskey(data, type)
                        push!(data[type], DataFile(file, ext, delim, rename))
                    else
                        push!(data, type=>[DataFile(file, ext, delim, rename)])
                    end
                end
            end
        end
    end

    return data
end

function dirs_in_folder(folder, keep_root)
    contents = readdir(folder, join=true)
    dirs = contents[isdir.(contents)]
    return keep_root ? dirs : basename.(dirs)
end

function files_with_val(datafiles, val, ext="")
    files = DataFile[]
    for file in datafiles
        if filevalue(file, ext) == val
            push!(files, file)
        end
    end

    return files
end

foldervalue(datafile) = splitpath(datafile.filename)[2]

function filevalues(datafiles)
    vals = Dict{String,Set{String}}()
    for file in datafiles
        fdv = foldervalue(file)
        fv = filevalue(file)
        if haskey(vals, fdv)
            push!(vals[fdv], fv)
        else
            push!(vals, fdv=>Set{String}([fv]))
        end
    end

    return vals
end

function common_values(data)
    fvs = filevalues(data)

    common = collect(fvs[first(keys(fvs))])

    for (k,v) in fvs
        intersect!(common, v)
    end

    return common
end

process_data(type::String, data; args...) = process_data(Val(Symbol(type)), data; args...)

value_index(datafile) = 3

function filevalue(datafile, ext="")
    filename = isempty(ext) ? datafile.filename : replace(datafile.filename, ext=>"")
    fn = basename(filename)
    parts = split(fn, '_')
    idx = value_index(datafile)
    replace(parts[idx], " "=>"")
end

function groupbyfolder(datafiles)
    data = Dict{String,Vector{DataFile}}()
    for df in datafiles
        f = foldervalue(df)
        if haskey(data, f)
            push!(data[f], df)
        else
            data[f] = [df]
        end
    end
    return data
end

function results(folder, type, parameter_val, file_val; processed=true, reduction=identity)
    if processed
        data = find_files(folder, exclude_with=r"!", select_with=".dat", rename=false)
        ext = ".dat"
    else
        data = find_files(folder)
        ext = ""
    end

    grouped = groupbyfolder(data[type])
    files = files_with_val(grouped[parameter_val], file_val, ext)
    datafile = only(reduction(files))

    if processed
        df = read_file(datafile, 3, false)
    else
        df = read_file(datafile)
    end

    return datafile, df
end

function clear(dir, to_delete)
    for (root,dirs,files) in walkdir(dir)
        for file in files
            if occursin(to_delete, file)
                rm(joinpath(root, file))
                @info "Deleted $(joinpath(root, file))"
            end
        end
    end
end


text:
module Plots

# using Pkg
using ..DISSEQT
using LinearAlgebra
using Gadfly
using DataFrames
using Colors
using Cairo
using Fontconfig

import PlotlyJS; const py=PlotlyJS # import to avoid clash with Gadfly
#using ORCA # Causes lot's of fontconfig warnings on some systems. If user wants to export PlotlyJS figures, then "import ORCA" can be run in the user script insted.

export
    groupcolors,
    linearregressionplot,
    pairwisescatterplot,
    talusplot,
    fitnesslandscapeplot,
    saveplot

include("colors.jl")
include("stopcodons.jl")
include("pairwise.jl")
include("talus.jl")
include("fitnesslandscape.jl")
include("save.jl")
include("saveplotly.jl")

end


text:
"""
    create_eps(A::PetscMat)

For a standard eigenvalue prolem.
"""
function create_eps(A::PetscMat; auto_setup = false)
    eps = EPSCreate()
    EPSSetOperators(eps, A)

    if (auto_setup)
        set_from_options!(eps)
        set_up!(eps)
    end

    return eps
end

"""
    create_eps(A::PetscMat, B::PetscMat)

For a generalized eigenvalue problem
"""
function create_eps(A::PetscMat, B::PetscMat; auto_setup = false)
    eps = EPSCreate()
    EPSSetOperators(eps, A, B)

    if (auto_setup)
        set_from_options!(eps)
        set_up!(eps)
    end

    return eps
end

PetscWrap.set_up!(eps::SlepcEPS) = EPSSetUp(eps)

PetscWrap.set_from_options!(eps::SlepcEPS) = EPSSetFromOptions(eps)

PetscWrap.solve!(eps::SlepcEPS) = EPSSolve(eps)

const neigs = EPSGetConverged

"""
    get_eig(eps::SlepcEPS, ieig::Integer)

Get the `ieig`-th eigenvalue as a complex number. Starts at 1 (Julia 1-based indexing)
"""
function get_eig(eps::SlepcEPS, ieig::Integer)
    eigr, eigi = EPSGetEigenvalue(eps, ieig - 1)
    return isa(eigr, Real) ? complex(eigr, eigi) : eigr
end

"""
    get_eig(eps::SlepcEPS, ieigs)

Get the `ieig`-th eigenvalues as a complex array. Starts at 1 (Julia 1-based indexing)
"""
function get_eig(eps::SlepcEPS, ieigs)
    eigs = zeros(Complex, size(ieigs))
    for i in ieigs
        eigs[i] = get_eig(eps, ieigs[i])
    end
    return eigs
end

const get_eigenvalue = get_eig

"""
    get_eigs(eps::SlepcEPS)

Returns all the converged eigenvalues as one complex array
"""
get_eigs(eps::SlepcEPS) = get_eig(eps, 1:neigs(eps))

const get_eigenvalues = get_eigs

const get_tolerances = EPSGetTolerances

get_eigenpair(eps::SlepcEPS, ieig) = EPSGetEigenpair(eps, ieig - 1)

PetscWrap.destroy!(eps::SlepcEPS) = EPSDestroy(eps)


text:
using Test
using .Stubs
using DandelionWebSockets

@testset "InProcessIO" begin
    @testset "Write 1-5 on endpoint 1; Can read 1-5 on endpoint 2" begin
        iopair = InProcessIOPair()
        readvalues = UInt8[]
        @sync begin
            @async begin
                for i = 1:5
                    write(iopair.endpoint1, UInt8(i))
                end
            end

            @async begin
                for i = 1:5
                    x = read(iopair.endpoint2, UInt8)
                    push!(readvalues, x)
                end
            end
        end

        @test readvalues == UInt8[1, 2, 3, 4, 5]
    end

    @testset "Write 1-5 on endpoint 1, and 100-110 on endpoint 2; Values can be read independently" begin
        iopair = InProcessIOPair()
        readon1 = UInt8[]
        readon2 = UInt8[]
        @sync begin
            @async begin
                for i = 1:5
                    write(iopair.endpoint1, UInt8(i))
                end
                close(iopair.endpoint1)
            end

            @async begin
                for i = 100:110
                    write(iopair.endpoint2, UInt8(i))
                end
                close(iopair.endpoint2)
            end

            @async begin
                try
                    while true
                        x = read(iopair.endpoint1, UInt8)
                        push!(readon1, x)
                    end
                catch ex
                    if typeof(ex) != EOFError
                        rethrow()
                    end
                end
            end

            @async begin
                try
                    while true
                        x = read(iopair.endpoint2, UInt8)
                        push!(readon2, x)
                    end
                catch ex
                    if typeof(ex) != EOFError
                        rethrow()
                    end
                end
            end
        end

        @test readon1 == UInt8[100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110]
        @test readon2 == UInt8[1, 2, 3, 4, 5]
    end

    @testset "Send 5 values, and then close endpoint 1; EOF is reported on endpoint 2 after 5 values" begin
        iopair = InProcessIOPair()
        readvalues = []
        for i = 1:5
            write(iopair.endpoint1, UInt8(i))
        end
        close(iopair.endpoint1)

        for i = 1:5
            x = read(iopair.endpoint2, UInt8)
            push!(readvalues, x)
        end

        @test readvalues == UInt8[1, 2, 3, 4, 5]
        @test eof(iopair.endpoint2)
    end

    @testset "Send 5 values, and then close endpoint 1; EOF is not reported on endpoint 2 until after 5 values have been read" begin
        iopair = InProcessIOPair()
        for i = 1:5
            write(iopair.endpoint1, UInt8(1))
        end
        close(iopair.endpoint1)

        for i = 1:5
            @test !eof(iopair.endpoint2)
            read(iopair.endpoint2, UInt8)
        end
    end

    @testset "Send a frame on endpoint 1; Received a frame on endpoint 2" begin
        iopair = InProcessIOPair()
        sentframe = Stubs.createserverframe(DandelionWebSockets.OPCODE_TEXT, b"Hello")
        write(iopair.endpoint1, sentframe)
        receivedframe = read(iopair.endpoint2, DandelionWebSockets.Frame)

        @test sentframe.payload == receivedframe.payload
    end

    @testset "Read from empty stream on endpoint 1 and close endpoint 2; Read throws exception" begin
        eofthrown = false
        iopair = InProcessIOPair()

        @sync begin
            @async begin
                try
                    read(iopair.endpoint1, UInt64)
                catch ex
                    if typeof(ex) == EOFError
                        eofthrown = true
                    else
                        rethrow()
                    end
                end
            end

            @async begin
                sleep(0.5)
                close(iopair.endpoint2)
            end
        end

        @test eofthrown
    end
end


text:
using StatsBase

function runall()
    nbuckets = 16
    unfinish = true
    spinner = ReentrantLock()
    buckets = rand(1:99, nbuckets)
    totaltrans = 0

    bucketsum() = sum(buckets)
    smallpause() = sleep(rand() / 2000)
    picktwo() = (samplepair(nbuckets)...)
    function equalizer()
        while unfinish
            smallpause()
            if trylock(spinner)
                i, j = picktwo()
                sm = buckets[i] + buckets[j]
                m = fld(sm + 1, 2)
                buckets[i], buckets[j] = m, sm - m
                totaltrans += 1
                unlock(spinner)
            end
        end
    end
    function redistributor()
        while unfinish
            smallpause()
            if trylock(spinner)
                i, j = picktwo()
                sm = buckets[i] + buckets[j]
                buckets[i] = rand(0:sm)
                buckets[j] = sm - buckets[i]
                totaltrans += 1
                unlock(spinner)
            end
        end
    end
    function accountant()
        count = 0
        while count < 16
            smallpause()
            if trylock(spinner)
                println("Current state of buckets: $buckets. Total in buckets: $(bucketsum())")
                unlock(spinner)
                count += 1
                sleep(1)
            end
        end
        unfinish = false
    end
    t = time()
    @async equalizer()
    @async redistributor()
    @async accountant()
    while unfinish sleep(0.25) end
    println("Total transactions: $totaltrans ($(round(Int, totaltrans / (time() - t))) unlocks per second).")
end

runall()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#Copyright Nathan Parrish 2018
#University of Colorado, Boulder
#Colorado Center for Astrodynamics Research
#Celestial and Spaceflight Mechanics Laboratory
#
#Helper functions for LowThrustOpt

#Simple function to calculate Jacobi constant
function jacobiConstant(state, MU, DU)
    r1 = sqrt.((state[1,:]+MU).^2 + state[2,:].^2 + state[3,:].^2)
    r2 = sqrt.((state[1,:]+MU-1).^2 + state[2,:].^2 + state[3,:].^2)
    v = norm_many(state[4:6,:])
    C = state[1,:].^2 + state[2,:].^2 + 2*(1-MU)./r1 + 2*MU./r2 - v.^2
end

#Spline interpolation for initial and final states as function of τ
function interpInitialStates(p1, X0_times, X0_states, MU)
    #check bounds on endpoint parameters (wrap to 0,1)
    while p1 > 1
    p1 -= 1
    end
    while p1 < 0
    p1 += 1
    end

    #build interpolation object:
    X0_itp = Interpolations.scale(interpolate(X0_states', (BSpline(Cubic(Natural())), NoInterp()), OnGrid()), X0_times, 1:6)

    #endpoints from fixed parameters given:
    state_0 = [X0_itp[p1, 1], X0_itp[p1, 2], X0_itp[p1, 3], X0_itp[p1, 4], X0_itp[p1, 5], X0_itp[p1, 6]]

    #outputs:
    state_0
end

#finds the value of τ that matches the X_states most closely to a new endpoint
function find_τ(X_times, X_states, state, MU)
    τ_trial = linspace(0, 1, 1001)
    d = zeros(length(τ_trial))
    for ind = 1:length(τ_trial)
        state_trial = interpInitialStates(τ_trial[ind], X_times, X_states, MU)
        d[ind] = norm(state_trial - state)
    end

    #return close approach τ:
    τ_trial[d .== minimum(d)][1]
end

#Add more intermediate points in a solution
function densify(XC_all, t_TU, params, n_desired)

    XC_dense = zeros(size(XC_all,1), 0)
    t_dense_db = zeros(0)
    t_dense = linspace(t_TU[1], t_TU[end], n_desired)
    for ind = 1:(size(XC_all, 2)-1)

        # t_save = linspace(t_TU[ind], t_TU[ind+1],
        #     Int32(round(n_desired / (size(XC_all, 2)-1))))

        #propagate and save points
        tspan = (t_TU[ind], t_TU[ind+1])
        prob = ODEProblem(CRTBP_stateCostate_deriv!, XC_all[:,ind], tspan, params)
        sol_forward = DifferentialEquations.solve(prob, Vern8(),
            reltol=1e-13, abstol=1e-13)

        # temp = t_dense[(t_dense .>= tspan[1]) .& (t_dense .< tspan[2])]
        # temp = 1
        # return temp

        if any((t_dense .>= tspan[1]) .& (t_dense .< tspan[2]))
            temp = sol_forward(t_dense[(t_dense .>= tspan[1]) .& (t_dense .< tspan[2])])[:,:]
        else
            temp = zeros(size(XC_dense,1), 0)
        end

        # if any(t_dense .== tspan[2])
        #     warn("argh")
        #     println("t overlap at $(tspan[2])")
        # end

        # display(tspan)
        # display(t_dense[(t_dense .>= tspan[1]) .& (t_dense .< tspan[2])])
        # println("---")

        # XC_dense = hcat(XC_dense, sol_forward[:,1:end-1])
        XC_dense = hcat(XC_dense, temp)
        t_dense_db = append!(t_dense_db, t_dense[(t_dense .>= tspan[1]) .& (t_dense .< tspan[2])])

        #add on the last one:
        if ind == (size(XC_all, 2)-1)
            XC_dense = hcat(XC_dense, sol_forward[:,end])
            t_dense_db = append!(t_dense_db, t_TU[end]) #debug
        end
    end

    (XC_dense, t_dense)
end

#Reduce rho algorithmically with successive indirect multiple shooting
function reduceFuel_indirect(XC_all, t_TU, MU, DU, TU, n_nodes, mass,
    thrustLimit, rho_current, rho_target)

    #The whole point is to reduce rho down to rho_target. Here, catch mistakes:
    if rho_target > rho_current
        warn("Non-sensical rho_target. Fixing...")
        rho_target = rho_current
    end

    #initial values. p is always 1, but rho gets decreased
    p = 1.
    rho_temp = rho_current
    maxIter = 10
    plot_yn = false
    flag_adjointsOnly = false

    #Check convergence on initial value of rho
    (XC_all_new, defect, status_flag) = multiShoot_CRTBP_indirect(XC_all, t_TU, MU, DU, TU, n_nodes,
        mass, thrustLimit, plot_yn, flag_adjointsOnly, maxIter, p, rho_temp)

    if (status_flag == 0) && (rho_current == rho_target)
        #We're golden!
        println("   multiShoot converged at initial rho = $rho_temp")

        return (XC_all_new, defect, status_flag)
    end

    #while not converged, increase rho
    while !(status_flag == 0) && (rho_temp < 1)
        #radius of convergence is larger for larger rho

        println("   Did not converge with rho = $rho_temp. Increasing rho...")
        rho_temp *= 5
        if rho_temp > 1.
            rho_temp = 1.
        end
        (XC_all_new, defect, status_flag) = multiShoot_CRTBP_indirect(XC_all, t_TU, MU, DU, TU, n_nodes,
            mass, thrustLimit, plot_yn, flag_adjointsOnly, maxIter, p, rho_temp)
    end

    if (rho_temp == 1) && !(status_flag == 0)
        warn("  Did not converge at rho = 1. Stopping...")
        return (XC_all_new, defect, status_flag)
    end
    if status_flag == 0
        println("   multiShoot converged at rho = $rho_temp")

        XC_all = copy(XC_all_new)
    end

    #Now, decrease rho again, if needed
    count = 0
    while (rho_temp > rho_target) || !(status_flag == 0)
        count += 1
        if count > 100
            status_flag = 3
            warn("  Did not converge after $count iterations!")

            # XC_all_new *= NaN
            return (XC_all_new, defect, status_flag)
        end

        #accept the new solution if it converged (need to check the defect magnitudes instead of this)
        if status_flag == 0
            XC_all = copy(XC_all_new)

            rho_temp /= 2

            if rho_temp < rho_target
                rho_temp = rho_target
            end
            println("   Reducing rho, new value:   rho = ", rho_temp)

        else
            rho_temp *= 3 * (1 + rand())
            println("   Increasing rho, new value: rho = ", rho_temp)
        end

        #Run single-shooting algorithm again:
        (XC_all_new, defect, status_flag) = multiShoot_CRTBP_indirect(XC_all, t_TU, MU, DU, TU, n_nodes,
            mass, thrustLimit, plot_yn, flag_adjointsOnly, maxIter, p, rho_temp)
    end

    #Outputs:
    (XC_all_new, defect, status_flag)
end

#Add time to the final endpoint
function addTimeFinal(XC_all, t_TU, Δt)

    #-propagate end state forward in time.
    XC_all[7:12,end] = 0.
        tspan_edge = (t_TU[end], t_TU[end] + Δt)
        prob = ODEProblem(CRTBP_stateCostate_deriv!, XC_all[:,end], tspan, params)
        sol = DifferentialEquations.solve(prob, Vern8(), reltol=1e-13, abstol=1e-13)

        XC_all = hcat(XC_all, sol[:, end])
        t_TU = append!(t_TU, tspan_edge[2])

    #-densify, then interpolate on the densified points to find the new XC_all
    n_desired = 200
        (XC_dense, t_dense) = densify(XC_all, t_TU, params, n_desired)

    #build interpolation object:
    XC_itp = Interpolations.scale(interpolate(XC_dense', (BSpline(Cubic(Natural())), NoInterp()), OnGrid()), t_dense, 1:12)

    #endpoints from fixed parameters given:
    #(This is also an opportunity to change the number / distribution of nodes.)
    # n_nodes = 10
    t_TU_new = collect(linspace(t_TU[1], t_TU[end], n_nodes))
    XC_all_new = zeros(12, n_nodes)
    for ind = 1:n_nodes
        XC_all_new[:,ind] = [XC_itp[t_TU_new[ind], 1], XC_itp[t_TU_new[ind], 2], XC_itp[t_TU_new[ind], 3],
            XC_itp[t_TU_new[ind], 4], XC_itp[t_TU_new[ind], 5], XC_itp[t_TU_new[ind], 6],
            XC_itp[t_TU_new[ind], 7], XC_itp[t_TU_new[ind], 8], XC_itp[t_TU_new[ind], 9],
            XC_itp[t_TU_new[ind], 10], XC_itp[t_TU_new[ind], 11], XC_itp[t_TU_new[ind], 12]]
    end

    #-call the function to find the value of τ that matches the Xf_states most closely to the new endpoint
    τ = find_τ(Xf_times, Xf_states, XC_all_new[1:6, end])

    #-Set XC_all[1:6, end] equal to that interpolated Xf_states value
    XC_all_new[1:6, end] = interpInitialStates(τ, Xf_times, Xf_states, MU)

    #-Re-solve indirect problem and repeat
    plot_yn = true
    # (XC_all_new, defect, status_flag) = reduceFuel_indirect(XC_all_new, t_TU, MU, DU, TU, n_nodes, mass,
    #     thrustLimit, rho, rho)
    (XC_all_new, defect, status_flag) = multishoot_CRTBP_indirect_fixedEnds(XC_all_new, t_TU, MU, DU, TU, n_nodes,
        mass, thrustLimit, false, false, maxIter, p, rho)

    if status_flag == 0
        println("Successfully added time!")
        XC_all = copy(XC_all_new)
        t_TU = copy(t_TU_new)
    else
        warn("Failed! Returning original values")
        XC_all = XC_all[:, 1:n_nodes]
        t_TU = t_TU[1:n_nodes]
    end

    (XC_all, t_TU)
end

max_stars_repo_name: travelingspaceman/LowThrustOpt

max_issues_repo_name: travelingspaceman/LowThrustOpt

max_forks_repo_name: travelingspaceman/LowThrustOpt


text:
# LSRTM of Marmousi using stochastic gradient descent
# Author: Philipp Witte, pwitte@eoas.ubc.ca
# Date: August 2017

using PyCall,PyPlot,HDF5,opesciSLIM.TimeModeling,opesciSLIM.LeastSquaresMigration,SeisIO

# Load migration velocity model
n,d,o,m0 = h5open("../data/marmousi_migration_velocity.h5","r") do file
	read(file, "n", "d", "o", "m0")
end

# Load data
block = segy_read("../data/marmousi_2D.segy")
dD = joData(block)

# Set up model structure
model0 = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m0)

# Set up wavelet
src_geometry = Geometry(block; key="source", segy_depth_key="SourceDepth")
wavelet = source(src_geometry.t[1],src_geometry.dt[1],0.03)	# 30 Hz wavelet
q = joData(src_geometry,wavelet)

# Set up info structure
ntComp = get_computational_nt(q.geometry,dD.geometry,model0)	# no. of computational time steps
info = Info(prod(model0.n),dD.nsrc,ntComp)

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

# Setup operators
F = joModeling(info,model0,q.geometry,dD.geometry)
J = joJacobian(F,q)

# Right-hand preconditioners (model topmute)
Mr = opTopmute(model0.n,42,10)	# mute up to grid point 42, with 10 point taper

# Stochastic gradient
x = zeros(Float32,info.n)
batchsize = 10
niter = 32
fval = zeros(Float32,niter)
for j=1:niter
	println("Iteration: ", j)

	# Select batch and set up left-hand preconditioner
	idx = randperm(dD.nsrc)[1:batchsize]
	Jsub = subsample(J,idx)
	dsub = subsample(dD,idx)
	Ml = opMarineTopmute(30,dsub.geometry)	# data topmute starting at nz=30

	# Compute residual and gradient
	r = Ml*Jsub*Mr*x - Ml*dsub
	g = Mr'*Jsub'*Ml'*r

	# Step size and update variable
	fval[j] = .5*norm(r)^2
	t = norm(r)^2/norm(g)^2
	x -= t*g
end

# Save result
h5open("marmousi_lsrtm.h5", "w") do file
	write(file, "x", reshape(x,model0.n))
end

figure(); imshow(reshape(x,model0.n)')

max_stars_repo_name: SINBADconsortium/SLIM-release-jlapps

max_issues_repo_name: SINBADconsortium/SLIM-release-jlapps

max_forks_repo_name: SINBADconsortium/SLIM-release-jlapps


text:
"""
Low-level binding to htslib.
"""
module htslib

using htslib_jll: libhts
import TranscodingStreams

# A macro that makes function binding definitions concise.
# For example, @defun foo(x::S, y::T)::U is expanded to:
#   foo(x, y) = ccall((:foo, libhts), U, (S, T), x, y)
macro defun(ex::Expr)
    @assert ex.head === :(::)
    call = ex.args[1]
    rettype = ex.args[2]
    eval(rettype)  # check that the type is defined
    funcname = call.args[1]
    argnames = []
    argtypes = []
    for i in 2:lastindex(call.args)
        arg = call.args[i]
        @assert arg.head === :(::)
        @assert arg.args[1] isa Symbol
        argname, argtype = arg.args
        push!(argnames, argname)
        push!(argtypes, argtype)
        eval(argtype)  # check that the type is defined
    end
    Expr(:(=),
         Expr(:call, funcname, argnames...),
         Expr(:call,
              :ccall,
              Expr(:tuple, QuoteNode(funcname), :libhts),
              rettype,
              Expr(:tuple, argtypes...),
              argnames...))
end

include("types.jl")
include("constants.jl")
include("functions.jl")

# The following files are machine-generated.  Do not edit them by hand.
include("htslib/bgzf.jl")
include("htslib/knetfile.jl")
include("htslib/kfunc.jl")
include("htslib/hfile.jl")
include("htslib/hts_os.jl")
include("htslib/hts_log.jl")
include("htslib/hts.jl")
include("htslib/faidx.jl")
include("htslib/tbx.jl")
include("htslib/regidx.jl")
include("htslib/sam.jl")

# Implementation of hFILE backend using julia's IO type
include("hfile_backend.jl")

end


text:
# Conductance types (conductance as in "g")
abstract type AbstractConductanceSystem <: AbstractTimeDependentSystem end

# Linear Ohmic/Nernst vs non-linear GHK
@enum IVCurvature Linear Rectifying

permeability(x::AbstractConductanceSystem) = getfield(x, :ion)
get_inputs(x::AbstractConductanceSystem) = getfield(x, :inputs)
get_output(x::AbstractConductanceSystem) = getfield(x, :output)
# Abstract types without parametrics
struct ConductanceSystem{S<:AbstractTimeDependentSystem} <: AbstractConductanceSystem
    output::Num # 'g' by default
    ion::IonSpecies # ion permeability
    gate_vars::Vector{<:AbstractGatingVariable}
    inputs::Set{Num} # required inputs
    sys::S
    linearity::IVCurvature
    transmatrix::Union{Matrix, Nothing}
    name::Symbol
end

Base.convert(::Type{ODESystem}, x::ConductanceSystem{ODESystem}) = getfield(x, :sys)

# Forward getters to internal system
MTK.get_systems(x::AbstractConductanceSystem) = get_systems(getfield(x, :sys))
MTK.get_eqs(x::AbstractConductanceSystem) = get_eqs(getfield(x, :sys))
MTK.get_dvs(x::AbstractConductanceSystem) = get_dvs(getfield(x, :sys))
MTK.has_ps(x::AbstractConductanceSystem) = MTK.has_ps(getfield(x, :sys))
MTK.get_ps(x::AbstractConductanceSystem) = get_ps(getfield(x, :sys))
MTK.get_defaults(x::AbstractConductanceSystem) = get_defaults(getfield(x, :sys))
MTK.get_states(x::AbstractConductanceSystem) = get_states(getfield(x, :sys))
MTK.get_ivs(x::AbstractConductanceSystem) = get_ivs(getfield(x, :sys))
MTK.get_iv(x::AbstractConductanceSystem) = get_iv(getfield(x, :sys))
MTK.independent_variables(x::AbstractConductanceSystem) = MTK.independent_variables(getfield(x, :sys))
MTK.get_observed(x::AbstractConductanceSystem) = MTK.get_observed(getfield(x, :sys))

function ConductanceSystem(g::Num, ion::IonSpecies, gate_vars::Vector{GatingVariable};
        gbar::Num, linearity::IVCurvature = Linear, extensions::Vector{ODESystem} = ODESystem[],
                           defaults = Dict(), name::Symbol = Base.gensym("Conductance"))

    eqs = Equation[]
    inputs = Set{Num}()
    gate_var_outputs = Set{Num}()
    embed_defaults = Dict()
    params = Set{Num}()

    isparameter(gbar) && push!(params, gbar)

    for var in gate_vars
        x, x∞, τₓ = output(var), steadystate(var), timeconstant(var)
        push!(gate_var_outputs, x)
        eq = D(x) ~ inv(τₓ)*(x∞ - x)
        get_variables!(inputs, eq)
        push!(eqs, eq)
    end

    for sym in inputs
        isparameter(sym) && push!(params, sym)
        hasdefault(sym) && push!(embed_defaults, sym => getdefault(sym))
    end

    # Remove parameters + generated states
    setdiff!(inputs, params, gate_var_outputs)
    if length(gate_vars) > 0
        push!(eqs, g ~ gbar * prod(hasexponent(x) ? output(x)^exponent(x) : output(x) for x in gate_vars))
    else
        push!(eqs, g ~ gbar)
    end
    sys = ODESystem(eqs, t, union(inputs, g, gate_var_outputs), params;
                    defaults = merge(embed_defaults, defaults), name = name)

    for ext in extensions
        sys = extend(sys, ext)
    end

    return ConductanceSystem(g, ion, gate_vars, inputs, sys, linearity, nothing, name)
end

function IonChannel(ion::IonSpecies,
                    gate_vars::Vector{GatingVariable} = GatingVariable[];
                    max_g::Union{Num, SpecificConductance} = 0mS/cm^2,
                    extensions::Vector{ODESystem} = ODESystem[],
                    name::Symbol = Base.gensym("IonChannel"),
                    linearity::IVCurvature = Linear, defaults = Dict())
    if max_g isa SpecificConductance
        gbar_val = ustrip(Float64, mS/cm^2, max_g)
        @parameters gbar
        push!(defaults, gbar => gbar_val)
    else
        gbar = max_g
        if hasdefault(gbar)
            gbar_val = getdefault(gbar)
            if gbar_val isa SpecificConductance
                gbar_val = ustrip(Float64, mS/cm^2, gbar_val)
                gbar = setdefault(gbar, gbar_val)
            end
        end
    end
    @variables g(t)
    g = setmetadata(g, ConductorUnits, mS/cm^2) # TODO: rework with MTK native unit system
    ConductanceSystem(g, ion, gate_vars;
                      gbar = gbar, name = name, defaults = defaults,
                      extensions = extensions, linearity = linearity)
end

function SynapticChannel(ion::IonSpecies,
                         gate_vars::Vector{GatingVariable} = GatingVariable[],
                         max_s::Union{Num, ElectricalConductance} = 0mS;
                         extensions::Vector{ODESystem} = ODESystem[],
                         name::Symbol = Base.gensym("SynapticChannel"),
                         linearity::IVCurvature = Linear, defaults = Dict())
    if max_s isa Electrical Conductance
        sbar_val = ustrip(Float64, mS, max_s)
        @parameters sbar
        push!(defaults, sbar => sbar_val)
    else
        sbar = max_s
        if hasdefault(sbar)
            sbar_val = getdefault(sbar)
            if sbar_val isa ElectricalConductance
                sbar_val = ustrip(Float64, mS, sbar_val)
                sbar = setdefault(sbar, sbar_val)
            end
        end
    end
    @variables s(t)
    s = setmetadata(s, ConductorUnits, mS) # TODO: rework iwth MTK native unit system
    ConductanceSystem(s, ion, gate_vars;
                      gbar = sbar, name = name, defaults = defaults,
                      extensions = extensions, linearity = linearity)
end

function (cond::AbstractConductanceSystem)(newgbar::Quantity)
    newcond = deepcopy(cond)
    g = get_output(newcond)
    outunits = getmetadata(g, ConductorUnits)
    if dimension(outunits) !== dimension(newgbar)
        @error "Input Dimensions do not match output of ConductanceSystem"
    end
    gbar_val = ustrip(Float64, outunits, newgbar)
    gbar_sym = getproperty(cond, :gbar, namespace=false)
    push!(get_defaults(newcond), gbar_sym => gbar_val)
    return newcond
end

macro ionchannel(chan::Expr, ex::Expr, gbar)
    gbar.args[1] != :gbar && throw("Please use `gbar` to define a channel conductance.")
    name, ion, gates = _make_ionchannel(chan, ex)
    IonChannel(eval(ion), gates; max_g = eval(gbar.args[2]), name = name)
end

macro ionchannel(chan::Expr, ex::Expr)
    name, ion, gates = _make_ionchannel(chan, ex)
    IonChannel(eval(ion), gates; name = name)
end

function _make_ionchannel(chan::Expr, ex::Expr)
    ex = MacroTools.striplines(ex)
    !@capture(chan, name_Symbol{ion_}) && throw("An ion type must be given `name{I<:Ion}`")
    gates = GatingVariable[]
    for gate in ex.args
      push!(gates, eval(gate))
    end
    name, ion, gates
end


text:
import Base.size

immutable NN{n}
    ws::NTuple{n, Matrix{Float64}}   # weight matrices
    bs::NTuple{n, Vector{Float64}}   # bias vectors
    ns::NTuple{n, NodeFunc}          # node function (sigmoid, tanh, etc.)
    ls::NTuple{Int}                  # size vector

    function NN(ws::NTuple{n, Matrix{Float64}},
                bs::NTuple{n, Vector{Float64}},
                ns::NTuple{n, NodeFunc})

        for i in 1:length(ws)-1
            @assert size(ws[i],1) == size(ws[i+1],2) "size of levels $i and $(i+1) do not fit"
        end

        for i in 1:length(ws)
            @assert size(ws[i],1) == size(bs[i],1)   "size of bias and weight matrix do not fit at level $i"
        end

        ls = tuple([ size(w,1) for w in ws ]...)
        new(ws, bs, ns, ls)
    end
end

# NN{n}(ws::NTuple{n, Matrix{FloatingPoint}}, ns::NTuple{n, NodeFunc}) = NN{n}(ws, ns)

size(nn::NN) = nn.ls
inputsize(nn::NN) = size(nn.ws[1],2)
depth{n}(nn::NN{n}) = n

function NN(sz::Vector{Int}, nf=Relu(); seed::Integer=0)
    @assert length(sz) > 1 "Net needs at least 2 layers"
    ws = Matrix{Float64}[]
    bs = Vector{Float64}[]

    (seed != 0) && srand(seed)
    for i in 1:length(sz)-1
        push!(ws, randn(sz[i+1], sz[i]) .* sqrt(6/(sz[i+1]+sz[i])))
        push!(bs, zeros(sz[i+1]))
    end

    if isa(nf, NodeFunc)
        ns = NodeFunc[ nf for i in 1:length(ws)]
    elseif isa(nf, Vector) && length(nf)==1
        ns = NodeFunc[ nf[1] for i in 1:length(ws)]
    elseif isa(nf, Vector) && length(nf)==2
        ns = NodeFunc[ [ nf[1] for i in 1:length(ws)-1] ; nf[2] ]
    elseif isa(nf, Vector) && length(nf)==length(ws)
        ns = NodeFunc[ nf[i] for i in 1:length(ws) ]
    else
        error("net depth and node function length do not match")
    end

    NN{length(ws)}(tuple(ws...), tuple(bs...), tuple(ns...))
end

function show(io::IO, nn::NN)
    print(io, "neural net : ")
    print(io, "layers sizes $( Int[ inputsize(nn) , size(nn)... ]), ")
    print(io, "node functions $( map(typeof, nn.ns) )")
    println(io,"")
end

max_stars_repo_name: fredo-dedup/SimpleNeuralNets

max_issues_repo_name: fredo-dedup/SimpleNeuralNets

max_forks_repo_name: fredo-dedup/SimpleNeuralNets


text:
module TestRandomDraws
using Test
using EcologicalNetworks
using LinearAlgebra

B = UnipartiteProbabilisticNetwork(Matrix{Float64}(I, (10,10)))
C = rand(B)
@test typeof(C) <:UnipartiteNetwork
@test C[1,1]
@test !C[1,2]

B = BipartiteProbabilisticNetwork(Matrix{Float64}(I, (10,10)))
C = rand(B)
@test typeof(C) <:BipartiteNetwork
@test C[1,1]
@test !C[1,2]

B = BipartiteProbabilisticNetwork(Matrix{Float64}(I, (10,10)))
C = rand(B, 10)
@test length(C) == 10
@test eltype(C) <: BipartiteNetwork

B = BipartiteProbabilisticNetwork(Matrix{Float64}(I, (10,10)))
D = rand(B, (3,3))
@test size(D) == (3,3)
@test eltype(D) <: BipartiteNetwork

end


text:
# functions to perform Global Sensitivity Anlaysis

struct GSAResult
  output_names::Vector{Symbol}
  params_names::Vector{Symbol}
  pearson::Matrix{Float64}
  partial::Matrix{Float64}
  standard::Matrix{Float64}
end

parameters(gsar::GSAResult) = gsar.params_names
output(gsar::GSAResult) = gsar.output_names
pearson(gsar::GSAResult) = gsar.pearson
partial(gsar::GSAResult) = gsar.partial
standard(gsar::GSAResult) = gsar.standard

function DataFrames.DataFrame(gsar::GSAResult, coef::Symbol)
  if coef == :pearson
    mat = pearson(gsar)
  elseif coef == :partial
    mat = partial(gsar)
  elseif coef == :standard
    mat = standard(gsar)
  else
    throw("Coeffitient type $coef not supported.")
  end
  df = DataFrame(mat, output(gsar))
  insertcols!(df, 1, :parameter => parameters(gsar))
  return df
end

"""
    gsa(mcr::MCResult, timepoint::Number)

Computes Pearson Correlation Coeffitients, Partial Regression Coeffietients and Standard Regression Coeffitients for
parameters vector and output at a given timepoint

Arguments:

- `mcr`: Monte-Carlo results of type `MCResult`
- `timepoint`: Time to compute coeffitients at
"""
function gsa(mcr::MCResult, timepoint::Number)

  params = parameters(mcr)
  params_mat = vecvec_to_mat(VectorOfArray(parameters(mcr)))'
  outvals_mat = vecvec_to_mat(VectorOfArray([mcr[i](timepoint) for i in 1:length(mcr)]))'

  output_names = observables(mcr)
  params_names = collect(keys(params[1]))

  pearson = _calculate_correlation_matrix(outvals_mat, params_mat)
  partial = _calculate_partial_correlation_coefficients(outvals_mat, params_mat)
  standard = _calculate_standard_regression_coefficients(outvals_mat, params_mat)

  return GSAResult(output_names, params_names, pearson, partial, standard)
end

# the following code was coppied from
# https://github.com/SciML/GlobalSensitivity.jl/blob/master/src/regression_sensitivity.jl

function _calculate_standard_regression_coefficients(X, Y)
  β̂ = X' \ Y'
  srcs = (β̂ .* std(X, dims = 2) ./ std(Y, dims = 2)')
  return Matrix(transpose(srcs))
end

function _calculate_correlation_matrix(X, Y)
  corr = cov(X, Y, dims = 2) ./ (std(X, dims = 2) .* std(Y, dims = 2)')
  return Matrix(transpose(corr))
end

function _calculate_partial_correlation_coefficients(X, Y)
  XY = vcat(X, Y)
  corr = cov(XY, dims = 2) ./ (std(XY, dims = 2) .* std(XY, dims = 2)')
  prec = pinv(corr) # precision matrix
  pcc_XY = -prec ./ sqrt.(diag(prec) .* diag(prec)')
  # return partial correlation matrix relating f: X -> Y model values
  return Matrix(transpose(pcc_XY[axes(X, 1), lastindex(X, 1) .+ axes(Y, 1)]))
end


text:
abstract DPMMInitializer
#abstract type DPMMInitializer end

immutable PriorProposalDPMMInitializer <: DPMMInitializer
    datum::Float64 # the first data point in the data set
    hypers::Hyperparameters
end

function add_first_datum!(init::PriorProposalDPMMInitializer, state::DPMMState)
    add_datum_prior!(state, 1, init.datum)
end

function remove_first_datum!(init::PriorProposalDPMMInitializer, state::DPMMState)
    remove_datum_prior!(state, 1, init.datum)
end

immutable OptimalProposalDPMMInitializer <: DPMMInitializer
    datum::Float64 # the first data point in the data set
    hypers::Hyperparameters
end

function add_first_datum!(init::OptimalProposalDPMMInitializer, state::DPMMState)
    add_datum_conditional!(state, 1, init.datum)
end

function remove_first_datum!(init::OptimalProposalDPMMInitializer, state::DPMMState)
    remove_datum_conditional!(state, 1, init.datum)
end

function SMC.forward(init::DPMMInitializer)
    alpha = rand(init.hypers.alpha_prior)
    nign_params = NIGNParams(init.hypers.nign_params_prior)
	state = DPMMState(alpha, CRPState(), OrderedDict{Int,NIGN}(),
                      nign_params, init.hypers)
    log_weight = add_first_datum!(init, state)
	@assert num_assignments(state.crp) == 1
	@assert length(state.components) == 1
    (state, log_weight)
end

function SMC.backward(init::DPMMInitializer, state::DPMMState)
	@assert num_assignments(state.crp) == 1
	@assert length(state.components) == 1
    # we don't actually want to remove the datum from the state
    # TODO implement a second function that does't remove it
    tmp_state = deepcopy(state)
    log_weight = remove_first_datum!(init, tmp_state)
	@assert num_assignments(tmp_state.crp) == 0
	@assert length(tmp_state.components) == 0
    # TODO it seems strange that we don't need to check the prior probs on parameters?
    log_weight
end

abstract DPMMIncrementer
#abstract type DPMMIncrementer end

immutable PriorProposalDPMMIncrementer <: DPMMIncrementer
    data::Vector{Float64}
	datum_index::Int
    rejuvenation_sweeps::Int
end

function add_datum!(incr::PriorProposalDPMMIncrementer, state::DPMMState)
    add_datum_prior!(state, incr.datum_index, incr.data[incr.datum_index])
end

function remove_datum!(incr::PriorProposalDPMMIncrementer, state::DPMMState)
    remove_datum_prior!(state, incr.datum_index, incr.data[incr.datum_index])
end

immutable OptimalProposalDPMMIncrementer <: DPMMIncrementer
    data::Vector{Float64}
	datum_index::Int
    rejuvenation_sweeps::Int
end

function add_datum!(incr::OptimalProposalDPMMIncrementer, state::DPMMState)
    add_datum_conditional!(state, incr.datum_index, incr.data[incr.datum_index])
end

function remove_datum!(incr::OptimalProposalDPMMIncrementer, state::DPMMState)
    remove_datum_conditional!(state, incr.datum_index, incr.data[incr.datum_index])
end

function SMC.forward(incr::DPMMIncrementer, state::DPMMState)
	# sample the new clustering from the conditional distribution,
    # and obtain the normalizing constant
    datum = incr.data[incr.datum_index]
	@assert num_assignments(state.crp) == incr.datum_index - 1
	new_state = deepcopy(state)

    # do a Rejuvenation sweep over the already incorporated data points
    for sweep=1:incr.rejuvenation_sweeps
        for i=1:incr.datum_index-1
            gibbs_step!(new_state, i, incr.data[i])
            alpha_mh_update!(new_state)
            nign_params_update!(new_state)
            nign_m_update!(new_state)
            nign_r_update!(new_state)
            nign_nu_update!(new_state)
            nign_s_update!(new_state)
        end
        nign_params_update!(new_state)
        nign_m_update!(new_state)
        nign_r_update!(new_state)
        nign_nu_update!(new_state)
        nign_s_update!(new_state)
    end

    log_weight = add_datum!(incr, new_state)
	@assert num_assignments(new_state.crp) == incr.datum_index
	(new_state, log_weight)
end

function SMC.backward(incr::DPMMIncrementer, new_state::DPMMState)
	# just get conditional sampling normalizing constant
    datum = incr.data[incr.datum_index]
	@assert num_assignments(new_state.crp) == incr.datum_index
	state = deepcopy(new_state)
    log_weight = remove_datum!(incr, state)
	@assert num_assignments(state.crp) == incr.datum_index - 1

    # Rejuvenation sweep (in reverse order)
    for sweep=1:incr.rejuvenation_sweeps
        nign_s_update!(state)
        nign_nu_update!(state)
        nign_r_update!(state)
        nign_m_update!(state)
        nign_params_update!(state)
        for i=incr.datum_index-1:1
            nign_s_update!(state)
            nign_nu_update!(state)
            nign_r_update!(state)
            nign_m_update!(state)
            nign_params_update!(state)
            alpha_mh_update!(state)
            gibbs_step!(state, i, incr.data[i])
        end
    end
	(state, log_weight)
end

function make_prior_proposal_smc_scheme(data::Vector{Float64}, hypers::Hyperparameters,
                                         rejuvenation_sweeps::Int, num_particles::Int)
    initializer = PriorProposalDPMMInitializer(data[1], hypers)
    incrementers = Vector{PriorProposalDPMMIncrementer}(length(data)-1)
    for i=2:length(data)
        incrementers[i-1] = PriorProposalDPMMIncrementer(data, i, rejuvenation_sweeps)
    end
    SMCScheme(initializer, incrementers, num_particles)
end

function make_prior_proposal_smc_sampler(data::Vector{Float64}, hypers::Hyperparameters,
                                         rejuvenation_sweeps::Int, num_particles::Int)
    scheme = make_prior_proposal_smc_scheme(data, hypers, rejuvenation_sweeps, num_particles)
    SMCSampler(scheme, log_joint_density)
end

function make_optimal_proposal_smc_scheme(data::Vector{Float64}, hypers::Hyperparameters,
                                           rejuvenation_sweeps::Int, num_particles::Int)
    initializer = OptimalProposalDPMMInitializer(data[1], hypers)
    incrementers = Vector{OptimalProposalDPMMIncrementer}(length(data)-1)
    for i=2:length(data)
        incrementers[i-1] = OptimalProposalDPMMIncrementer(data, i, rejuvenation_sweeps)
    end
    SMCScheme(initializer, incrementers, num_particles)
end

function make_optimal_proposal_smc_sampler(data::Vector{Float64}, hypers::Hyperparameters,
                                           rejuvenation_sweeps::Int, num_particles::Int)
    scheme = make_optimal_proposal_smc_scheme(data, hypers, rejuvenation_sweeps, num_particles)
    SMCSampler(scheme, log_joint_density)
end


text:
include("miptests.jl")

testset = [
    (test1, "JG", "Example"),

    (test2, "GB", "MIP-p"),
    (test_P1_Guilherme, "GB", "MIP-g1"),
    (test_TSPbin7_Guilherme, "GB", "MIP-g1-TSP-bin"),
    (test_TSPmip7_Guilherme, "GB", "MIPg1-TSP-mip"),
    (test_TSPmip20_Guilherme, "GB", "MIP-TSP20-mip"),
    (test_TSPmip25_Guilherme, "GB", "MIP-TSP25-mip"),
    (test_TSPmip30_Guilherme, "GB", "MIP-TSP30-mip"),
    (test_TSPmip40_Guilherme, "GB", "MIP-TSP40-mip"),
    (test_TSPmip50_Guilherme, "GB", "MIP-TSP50-mip"),
    (test_PL_Infeasible_Guilherme, "GB", "MIP-inf-p"),
    (test_MIP_medio_Guilherme, "GB", "MIP-m"),
    (test_MIP_Grande_Guilherme, "GB", "MIP-g"),
    (test_PL_Unbounded_Guilherme, "GB", "LP-unb"),
    (test_MIP_Unbounded_Guilherme, "GB", "MIP-unb"),
    (test_MIP_Infeasible_Minimal_Guilherme, "GB", "MIP-inf-m"),
    (test_MIP_Infeasible_Pequeno_Guilherme, "GB", "MIP-inf-m"),
    (test_PL_Guilherme, "GB", "LP-p"),

    (test_P1_Brito, "EB", "MIP-g1"),
    (test_PL_Simples_Brito, "EB", "LP-p"),
    (test_PL_Infeasible_Brito, "EB", "LP-inf-p"),
    (test_PL_Unbounded_Brito, "EB", "LP-unb"),
    (test_MIP_Minimal_Brito, "EB", "MIP-pp"),
    (test_MIP_Pequeno_Brito, "EB", "MIP-p"),

    (testSudoku, "RS", "MIP-m"),
    (testInfeasibleKnapsack, "RS", "MIP-inf-p"),
    (testRobustCCUC, "RS", "MIP-g1"),
    (testUnboundedKnapsack, "RS", "MIP-unb"),
    (testInfeasibleUC, "RS", "MIP-inf"),
    (test_PL_Simples_Raphael, "RS", "LP-opt"),
    (test_PL_Infeasible_Raphael, "RS", "LP-inf"),
    (test_Minimal_UC, "RS", "MIP-pp"),
    (testSudoku4x4, "RS", "MPI-p"),

    (testCaminho, "CG", "MIP-p"),
    (test_optimal_dispatch, "CG", "MIP-g"),

    (test3, "AR", "MIP-p"),
    (test3_2, "AR", "MIP-unb"),
    (test3_3, "AR", "MPI-inf"),
    (test_feature_selection_pequeno_viavel, "AR", "MIP-p"),
    (test_feature_selection_medio, "AR", "MIP-m"),
    (test_feature_selection_grande, "AR", "MIP-g"),
    (test_feature_selection_pequeno_inviavel, "AR", "MIP-inf-m"),
    (teste_PL_andrew_unbounded, "AR", "LP-unb"),
    (teste_PL_andrew_viavel, "AR", "LP-opt"),
    (teste_PL_andrew_inviavel, "AR", "LP-inf"),

    (test_P1_Andrew_Bianca_viavel, "AR", "MIP-g1"),

    (test_rv_1, "RV", "MIP-pp"),
    (test_rv_3, "RV", "MIP-p"),
    (test_rv_5, "RV", "MIP-m"),
    (test_rv_2, "RV", "MIP-inf-p"),
    (test_rv_4, "RV", "LP-inf"),
    (test_rv_7, "RV", "LP-opt"),
    (test_rv_8, "RV", "LP-inf"),
    (test_rv_p1, "RV", "MIP-p1"),

]

function runtests(solveMIP::Function, solver::MathProgBase.AbstractMathProgSolver = JuMP.UnsetSolver();
    author = String[], kind = String[], ignore = String[], name = "")

    table = MIPSolution[]
    @testset "Main" begin
        for teste in testset
            if !("$(teste[1])" in ignore) && (isempty(author) || teste[2] in author) && (isempty(kind) || teste[3] in kind)
                println("$(teste[1])")
                line = teste[1](solveMIP, solver)
                line.name = "$(teste[1]) - $(teste[2])"
                push!(table, line)
            end
        end
        mat = table2mat(table)
        if name == ""
            name = "result_$(solveMIP)_$(solver)"
        end
        writecsv(name*".out",mat)
    end

    return table
end

function table2mat(table)
    out = Any["Name" "Pass" "Objective" "BestBound" "Time" "Nodes" "IntSols" "Status"]
    for line in table
        out = vcat(out, Any[line.name line.pass line.objective line.bestbound line.time line.nodes line.intsols line.status])
    end
    return out
end

function jumpsolve(m)
    tic()
    sol = solve(m)
    time = toq()
    m.ext[:status] = sol
    m.ext[:time] = time
end


text:
using OOPMacro
import OOPMacro: findFnCall, findFnSelfArgNameSymbol, setFnSelf!, deleteFnSelf!, setFnName!, getFnName, setFnSelfArgType!

include("fnUtil.jl")
include("basic.jl")
include("constructor.jl")
include("inheritence.jl")


text:
module JetPackDSP

using DSP, FFTW, Jets

include("jop_convolve.jl")
include("jop_envelope.jl")
include("jop_filter.jl")

end