text:
export samplecdf,normalizedcumsum!

##bisection inverse

bisectioninv(f::Fun{S,T},x::Real;opts...) where {S,T} = first(bisectioninv(f,[x];opts...))

# gives a generalization of midpoint for when `a` or `b` is infinite
function genmidpoint(a::T,b::T) where T
    if isinf(a) && isinf(b)
        zero(T)
    elseif isinf(a)
        b - 100
    elseif isinf(b)
        a + 100
    else
        (a+b)/2
    end
end

function bisectioninv(f::Fun{S,T},x::Float64;numits::Int=47) where {S,T}
    d=domain(f)
    a = minimum(d);b = maximum(d)

    for k=1:numits  #TODO: decide 47
        m= genmidpoint(a,b)
        val = f(m)

            (val<= x) ? (a = m) : (b = m)
    end
    .5*(a+b)
end

bisectioninv(f::Fun{S,T},x::AbstractVector;opts...) where {S,T} = Float64[bisectioninv(f,xx;opts...) for xx in x]

## Clenshaw bisection

function chebbisectioninv(c::AbstractVector{Float64},x::Float64;numits::Int=47)
    a = -1.;b = 1.

    C=Chebyshev()
    for k=1:numits  #TODO: decide 47
        m=.5*(a+b)
        val = clenshaw(C,c,m)

            (val<= x) ? (a = m) : (b = m)
    end
    .5*(a+b)
end

chebbisectioninv(c::AbstractVector{Float64},xl::AbstractVector{Float64}) =
    (n=length(xl);chebbisectioninv(c,xl,ClenshawPlan(Float64,Chebyshev(),length(c),n)))
function chebbisectioninv(c::AbstractVector{Float64},xl::AbstractVector{Float64},plan::ClenshawPlan{Chebyshev{D,R},Float64}) where {D<:Domain,R}
    n = length(xl)
    a = -fill(1.0,n)
    b = fill(1.0,n)

    for k=1:47  #TODO: decide 47
        m=.5*(a+b);
        vals = clenshaw(c,m,plan)

        for j = 1:n
            (vals[j] <= xl[j]) ? (a[j] = m[j]) : (b[j] = m[j])
        end
    end
    m=.5*(a+b)
end

#here, xl is vector w/ length == #cols of c
chebbisectioninv(c::AbstractMatrix{Float64},xl::AbstractVector{Float64}) =
    (n=length(xl);chebbisectioninv(c,xl,ClenshawPlan(Float64,Chebyshev(),size(c,1),n)))
function chebbisectioninv(c::AbstractMatrix{Float64},xl::AbstractVector{Float64},plan::ClenshawPlan{Chebyshev{D,R},Float64}) where {D<:Domain,R}
    @assert size(c)[2] == length(xl)

    n = length(xl)
    a = -fill(1.0,n)
    b = fill(1.0,n)

    for k=1:47  #TODO: decide 47
        m=.5*(a+b);
        vals = clenshaw(c,m,plan)

        for j = 1:n
            (vals[j] <= xl[j]) ? (a[j] = m[j]) : (b[j] = m[j])
        end
    end
    m=.5*(a+b)
end

for TYP in (:Vector,:Float64)
    @eval begin
        bisectioninv(cf::Fun{SP,Float64},x::$TYP;opts...) where {SP<:Chebyshev} =
            fromcanonical.(Ref(space(cf)),chebbisectioninv(coefficients(cf),x;opts...))
#        bisectioninv{SP<:LineSpace}(cf::Fun{SP,Float64},x::$TYP;opts...)=fromcanonical(cf,chebbisectioninv(coefficients(cf),x;opts...))
    end
end

##normalized cumsum

function normalizedcumsum(f::Fun)
    cf = cumsum(f)
    cf = cf/last(cf)

    cf
end

function subtract_zeroatleft!(f::AbstractVector{Float64})
    for k=2:length(f)
        @inbounds f[1] += (-1.)^k.*f[k]
    end

    f
end

function subtract_zeroatleft!(f::AbstractMatrix{Float64})
    for k=2:size(f)[1],j=1:size(f)[2]
        @inbounds f[1,j] += (-1.)^k.*f[k,j]
    end

    f
end

function multiply_oneatright!(f::AbstractVector{Float64})
    val=0.0
    for k=1:length(f)
        val += f[k]
    end

    val = 1/val

    for k=1:length(f)
        @inbounds f[k] *= val
    end

    f
end

function multiply_oneatright!(f::AbstractMatrix{Float64})

    for j=1:size(f)[2]
        val=0.0
        for k=1:size(f)[1]
            val += f[k,j]
        end

        val = 1/val

        for k=1:size(f)[1]
            @inbounds f[k,j] *= val
        end
    end

    f
end

function chebnormalizedcumsum!(f)
    ultraconversion!(f)
    ultraint!(f)
    subtract_zeroatleft!(f)
    multiply_oneatright!(f)
end

# For RandomMatrices compatibility
normalizedcumsum!(f::AbstractVector{Float64})=chebnormalizedcumsum!(f)

## Sampling

sample(f::Fun,n::Integer) = samplecdf(normalizedcumsum(f),n)

samplecdf(cf::Fun,n::Integer) = bisectioninv(cf,rand(n))

sample(f::Fun) = sample(f,1)[1]
samplecdf(f::Fun) = samplecdf(f,1)[1]
samplecdf(v::AbstractVector) = chebbisectioninv(v,rand())

##2D sample

sample(f::Fun{TS},k::Integer) where {TS<:AbstractProductSpace}  =sample(ProductFun(f),k)

function sample(f::LowRankFun,n::Integer)
    rx=sample(sum(f,2),n)
    fA=evaluate.(f,rx,:)
    ry=map(sample,fA)
    [rx ry]
end

function sample(f::LowRankFun{C,C,TensorSpace{Tuple{C,C},DD,RR},Float64},n::Integer) where {C<:Chebyshev,DD<:Domain2d,RR<:Real}
    ry=sample(sum(f,1),n)
    fA=evaluate.(f.A,ry')
    CB=coefficientmatrix(f.B)
    AB=CB*fA
    chebnormalizedcumsum!(AB)
    rx=chebbisectioninv(AB,rand(n))
    [fromcanonical.(Ref(domain(f,1)),rx) ry]
end

sample(f::ProductFun,n)=sample(LowRankFun(f),n)
sample(f::MultivariateFun)=sample(f,1)[1,:]

## Special spaces

# Rays may be schwartz at right endpoint so we project
function sample(f::Fun{JacobiWeight{SS,DD,RR,TT},Float64}, n::Integer) where {SS,DD<:Ray,RR,TT}
    if space(f).α == 0
        samplecdf(normalizedcumsum(f),n)
    else
        sample(default_Fun(f,JacobiWeight(1,space(f).α,domain(f))),n)
    end
end


text:
function erwkpath{T}(g::AbstractGraph{T}, κ::Int, ρ::Int)

    @graph_requires g edge_map incidence_list vertex_list vertex_map

    V = vertices(g)
    N = num_vertices(g)
    weights = ones(num_edges(g))

    for i=1:ρ
        visited_edges = Set{Edge{T}}()
        v = V[rand(1:N)]
        candidate_edges = out_edges(v, g)
        j = 0
        while j < κ && !isempty(candidate_edges)
            e = candidate_edges[rand(1:length(candidate_edges))]
            push!(visited_edges, e)
            weights[edge_index(e, g)] += 1.0
            w = target(e, g)
            E = out_edges(w, g)
            candidate_edges = collect(setdiff(Set(E), visited_edges))
            j += 1
        end
    end
    weights /= ρ
end

function erwkpath1{T}(g::AbstractGraph{T}, κ::Int, ρ::Int)

    @graph_requires g edge_map incidence_list vertex_list vertex_map

    V = vertices(g)
    N = num_vertices(g)
    weights = ones(num_edges(g))

    for i=1:ρ
        visited_edges = Set{Int}()
        u = V[rand(1:N)]
        candidate_edges = Int[edge_index(e, g) for e in out_edges(u, g)]
        j = 0
        while j < κ && !isempty(candidate_edges)
            em_idx = candidate_edges[rand(1:length(candidate_edges))]
            v = target(edges(g)[em_idx], g)
            weights[em_idx] += 1
            push!(visited_edges, em_idx)
            E = Int[edge_index(e, g) for e in out_edges(v, g)]
            candidate_edges = collect(setdiff(Set(E), visited_edges))
            j += 1
        end
    end
    weights /= (ρ+1)
end

ext: jl


text:
#================================
    Service
================================#
"""
    Service

Information about services on open ports
"""
Base.@kwdef mutable struct Service <:Marsh
    devicetype::String      = ""
	extrainfo::String       = ""
	highversion::String     = ""
	hostname::String        = ""
	lowversion::String      = ""
	method::String          = ""
	name::String            = ""
	ostype::String          = ""
	product::String         = ""
	proto::String           = ""
	rcpum::String           = ""
	servicefp::String       = ""
	tunnel::String          = ""
	version::String         = ""
	confidence::Int         = 0
	cpe::Vector{String}     = Vector{String}()
end
StructTypes.names(::Type{Service}) = ((:confidence, :conf),)

confidence(service::Service)  = service.confidence
fingerprint(service::Service) = service.servicefp
name(service::Service) = service.name
hostname(service::Service) = hostname(service.hostname)
product(service::Service)  = hostname.product

Base.string(service::Service) = name(service)

"""
    ServiceState

Information on a port's status, e.g. open, closed, reason, etc
"""
Base.@kwdef mutable struct ServiceState <:Leaf
    state::String        = ""
    reason::String       = ""
    reason_ip::String    = ""
    reason_ttl::String   = ""
end
state(ss::ServiceState) = ss.state
Base.string(ss::ServiceState) = state(ss)

#=======================
    END Service
=======================#

"""
    Owner

Represents a port owner
"""
Base.@kwdef mutable struct Owner <:Leaf
    name::String      = ""
end
name(o::Owner) = o.name
Base.string(o::Owner) = name(o)

"""
    Reason

Reason why a port is closed or filtered

Available when `withReason` is used
"""
Base.@kwdef mutable struct Reason <:Marsh
    reason::String  = ""
    count::Int      = 0
    proto::String   = ""
    ports::String   = ""
end

"""
    ExtraPort

Info on closed and filtered ports
"""
Base.@kwdef mutable struct ExtraPort <:Marsh
    state::String           = ""
    count::Int              = 0
    reasons::Vector{Reason} = Vector{Reason}()
end
StructTypes.names(::Type{ExtraPort}) = ((:reasons, :extrareasons),)

"""
    Port

Details on a scanned port
"""
mutable struct Port <: Marsh
    id::Int64
    protocol::String
    owner::Owner
    state::ServiceState
    service::Service
    scripts::Vector{Script}

    function Port()
        return new(
            0, "", Owner(), ServiceState(), Service(), Vector()
        )
    end
end
StructTypes.names(::Type{Port}) = ((:id, :portid), (:state, :state), (:scripts, :script))
id(port::Port)       = port.id
name(port::Port)     = name(port.service)
owner(port::Port)    = port.owner
protocol(port::Port) = port.protocol
scripts(port::Port)  = port.scripts
service(port::Port)  = port.service
state(port::Port)    = port.state


text:
include("../src/interface.jl")
using DifferentialEquations
using LinearAlgebra
using Plots
using Random

function getmindif(tab::Vector{Vector{BigFloat}})
    nmmin = Inf
    ret = (0, 0)
    for i = 1:size(tab, 1)
        for j = (i+1):size(tab, 1)
            nm = norm(tab[i] - tab[j], Inf)
            if nm < nmmin
                nmmin = nm
                ret = i, j
            end
        end
    end
    return ret, nmmin
end
# henon_heiles
function fct(u, p, t)
    return [0, u[4], -2u[1] * u[2], -u[2] - u[1]^2 + u[2]^2]
end

function fctmain(n_tau, prec)
    Random.seed!(19999909)
    setprecision(512)
    u0 = 2rand(BigFloat, 4) - ones(BigFloat, 4)
    setprecision(prec)
    u0 = BigFloat.(u0)
    println("u0=$u0")
    t_max = big"1.0"
    epsilon = big"1e-7"
    println("epsilon=$epsilon")
    nbmaxtest = 15
    ordmax = 15
    debord = 3
    pasord = 1
    y = ones(Float64, nbmaxtest, div(ordmax - debord, pasord) + 1)
    res_gen = Array{Array{BigFloat,1},2}(undef, nbmaxtest, div(ordmax - debord, pasord) + 1)
    x = zeros(Float64, nbmaxtest)
    A = [0 0 1 0; 0 0 0 0; -1 0 0 0; 0 0 0 0]
    t_0 = big"0.0"
    t_max = big"1.0"
    prob = HiOscODEProblem(fct, u0, (t_0, t_max), missing, A, epsilon)
    nm = NaN
    ordmax += 1
    ordprep = undef
    nb = 10 * 2^(nbmaxtest)
    solref = undef
    while isnan(nm)
        ordmax -= 1
        @time sol = solve(
            prob,
            nb_t = nb,
            order = ordmax,
            getprecision = false,
            nb_tau = n_tau,
            dense = false,
        )
        solref = sol[end]
        nm = norm(solref, Inf)
    end

    tabsol = Array{Array{BigFloat,1},1}(undef, 1)

    tabsol[1] = solref

    println("ordmax=$ordmax")

    solref_gen = solref
    indref = 1

    ind = 1
    for order = debord:pasord:ordmax
        ordprep = order + 2
        println("eps=$epsilon solRef=$solref order=$order")
        nb = 10
        indc = 1
        labels = Array{String,2}(undef, 1, order - debord + 1)
        sol = undef
        par_u0 = missing
        println("preparation ordre $order + 2")
        while indc <= nbmaxtest
            @time solall = solve(
                prob,
                nb_t = nb,
                order = order,
                getprecision = false,
                nb_tau = n_tau,
                par_u0 = par_u0,
                dense = false,
            )
            par_u0 = solall.par_u0
            sol = solall[end]
            push!(tabsol, sol)
            res_gen[indc, ind] = sol
            (a, b), nm = getmindif(tabsol)
            if a != indref
                println("New solref !!!! a=$a, b=$b nm=$nm")
                indref = a
                solref = tabsol[a]
                for i = 1:ind
                    borne = (i < ind) ? size(res_gen, 1) : indc
                    for j = 1:borne
                        nm2 = min(norm(res_gen[j, i] - solref, Inf), 1.1)
                        y[j, i] = nm2 == 0 ? nm : nm2
                    end
                end
            else
                diff = solref - sol
                y[indc, ind] = min(norm(diff, Inf), 1.1)
            end
            x[indc] = 1.0 / nb
            println("result=$y")
            println("log2(y)=$(log2.(y))")
            nb *= 2
            indc += 1
        end
        for i = debord:pasord:order
            labels[1, (i-debord)÷pasord+1] = " eps,order=$(convert(Float32,epsilon)),$i "
        end
        p = Plots.plot(
            x,
            view(y, :, 1:ind),
            xlabel = "delta t",
            xaxis = :log,
            ylabel = "error",
            yaxis = :log,
            legend = :bottomright,
            label = labels,
            marker = 2,
        )
        prec_v = precision(BigFloat)
        eps_v = convert(Float32, epsilon)
        Plots.savefig(p, "out/r4_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_henon_heiles.pdf")
        ind += 1
    end
end
fctmain(32, 512)


text:
using PyPlot
ip = (x, y) -> sum([x_i * y_i for (x_i, y_i) in zip(x, y)])
ct = length(filter(x -> x % 2 == 0, [0:99;]))
pairs = [(2, 5), (4, 2), (9, 8), (12, 10)]
ct2 = length(filter(x -> x[1] % 2 == 0 && x[2] % 2 == 0, pairs))
ct3 = sum([(x % 2 == 0) & (y % 2 == 0) for (x, y) in pairs])
p = (x, coeff) -> sum([a_i * x^(i-1) for (i, a_i) in enumerate(coeff)])
ct4 = str -> sum([a == A for (a, A) in zip(str, uppercase(str))])
isin = (seq_a, seq_b) -> length(seq_a) == sum([in(a, seq_b) for a in seq_a])

function linapprox(f, a, b, n, x)
  delta = (b-a)/n
  points = [a + delta*i for i in 0:n]
  fs = map(f, points)
  m = round(Int, div((x-a), delta))

  plot(points, fs)

  return fs[m] + (fs[m+1]-fs[m])/delta * (x - delta*m)
end

println(ip([1,2], [2,1]))
println(ct)
println(ct2)
println(ct3)
println(p(1, [1, 1, 0]))
println(ct4("Hello"))
println(isin([2:10;], [0:10;]))
println(linapprox(x -> x^2, 0, 1, 5, 0.39))


text:
struct RelandscapeAlg <: Algorithm end
const Relandscape = RelandscapeAlg()
using Rotations
using Flux: Params
using ZenUtils: @grab

# Outstanding
# Need to make sure omegas have same dimensions
# Implement optimize

# wow(x, y) = sin(x*100) * sin(y*100) * exp(sqrt(x^2 + y^2) /  π)
# wow2(x, y) = sin(sqrt(x^2 + y^2) * 10π)
# wow3(x, y) = (sin(10π * x) * sin(10π * y))^2

# sin2d(x, y) = cos(2π * x) = cos(2π * y)

# vectorize(ω) = [x.data for x in values(ω)]

# So this is a little janky but it might work
# Optimize in constrained space

"ω in argmax_Ω(f) initialised at ωinit"
function argmax(f, ω::Ω)
  f(ω) # Init
  # @grab f
  # @grab ω
  ωvals = collect(values((ω)))
  # @grab ωvals
  for i = 1:10000
    # println("ω", ωvals)
    # @show i
    gs = Tracker.gradient(() -> f(ω), Params(ωvals))
    # @grab ωvals
    # @grab gs
    for ωvec in ωvals
      Δ = gs[ωvec]
      l = 0.0001
      # @show Δ, Δ * l
      Flux.Tracker.update!(ωvec, Δ * l)
    end
    # @show f(ω)
  end
  @show f(ω)
  @show exp(f(ω))
  ω
end

"Euclidean Distance between `ω1` and `ω2`"
function ωdist(ω1, ω2)
  # @show keys(ω1)
  # @show keys(ω2)
  # @pre keys(ω1) == keys(ω2)
  total = 0.0
  for i in keys(ω1)
    # FIXME: Should I be able to index ω by vector? i..e ω[i]?
    total += sum((ω1.vals[i] .- ω2.vals[i]).^2)
  end
  sqrt(total)
end

const ϵ = 0.0001
@inline kernel(x, y; beta = 1.0) = x^(beta * y)
@inline logkernel(logx, y; beta = 1.0) = logx > -ϵ ? logx + log((1 + y / exp(logx))) : logx

"Sample from `x | y == true` with Single Site Metropolis Hasting"
function Base.rand(x::RandVar,
                   alg::RelandscapeAlg,
                   ΩT::Type{OT};
                   cb = donothing) where {OT <: Ω}
  # Random phase
  unifs = rand(2)

  # Period, adjust per problem
  period = 10

  # Random rotation
  r = rand(RotMatrix{2})

  # To bound betwee [0, 1]
  sawtooth(x) = x - floor(Flux.data(x))
  function wave(ω)
    res = 1.0

    # Each arrow only has one element
    vs = [x[1] for x in values(ω)]
    vs = sawtooth.(r * vs)
    for (i, v) in enumerate(vs)
      res *= sin((v + unifs[i]) * 2 * pi * period * 1/5)
    end
    log((res)^2) # Work in log scale
  end
  ω = ΩT()
  # U(ω) = logerr(indomainₛ(x, ω)) + wave(ω)
  U(ω) = min(logerr(indomainₛ(x, ω)),  wave(ω))

  @grab U
  # @grab x
  @grab wave
  U(ω) # Init
  argmax(U, ω)
end

function Base.rand(x::RandVar,
                   n::Int,
                   alg::RelandscapeAlg,
                   ΩT::Type{OT};
                   cb = donothing) where {OT <: Ω}
  samples = []
  i = 1
  while i < n
    # try
      sample = rand(x, alg, ΩT, cb = cb)
      push!(samples, sample)
      i += 1
      # @assert false
    # catch e
    #   @show e
    #   @assert false
    # end
  end
  # [rand(x, alg, ΩT, cb = cb) for i = 1:n]
  samples
end

## Example


text:
__precompile__()

module PiecewiseLinearOpt

using JuMP
using LinearAlgebra
using Random

export PWLFunction, UnivariatePWLFunction, BivariatePWLFunction, piecewiselinear

include("types.jl")
include("jump.jl")

end # module


text:
# blind deconvolution
module FocusedBlindDecon

using Misfits
using ProgressMeter
using Printf
using Conv
using Optim
using LineSearches
using RecipesBase
using DataFrames
using StatsBase
using IterativeSolvers
using LinearMaps
using TimerOutputs
using CSV
using Random
using DSP
using FFTW
using LinearAlgebra
using Dates
using JLD

global const to = TimerOutput(); # create a timer object

global STF_FLAG=false
global SDP_FLAG=false
global FILT_FLAG=false

struct UseOptim end
struct UseIterativeSolvers end

global optG
global optS

export P_fbd, fbd!, fibd!, fpr!, lsbd!, ibd!

function __init__(;stf=false, filt=false, sdp=false, optg="iterativesolvers",
		  opts=optg)
	global optG, optS, STF_FLAG, SDP_FLAG
	if(isequal(optg,"iterativesolvers"))
		@info "FocusBD using IterativeSolvers.jl for G"
		optG=UseIterativeSolvers()
	else
		@info "FocusBD using Optim.jl for G"
		optG=UseOptim()
	end
	if(isequal(opts,"iterativesolvers"))
		@info "FocusBD using IterativeSolvers.jl for S"
		optS=UseIterativeSolvers()
	else
		@info "FocusBD using Optim.jl for S"
		optS=UseOptim()
	end

	if(stf)
		@info "S will always be positive"
		STF_FLAG=true
	end
	if(sdp)
		@info "positive semidefinite programming for S"
		SDP_FLAG=true
	end
	if(filt)
		FILT_FLAG=true
	end
	return nothing
end

const FBD=FocusedBlindDecon
export FBD

struct Sxparam
	n::Int64
	attrib::Symbol
end

include("inversion.jl")
include("toygreen.jl")
include("types.jl")
include("toeplitz.jl")
include("fpr.jl")
include("bd.jl")
include("common.jl")
include("ibd.jl")
include("fbd.jl")
include("misfits.jl")
include("ddecon.jl")
include("operators.jl")
include("updates.jl")
include("save.jl")
include("misc.jl")
include("plots.jl")
include("getprop.jl")

end # module

#=

function update_func_grad!(pa; goptim=nothing, soptim=nothing, gαvec=nothing, sαvec=nothing)
	# they will be changed in this program, so make a copy
	ssave=copy(pa.optm.cal.s);
	gsave=copy(pa.optm.cal.g);
	dcalsave=copy(pa.optm.cal.d);

	(goptim===nothing) && (goptim=[:ls])
	(gαvec===nothing) && (gαvec=ones(length(goptim)))

	(soptim===nothing) && (soptim=[:ls])
	(sαvec===nothing) && (sαvec=ones(length(soptim)))

	# dfg for optimization functions
	optim_funcg=Vector{Function}(length(goptim))
	optim_gradg=Vector{Function}(length(goptim))
	for iop in 1:length(goptim)
		if (goptim[iop]==:ls)
			optim_funcg[iop]= x->func_grad!(nothing, x,  pa)
			optim_gradg[iop]=(storage, x)->func_grad!(storage, x,  pa)
		elseif(goptim[iop]==:weights)
			optim_funcg[iop]= x -> func_grad_g_weights!(nothing, x, pa)
			optim_gradg[iop]= (storage, x) -> func_grad_g_weights!(storage, x, pa)
		elseif(goptim[iop]==:acorr_weights)
			optim_funcg[iop]= x -> func_grad_g_acorr_weights!(nothing, x, pa)
			optim_gradg[iop]= (storage, x) -> func_grad_g_acorr_weights!(storage, x, pa)
		else
			error("invalid optim_funcg")
		end
	end
	pa.attrib_inv=:g
	# multi-objective framework
	paMOg=ParamMO(noptim=length(goptim), ninv=length(pa.gx.x), αvec=gαvec,
			    		optim_func=optim_funcg,optim_grad=optim_gradg,
					x_init=randn(length(pa.gx.x),10))
	# create dfg for optimization
	pa.gx.func = x -> paMOg.func(x, paMOg)
	pa.gx.grad! = (storage, x) -> paMOg.grad!(storage, x, paMOg)
#	pa.dfg = OnceDifferentiable(x -> paMOg.func(x, paMOg),
#			    (storage, x) -> paMOg.grad!(storage, x, paMOg), )

	# dfs for optimization functions
	optim_funcs=Vector{Function}(length(soptim))
	optim_grads=Vector{Function}(length(soptim))
	for iop in 1:length(soptim)
		if (soptim[iop]==:ls)
			optim_funcs[iop]=x->func_grad!(nothing, x,  pa)
			optim_grads[iop]=(storage, x)->func_grad!(storage, x,  pa)
		else
			error("invalid optim_funcs")
		end
	end

	pa.attrib_inv=:s
	# multi-objective framework
	paMOs=ParamMO(noptim=length(soptim), ninv=length(pa.sx.x), αvec=sαvec,
			    		optim_func=optim_funcs,optim_grad=optim_grads,
					x_init=vcat(ones(1,10),randn(length(pa.sx.x)-1,10)))
#	pa.dfs = OnceDifferentiable(x -> paMOs.func(x, paMOs),
#			    (storage, x) -> paMOs.grad!(storage, x, paMOs))
	pa.sx.func = x -> paMOs.func(x, paMOs)
	pa.sx.grad! =  (storage, x) -> paMOs.grad!(storage, x, paMOs)

	copyto!(pa.optm.cal.s, ssave)
	copyto!(pa.optm.cal.g, gsave)
	copyto!(pa.optm.cal.d,dcalsave)

	return pa
end
=#


text:
# The main Shapley sampling algorithm.
function _shap_sample(explain::DataFrame,
                      reference::DataFrame,
                      n_instances::Int64,
                      n_instances_explain::Int64,
                      n_features::Int64,
                      n_target_features::Int64,
                      target_features::Array{String,1},
                      feature_names::Array{String,1},
                      feature_names_symbol::Array{Symbol,1},
                      sample_size::Int64,
                      parallel::Symbol,
                      seeds::Union{Array, Integer},
                      chunk::Bool,
                      model,
                      predict_function::Function,
                      precision::Union{Integer, Nothing}
                      )

    if any(parallel .== [:none, :features])
        data_sample = Array{DataFrame}(undef, sample_size)
    elseif any(parallel .== [:samples, :both])
        data_sample = Array{DataFrame}(undef, 1)  # Parallel with pmap().
    end

    for i in if any(parallel .== [:none, :features])
        1:sample_size  # Non-parallel loop over Monte Carlo samples.
    elseif any(parallel .== [:samples, :both])
        1  # Parallel with pmap().
    end

        # Shuffle the column indices, keeping all column indices.
        Random.seed!(seeds[i])
        feature_indices_random = Random.randperm(n_features)

        feature_names_random = feature_names[feature_indices_random]

        # Select a reference instance that all instances in explain will be compared to in
        # this Monte Carlo iteration.
        Random.seed!(seeds[i])
        reference_index = rand(1:n_instances)

        # Shuffle the column order for the randomly selected instance.
        reference_instance = reference[[reference_index], feature_indices_random]
        # Make the reference instance have the same row dimension as "explain" for DataFrame concatenation.
        reference_instance = repeat(reference_instance, n_instances_explain)

        # For the instance(s) to be explained, shuffle the columns to match the randomly selected and shuffled instance.
        explain_instances = explain[!, feature_indices_random]

        target_feature_indices_shuffled = map(x -> (1:n_features)[x .== feature_names_random][1], target_features)
        #----------------------------------------------------------------------
        # Inner loop sampling over target features, creating Frankenstein instances.
        data_sample_feature = Array{DataFrame}(undef, n_target_features)

        if any(parallel .== [:none, :samples])  # Single threaded based on user input.

            for j in 1:n_target_features  # Loop over model features in target_features.

                data_sample_feature[j] = _shap_sample_features(explain_instances,
                                                               reference_instance,
                                                               n_instances_explain,
                                                               n_features,
                                                               target_feature_indices_shuffled[j]
                                                               )
            end  # End single-threaded loop over target features.

        elseif any(parallel .== [:features, :both])  # Multi-threaded based on user input

            Base.Threads.@threads for j in 1:n_target_features  # Loop over model features in target_features.

                data_sample_feature[j] = _shap_sample_features(explain_instances,
                                                               reference_instance,
                                                               n_instances_explain,
                                                               n_features,
                                                               target_feature_indices_shuffled[j]
                                                               )
            end  # End multi-threaded loop over target features.
        end  # End inner loop with feature shuffing for each Monte Carlo sample.
        #----------------------------------------------------------------------
        data_sample[i] = vcat(data_sample_feature...)

        # Two Frankenstein instances per explained instance per target feature.
        # The "* 2" multiplier is because each instance has two Frankenstein instances.

        # Re-order columns for the user-defined predict() function
        data_sample[i] = data_sample[i][!, feature_names_symbol]
        data_sample[i].index = repeat(repeat(1:n_instances_explain, outer = 2), outer = n_target_features)
        data_sample[i].feature_group = repeat(repeat(["real_target", "fake_target"], inner = n_instances_explain), outer = n_target_features)
        data_sample[i].feature_name = repeat(target_features, inner = n_instances_explain * 2)
        data_sample[i].sample = repeat(repeat([i], n_instances_explain * 2), outer = n_target_features)

        # This block of code with 'chunk == true' calls predict() once per Monte Carlo
        # sample which can increase speed and reduce memory on large datasets with many instances and/or features.
        if chunk
            # User-defined predict() function.
            data_predicted = predict_function(model, data_sample[i][:, 1:n_features])
            data_sample[i] = hcat(data_sample[i][:, (n_features + 1):size(data_sample[i], 2)], data_predicted, copycols = false)
        end

    end  # End 'i' loop for data_sample.

    return data_sample
end  # End _shap_sample().
#------------------------------------------------------------------------------

# The inner loop function that shuffles features to create Frankenstein instances.
function _shap_sample_features(explain_instances::DataFrame,
                               reference_instance::DataFrame,
                               n_instances_explain::Int64,
                               n_features::Int64,
                               target_feature_index_shuffled::Int
                               )

        # Create the Frankenstein instances: a combination of the instance to be explained with the
        # reference instance to create a new instance that [likely] does not exist in the dataset.

        # These instances have the real target feature and all features to the right of the shuffled
        # target feature index are from the random reference instance.

        # Then, the marginal feature effect, or stochastic Shapley value approximation,
        # is the difference in predicted values between 1 Frankenstein instance
        # that also replaces the target feature from the reference group and 1 Frankenstein
        # instance where the target feature remains unchanged from its value in explain.

        # Only create a Frankenstein instance if the target is not the last feature and there is actually
        # one or more features to the right of the target to replace with the reference.
        if target_feature_index_shuffled < n_features
          explain_instances = hcat(explain_instances[!, 1:target_feature_index_shuffled],
                                   reference_instance[!, (target_feature_index_shuffled + 1):n_features], copycols = false)
        end

        # These instances are otherwise the same as the Frankenstein instance created above with the
        # exception that the target feature is now replaced with the target feature in the random reference
        # instance. The difference in model predictions between these two Frankenstein instances is
        # what gives us the stochastic Shapley value approximation.
        explain_instances_fake_target = copy(explain_instances)
        explain_instances_fake_target[!, target_feature_index_shuffled] = reference_instance[!, target_feature_index_shuffled]
        #------------------------------------------------------------------
        explain_instances = vcat(explain_instances, explain_instances_fake_target)

        return explain_instances
end  # End _shap_sample_features

max_stars_repo_name: nredell/shapML

max_issues_repo_name: nredell/shapML

max_forks_repo_name: nredell/shapML


text:
function launch()

    @static if Sys.isapple()
        # OpenGL 3.2 + GLSL 150
        glsl_version = 150
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
        GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
        GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
    else
        # OpenGL 3.0 + GLSL 130
        glsl_version = 130
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
        # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
        # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
    end

    # setup GLFW error callback
    error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
    GLFW.SetErrorCallback(error_callback)

    # create window
    window = GLFW.CreateWindow(1280, 720, "ElectroDermal Activity Analysis")
    @assert window != C_NULL
    GLFW.MakeContextCurrent(window)
    GLFW.SwapInterval(1)  # enable vsync

    # setup Dear ImGui context
    ctx = CImGui.CreateContext()

    # setup Dear ImGui style
    # CImGui.StyleColorsDark()
    # CImGui.StyleColorsClassic()
    CImGui.StyleColorsLight()

    # load Fonts
    # - If no fonts are loaded, dear imgui will use the default font. You can also load multiple fonts and use `CImGui.PushFont/PopFont` to select them.
    # - `CImGui.AddFontFromFileTTF` will return the `Ptr{ImFont}` so you can store it if you need to select the font among multiple.
    # - If the file cannot be loaded, the function will return C_NULL. Please handle those errors in your application (e.g. use an assertion, or display an error and quit).
    # - The fonts will be rasterized at a given size (w/ oversampling) and stored into a texture when calling `CImGui.Build()`/`GetTexDataAsXXXX()``, which `ImGui_ImplXXXX_NewFrame` below will call.
    # - Read 'fonts/README.txt' for more instructions and details.
    # fonts_dir = joinpath(@__DIR__, "..", "fonts")
    # fonts = CImGui.GetIO().Fonts
    # default_font = CImGui.AddFontDefault(fonts)
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Cousine-Regular.ttf"), 15)
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "DroidSans.ttf"), 16)
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Karla-Regular.ttf"), 10)
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "ProggyTiny.ttf"), 10)
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Roboto-Medium.ttf"), 16)
    # @assert default_font != C_NULL

    # setup Platform/Renderer bindings
    ImGui_ImplGlfw_InitForOpenGL(window, true)
    ImGui_ImplOpenGL3_Init(glsl_version)
    clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]
    Open_files = false
    while !GLFW.WindowShouldClose(window)
        GLFW.PollEvents()
        # start the Dear ImGui frame
        ImGui_ImplOpenGL3_NewFrame()
        ImGui_ImplGlfw_NewFrame()
        CImGui.NewFrame()
        begin
            CImGui.Begin("Menu")
            CImGui.Text("Please select files that you want to Analysis.")
            @c CImGui.Checkbox("Open files", &Open_files)
            if Open_files
                #df = CSV.read("F:\\julia\\CSV\\EDA.csv", header = ["EDA"])
                df = CSV.read("F:\\julia\\CSV\\HR.csv", header = ["HR"])
                df_data =sort(Cfloat.(df[3:end,1]))
                st = df[1,1]
                freq = df[2,1]
                max_value = df_data[end,1]
                #max = Cfloat(max_value)
                #length = size(df,1)-2
                len = size(df,1)-2
                #@show typeof(len)
                start_time, end_time = @cstatic start_time=Cint(1) end_time=Cint(2) begin
                    @c CImGui.SliderInt("Start Time", &start_time, 1,len)
                    CImGui.SameLine(0.0, CImGui.GetStyle().ItemInnerSpacing.x)
                    time1 = Dates.unix2datetime(st+ (start_time-1)/freq)
                    CImGui.Text(string(Time(time1)))
                    @c CImGui.SliderInt("End Time", &end_time, 2,len)
                    CImGui.SameLine(0.0, CImGui.GetStyle().ItemInnerSpacing.x)
                    time2 = Dates.unix2datetime(st+ (end_time-1)/freq)
                    CImGui.Text(string(Time(time2)))
                    if start_time > end_time -1
                        start_time = end_time - Cint(1)
                    elseif start_time < end_time - 2000*freq
                        start_time = end_time - Cint(2000*freq)
                    end
                    #if end_time < start_time + 1
                    #    end_time =  start_time + Cint(1)
                    #elseif end_time > start_time + 7000
                    #    end_time = start_time + Cint(7000)
                    #end

                end
                #@cstatic length = Cint(size(df, 1) - 2)
                #@show typeof(length)
                #@cstatic start_time = Cint(1) end_time = Cint(70) begin
                #    @c CImGui.SliderInt("Start Time", &start_time, 1, length)
                #    @c CImGui.SliderInt("End Time", &end_time, 1, length)
                #end
                #data = Cfloat.(df[3:end, 1])

                data = Cfloat.(df[(start_time+2):(end_time+2),1])
                CImGui.Text("Primitives")
                    sz, thickness, col = @cstatic sz=Cfloat(36.0) thickness=Cfloat(4.0) col=Cfloat[1.0,0.0,0.4,0.2] begin
                        CImGui.ColorEdit4("Color", col)
                    end
                    p = CImGui.GetCursorScreenPos()
                    col32 = CImGui.ColorConvertFloat4ToU32(ImVec4(col...))

                    begin
                        width = 1200
                        height = 200

                        CImGui.PlotLines("EDA Measurements", data, length(data), 0 , "EDA", 0, Cfloat(max_value*1.2), (width,height))
                        draw_list = CImGui.GetWindowDrawList()
                        x::Cfloat = p.x
                        y::Cfloat = p.y
                        spacing = 8.0
                        # Draws (almost transparent) horizontal bars
                        for yₙ in range(y, step = 40, stop = y + height - 40)
                            CImGui.AddRectFilled(draw_list, ImVec2(x, yₙ), ImVec2(x+width, yₙ+20), col32);
                        end
                    end

                    #draw the x axis
                    p = CImGui.GetCursorScreenPos()
                    begin
                        width = 1200
                        col = Cfloat[0.0,0.0,0.0,1.0]
                        col32 = CImGui.ColorConvertFloat4ToU32(ImVec4(col...))
                        draw_list = CImGui.GetWindowDrawList()
                        x = p.x
                        y = p.y
                        #CImGui.display(typeof(p))

                        time = Dates.unix2datetime(st+ (start_time-1)/freq)
                        f = Cfloat(1 / freq)
                        CImGui.AddLine(draw_list, ImVec2(x, y), ImVec2(x+width, y), col32, Cfloat(1));
                        for xₙ in range(x, step = 1200/ (end_time-start_time+1), stop = x + width)
                            hou = hour(time);
                            min = minute(time);
                            sec = second(time);
                            mil = millisecond(time);
                            if sec ==0 && mil ==0
                                CImGui.AddLine(draw_list, ImVec2(xₙ, y), ImVec2(xₙ, y-5), col32, Cfloat(1));
                            end
                            #if mil ==0
                            if sec ==0 && mil ==0
                                #CImGui.AddText(draw_list, ImVec2(xₙ, y), col32, string(string(min), ":", string(sec)));
                                CImGui.AddText(draw_list, ImVec2(xₙ, y), col32, string(string(hou), ":", string(min)));
                            end
                            next_time = time + Dates.Millisecond(f*1000);
                            time = next_time;
                        end
                    end

                    #date = Dates.unix2datetime(ts)

                    #day(date)
                    #month(date)
                    #year(date)
                    #hour(date)
                    #minute(date)
                    #second(date)
                    #millisecond(date)
                #x = range(0, stop = 2pi, length = 80)
                #f1(x) = sin.(x)
                #f2(x) = exp.(-x) .* cos.(2pi*x)
                #y1 = f1(x)
                #y2 = f2(x)

                #animate, _ = @cstatic animate=true arr=Cfloat[0.6, 0.1, 1.0, 0.5, 0.92, 0.1, 0.2] begin
                #    @c CImGui.Checkbox("Animate", &animate)
                #     data = view(df, 1)
                #     data = convert(Matrix, df)
                #     data = transpose(Matrix(df))
                #     data = convert(Array,data)
                #     CImGui.PlotLines("Result", data, length(data))

                    #@show typeof(data)
                    #@show typeof(Cfloat)
                #     convert(Matrix, df)
                    # create a dummy array of contiguous float values to plot
                    # Tip: If your float aren't contiguous but part of a structure, you can pass a pointer to your first float and the sizeof() of your structure in the Stride parameter.
                #    @cstatic values=fill(Cfloat(0),90) values_offset=Cint(0) refresh_time=Cdouble(0) begin
                #        (!animate || refresh_time == 0.0) && (refresh_time = CImGui.GetTime();)

                #        while refresh_time < CImGui.GetTime() # create dummy data at fixed 60 hz rate for the demo
                #            @cstatic phase=Cfloat(0) begin
                #                values[values_offset+1] = cos(phase)
                #                values_offset = (values_offset+1) % length(values)
                #                phase += 0.10*values_offset
                #                refresh_time += 1.0/60.0
                #            end
                #        end
                #        CImGui.PlotLines("Lines", values, length(values), values_offset, "avg 0.0", -1.0, 1.0, (0,80))
                #        CImGui.PlotHistogram("Histogram", arr, length(arr), 0, C_NULL, 0.0, 1.0, (0,80))
                #    end
                #end # @cstatic
                # use functions to generate output
                # FIXME: This is rather awkward because current plot API only pass in indices. We probably want an API passing floats and user provide sample rate/count.

                #CImGui.Text("")
                #Sin(::Ptr{Cvoid}, i::Cint) = Cfloat(sin(i * 0.1))
                #Saw(::Ptr{Cvoid}, i::Cint) = Cfloat((i & 1) != 0 ? 1.0 : -1.0)
                #Sin_ptr = @cfunction($Sin, Cfloat, (Ptr{Cvoid}, Cint))
                #Saw_ptr = @cfunction($Saw, Cfloat, (Ptr{Cvoid}, Cint))
                #@cstatic func_type=Cint(0) display_count=Cint(70) begin
                #    CImGui.Separator()
                #    CImGui.PushItemWidth(100)
                #    @c CImGui.Combo("func", &func_type, "Sin\0Saw\0")
                #    CImGui.PopItemWidth()
                #    CImGui.SameLine()
                #    @c CImGui.SliderInt("Sample count", &display_count, 1, 400)
                #    func = func_type == 0 ? Sin_ptr : Saw_ptr
                #    CImGui.PlotLines("Lines", func, C_NULL, display_count, 0, C_NULL, -1.0, 1.0, (0,80))
                #    CImGui.PlotHistogram("Histogram", func, C_NULL, display_count, 0, C_NULL, -1.0, 1.0, (0,80))
                #    CImGui.Separator()
                #end

                # animate a simple progress bar
                # @cstatic progress=Cfloat(0) progress_dir=Cfloat(1) begin
                    # if animate
                    #    progress += progress_dir * 0.4 * CImGui.GetIO().DeltaTime
                    #    progress ≥ 1.1 && (progress = 1.1; progress_dir *= -1.0;)
                    #    progress ≤ -0.1 && (progress = -0.1; progress_dir *= -1.0;)
                    # end

                    # typically we would use ImVec2(-1.0,0.0) to use all available width, or ImVec2(width,0.0) for a specified width. ImVec2(0.0,0.0) uses ItemWidth.
                    # CImGui.ProgressBar(progress, ImVec2(0.0,0.0))
                    # CImGui.SameLine(0.0, CImGui.GetStyle().ItemInnerSpacing.x)
                    # CImGui.Text("Progress Bar")

                    # progress_saturated = (progress < 0.0) ? 0.0 : (progress > 1.0) ? 1.0 : progress
                    # buf = @sprintf("%d/%d", progress_saturated*1753, 1753)
                    # CImGui.ProgressBar(progress, ImVec2(0,0), buf)
                # end
            end
            CImGui.End()
        end

        # show a simple window that we create ourselves.
        # we use a Begin/End pair to created a named window.
        # @cstatic f=Cfloat(0.0) counter=Cint(0) begin
        #    CImGui.Begin("Hello, world!")  # create a window called "Hello, world!" and append into it.
        #    CImGui.Text("This is some useful text.")  # display some text
        #    @c CImGui.Checkbox("Demo Window", &show_demo_window)  # edit bools storing our window open/close state
        #    @c CImGui.Checkbox("Another Window", &show_another_window)

        #    @c CImGui.SliderFloat("float", &f, 0, 1)  # edit 1 float using a slider from 0 to 1
        #    CImGui.ColorEdit3("clear color", clear_color)  # edit 3 floats representing a color
        #    CImGui.Button("Button") && (counter += 1)

        #    CImGui.SameLine()
        #    CImGui.Text("counter = $counter")
        #    CImGui.Text(@sprintf("Application average %.3f ms/frame (%.1f FPS)", 1000 / CImGui.GetIO().Framerate, CImGui.GetIO().Framerate))

        #    CImGui.End()
        # end

        # show another simple window.
        # if show_another_window
        #    @c CImGui.Begin("Another Window", &show_another_window)  # pass a pointer to our bool variable (the window will have a closing button that will clear the bool when clicked)
        #    CImGui.Text("Hello from another window!")
        #    CImGui.Button("Close Me") && (show_another_window = false;)
        #    CImGui.End()
        # end

        # rendering
        CImGui.Render()
        GLFW.MakeContextCurrent(window)
        display_w, display_h = GLFW.GetFramebufferSize(window)
        glViewport(0, 0, display_w, display_h)
        glClearColor(clear_color...)
        glClear(GL_COLOR_BUFFER_BIT)
        ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

        GLFW.MakeContextCurrent(window)
        GLFW.SwapBuffers(window)
    end

    # cleanup
    ImGui_ImplOpenGL3_Shutdown()
    ImGui_ImplGlfw_Shutdown()
    CImGui.DestroyContext(ctx)

    GLFW.DestroyWindow(window)
end


text:
module TeukEvolution

include("Fields.jl")
include("Io.jl")
include("Radial.jl")
include("Sphere.jl")
include("Id.jl")
include("Evolution.jl")
include("GHP.jl")
include("BackgroundNP.jl")
include("LinearEvolution.jl")

using .Fields: Field, Initialize_Field
import .Io
import .Radial
import .Sphere
import .Id
import .BackgroundNP
using .GHP: GHP_ops, Initialize_GHP_ops
using .Evolution: Evo_lin_f, Initialize_Evo_lin_f, Evolve_lin_f!
using .LinearEvolution: Linear_evolution!

import TOML

#import MPI
#MPI.Init()
#const comm = MPI.COMM_WORLD
#const nm = MPI.Comm_size(comm) ## number of m angular numbers

#const mi = 1 + MPI.Comm_rank(comm) ## m index
#const m_val = params["m_vals"][mi] ## value of m angular number

function launch(paramfile::String)::Nothing
   println("Launching run, params=$paramfile")
   params = TOML.parsefile(paramfile)

   nx = convert(Int64,params["nx"])
   ny = convert(Int64,params["ny"])
   nt = convert(Int64,params["nt"])
   ts = convert(Int64,params["ts"])

   psi_spin    = convert(Int64,params["psi_spin"])
   psi_falloff = convert(Int64,params["psi_falloff"])

   cl  = convert(Float64,params["cl"])
   cfl = convert(Float64,params["cfl"])
   bhs = convert(Float64,params["bhs"])
   bhm = convert(Float64,params["bhm"])

   runtype = convert(String,params["runtype"])

   outdir = params["outdir"]
   ##===================
   ## Derived parameters
   ##===================
   nm   = length(params["m_vals"])
   minr = bhm*(
      1.0 + sqrt(1.0+(bhs/bhm))*sqrt(1.0-(bhs/bhm))
     ) # horizon (uncompactified)
   maxR = (cl^2)/minr
   dr   = maxR/(nx-1.0)
   dt   = min(cfl*dr,6.0/ny^2)

   println("Number of threads: $(Threads.nthreads())")

   println("Setting up output directory")
   if !isdir(outdir)
      mkdir(outdir)
   else
      rm(outdir,recursive=true)
      mkdir(outdir)
   end
   println("Initializing constant fields")
   Rv = Radial.R_vals(nx, dr)
   Yv = Sphere.Y_vals(ny)
   Cv = Sphere.cos_vals(ny)
   Sv = Sphere.sin_vals(ny)
   Mv = params["m_vals"]
   time = 0.0

   ##=================
   ## Dynamical fields
   ##=================
   println("Initializing linear psi4")
   lin_f = Initialize_Field(name="lin_f",spin=psi_spin,boost=psi_spin,falloff=psi_falloff,Mvals=Mv,nx=nx,ny=ny)
   lin_p = Initialize_Field(name="lin_p",spin=psi_spin,boost=psi_spin,falloff=psi_falloff,Mvals=Mv,nx=nx,ny=ny)

   println("Initializing metric reconstruction fields")
   psi3_f = Initialize_Field(name="psi3",spin=-1,boost=-1,falloff=2,Mvals=Mv,nx=nx,ny=ny)
   psi2_f = Initialize_Field(name="psi2",spin= 0,boost= 0,falloff=3,Mvals=Mv,nx=nx,ny=ny)

   lam_f = Initialize_Field(name="lam",spin=-2,boost=-1,falloff=1,Mvals=Mv,nx=nx,ny=ny)
   pi_f  = Initialize_Field(name="pi", spin=-1,boost= 0,falloff=2,Mvals=Mv,nx=nx,ny=ny)

   muhll_f = Initialize_Field(name="muhll",spin= 0,boost=1,falloff=3,Mvals=Mv,nx=nx,ny=ny)
   hlmb_f  = Initialize_Field(name="hlmb" ,spin=-1,boost=1,falloff=2,Mvals=Mv,nx=nx,ny=ny)
   hmbmb_f = Initialize_Field(name="hmbmb",spin=-2,boost=0,falloff=1,Mvals=Mv,nx=nx,ny=ny)

   println("Initializing independent residuals")
   res_bianchi3_f = Initialize_Field(name="res_bianchi3",spin=-2,boost=-1,falloff=2,Mvals=Mv,nx=nx,ny=ny)
   res_bianchi2_f = Initialize_Field(name="res_bianchi2",spin=-1,boost= 0,falloff=2,Mvals=Mv,nx=nx,ny=ny)
   res_hll_f      = Initialize_Field(name="res_hll",     spin= 0,boost= 2,falloff=2,Mvals=Mv,nx=nx,ny=ny)

   println("Initializing 2nd order psi4")
   scd_f = Initialize_Field(name="scd_f",spin=psi_spin,boost=psi_spin,falloff=psi_falloff,Mvals=Mv,nx=nx,ny=ny)
   scd_p = Initialize_Field(name="scd_p",spin=psi_spin,boost=psi_spin,falloff=psi_falloff,Mvals=Mv,nx=nx,ny=ny)

   ##=======================================
   ## Fixed fields (for evolution equations)
   ##=======================================
   println("Initializing psi4 evolution operators")
   evo_psi4 = Initialize_Evo_lin_f(Rvals=Rv,Cvals=Cv,Svals=Sv,Mvals=Mv,bhm=bhm,bhs=bhs,cl=cl,spin=psi_spin)

   println("Initializing GHP operators")
   ghp = Initialize_GHP_ops(Rvals=Rv,Cvals=Cv,Svals=Sv,Mvals=Mv,bhm=bhm,bhs=bhs,cl=cl)

   println("Initializing Background NP operators")
   bkgrd_np = BackgroundNP.NP_0(Rvals=Rv,Yvals=Yv,Cvals=Cv,Svals=Sv,bhm=bhm,bhs=bhs,cl=cl)
   ##=============
   ## Initial data
   ##=============
   println("Initial data")

   if params["id_kind"]=="gaussian"
      for (mi,mv) in enumerate(Mv)
         Id.set_gaussian!(lin_f[mv], lin_p[mv],
            psi_spin,
            mv,
            params["id_l_ang"][mi],
            params["id_ru"][mi],
            params["id_rl"][mi],
            params["id_width"][mi],
            params["id_amp"][mi][1] + params["id_amp"][mi][2]*im,
            cl, Rv, Yv
         )
         Io.save_csv(tc=0,mv=mv,Rv=Rv,Yv=Yv,outdir=outdir,f=lin_f[mv])
         if runtype=="reconstruction"
            Io.save_csv(tc=0,mv=mv,Rv=Rv,Yv=Yv,outdir=outdir,f=lam_f[mv])
         end
      end
   elseif params["id_kind"]=="qnm"
      Id.set_qnm!()
   else
      throw(DomainError(params["id_kind"],"Unsupported `id_kind` in parameter file"))
   end

   ##===================
   ## Time evolution
   ##===================
   println("Beginning evolution")

   for tc=1:nt
      if runtype=="reconstruction"
         for mv in Mv
            if mv>=0
               Linear_evolution!(
                  psi4_f_pm=lin_f[mv],  psi4_f_nm=lin_f[-mv],
                  psi4_p_pm=lin_p[mv],  psi4_p_nm=lin_p[-mv],
                  psi3_pm=psi3_f[mv],   psi3_nm=psi3_f[-mv],
                  psi2_pm=psi2_f[mv],   psi2_nm=psi2_f[-mv],
                  lam_pm=lam_f[mv],     lam_nm=lam_f[-mv],
                  pi_pm=pi_f[mv],       pi_nm=pi_f[-mv],
                  hmbmb_pm=hmbmb_f[mv], hmbmb_nm=hmbmb_f[-mv],
                  hlmb_pm=hlmb_f[mv],   hlmb_nm=hlmb_f[-mv],
                  muhll_pm=muhll_f[mv], muhll_nm=muhll_f[-mv],
                  Evo_pm=evo_psi4[mv],  Evo_nm=evo_psi4[-mv],
                  Op_pm=ghp[mv],        Op_nm=ghp[-mv],
                  NP=bkgrd_np,
                  R=Rv,
                  m_ang=mv,
                  bhm=bhm,
                  cl=cl, dr=dr, dt=dt
               )
            end
         end
         for mv in Mv
            lin_f_n = lin_f[mv].n;   lin_f_np1 = lin_f[mv].np1
            lin_p_n = lin_p[mv].n;   lin_p_np1 = lin_p[mv].np1
            psi3_n  = psi3_f[mv].n;  psi3_np1  = psi3_f[mv].np1
            psi2_n  = psi2_f[mv].n;  psi2_np1  = psi2_f[mv].np1
            lam_n   = lam_f[mv].n;   lam_np1   = lam_f[mv].np1
            pi_n    = pi_f[mv].n;    pi_np1    = pi_f[mv].np1
            hmbmb_n = hmbmb_f[mv].n; hmbmb_np1 = hmbmb_f[mv].np1
            hlmb_n  = hlmb_f[mv].n;  hlmb_np1  = hlmb_f[mv].np1
            muhll_n = muhll_f[mv].n; muhll_np1 = muhll_f[mv].np1

            for j=1:ny
               for i=1:nx
                  lin_f_n[i,j] = lin_f_np1[i,j]
                  lin_p_n[i,j] = lin_p_np1[i,j]
                  psi3_n[i,j]  = psi3_np1[i,j]
                  psi2_n[i,j]  = psi2_np1[i,j]
                  lam_n[i,j]   = lam_np1[i,j]
                  pi_n[i,j]    = pi_np1[i,j]
                  hmbmb_n[i,j] = hmbmb_np1[i,j]
                  hlmb_n[i,j]  = hlmb_np1[i,j]
                  muhll_n[i,j] = muhll_np1[i,j]
               end
            end
         end
      elseif runtype=="linear_field"
         Threads.@threads for mv in Mv
            Evolve_lin_f!(lin_f[mv],lin_p[mv],evo_psi4[mv],dr,dt)

            lin_f_n   = lin_f[mv].n
            lin_p_n   = lin_p[mv].n
            lin_f_np1 = lin_f[mv].np1
            lin_p_np1 = lin_p[mv].np1

            for j=1:ny
               for i=1:nx
                  lin_f_n[i,j] = lin_f_np1[i,j]
                  lin_p_n[i,j] = lin_p_np1[i,j]
               end
            end
         end
      else
         throw(DomainError(runtype,"Unsupported `runtype` in parameter file"))
      end
      if tc%ts==0
         println("time/bhm ", tc*dt/bhm)
         Threads.@threads for mv in Mv
            Io.save_csv(tc=tc,mv=mv,Rv=Rv,Yv=Yv,outdir=outdir,f=lin_f[mv])
            if runtype=="reconstruction"
               Io.save_csv(tc=tc,mv=mv,Rv=Rv,Yv=Yv,outdir=outdir,f=lam_f[mv])
            end
         end
      end
   end
   println("Finished evolution")
   return nothing
end

end


text:
using FinancialToolbox
if !(VERSION.major == 0 && VERSION.minor <= 6)
    using Test, Dates
else
    using Base.Test
end
function yearFractionTester(StartDate, EndDate, MatlabResults, TestToll = 1e-14)
    for i = 0:FinancialToolbox.currMaxImplemented
        @test(abs(MatlabResults[i+1] - yearfrac(StartDate, EndDate, i)) < TestToll)
    end
    return
end

StartDate = Date(1992, 12, 14);
EndDate = Date(1996, 2, 28);
@test(-yearfrac(StartDate, StartDate, 0) == 0)
@test(-yearfrac(StartDate, EndDate, 0) == yearfrac(EndDate, StartDate, 0))
MatlabResults = [3.20821917808219; 3.20555555555556; 3.25277777777778; 3.20821917808219; 3.20555555555556; 3.20555555555556; 3.20555555555556; 3.20821917808219; 3.20821917808219; 3.25277777777778; 3.20821917808219; 3.20555555555556; 3.20765027322409];
yearFractionTester(StartDate, EndDate, MatlabResults);

EndDate2 = Date(1996, 2, 29);
MatlabResults2 = [3.21095890410959; 3.20833333333333; 3.25555555555556; 3.21095890410959; 3.20833333333333; 3.20833333333333; 3.20833333333333; 3.21095890410959; 3.21095890410959; 3.25555555555556; 3.21095890410959; 3.20833333333333; 3.21038251366122];
yearFractionTester(StartDate, EndDate2, MatlabResults2);

EndDate3 = Date(2026, 2, 28);
MatlabResults3 = [33.2301369863014; 33.2055555555556; 33.6916666666667; 33.2301369863014; 33.2055555555556; 33.2055555555556; 33.2055555555556; 33.2082191780822; 33.2301369863014; 33.6916666666667; 33.2301369863014; 33.2055555555556; 33.2080844374580];
TestToll = 1e-13;
yearFractionTester(StartDate, EndDate3, MatlabResults3, TestToll);

StartDate2 = Date(1992, 2, 29);
EndDate4 = Date(1996, 2, 29);
MatlabResults4 = [3.99180327868852; 4; 4.05833333333333; 4.00273972602740; 3.99722222222222; 4; 4; 4; 3.99180327868852; 4.05833333333333; 4.00273972602740; 4; 4]
yearFractionTester(StartDate2, EndDate4, MatlabResults4);

StartDate3 = Date(1992, 3, 30);
EndDate5 = Date(1996, 12, 31);
MatlabResults5 = [4.75890410958904; 4.75000000000000; 4.82500000000000; 4.75890410958904; 4.75000000000000; 4.75000000000000; 4.75000000000000; 4.75616438356164; 4.75890410958904; 4.82500000000000; 4.75890410958904; 4.75000000000000; 4.75409836065569]
yearFractionTester(StartDate3, EndDate5, MatlabResults5);

StartDate4 = Date(1992, 3, 31);
EndDate6 = Date(1996, 12, 30);
MatlabResults6 = [4.75342465753425; 4.75000000000000; 4.81944444444445; 4.75342465753425; 4.75000000000000; 4.75000000000000; 4.75000000000000; 4.75068493150685; 4.75342465753425; 4.81944444444445; 4.75342465753425; 4.75000000000000; 4.74863387978144]
yearFractionTester(StartDate4, EndDate6, MatlabResults6);

@test_throws(ErrorException, yearfrac(StartDate, EndDate3, FinancialToolbox.currMaxImplemented + 1));
@test_throws(ErrorException, yearfrac(StartDate, EndDate3, -1));

@test(daysact(Date(1992, 12, 14), Date(1992, 12, 13)) == -1)

##Excel test
@test(fromExcelNumberToDate(33952) == Date(1992, 12, 14))

println("Test Dates Passed")


text:
export reconstruction_simple, reconstruction_multiEcho, reconstruction_multiCoil, reconstruction_multiCoilMultiEcho, reconstruction_lowRank, RecoParameters

"""
Performs iterative image reconstruction independently for the data of all coils,
contrasts and slices

# Arguments
* `acqData::AcquisitionData`            - AcquisitionData object
* `reconSize::NTuple{2,Int64}`              - size of image to reconstruct
* `reg::Regularization`                 - Regularization to be used
* `sparseTrafo::AbstractLinearOperator` - sparsifying transformation
* `weights::Vector{Vector{ComplexF64}}` - sampling density of the trajectories in acqData
* `solvername::String`                  - name of the solver to use
* (`correctionMap::Array{ComplexF64}`)  - fieldmap for the correction of off-resonance effects
* (`method::String="nfft"`)             - method to use for time-segmentation when correctio field inhomogeneities
* (`normalize::Bool=false`)             - adjust regularization parameter according to the size of k-space data
* (`params::Dict{Symbol,Any}`)          - Dict with additional parameters
"""
function reconstruction_simple( acqData::AcquisitionData
                              , reconSize::NTuple{2,Int64}
                              , reg::Regularization
                              , sparseTrafo::AbstractLinearOperator
                              , weights::Vector{Vector{ComplexF64}}
                              , solvername::String
                              , correctionMap::Array{ComplexF64}=ComplexF64[]
                              , method::String="nfft"
                              , normalize::Bool=false
                              , params::Dict{Symbol,Any}=Dict{Symbol,Any}())

  numContr, numChan, numSl = numContrasts(acqData), numChannels(acqData), numSlices(acqData)

  # set sparse trafo in reg
  reg.params[:sparseTrafo] = sparseTrafo

  # reconstruction
  Ireco = zeros(ComplexF64, prod(reconSize), numSl, numContr, numChan)
  for k = 1:numSl
    F = encodingOps2d_simple(acqData, reconSize, slice=k, correctionMap=correctionMap, method=method)
    for j = 1:numContr
      W = WeightingOp(weights[j])
      for i = 1:numChan
        kdata = kData(acqData,j,i,k).* weights[j]

        reg2 = deepcopy(reg)
        if normalize
          RegularizedLeastSquares.normalize!(reg2, kdata)
        end
        solver = createLinearSolver(solvername, W*F[j]; reg=reg2, params...)

        I = solve(solver, kdata, startVector=get(params,:startVector,ComplexF64[]),
                                 solverInfo=get(params,:solverInfo,nothing))

        if isCircular( trajectory(acqData, j) )
          circularShutter!(reshape(I, reconSize), 1.0)
        end
        Ireco[:,k,j,i] = I
      end
    end
  end

  Ireco = reshape(Ireco, reconSize..., numSl, numContr, numChan)
  return makeAxisArray(Ireco, acqData)
end

"""
Performs a iterative image reconstruction jointly for all contrasts. Different slices and coil images
are reconstructed independently.

# Arguments
* `acqData::AcquisitionData`            - AcquisitionData object
* `reconSize::NTuple{2,Int64}`              - size of image to reconstruct
* `reg::Regularization`                 - Regularization to be used
* `sparseTrafo::AbstractLinearOperator` - sparsifying transformation
* `weights::Vector{Vector{ComplexF64}}` - sampling density of the trajectories in acqData
* `solvername::String`                  - name of the solver to use
* (`correctionMap::Array{ComplexF64})`  - fieldmap for the correction of off-resonance effects
* (`method::String="nfft"`)             - method to use for time-segmentation when correctio field inhomogeneities
* (`normalize::Bool=false`)             - adjust regularization parameter according to the size of k-space data
* (`params::Dict{Symbol,Any}`)          - Dict with additional parameters
"""
function reconstruction_multiEcho(acqData::AcquisitionData
                              , reconSize::NTuple{2,Int64}
                              , reg::Regularization
                              , sparseTrafo::AbstractLinearOperator
                              , weights::Vector{Vector{ComplexF64}}
                              , solvername::String
                              , correctionMap::Array{ComplexF64}=ComplexF64[]
                              , method::String="nfft"
                              , normalize::Bool=false
                              , params::Dict{Symbol,Any}=Dict{Symbol,Any}())

  numContr, numChan, numSl = numContrasts(acqData), numChannels(acqData), numSlices(acqData)

  # set sparse trafo in reg
  reg.params[:sparseTrafo] = diagOp( repeat([sparseTrafo],numContr) )

  W = WeightingOp( vcat(weights...) )

  # reconstruction
  Ireco = zeros(ComplexF64, prod(reconSize)*numContr, numChan, numSl)
  for i = 1:numSl
    F = encodingOps2d_multiEcho(acqData, reconSize, slice=k, correctionMap=correctionMap, method=method)
    for j = 1:numChan
      kdata = multiEchoData(acqData, j, i) .* weights

      reg2 = deepcopy(reg)
      if normalize
        RegularizedLeastSquares.normalize!(reg2, kdata)
      end
      solver = createLinearSolver(solvername, W*F; reg=reg2, params...)

      Ireco[:,j,i] = solve(solver,kdata)
      # TODO circular shutter
    end
  end

  Ireco = reshape(Ireco, reconSize..., numContr, numChan, numSl)
  return makeAxisArray(permutedims(Ireco,[1,2,5,3,4]), acqData)
end

"""
Performs a SENSE-type iterative image reconstruction. Different slices and contrasts images
are reconstructed independently.

# Arguments
* `acqData::AcquisitionData`            - AcquisitionData object
* `reconSize::NTuple{2,Int64}`              - size of image to reconstruct
* `reg::Regularization`                 - Regularization to be used
* `sparseTrafo::AbstractLinearOperator` - sparsifying transformation
* `weights::Vector{Vector{ComplexF64}}` - sampling density of the trajectories in acqData
* `solvername::String`                  - name of the solver to use
* `senseMaps::Array{ComplexF64}`        - coil sensitivities
* (`correctionMap::Array{ComplexF64}`)  - fieldmap for the correction of off-resonance effects
* (`method::String="nfft"`)             - method to use for time-segmentation when correctio field inhomogeneities
* (`normalize::Bool=false`)             - adjust regularization parameter according to the size of k-space data
* (`params::Dict{Symbol,Any}`)          - Dict with additional parameters
"""
function reconstruction_multiCoil(acqData::AcquisitionData
                              , reconSize::NTuple{2,Int64}
                              , reg::Regularization
                              , sparseTrafo::AbstractLinearOperator
                              , weights::Vector{Vector{ComplexF64}}
                              , solvername::String
                              , senseMaps::Array{ComplexF64}
                              , correctionMap::Array{ComplexF64}=ComplexF64[]
                              , method::String="nfft"
                              , normalize::Bool=false
                              , params::Dict{Symbol,Any}=Dict{Symbol,Any}())

  numContr, numChan, numSl = numContrasts(acqData), numChannels(acqData), numSlices(acqData)

  # set sparse trafo in reg
  reg.params[:sparseTrafo] = sparseTrafo

  # solve optimization problem
  Ireco = zeros(ComplexF64, prod(reconSize), numSl, numContr, 1)
  for k = 1:numSl
    E = encodingOps2d_parallel(acqData, reconSize, senseMaps, slice=k, correctionMap=correctionMap, method=method)
    for j = 1:numContr
      W = WeightingOp(weights[j],numChan)
      kdata = multiCoilData(acqData, j, k) .* repeat(weights[j], numChan)

      reg2 = deepcopy(reg)
      if normalize
        RegularizedLeastSquares.normalize!(reg2, kdata)
      end

      solver = createLinearSolver(solvername, W*E[j]; reg=reg2, params...)

      I = solve(solver, kdata)

      if isCircular( trajectory(acqData, j) )
        circularShutter!(reshape(I, reconSize), 1.0)
      end
      Ireco[:,k,j] = I
    end
  end

  Ireco = reshape(Ireco, reconSize..., numSl, numContr, 1)
  return makeAxisArray(Ireco, acqData)
end

"""
Performs a SENSE-type iterative image reconstruction which reconstructs all contrasts jointly.
Different slices are reconstructed independently.

# Arguments
* `acqData::AcquisitionData`            - AcquisitionData object
* `reconSize::NTuple{2,Int64}`              - size of image to reconstruct
* `reg::Regularization`                 - Regularization to be used
* `sparseTrafo::AbstractLinearOperator` - sparsifying transformation
* `weights::Vector{Vector{ComplexF64}}` - sampling density of the trajectories in acqData
* `solvername::String`                  - name of the solver to use
* `senseMaps::Array{ComplexF64}`        - coil sensitivities
* (`correctionMap::Array{ComplexF64}`)  - fieldmap for the correction of off-resonance effects
* (`method::String="nfft"`)             - method to use for time-segmentation when correctio field inhomogeneities
* (`normalize::Bool=false`)             - adjust regularization parameter according to the size of k-space data
* (`params::Dict{Symbol,Any}`)          - Dict with additional parameters
"""
function reconstruction_multiCoilMultiEcho(acqData::AcquisitionData
                              , reconSize::NTuple{2,Int64}
                              , reg::Regularization
                              , sparseTrafo::AbstractLinearOperator
                              , weights::Vector{Vector{ComplexF64}}
                              , solvername::String
                              , senseMaps::Array{ComplexF64}
                              , correctionMap::Array{ComplexF64}=ComplexF64[]
                              , method::String="nfft"
                              , normalize::Bool=false
                              , params::Dict{Symbol,Any}=Dict{Symbol,Any}())

  numContr, numChan, numSl = numContrasts(acqData), numChannels(acqData), numSlices(acqData)

  # set sparse trafo in reg
  reg.params[:sparseTrafo] = diagOp( repeat([sparseTrafo],numContr) )

  W = WeightingOp( vcat(weights)..., numChan )

  Ireco = zeros(ComplexF64, prod(reconSize), numContr, numSl)
  for i = 1:numSl
    E = encodingOp_2d_multiEcho_parallel(acqData, reconSize, senseMaps, slice=k, correctionMap=correctionMap, method=method)

    kdata = multiCoilMultiEchoData(acqData, i) .* repeat(weights, numChan)

    reg2 = deepcopy(reg)
    if normalize
      RegularizedLeastSquares.normalize!(reg2, acqData.kdata)
    end
    solver = createLinearSolver(solvername, W*E; reg=reg2, params...)

    Ireco[:,:,i] = solve(solver, kdata)
  end

  Ireco = reshape( permutedims(Ireco, [1,3,2]), recoParams[:reconSize]..., numSl, numContr )
end

"""
Auxilary struct that holds parameters relevant for image reconstruction

# Fields
* `reconSize::NTuple{2,Int64}`              - size of image to reconstruct
* `weights::Vector{Vector{ComplexF64}}` - sampling density of the trajectories in acqData
* `sparseTrafo::AbstractLinearOperator` - sparsifying transformation
* `reg::Regularization`                 - Regularization to be used
* `normalize::Bool`                     - adjust regularization parameter according to the size of k-space data
* `solvername::String`                  - name of the solver to use
* `senseMaps::Array{ComplexF64}`        - coil sensitivities
* `correctionMap::Array{ComplexF64}`    - fieldmap for the correction of off-resonance effects
* `method::String="nfft"`               - method to use for time-segmentation when correctio field inhomogeneities
"""
mutable struct RecoParameters{N}
  reconSize::NTuple{N,Int64}
  weights::Vector{Vector{ComplexF64}}
  sparseTrafo::AbstractLinearOperator
  reg::Regularization
  normalize::Bool
  solvername::String
  senseMaps::Array{ComplexF64}
  correctionMap::Array{ComplexF64}
  method::String
end

"""
    setupIterativeReco(acqData::AcquisitionData, recoParams::Dict)

builds relevant parameters and operators from the entries in `recoParams`

# relevant parameters
* `reconSize::NTuple{2,Int64}`              - size of image to reconstruct
* `weights::Vector{Vector{ComplexF64}}` - sampling density of the trajectories in acqData
* `sparseTrafo::AbstractLinearOperator` - sparsifying transformation
* `reg::Regularization`                 - Regularization to be used
* `normalize::Bool`                     - adjust regularization parameter according to the size of k-space data
* `solvername::String`                  - name of the solver to use
* `senseMaps::Array{ComplexF64}`        - coil sensitivities
* `correctionMap::Array{ComplexF64}`    - fieldmap for the correction of off-resonance effects
* `method::String="nfft"`               - method to use for time-segmentation when correctio field inhomogeneities

`sparseTrafo` and `reg` can also be speficied using their names in form of a string.
"""
function setupIterativeReco(acqData::AcquisitionData, recoParams::Dict)

  red3d = dims(trajectory(acqData,1))==2 && length(recoParams[:reconSize])==3
  if red3d  # acqData is 3d data converted to 2d
    reconSize = (recoParams[:reconSize][2], recoParams[:reconSize][3])
    recoParams[:reconSize] = reconSize
  else  # regular case
    reconSize = recoParams[:reconSize]
  end

  # density weights
  densityWeighting = get(recoParams,:densityWeighting,true)
  if densityWeighting
    weights = samplingDensity(acqData,reconSize)
  else
    numContr = numContrasts(acqData)
    weights = Array{Vector{ComplexF64}}(undef,numContr)
    for contr=1:numContr
      numNodes = size(acqData.kdata[contr],1)
      weights[contr] = [1.0/sqrt(prod(reconSize)) for node=1:numNodes]
    end
  end

  # sparsifying transform
  if haskey(recoParams,:sparseTrafo) && typeof(recoParams[:sparseTrafo]) != String
    sparseTrafo = recoParams[:sparseTrafo]
  else
    sparseTrafoName = get(recoParams, :sparseTrafo, "nothing")
    sparseTrafo = SparseOp(sparseTrafoName, reconSize; recoParams...)
  end

  # bare regularization (without sparsifying transform)
  regName = get(recoParams, :regularization, "L1")
  λ = get(recoParams,:λ,0.0)
  reg = Regularization(regName, λ; shape=reconSize, recoParams...)

  # normalize regularizer ?
  normalize = get(recoParams, :normalizeReg, false)

  # solvername
  solvername = get(recoParams, :solver, "fista")

  # sensitivity maps
  senseMaps = get(recoParams, :senseMaps, ComplexF64[])
  if red3d && !isempty(senseMaps) # make sure the dimensions match the trajectory dimensions
    senseMaps = permutedims(senseMaps,[2,3,1,4])
  end

  # field map
  cmap = get(recoParams, :correctionMap, ComplexF64[])
  method = get(recoParams, :method, "nfft")

  return RecoParameters(reconSize, weights, sparseTrafo, reg, normalize, solvername, senseMaps, cmap, method)
end


text:
using BenchmarkTools, Base.Cartesian
using Interpolations

const suite = BenchmarkGroup()

suite["bsplines"] = BenchmarkGroup()
for bspline in ["constant", "linear", "quadratic", "cubic"]
    suite["bsplines"][bspline] = BenchmarkGroup()
end

# To evaluate at fractional positions without any "unnecessary"
# overhead, safer to use Base.Cartesian
@generated function sumvalues(itp::AbstractInterpolation{T,N}, inds) where {T,N}
    quote
        @nexprs $N d->inds_d = inds[d]
        s = zero(eltype(itp))
        @inbounds @nloops $N i d->inds_d begin
            s += @ncall($N, itp, i)
        end
        s
    end
end

function sumvalues_indices(itp)
    inds = axes(itp)
    n = Int(round(10^(3/ndims(itp))))
    ntuple(d->collect(range(first(inds[d])+0.001, stop=last(inds[d])-0.001, length=n)), ndims(itp))
end

strip_prefix(str::AbstractString) = replace(str, "Interpolations."=>"")
benchstr(::Type{T}) where {T<:Interpolations.GridType} = strip_prefix(string(T))

benchstr(::Type{Constant}) = "Constant()"
benchstr(::Type{Linear}) = "Linear()"
benchstr(::Type{Quadratic{BC}}, ::Type{GT}) where {BC<:Interpolations.BoundaryCondition,GT<:Interpolations.GridType} =
    string("Quadratic(", strip_prefix(string(BC)), "(", strip_prefix(string(GT)), "()))")
benchstr(::Type{Cubic{BC}}, ::Type{GT}) where {BC<:Interpolations.BoundaryCondition,GT<:Interpolations.GridType} =
    string("Cubic(", strip_prefix(string(BC)), "(", strip_prefix(string(GT)), "()))")

groupstr(::Type{Constant}) = "constant"
groupstr(::Type{Linear}) = "linear"
groupstr(::Type{Quadratic}) = "quadratic"
groupstr(::Type{Cubic}) = "cubic"

for A in (collect(Float64, 1:3),
          reshape(collect(Float64, 1:9), 3, 3),
          reshape(collect(Float64, 1:27), 3, 3, 3))
    # Constant & Linear
    for D in (Constant, Linear)
        gstr = groupstr(D)
        Ac = copy(A)
        idstr = string(ndims(A), "d_", benchstr(D), '_', benchstr(OnGrid))
        suite["bsplines"][gstr][string(idstr, "_construct")] =
            @benchmarkable interpolate($Ac, BSpline($D()))
        itp = interpolate(copy(A), BSpline(D()))
        inds = sumvalues_indices(itp)
        suite["bsplines"][gstr][string(idstr, "_use")] =
            @benchmarkable sumvalues($itp, $inds)
    end
    # Quadratic
    gstr = groupstr(Quadratic)
    for BC in (Flat,Line,Free,Periodic,Reflect,Natural), GT in (OnGrid, OnCell)
        Ac = copy(A)
        idstr = string(ndims(A), "d_", benchstr(Quadratic{BC}, GT))
        suite["bsplines"][gstr][string(idstr, "_construct")] =
            @benchmarkable interpolate($Ac, BSpline(Quadratic($BC($GT()))))
        itp = interpolate(copy(A), BSpline(Quadratic(BC(GT()))))
        inds = sumvalues_indices(itp)
        suite["bsplines"][gstr][string(idstr, "_use")] =
            @benchmarkable sumvalues($itp, $inds)
    end
    for BC in (InPlace,InPlaceQ)
        Ac = copy(A)
        idstr = string(ndims(A), "d_", benchstr(Quadratic{BC}, OnCell))
        suite["bsplines"][gstr][string(idstr, "_construct")] =
            @benchmarkable interpolate!($Ac, BSpline(Quadratic($BC(OnCell()))))
        itp = interpolate!(copy(A), BSpline(Quadratic(BC(OnCell()))))
        inds = sumvalues_indices(itp)
        suite["bsplines"][gstr][string(idstr, "_use")] =
            @benchmarkable sumvalues($itp, $inds)
    end
    # Cubic
    gstr = groupstr(Cubic)
    for BC in (Flat,Line,Free,Periodic), GT in (OnGrid, OnCell)
        Ac = copy(A)
        idstr = string(ndims(A), "d_", benchstr(Cubic{BC}, GT))
        suite["bsplines"][gstr][string(idstr, "_construct")] =
            @benchmarkable interpolate($Ac, BSpline(Cubic($BC($GT()))))
        itp = interpolate(copy(A), BSpline(Cubic(BC(GT()))))
        inds = sumvalues_indices(itp)
        suite["bsplines"][gstr][string(idstr, "_use")] =
            @benchmarkable sumvalues($itp, $inds)
    end
end

paramspath = joinpath(dirname(@__FILE__), "params.json")

if isfile(paramspath)
    loadparams!(suite, BenchmarkTools.load(paramspath)[1], :evals);
else
    @info "Tuning suite (this may take a while)"
    tune!(suite)
    BenchmarkTools.save(paramspath, params(suite));
end

# To run the benchmarks:
# results = run(suite, verbose = true, seconds = 1)
# BenchmarkTools.save(filename, results)


text:
module TestLatinHypercube

using Test
using MLJBase
using MLJTuning
using LatinHypercubeSampling
import Distributions
import Random
using StableRNGs
using ..Models

const Dist = Distributions

rng = StableRNGs.StableRNG(1234)

x1 = rand(rng, 100);
x2 = rand(rng, 100);
x3 = rand(rng, 100)
X = (x1=x1, x2=x2, x3=x3);
y = 2*x1 .+ 5*x2 .- 3*x3 .+ 0.2*rand(rng, 100);

mutable struct DummyModel <: Deterministic
    lambda::Int
    alpha::Int
    kernel::Char
end

mutable struct SuperModel <: Deterministic
    K::Int64
    model1::DummyModel
    model2::DummyModel
end

MLJBase.fit(::DummyModel, verbosity::Int, X, y) = mean(y), nothing, nothing
MLJBase.predict(::DummyModel, fitresult, Xnew) =
    fill(fitresult, nrows(Xnew))
MLJBase.fit(::SuperModel, verbosity::Int, X, y) = mean(y), nothing, nothing
MLJBase.predict(::SuperModel, fitresult, Xnew) =
    fill(fitresult, nrows(Xnew))

dummy_model = DummyModel(1, 9, 'k')
super_model = SuperModel(4, dummy_model, deepcopy(dummy_model))

@testset "Two ranges with scale" begin
    #ok this works
    model = DummyModel(1,1,'k')
    r1 = range(model, :lambda, lower=1, upper=9);
    r2 = range(model, :alpha, lower=0.4, upper=1.0, scale=:log);
    my_latin = LatinHypercube(gens=2, popsize= 120, rng = rng)
    self_tuning_model = TunedModel(model=model,
                                   tuning=my_latin,
                                   resampling=CV(nfolds=6),
                                   range=[r1, r2],
                                   measure=rms);
end

@testset "Range with infinity" begin
    #ok this works
    model = DummyModel(1, 9, 'k')
    r1 = range(model, :lambda, lower=1, upper=9);
    r2 = range(model, :alpha, lower=0, upper=Inf, origin=2,
               unit=3, scale = :log);
    my_latin = LatinHypercube(gens=4, popsize=100,
                              ntour=2, ptour=0.3,
                              interSampleWeight=0.8, ae_power=2,
                              periodic_ae=false, rng=rng)
    self_tuning_model = TunedModel(model=model,
                                   tuning=my_latin,
                                   resampling=CV(nfolds=6),
                                   range=[r1, r2],
                                   measure=rms);
end

@testset "Full features of latin hypercube" begin
    #ok this works
    model = DummyModel(1, 9, 'k')
    supermodel = SuperModel(4, model, deepcopy(model))

    r1 = range(supermodel, :(model1.lambda), lower=1, upper=9);
    r2 = range(supermodel, :K, lower=0.4, upper=1.5);
    my_latin = LatinHypercube(gens=4, popsize=100,
                              ntour=2, ptour=0.3,
                              interSampleWeight=0.5, ae_power=1.7,
                              periodic_ae=true, rng=rng)

    self_tuning_forest_model = TunedModel(model=supermodel,
                                          tuning=my_latin,
                                          resampling=CV(nfolds=6),
                                          range=[r1, r2],
                                          measure=rms);
end

@testset "_transform" begin
    r_cat = range(Char, :cat, values=['a', 'b', 'c'])
    r_float = range(Float64, :float, lower=1, upper =100, scale=:log10)
    r_int = range(Int32, :int, lower=1, upper=100, scale=:log10)
    @test MLJTuning._transform(r_cat, 2) == 'b'
    @test MLJTuning._transform(r_float, 2.0) == 100.0
    @test MLJTuning._transform(r_int, 2.5) == round(Int32, 10^2.5)
end

@testset "setup" begin
    model = DummyModel(1, 9, 'k')
    r1 = range(model, :lambda, lower=1, upper=9)
    r2 = range(model, :alpha, lower=1, upper=100, scale=:log10);
    my_latin = LatinHypercube(gens=4, popsize=100,
                              ntour=2, ptour=0.3,
                              interSampleWeight=0.5, ae_power=1.7,
                              periodic_ae=true, rng=rng)
    state = MLJTuning.setup(my_latin, model, [r1,r2], 11, 1)
    models = state.models
    @test length(models) == 11
    @test all(models) do model
        model.lambda in 1:9 && model.alpha in 1:100
    end
    @test state.fields == [:lambda, :alpha]
    @test state.parameter_scales == [:linear, :log10]
end

@testset "Scale not a symbol" begin
    #ok this works
    model = DummyModel(1,9,'k')
    r1 = range(model, :lambda, lower=1, upper=9,scale = x->10^x)
    r2 = range(model, :alpha, lower=0.4, upper=1.5);
    @test_throws ArgumentError MLJTuning._create_bounds_and_dims_type(2,[r1,r2])
end

@testset "Return value for ranges" begin
    #ok this works
    model = DummyModel(1,9,'k')
    d = 2
    r1 = range(model, :alpha, lower=1., upper=9.);
    r2 = range(model, :lambda, lower=0.4, upper=1.5);
    bounds, dims = MLJTuning._create_bounds_and_dims_type(d,[r1,r2])
    @test all(bounds[1] .≈ (1.0,9.0))
    @test all(bounds[2] .≈ (0.4,1.5))
    @test all(dims .== [LatinHypercubeSampling.Continuous(),
                        LatinHypercubeSampling.Continuous()])

    r3 = range(model, :lambda, lower=1., upper=9,scale =:log);
    r4 = range(model, :alpha, lower=0.4, upper=1.5);
    bounds, dims = MLJTuning._create_bounds_and_dims_type(d,[r3,r4])
    @test all(bounds[1] .≈ (0.0, log(9)))
    @test all(bounds[2] .≈ (0.4,1.5))
    @test all(dims .== [LatinHypercubeSampling.Continuous(),
                        LatinHypercubeSampling.Continuous()])

    r5 = range(model, :kernel, values=collect("abc"))
    r6 = range(model, :lambda, lower=-Inf, upper=+Inf,
               origin = 0., unit = 3.)
    bounds, dims = MLJTuning._create_bounds_and_dims_type(d,[r5,r6])
    @test all(bounds[1] .≈ (1,3))
    @test all(bounds[2] .≈ (-3.0, 3.0))
    @test all(dims .== [LatinHypercubeSampling.Categorical(3,1.0),
                       LatinHypercubeSampling.Continuous()])

    r7 = range(model, :lambda, lower=0., upper=1.0)
    r8 = range(model, :alpha, lower = -Inf, upper = 15.,
               origin = 4., unit = 10., scale =:linear)
    bounds, dims = MLJTuning._create_bounds_and_dims_type(d,[r7,r8])
    @test all(bounds[1] .≈ (0.0, 1.0))
    @test all(bounds[2] .≈ (-5.,15.0))
    @test all(dims .== [LatinHypercubeSampling.Continuous(),
                        LatinHypercubeSampling.Continuous()])

    r9 = range(model, :alpha, lower=-Inf, upper=+Inf,
               origin = 5.0, unit = 1.0, scale = :log2)
    r10 = range(model, :lambda, lower = 5.0, upper = +Inf,
                origin = 10., unit = 10., scale =:linear)
    bounds, dims = MLJTuning._create_bounds_and_dims_type(d,[r9,r10])
    @test all(bounds[1] .≈ (2.0, log(2, 6)))
    @test all(bounds[2] .≈ (5.0, 25.0))
    @test all(dims .== [LatinHypercubeSampling.Continuous(),
                        LatinHypercubeSampling.Continuous()])

end

@testset "integration test 1 - dummy model" begin

    model = DummyModel(1, 9, 'k')
    supermodel = SuperModel(4, model, deepcopy(model))

    r2 = range(supermodel, :K, lower=10, upper=100, scale=:log10)
    r3 = range(supermodel, :(model2.kernel), values = collect("abcdefghijk"))
    tuning = LatinHypercube()
    r = [r2, r3]
    holdout = Holdout(fraction_train=0.8)
    tuned_model = TunedModel(n =5,
                             model=supermodel,
                             tuning=tuning,
                             resampling=holdout,
                             range=r,
                             measure=rms);

    tuned = machine(tuned_model, X, y)

    fit!(tuned, verbosity=0)
    rep = MLJBase.report(tuned)
    models = map(first, rep.history)
    @test length(models) == 5
    @test all(models) do model
        model.K in 1:100 && model.model2.kernel in collect("abcdefghijk")
    end

end

@testset "integration test 2 - compare with grid search in real case" begin

    X, y = @load_boston();
    model = DecisionTreeRegressor()

    r = range(model, :min_samples_split, lower=2, upper=20)

    # latin:
    tuning = LatinHypercube(rng=rng)
    tuned_model = TunedModel(model=model,
                             tuning=tuning,
                             resampling=Holdout(),
                             range=r,
                             measure=rms,
                             n=100);
    tuned = machine(tuned_model, X, y)
    fit!(tuned, verbosity=0)
    best_latin = fitted_params(tuned).best_model.min_samples_split
    best_rms_latin = report(tuned).best_history_entry.measurement

    # grid:
    tuning = Grid(goal=100, rng=rng)
    tuned_model = TunedModel(model=model,
                             tuning=tuning,
                             resampling=Holdout(),
                             range=r,
                             measure=rms,
                             n=100);
    tuned = machine(tuned_model, X, y)
    fit!(tuned, verbosity=0)
    best_grid= fitted_params(tuned).best_model.min_samples_split
    best_rms_grid = report(tuned).best_history_entry.measurement

    @test best_latin == best_grid
    @test best_rms_latin ≈ best_rms_grid

end

end # module

true

max_stars_repo_name: alan-turing-institute/MLJTuning

max_issues_repo_name: alan-turing-institute/MLJTuning

max_forks_repo_name: alan-turing-institute/MLJTuning


text:
"""
    train!(model::AbstractGPModel, X::AbstractMatrix, y::AbstractArray; obsdim = 1, iterations::Int=10,callback=nothing,conv=0)
    train!(model::AbstractGPModel, X::AbstractVector, y::AbstractArray; iterations::Int=20,callback=nothing,conv=0)

Function to train the given GP `model`.

## Keyword Arguments

there are options to change the number of max iterations,
- `iterations::Int` : Number of iterations (not necessarily epochs!)for training
- `callback::Function` : Callback function called at every iteration. Should be of type `function(model,iter) ...  end`
- `conv::Function` : Convergence function to be called every iteration, should return a scalar and take the same arguments as `callback`
"""
function train!(
    m::OnlineSVGP,
    X::AbstractMatrix,
    y::AbstractArray,
    state=nothing;
    iterations::Int=20,
    callback=nothing,
    convergence=nothing,
    obsdim::Int=1,
)
    return train!(
        m,
        KernelFunctions.vec_of_vecs(X; obsdim=obsdim),
        y,
        state;
        iterations=iterations,
        callback=callback,
    )
end

function train!(
    m::OnlineSVGP,
    X::AbstractVector,
    y::AbstractArray,
    state=nothing;
    iterations::Int=20,
    callback::Union{Nothing,Function}=nothing,
    conv::Union{Nothing,Function}=nothing,
)
    iterations > 0 || error("Number of iterations should be positive")
    X, _ = wrap_X(X)
    y = check_data!(y, likelihood(m))

    data = wrap_data(X, y) # Set the data in the model
    if is_stochastic(m)
        0 < batchsize(inference(m)) <= n_sample(data) || error(
            "The size of mini-batch $(batchsize(inference(m))) is incorrect (negative or bigger than number of samples), please set `batchsize` correctly in the inference object",
        )
        set_ρ!(model, n_sample(data) / batchsize(inference))
    else
        set_batchsize!(inference(m), n_sample(data))
    end

    if n_iter(m) == 0 # The first time data is seen, initialize all parameters
        init_online_model(m, X)
    else
        state = save_old_parameters!(m, state)
        updateZs!(m, X)
    end
    state = isnothing(state) ? init_state(m) : state

    # model.evol_conv = [] # Array to check on the evolution of convergence
    local_iter = 1
    conv = Inf

    while true # Loop until one condition is matched
        try # Allow for keyboard interruption without losing the model
            if local_iter == 1 # This is needed for dealing with change of sizes
                if n_iter(m) == 0 # First iteration Zₐ is empty
                    state = compute_kernel_matrices(m, state, X, true)
                else
                    state = compute_old_matrices(m, state, X)
                end
                local_vars = local_updates!(
                    state.local_vars,
                    likelihood(m),
                    y,
                    mean_f(m, state.kernel_matrices),
                    var_f(m, state.kernel_matrices),
                )
                # We compute the updates given current data but with
                # the previous inducing points to get directly
                # the right size for the variational parameters
                ∇E_μs = ∇E_μ(m, y, local_vars)
                ∇E_Σs = ∇E_Σ(m, y, local_vars)
                state = compute_kernel_matrices(m, state, X, true)
                natural_gradient!.(
                    m.f,
                    ∇E_μs,
                    ∇E_Σs,
                    ρ(m),
                    opt(inference(m)),
                    Zviews(m),
                    state.kernel_matrices,
                    state.opt_state,
                )
                state = global_update!(m, state)

            else
                state = update_parameters!(m, state, X, y) # Update all the variational parameters
            end
            set_trained!(m, true)
            if !isnothing(callback)
                callback(m, state, n_iter(m)) # Use a callback method if given by user
            end
            if (n_iter(m) % m.atfrequency == 0) && n_iter(m) >= 3
                state = update_hyperparameters!(m, state, X, y) # Update the hyperparameters
            end
            if verbose(m) > 2 || (verbose(m) > 1 && local_iter % 10 == 0)
                print("Iteration : $(n_iter(m)), ")
                print("ELBO is : $(objective(m, state, y))")
                print("\n")
                println("number of points : $(dim(m[1]))")
            end
            ### Print out informations about the convergence
            local_iter += 1
            m.inference.n_iter += 1
            local_iter <= iterations || break # Check if the number of maximum iterations has been reached
        # (iter < model.nEpochs && conv > model.ϵ) || break; # Check if any condition has been broken
        catch e
            if e isa InterruptException
                @warn "Training interrupted by user at iteration $local_iter"
                break
            else
                rethrow(e)
            end
        end
    end
    if verbose(m) > 0
        println(
            "Training ended after",
            local_iter,
            " iterations. Total number of iterations ",
            n_iter(model),
        )
    end
    state = compute_kernel_matrices(m, state, X, true) #Compute final version of the matrices for prediction
    set_trained!(m, true)
    return m, state
end

# Update all variational parameters of the online sparse
# variational GP Model
function update_parameters!(model::OnlineSVGP, state, X, y)
    state = compute_kernel_matrices(model, state, X) #Recompute the matrices if necessary (always for the stochastic case, or when hyperparameters have been updated)
    state = variational_updates(model, state, y)
    return state
end

function updateZs!(m::OnlineSVGP, x)
    for gp in m.f
        gp.Z = InducingPoints.updateZ(gp.Z, gp.Zalg, x; kernel=kernel(gp))
        gp.post.dim = length(Zview(gp))
    end
    setHPupdated!(inference(m), true)
    return nothing
end

function save_old_parameters!(m::OnlineSVGP, state)
    opt_state =
        map(m.f, state.opt_state, state.kernel_matrices) do gp, opt_state, kernel_matrices
            previous_gp = save_old_gp!(gp, opt_state.previous_gp, kernel_matrices.K)
            merge(opt_state, (; previous_gp))
        end
    return merge(state, (; opt_state))
end

function save_old_gp!(gp::OnlineVarLatent{T}, previous_gp, K) where {T}
    gp.Zₐ = deepcopy(gp.Z)
    gp.Z = InducingPoints.remove_point(Random.GLOBAL_RNG, gp.Z, gp.Zalg, Matrix(K))# Matrix(pr_cov(gp)))
    invDₐ = Symmetric(-2.0 * nat2(gp) - inv(K)) # Compute Σ⁻¹ₐ - K⁻¹ₐ
    prevη₁ = nat1(gp)
    prev𝓛ₐ = (-logdet(cov(gp)) + logdet(K) - dot(mean(gp), nat1(gp))) / 2
    return merge(previous_gp, (; invDₐ, prevη₁, prev𝓛ₐ))
end

function init_online_model(m::OnlineSVGP{T}, x) where {T<:Real}
    m.f = ntuple(n_latent(m)) do i
        init_online_gp!(m.f[i], x)
    end
    return setHPupdated!(inference(m), false)
end

function init_online_gp!(gp::OnlineVarLatent{T}, x, jitt::T=T(jitt)) where {T}
    Z = InducingPoints.initZ(gp.Zalg, x; kernel=kernel(gp))
    k = length(Z)
    post = OnlineVarPosterior{T}(k)
    prior = GPPrior(kernel(gp), pr_mean(gp))
    return OnlineVarLatent(
        prior, post, Z, typeof(Z)([]), gp.Zalg, gp.Zupdated, gp.opt, gp.Zopt
    )
end

function compute_old_matrices(m::OnlineSVGP{T}, state, x) where {T}
    kernel_matrices = if haskey(state, :kernel_matrices)
        state.kernel_matrices
    else
        ntuple(x -> (;), n_latent(m))
    end
    kernel_matrices = map(m.f, kernel_matrices) do gp, k_mat
        return compute_old_matrices(gp, k_mat, x, T(jitt))
    end
    return merge(state, (; kernel_matrices))
end

function compute_old_matrices(gp::OnlineVarLatent, state, X::AbstractVector, jitt::Real)
    K = cholesky(kernelmatrix(kernel(gp), gp.Zₐ) + jitt * I)
    Knm = kernelmatrix(kernel(gp), X, gp.Zₐ)
    κ = Knm / K
    K̃ = kernelmatrix_diag(kernel(gp), X) .+ jitt - diag_ABt(κ, Knm)
    all(>(0), K̃) || error("K̃ has negative values")
    return merge(state, (; K, Knm, κ, K̃))
end