text:
#!/usr/bin/env julia

#=
-----------------------
Preprocessing functions
-----------------------
=#
"""
    resample(self::RheoTimeData, elperiods::Union{Vector{K}, K}; time_boundaries::Union{Nothing, Vector{T}} = nothing)

Resample data with new sample rate(s).

`resample` can downsample or upsample data. If the number of `elperiods` is negative it is going to reduce the number of samples,
vice versa if it is positive. If `time_boundaries` are not specified, resampling is applied to the whole set of data.
If number of elements per period (`elperiods`) is `1` or `-1` it returns the original `RheoTimeData`, whilst `0` is not accepted as a valid
argument for `elperiods`.

The last element may or may not be included with a normal downsampling. By default the last element is forced to be included
but this can be negated by providing the keyword argument `includelastel=false`.
"""
# function resample(self::RheoTimeData, elperiods::Union{Vector{K}, K}; time_boundaries::Union{Nothing, Vector{T}} = nothing, includelast=true) where {K<:Integer,T<:Real}
#
#     @assert (count(iszero, elperiods)==0) "Number of elements cannot be zero"
#     # convert boundaries from times to element indices
#     if isnothing(time_boundaries)
#         boundaries = [1,length(self.t)];
#     else
#         boundaries = closestindices(self.t, time_boundaries)
#     end
#
#     check = RheoTimeDataType(self)
#     if (check == strain_only)
#         (time, epsilon) = fixed_resample(self.t, self.ϵ, boundaries, elperiods; includelastel=includelast)
#         sigma = [];
#     elseif (check == stress_only)
#         (time, sigma) = fixed_resample(self.t, self.σ, boundaries, elperiods; includelastel=includelast)
#         epsilon = [];
#     elseif (check == strain_and_stress)
#         (time, sigma) = fixed_resample(self.t, self.σ, boundaries, elperiods; includelastel=includelast)
#         (time, epsilon) = fixed_resample(self.t, self.ϵ, boundaries, elperiods; includelastel=includelast)
#     end
#
#     log = self.log == nothing ? nothing : [self.log; RheoLogItem( (type=:process, funct=:resample, params=(elperiods = elperiods, ), keywords=(time_boundaries = time_boundaries, ) ),
#                                         (comment="Resample the data",) ) ]
#
#     return RheoTimeData(sigma, epsilon, time, log)
#
# end

# function resample(d::RheoTimeData, t::Vector{T}) where T<:Real
#
#     @assert hastime(d) "Data without time information cannot be resampled."
#     σr = hasstress(d) ? Spline1D(d.t,d.σ)(t) : d.σ
#     ϵr = hasstress(d) ? Spline1D(d.t,d.ϵ)(t) : d.ϵ
#     log = d.log == nothing ? nothing : [d.log; RheoLogItem( (type=:process, funct=:resample, params=(t = t, ), keywords=NamedTuple() ),
#                                     (comment="Resample the data",) ) ]
#     return RheoTimeData(σr, ϵr, t, log)
#
# end
"""
    resample(d::RheoTimeData [, t = array of time value or range, dt = time step required, scale = multiplicator to apply to existing sampling rate] )

Resample data.

Usage:

* 'resample(d)' would keep the number of sampling points the same but interpolate to set a uniform time step.
* 'resample(d, t=-1:0.1:10)' would resample by interpolation to generate a new dataset with time points given by the range or array values passed with keyword 't'.
* 'resample(d, dt=0.1)' would resample by interpolation to generate a new dataset with time step 'dt'.
* 'resample(d, scale=2)' would resample by multiplying the sampling rate by 'scale'. This could down-sample (scale<1) or upsample (scale>1). If timesteps are non uniform, it would interpolate values accordingly.
"""
function resample(d::RheoTimeData; t::Union{Vector{T},R}=RheoFloat[], scale::T1=0, dt::T2=0) where {T<:Real, R <: AbstractRange, T1<:Real, T2<:Real}

    @assert hastime(d) "Data without time information cannot be resampled."
    if length(t)==0
        if scale > 0
            δt = (typeof(scale)<:Union{Rational,Int}) ? 1//scale : 1. /scale
            idxt=Array(1:length(d.t))
            newidxt=Array(1:δt:length(d.t))
            t=Spline1D(idxt,d.t)(newidxt)
            keywords=(scale = scale,)
        else
            mn, mx = extrema(d.t)
            dt==0 ? δt=(mx-mn)/(length(d.t)-1) : δt=dt
            t=Array{RheoFloat}(mn:δt:mx)
            keywords=(dt = dt,)
        end
    else
        keywords=(t = t,)
        if typeof(t) <: AbstractRange
            t = Vector{RheoFloat}(t)
        end
    end
    σr = hasstress(d) ? Spline1D(d.t,d.σ)(t) : d.σ
    ϵr = hasstrain(d) ? Spline1D(d.t,d.ϵ)(t) : d.ϵ
    log = d.log == nothing ? nothing : [d.log; RheoLogItem( (type=:process, funct=:resample, params=NamedTuple(), keywords=keywords ),
                                    (comment="Resample the data",) ) ]
    return RheoTimeData(σr, ϵr, t, log)

end

"""
    indexweight(self::RheoTimeData, elperiods::Union{Vector{K}, K}; time_boundaries::Union{Nothing, Vector{T}} = nothing, includelast=true) where {K<:Integer,T<:Real}

This function returns array indices (i.e. an array of integers) which can be sent to the `modelfit`, `modelstepfit` or `modeldiffeqfit` functions
to provide a weighted fitting whilst maintaining constant sample-rate.

`indexweight` can underweight indices or overweight them. If `elperiods` in a given boundary is negative, every `abs(n)` index will be used, where `n` is the `elperiod`
corresponding to a given boundary. If `n` is positive, then indicies will be duplicated `n` times such that they are given a higher weighting during the fitting procedure.
If `time_boundaries` are not specified, resampling is applied to the whole set of data. If number of elements per period (`elperiods`) is `1` or `-1` it returns the original
indicies for that boundary, whilst `0` is not accepted as a valid argument for `elperiods`.

The last element may or may not be included. By default the last element is forced to be included
but this can be negated by providing the keyword argument `includelast=false`.
"""
function indexweight(self::RheoTimeData, elperiods::Union{Vector{K}, K}; time_boundaries::Union{Nothing, Vector{T}} = nothing, includelast=true) where {K<:Integer,T<:Real}
    # get element-wise boundaries
    if isnothing(time_boundaries)
        boundaries = [1,length(self.t)]
    else
        boundaries = closestindices(self.t, time_boundaries)
    end

    # assert correct function signature
    @assert length(elperiods)==length(boundaries)-1 "Number of different sample periods must be 1 less than boundaries provided"
    @assert (count(iszero, elperiods)==0) "Number of elements cannot be zero"

    # initialise indices array
    indices = Integer[]

    # loop through boundaries
    for i in 1:(length(boundaries)-1)
        # upsampling, starts at each element then intepolates up N times
        if !signbit(elperiods[i])
            for j in boundaries[i]:(boundaries[i+1]-1)
                indices = vcat(indices, j*ones(Integer, elperiods[i]))
            end

        # downsampling, simply takes every N element as in downsample function
        elseif signbit(elperiods[i])
            append!(indices, collect(boundaries[i]:abs(elperiods[i]):(boundaries[i+1]-1)))
        end
    end

    # last element may be missed, so check and add if needed
    lastel = boundaries[end]
    numlastel = count(i->i==lastel, indices)
    countlastel = elperiods[end]
    if includelast && !signbit(countlastel) && numlastel<countlastel
        # up-sample case
        indices = vcat(indices, lastel*ones(Integer, countlastel-numlastel))

    elseif includelast && signbit(countlastel) && indices[end]<lastel
        append!(indices, lastel)
    end

    return indices
end

"""
    cutting(self::RheoTimeData, time_on::T1, time_off::T2) where {T1<:Number, T2<:Number}

Remove the data outside a specified time interval.

By specifing a time interval (`time_on`, `time_off`), a new `RheoTimeData` is returned without the data
lying outside time interval.
"""
function cutting(self::RheoTimeData, time_on::T1, time_off::T2) where {T1<:Number, T2<:Number}

    @assert isreal(time_on) && isreal(time_off) "Boundaries cannot be complex numbers"

    boundary_on = closestindex(self.t, time_on);
    boundary_off = closestindex(self.t, time_off);
    time = self.t[boundary_on:boundary_off]

    check = RheoTimeDataType(self)
    if (check == strain_only)
        ϵ = self.ϵ[boundary_on:boundary_off]
        σ = [];
    elseif (check == stress_only)
        σ = self.σ[boundary_on:boundary_off]
        ϵ = [];
    elseif (check == strain_and_stress)
        ϵ = self.ϵ[boundary_on:boundary_off]
        σ = self.σ[boundary_on:boundary_off]
    end

    log = self.log == nothing ? nothing : [self.log; RheoLogItem( (type=:process, funct=:cutting, params=(time_on = time_on, time_off = time_off), keywords=() ),
                                        (comment="Cut section of the data between $time_on and $time_off",) ) ]

    return RheoTimeData(σ, ϵ, time, log)

end

"""
    smooth(self::RheoTimeData, τ::Real; pad::String="reflect")

Smooth data using a Gaussian Kernel to time scale `τ` (approximately half power).

Smooths both `σ` and `ϵ`. Sampling frequency must be constant as it is based on FFT convolution. Essentially a
low pass filter with frequencies of 1/τ being cut to approximately half power. For other pad types available
see ImageFiltering documentation. As of doc writing, pad options are: `"replicate"` (repeat edge values to
infinity), `"circular"` (image edges "wrap around"), `"symmetric"` (the image reflects relative to a position
between pixels), `"reflect"` (the image reflects relative to the edge itself).
"""
function smooth(self::RheoTimeData, τ::Real; pad::String="reflect")

    @eval import ImageFiltering: imfilter, Kernel

    # get sample-rate and Gaussian kernel (std. dev)
    samplerate = 1.0/getsampleperiod(self.t)
    Σ = getsigma(τ, samplerate)

    # smooth signal and return
    check = RheoTimeDataType(self)
    if (check == strain_only)
        epsilon = Base.invokelatest(imfilter, self.ϵ, Base.invokelatest(Kernel.reflect, Base.invokelatest(Kernel.gaussian, (Σ,))), pad)
        sigma = [];
    elseif (check == stress_only)
        sigma = Base.invokelatest(imfilter, self.σ, Base.invokelatest(Kernel.reflect, Base.invokelatest(Kernel.gaussian, (Σ,))), pad)
        epsilon = [];
    elseif (check == strain_and_stress)
        sigma = Base.invokelatest(imfilter, self.σ, Base.invokelatest(Kernel.reflect, Base.invokelatest(Kernel.gaussian, (Σ,))), pad)
        epsilon = Base.invokelatest(imfilter, self.ϵ, Base.invokelatest(Kernel.reflect, Base.invokelatest(Kernel.gaussian, (Σ,))), pad)
    end

    log = self.log == nothing ? nothing : [self.log; RheoLogItem( (type=:process, funct=:smooth, params=(τ = τ,), keywords=(pad = pad,) ),
                                        (comment="Smooth data with timescale $τ",) ) ]

    RheoTimeData(sigma, epsilon, self.t, log)

end

"""
    extract(self::Union{RheoTimeData,RheoFreqData}, type::Union{TimeDataType,FreqDataType,Integer})

Extract specific fields form `RheoTimeData` or `RheoFreqData`.

Extract can copy one or more fields from a given `RheoXData` variable into a new `RheoXData` one. The fields that
are copied are identified by the specified type of data.
If self is a `RheoTimeData`, the type that can be extracted is `time_only` (or `0`), `stress_only` (or `1`), `strain_only` (or `2`).
Note that `strain_and_stress` (or `3`) is not allowed.
If self is a `RheoFreqData`, the type that can be extracted is `freq_only` (or `0`).
"""
function extract(self::RheoTimeData, type::Union{TimeDataType,Integer})

        type = typeof(type)==Int ? TimeDataType(type) : type
        @assert (typeof(type)==TimeDataType) || (typeof(type)==Int) "Cannot extract frequency data from RheoTimeData"
        @assert (type!= strain_and_stress) && (type!= invalid_time_data) "Cannot extract both stress and strain"
        @assert (type!= invalid_time_data) "Cannot extract information from invalid time data"
        check = RheoTimeDataType(self)

        log = self.log == nothing ? nothing : [self.log; RheoLogItem( (type=:process, funct=:extract, params=(type=type,), keywords=() ),
                                            (comment="Data field extraction, $type from $check",) ) ]

        if type == time_only
            @assert check!= invalid_time_data "Time not available"
            return RheoTimeData([], [], self.t,log)
        elseif type == strain_only
            @assert (check == strain_and_stress) || (check == strain_only) "Strain not available"
            return RheoTimeData([], self.ϵ, self.t,log)
        elseif type == stress_only
            @assert (check == strain_and_stress) || (check == stress_only) "Stress not available"
            return RheoTimeData(self.σ, [], self.t,log)
        end

end

function extract(self::RheoFreqData, type::Union{FreqDataType,Integer})

        type = typeof(type)==Int ? FreqDataType(type) : type
        @assert (type!= with_modulus) "Cannot extract frequency with moduli"
        @assert (type!= invalid_freq_data) "Cannot extract information from invalid frequency data"
        check = RheoFreqDataType(self)
        @assert (check == with_modulus) "Frequency and modulii required"

        log = self.log == nothing ? nothing : [self.log; RheoLogItem( (type=:process, funct=:extract, params=(type=type,), keywords=() ),
                                            (comment="Frequency field extraction",) ) ]
        return RheoFreqData([], [], self.ω,log)
end

#=
--------------------------------
Fitting and predicting functions
--------------------------------
=#
function fill_init_params(model::RheoModelClass, p0::Union{NamedTuple,Nothing})

    springpotnumber = length(findall(x->(x==:a)||(x==:β), model.params))

    if isnothing(p0)
        # fill all initial guesses as 0.5, if two springpots then
        # make sure dissipation coefficients conform to constraints
        p0a = ones(RheoFloat, length(model.params))/2
        if springpotnumber==2
            index = findfirst(model.params.==:a)
            p0a[index] = RheoFloat(0.8)
        end
        @warn "Initial values for model parameters are set to $p0a by default"

    elseif length(p0)<length(model.params)
        # check all provided symbols are part of the model, fill all parameters not provided
        for i in keys(p0)
            @assert (i in model.params) "Incorrect parameter name provided for model initial conditions:" * i
        end
        # propagate list
        p0a = Vector{RheoFloat}(undef, length(model.params))
        for (i,p) in enumerate(model.params)
            p0a[i] = (p in keys(p0)) ? p0[p] : RheoFloat(1/2)
        end
        # check for second springpot specification
        indexa = findfirst(model.params.==:a)
        indexβ = findfirst(model.params.==:β)
        if !(:a in keys(p0)) && !(:β in keys(p0))
            p0a[indexa] = RheoFloat(0.8)
        elseif !(:β in keys(p0))
            p0a[indexβ] = p0[:a]/2
        elseif !(:a in keys(p0))
            p0a[indexa] = minimum((p0[:β]*2, RheoFloat(0.999)))
        end

        @warn "Unspecified initial guesses filled by default values: $p0a"

    else
        p0a = check_and_reorder_parameters(p0, model.params,  err_string="initial guess")
    end

    @assert model.constraint(p0a)  "Initial guess not feasible"

    return p0a
end

function fill_lower_bounds(model::RheoModelClass, lo::Union{NamedTuple,Nothing})
    if isnothing(lo)
        loa = nothing

    elseif length(lo)<length(model.params)
        loa = Vector{RheoFloat}(undef, length(model.params))
        for (i,p) in enumerate(model.params)
            loa[i] = (p in keys(lo)) ? lo[p] : zero(RheoFloat)
        end

    else
        loa = check_and_reorder_parameters(lo, model.params,  err_string="low bounds")
    end

    return loa
end

function fill_upper_bounds(model::RheoModelClass, hi::Union{NamedTuple,Nothing})
    if isnothing(hi)
        hia = nothing

    elseif length(hi)<length(model.params)
        hia = Vector{RheoFloat}(undef, length(model.params))
        for (i,p) in enumerate(model.params)
            hia[i] = (p in keys(hi)) ? hi[p] : RheoFloat(Inf)
        end

    else
        hia = check_and_reorder_parameters(hi, model.params, err_string="high bounds")
    end

    return hia
end

"""
    modelfit(data::RheoTimeData, model::RheoModelClass, modloading::Symbol; p0::Union{NamedTuple,Nothing} = nothing, lo::Union{NamedTuple,Nothing} = nothing, hi::Union{NamedTuple,Nothing} = nothing, verbose::Bool = false, rel_tol = 1e-4, diff_method="BD", weights::Union{Vector{Integer},Nothing}=nothing)

Fit `RheologyData` struct to model and return a fitted model as a `RheologyModel` object.
For the fitting process RHEOS relies on the optimistion package NLopt.jl (https://nlopt.readthedocs.io/en/latest/).
RHEOS makes use of a local derivative free algorithm, specifically the Tom Rowan's "Subplex"

# Arguments

- `data`: `RheoTimeData` struct containing all data
- `model`: `RheoModelClass` containing moduli functions and named tuple parameters
- `modloading`: `strain_imposed` or `1`, `stress_imposed` or `2`
- `p0`: Initial parameters to use in fit (uses 0.5 for all parameters if not defined)
- `lo`: Lower bounds for parameters
- `hi`: Upper bounds for parameters
- `verbose`: If true, prints parameters on each optimisation iteration
- `rel_tol`: Relative tolerance of optimization, see NLOpt docs for more details
- `diff_method`: Set finite difference formula to use for derivative, currently `"BD"` or `"CD"`
- `weights`: Vector of indices weighted according to importance, can be generated by `indexweight` function
"""
function modelfit(data::RheoTimeData,
                  model::RheoModelClass,
                  modloading::LoadingType;
                  p0::Union{NamedTuple,Nothing} = nothing,
                  lo::Union{NamedTuple,Nothing} = nothing,
                  hi::Union{NamedTuple,Nothing} = nothing,
                  verbose::Bool = false,
                  rel_tol = 1e-4,
                  diff_method="BD",
                  weights::Union{Vector{Integer},Nothing} = nothing)

    p0a = fill_init_params(model, p0)
    loa = fill_lower_bounds(model, lo)
    hia = fill_upper_bounds(model, hi)

    rel_tol = convert(RheoFloat,rel_tol)

    check = RheoTimeDataType(data)
    @assert (check == strain_and_stress) "Both stress and strain are required"

    # check provided weights are all valid
    if !isnothing(weights)
        @assert isempty(weights[weights.<1]) "Invalid weighting indices provided"
    end

    # use correct method for derivative
    if diff_method=="BD"
        deriv = derivBD
    elseif diff_method=="CD"
        deriv = derivCD
    end

    # get modulus function and derivative
    if modloading == stress_imposed
        dcontrolled = deriv(data.σ, data.t)
        measured = data.ϵ
        modulus = _Ja(model)
        modsing = (t->model._J(t,p0a))
        modused = "J"
    elseif modloading == strain_imposed
        dcontrolled = deriv(data.ϵ, data.t)
        measured = data.σ
        modulus = _Ga(model)
        modsing = (t->model._G(t,p0a))
        modused ="G"
    end

    # check necessary modulus is defined
    @assert modulusexists(modsing) "Modulus to use not defined"
    # get singularity presence
    sing = singularitytest(modsing)
    # get time step (only needed for convolution, which requires constant dt so t[2]-t[1] is sufficient)
    dt = data.t[2] - data.t[1]

    # time must start at 0 for convolution to work properly!
    t_zeroed = data.t .- minimum(data.t)

    # fit
    is_constant = constantcheck(data.t)

    # indices weighting only used for constant sample-rate data
    if !is_constant && !isnothing(weights)
        @warn "Indices weighting not used as variable sample-rate data has been provided"
    end

    (minf, minx, ret), timetaken, bytes, gctime, memalloc =
                    @timed leastsquares_init(   p0a,
                                                loa,
                                                hia,
                                                modulus,
                                                t_zeroed,
                                                dt,
                                                dcontrolled,
                                                measured;
                                                insight = verbose,
                                                constant_sampling = is_constant,
                                                singularity = sing,
                                                _rel_tol = rel_tol,
                                                indweights = weights)

    println("Time: $timetaken s, Why: $ret, Parameters: $minx, Error: $minf")

    nt = NamedTuple{Tuple(model.params)}(minx)

    if data.log != nothing
        # Preparation of data for log item
        info=(comment="Fiting rheological model to data", model_name=model.name, model_params=nt, time_taken=timetaken, stop_reason=ret, error=minf)
        params=(model=model, modloading=modloading)
        keywords=(p0=p0, lo=lo, hi=hi, rel_tol=rel_tol, diff_method=diff_method)
        # Add data to the log
        push!(data.log, RheoLogItem( (type=:analysis, funct=:modelfit, params=params, keywords=keywords), info))
    end

    return RheoModel(model, nt);

end

"""
    modelpredict(data::RheoTimeData, model::RheoModel; diff_method="BD")

Given an incomplete data set (only either stress or strain missing) and model with values substituted into
parameters (`RheoModel`), return a new dataset based on the model.
If data is type of `stress_only`, then creep modulus (`:J`) is used; if data type is `strain_only` relaxation modulus (`:G`).
A complete `RheoTimeData` of type `strain_and_stress` is returned.
`diff_method` sets finite difference for calculating the derivative used in the hereditary integral and
can be either backwards difference (`"BD"`) or central difference (`"CD"`).
"""
function modelpredict(data::RheoTimeData, model::RheoModel; diff_method="BD")

    # use correct method for derivative
    if diff_method=="BD"
        deriv = derivBD
    elseif diff_method=="CD"
        deriv = derivCD
    end

    check = RheoTimeDataType(data)
    @assert (check == strain_only)||(check == stress_only) "Need either strain only or stress only data. Data provide: " * string(check)

    if (check == strain_only)
        modulus = _Ga(model)
        modsing = model._G
        dcontrolled = deriv(data.ϵ, data.t)
    elseif (check == stress_only)
        modulus = _Ja(model)
        modsing = model._J
        dcontrolled = deriv(data.σ, data.t)
    end

    # get singularity presence
    sing = singularitytest(modsing)

    # get time step (only needed for convolution, which requires constant so t[2]-t[1] is sufficient)
    dt = data.t[2] - data.t[1]

    # time must start at 0 for convolution to work properly
    t_zeroed = data.t .- minimum(data.t)

    # get convolution
    if !sing && constantcheck(data.t)
        convolved = boltzconvolve(modulus, t_zeroed, dt, dcontrolled)

    elseif sing && constantcheck(data.t)
        t_zeroed[1] = 0.0 + (t_zeroed[2] - t_zeroed[1])/singularity_offset
        convolved = boltzconvolve(modulus, t_zeroed, dt, dcontrolled)

    elseif !sing && !constantcheck(data.t)
        convolved = boltzintegral_nonsing(modulus, t_zeroed, dcontrolled)

    elseif sing && !constantcheck(data.t)
        convolved = boltzintegral_sing(modulus, t_zeroed, dcontrolled)

    end

    if (check == stress_only)
        sigma = data.σ
        epsilon = convolved
        pred_mod = "creep function J"
    elseif (check == strain_only)
        sigma = convolved
        epsilon = data.ϵ
        pred_mod = "relaxation function G"
    end

    log = data.log == nothing ? nothing : [ data.log;
            RheoLogItem( (type=:process, funct=:modelpredict, params=(model::RheoModel,), keywords=(diff_method = diff_method,)),
                         (comment="Predicted data - modulus: $pred_mod, parameters:$(model.params)",) ) ]

    return RheoTimeData(sigma, epsilon, data.t, log)

end

"""
    modelstepfit(data::RheoTimeData, model::RheoModelClass, modloading::Union{LoadingType,Integer}; step=nothing, p0::Union{NamedTuple,Nothing} = nothing, lo::Union{NamedTuple,Nothing} = nothing, hi::Union{NamedTuple,Nothing} = nothing, verbose::Bool = false, rel_tol = 1e-4, weights::Union{Vector{Integer},Nothing} = nothing)

Same as `modelfit` except assumes a step loading. If this assumption is appropriate for the data
then fitting can be sped up greatly by use of this function. If modloading is `strain_imposed`, relaxation modulus is used,
then the element in the middle of the strain is assumed to be the amplitude of the step. If modloading is `stress_imposed`,
the creep modulus is used, then the middle element of the stress is assumed to be the amplitude of the step.
Alternatively, it is possible to define the value of the step by defining the optional `step` parameter.

# Arguments

- `data`: `RheoTimeData` struct containing all data
- `model`: `RheoModelClass` containing moduli and parameters tuples
- `modloading`: `strain_imposed` for relaxation modulus, `stress_imposed` for creep modulus
- `step`: Optional amplitude for step
- `p0`: Named tuple of initial parameters to use in fit (uses 0.5 for all parameters if none given)
- `lo`: Named tuple of lower bounds for parameters
- `hi`: Named tuple of upper bounds for parameters
- `verbose`: If true, prints parameters on each optimisation iteration
- `rel_tol`: Relative tolerance of optimization, see NLOpt docs for more details
- `weights`: Vector of indices weighted according to importance, can be generated by `indexweight` function
"""
function modelstepfit(data::RheoTimeData,
                  model::RheoModelClass,
                  modloading::Union{LoadingType,Integer};
                  step = nothing,
                  p0::Union{NamedTuple,Nothing} = nothing,
                  lo::Union{NamedTuple,Nothing} = nothing,
                  hi::Union{NamedTuple,Nothing} = nothing,
                  verbose::Bool = false,
                  rel_tol = 1e-4,
                  weights::Union{Vector{Integer},Nothing} = nothing) where {T1<:Real, T2<:Real, T3<:Real}

    p0a = fill_init_params(model, p0)
    loa = fill_lower_bounds(model, lo)
    hia = fill_upper_bounds(model, hi)

    rel_tol = convert(RheoFloat, rel_tol)

    # get data type provided
    check = RheoTimeDataType(data)
    # if step amplitude not provided, use loading data specified
    if isnothing(step)
        @assert (check == strain_and_stress) "Both stress and strain are required if step amplitude not provided."
        if (modloading == stress_imposed)
            # step amplitude is set to the middle value of the 'loading' data array
            controlled = data.σ[round(Integer, length(data.σ)/2)]
            measured = data.ϵ
            modulus = _Ja(model)
            modsing = (t->model._J(t, p0a))
            modused = "J"
        elseif (modloading == strain_imposed)
            # step amplitude is set to the middle value of the 'loading' data array
            controlled = data.ϵ[round(Integer, length(data.σ)/2)]
            measured = data.σ
            modulus = _Ga(model)
            modsing = (t->model._G(t, p0a))
            modused = "G"
        end
    # if step provided then can use that value
    elseif !isnothing(step)
        if (modloading == stress_imposed)
            @assert (check == strain_only)||(check == strain_and_stress) "Strain required"
            modulus = _Ja(model)
            modsing = (t->model._J(t, p0a))
            controlled = convert(RheoFloat, step);
            measured = data.ϵ
            modused = "J"
        elseif (modloading == strain_imposed)
            @assert (check == stress_only)||(check == strain_and_stress) "Stress required"
            measured = data.σ
            modulus = _Ga(model)
            modsing = (t->model._G(t, p0a))
            controlled =convert(RheoFloat, step);
            modused = "G"
        end
    end

    # step is assumed to be at 0, warn if data does not start at 0
    if data.t[1] != 0
        @warn "Step fitting assumes that step occurs at time 0.0s but data does not start at 0.0s. If this mismatch is unintended, this fitting may return erroneous results."
    end

    # check necessary modulus is defined
    @assert modulusexists(modsing) "Modulus to use not defined"
    # get singularity presence
    sing = singularitytest(modsing)

    # time must start at 0 for modulus to be defined properly
    t_zeroed = data.t .- minimum(data.t)

    # start fit
    (minf, minx, ret), timetaken, bytes, gctime, memalloc =
                        @timed leastsquares_stepinit(p0a,
                                                    loa,
                                                    hia,
                                                    modulus,
                                                    t_zeroed,
                                                    controlled,
                                                    measured;
                                                    insight = verbose,
                                                    singularity = sing,
                                                    _rel_tol = rel_tol,
                                                    indweights = weights)

    nt = NamedTuple{Tuple(model.params)}(minx)

    if data.log != nothing
        # Preparation of data for log item
        info=(comment="Fiting rheological model to data (step input assumed)", model_name=model.name, model_params=nt, time_taken=timetaken, stop_reason=ret, error=minf)
        params = (model = model, modloading = modloading)
        keywords = (step = step, p0 = p0, lo = lo, hi = hi, rel_tol = rel_tol)
        # Add data to the log
        push!(data.log, RheoLogItem( (type=:analysis, funct=:modelstepfit, params=params, keywords=keywords), info))
    end

    print("Time: $timetaken s, Why: $ret, Parameters: $minx, Error: $minf")

    return RheoModel(model, nt)
end

"""
    modelsteppredict(data::RheoTimeData, model::RheoModel; step_on::Real = 0.0)

Same as `modelpredict` but assumes a step loading with step starting at `step_on` or closest actual value to that
specified. If the loading data is variable, the magnitude in the middle of the array is used.
"""
function modelsteppredict(data::RheoTimeData, model::RheoModel; step_on::Real = 0.0)

    check = RheoTimeDataType(data)
    @assert (check == strain_only)||(check == stress_only) "Need either strain only or stress only data. Data provide: " * string(check)

    if (check == strain_only)
        modulus = _Ga(model)
        modsing = model._G
        controlled = data.ϵ[round(Integer, length(data.ϵ)/2)]
    elseif (check == stress_only)
        modulus = _Ja(model)
        modsing = model._J
        controlled = data.σ[round(Integer, length(data.σ)/2)]
    end

    # get closest index and time to desired step-on time point
    stepon_el = closestindex(data.t, step_on)
    stepon_closest = data.t[stepon_el]

    # check singularity presence at time closest to step
    sing = singularitytest(modsing)

    # get predicted
    predicted = zeros(RheoFloat, length(data.t))
    if !sing
        predicted[stepon_el:end] = controlled*modulus(data.t[stepon_el:end] .- stepon_closest)
    elseif sing
        predicted[stepon_el:end] = controlled*modulus(data.t[stepon_el:end] .- stepon_closest)
        @warn "Singularity (Inf value) is present in modulus used for modelsteppredict"
    end

    if check == stress_only
        σ = data.σ
        ϵ = predicted
    elseif check == strain_only
        σ = predicted
        ϵ = data.ϵ
    end

    # store operation
    if (check == stress_only)
        sigma = data.σ
        epsilon = predicted
        pred_mod = "creep function J"
    elseif (check == strain_only)
        sigma = predicted
        epsilon = data.ϵ
        pred_mod = "relaxation function G"
    end
    time = data.t

    log = data.log == nothing ? nothing : [ data.log;
            RheoLogItem( (type=:process, funct=:modelsteppredict, params=(model::RheoModel,), keywords=(step_on = step_on,)),
                         (comment="Predicted step response - modulus: $pred_mod, parameters:$(model.params)",) ) ]

    return RheoTimeData(sigma,epsilon,time, log)
end

function obj_dynamic(params, grad, ω, dataGp, dataGpp, modelGp, modelGpp; _insight::Bool = false)

    if _insight
        println("Current Parameters: ", params)
    end

    costGp = sum((dataGp - modelGp(ω, params)).^2)
    costGpp = sum((dataGpp - modelGpp(ω, params)).^2)

    cost = costGp + costGpp
end

function obj_dynamic_mean(params, grad, ω, dataGp, dataGpp, modelGp, modelGpp, meanGp, meanGpp; _insight::Bool = false)

    if _insight
        println("Current Parameters: ", params)
    end

    costGp = sum((dataGp/meanGp - modelGp(ω, params)/meanGp).^2)
    costGpp = sum((dataGpp/meanGpp - modelGpp(ω, params)/meanGpp).^2)

    cost = costGp + costGpp

end

function obj_dynamic_log(params, grad, ω, dataGp, dataGpp, modelGp, modelGpp; _insight::Bool = false)

    if _insight
        println("Current Parameters: ", params)
    end

    costGp = sum((log.(dataGp) - log.(modelGp(ω, params))).^2)
    costGpp = sum((log.(dataGpp) - log.(modelGpp(ω, params))).^2)

    cost = costGp + costGpp

end

function obj_dynamic_local(params, grad, ω, dataGp, dataGpp, modelGp, modelGpp; _insight::Bool = false)

    if _insight
        println("Current Parameters: ", params)
    end

    costGp = sum(((dataGp - modelGp(ω, params))./dataGp).^2)
    costGpp = sum(((dataGpp - modelGpp(ω, params))./dataGpp).^2)

    cost = costGp + costGpp

end

function obj_dynamic_manual(params, grad, ω, dataGp, dataGpp, modelGp, modelGpp, weights; _insight::Bool = false)

    if _insight
        println("Current Parameters: ", params)
    end

    costGp = sum((dataGp - modelGp(ω, params)).^2)
    costGpp = sum((dataGpp - modelGpp(ω, params)).^2)

    cost = weights[1]*costGp + weights[2]*costGpp
end

"""
    dynamicmodelfit(data::RheoFreqData, model::RheoModelClass; p0::Union{NamedTuple,Nothing} = nothing, lo::Union{NamedTuple,Nothing} = nothing, hi::Union{NamedTuple,Nothing} = nothing, verbose::Bool = false, rel_tol = 1e-4) where T<:Real

Fits model to the frequency/loss+storage moduli data.

All arguments are as described below.
As this fitting procedure is fitting two functions simultaneously (the storage
and loss moduli), if left untransformed the fit would tend to favour the
modulus which is larger in magnitude and not fit the other modulus well. To avoid this,
RHEOS offers a number of data transforms which can be used by changing `weights` argument.

# Arguments

- `data`: `RheoFreqData` struct containing all data
- `model`: `RheoModelClass` containing moduli and symbols of parameters
- `p0`: Initial parameters to use in fit (uses 0.5 for all parameters if none given)
- `lo`: Lower bounds for parameters
- `hi`: Upper bounds for parameters
- `verbose`: If true, prints parameters on each optimisation iteration
- `rel_tol`: Relative tolerance of optimization, see NLOpt docs for more details
- `weights`: Weighting mode for storage and loss modulus (`"none"`, `"mean"`, `"log"`, `"local"` or manually specified)
"""
function dynamicmodelfit(data::RheoFreqData,
                model::RheoModelClass;
                p0::Union{NamedTuple,Nothing} = nothing,
                lo::Union{NamedTuple,Nothing} = nothing,
                hi::Union{NamedTuple,Nothing} = nothing,
                verbose::Bool = false,
                rel_tol::T = 1e-4,
                weights::Union{String, Vector{T}}="local") where T<:Real

    p0a = fill_init_params(model, p0)
    loa = fill_lower_bounds(model, lo)
    hia = fill_upper_bounds(model, hi)

    # check necessary moduli are defined
    modsingGp = (ω->model._Gp(ω,p0a))
    modsingGpp = (ω->model._Gpp(ω,p0a))
    @assert modulusexists(modsingGp) "Storage modulus not defined for this model"
    @assert modulusexists(modsingGpp) "Loss modulus not defined for this model"

    # initialise NLOpt.Opt object with :LN_SBPLX Subplex algorithm
    opt = Opt(:LN_SBPLX, length(p0a))

    if !isnothing(loa)
       lower_bounds!(opt, loa)
    end

    if !isnothing(hia)
       upper_bounds!(opt, hia)
    end

    # set relative tolerance
    rel_tol = convert(RheoFloat, rel_tol)
    xtol_rel!(opt, rel_tol)

    # set objective/cost function
    if weights=="none"
        min_objective!(opt, (params, grad) -> obj_dynamic(params, grad, data.ω, data.Gp, data.Gpp, model._Gpa, model._Gppa; _insight = verbose))

    elseif weights=="mean"
        meanGp = sum(data.Gp)/length(data.Gp)
        meanGpp = sum(data.Gp)/length(data.Gp)
        min_objective!(opt, (params, grad) -> obj_dynamic_mean(params, grad, data.ω, data.Gp, data.Gpp, model._Gpa, model._Gppa, meanGp, meanGpp; _insight = verbose))

    elseif weights=="log"
        @warn "Note that a logarithmic rescaling will fail if Gp or Gpp data contain 0.0 values as it will result in -Inf cost. Trying a different rescaling scheme, or not fitting around ω≈0.0 may alleviate the issue."
        min_objective!(opt, (params, grad) -> obj_dynamic_log(params, grad, data.ω, data.Gp, data.Gpp, model._Gpa, model._Gppa; _insight = verbose))

    elseif weights=="local"
        @warn "Note that a local rescaling will fail if Gp or Gpp data contain 0.0 values as it will result in division by 0.0. Trying a different rescaling scheme, or not fitting around ω≈0.0 may alleviate the issue."
        min_objective!(opt, (params, grad) -> obj_dynamic_local(params, grad, data.ω, data.Gp, data.Gpp, model._Gpa, model._Gppa; _insight = verbose))

    elseif typeof(weights)==Vector{T} && length(weights)==2
        min_objective!(opt, (params, grad) -> obj_dynamic_manual(params, grad, data.ω, data.Gp, data.Gpp, model._Gpa, model._Gppa, weights; _insight = verbose))

    end

    # timed fitting
    (minf, minx, ret), timetaken, bytes, gctime, memalloc = @timed NLopt.optimize(opt, p0a)

    nt = NamedTuple{Tuple(model.params)}(minx)

    if data.log != nothing
        # Preparation of data for log item
        info = (comment="Fiting rheological model to frequency spectrum", model_name=model.name, model_params=nt, time_taken=timetaken, stop_reason=ret, error=minf)
        params = (model=model, )
        keywords = (p0=p0, lo=lo, hi=hi, rel_tol=rel_tol, weights=weights)
        # Add data to the log
        push!(data.log, RheoLogItem( (type=:analysis, funct=:dynamicmodelfit, params=params, keywords=keywords), info))
    end

    print("Time: $timetaken s, Why: $ret, Parameters: $minx, Error: $minf")

    return RheoModel(model, nt)
end

"""
    dynamicmodelpredict(data::RheoFreqData, model::RheoModel)

Given dynamic rheology data with only frequency and model where parameters have been substituted.
Returns another `RheoFreqData` instance with the predicted `Gp` and `Gpp` based on the
frequencies and model given as arguments.
"""
function dynamicmodelpredict(data::RheoFreqData, model::RheoModel)

    predGp = model._Gpa(data.ω)
    predGpp = model._Gppa(data.ω)

    log = data.log == nothing ? nothing : [ data.log;
            RheoLogItem( (type=:process, funct=:dynamicmodelpredict, params=(model::RheoModel,), keywords=() ),
                         (comment="Calculated frequency spectrum - parameters:$(model.params)",) ) ]

    return RheoFreqData(predGp, predGpp, data.ω, log)

end


text:
using Pkg
Pkg.add(PackageSpec(name = "RDatasets", version = "0.7.5"))
Pkg.add(PackageSpec(name = "Flux", version = "0.12"))

# Import packages and set a seed
import Flux
using Random
Random.seed!(1);

# Load data and split it into training and test data
function rand_split_data(x, labels)
   nsamples = size(x, 2)
   testidxs = randperm(nsamples)[1:(round(Int, nsamples*0.3))]
   trainidxs = setdiff(1:nsamples, testidxs)
   x_train = x[:, trainidxs]
   x_test = x[:, testidxs]
   labels_train = labels[trainidxs]
   labels_test = labels[testidxs]
   y = Flux.onehotbatch(labels, unique(labels))
   y_train = y[:, trainidxs]
   y_test = y[:, testidxs]
   (training = (x = x_train, y = y_train),
         test = (x = x_test, y = y_test))
end

using RDatasets
import Tables
iris = dataset("datasets", "iris")
x = Tables.matrix(iris[:, 1:4])'
data = rand_split_data(x, iris[:, :Species])
trainingdata = data.training
testdata = data.test

# Load necessary features
import Flux
using Statistics

# Train the Flux model
model = Flux.Chain(
      Flux.Dense(4, 4, Flux.relu),
      Flux.Dense(4, 4, Flux.relu),
      Flux.Dense(4, 3),
      Flux.softmax)

loss(model, x, y) = Flux.crossentropy(model(x), y)
loss(model, data::NamedTuple) = loss(model, data.x, data.y)

function train_network!(model, x, y; epochs, callback)
   opt = Flux.ADAM()
   loss_(x, y) = loss(model, x, y)
   for i in 1:epochs
      Flux.train!(loss_, Flux.params(model), [(x, y)], opt)
      callback(i)
   end
end

epochs = 2500
train_losses = Vector{Float32}(undef, epochs);
test_losses = Vector{Float32}(undef, epochs);
train_network!(model, trainingdata.x, trainingdata.y, epochs = epochs,
      callback = (i) -> begin
         train_losses[i] = loss(model, trainingdata)
         test_losses[i] = loss(model, testdata)
      end)

# (Could plot losses here)

# Evaluate model
accuracy(model, data) =
      mean(Flux.onecold(model(data.x)) .== Flux.onecold(data.y))
accuracy(model, trainingdata)
accuracy(model, testdata)

ext: jl


text:
using Documenter, OpenML, DataFrames

makedocs(
    modules = [OpenML,],
    sitename = "OpenML.jl",
)

deploydocs(
    repo = "github.com/JuliaAI/OpenML.jl",
    push_preview = true
)


text:
"""
Remove spike `s` from event at index `k`.
"""
function remove_datapoint!(
    model::SeqModel,
    s::Spike,
    k::Int64
)

    # Nothing to do if spike is already in background.
    (k == -1) && return

    # Remove the contributions of datapoint i to the
    # sufficient statistics of sequence event k.
    t = s.timestamp
    n = s.neuron
    event = model.sequence_events[k]

    log_p_neuron = model.globals.neuron_response_log_proportions
    offsets = model.globals.neuron_response_offsets
    widths = model.globals.neuron_response_widths
    warps = model.priors.warp_values

    # If this is the last spike in the event, we can return early.
    (event.spike_count == 1) && (return remove_event!(model.sequence_events, k))

    # Otherwise, subtract off the sufficient statistics.
    event.spike_count -= 1

    for r = 1:num_sequence_types(model)
        for w in 1:num_warp_values(model)
            # Compute 1D precision and potential.
            v = 1 / (widths[n, r] * warps[w]^2)
            m = (t - offsets[n, r] * warps[w]) * v
            sv = event.summed_precisions[r, w]

            # Subtract sufficient statistics.
            event.summed_potentials[r, w] -= m
            event.summed_precisions[r, w] = max(0, sv - v) # prevent negative precision.
            event.summed_logZ[r, w] -= (log_p_neuron[n, r] - gauss_info_logZ(m, v))
            event.seq_type_posterior[r, w] = (
                model.globals.seq_type_log_proportions[r]
                + model.priors.warp_log_proportions[w]
                + event.summed_logZ[r, w]
                + gauss_info_logZ(
                    event.summed_potentials[r, w],
                    event.summed_precisions[r, w])
            )
        end
    end

    event.seq_type_posterior .-= logsumexp(event.seq_type_posterior)

end

"""
Add spike `s` to event at index `k`. Return assignment index `k`.
"""
function add_datapoint!(
        model::SeqModel,
        s::Spike,
        k::Int64;
        recompute_posterior::Bool=true
    )

    t = s.timestamp
    n = s.neuron
    event = model.sequence_events[k]

    log_p_neuron = model.globals.neuron_response_log_proportions
    offsets = model.globals.neuron_response_offsets
    widths = model.globals.neuron_response_widths
    warps = model.priors.warp_values

    event.spike_count += 1

    for r = 1:num_sequence_types(model)
        for w in 1:num_warp_values(model)
            # Compute 1D precision and potential.
            v = 1 / (widths[n, r] * warps[w]^2)
            m = (t - offsets[n, r] * warps[w]) * v

            # Add sufficient statistics.
            event.summed_potentials[r, w] += m
            event.summed_precisions[r, w] += v
            event.summed_logZ[r, w] += (log_p_neuron[n, r] - gauss_info_logZ(m, v))
        end
    end

    if recompute_posterior
        set_posterior!(model, k)
    end

    return k
end

function set_posterior!(
    model::SeqModel,
    k::Int64;
)
    event = model.sequence_events[k]

    for r = 1:num_sequence_types(model)
        for w in 1:num_warp_values(model)
            event.seq_type_posterior[r, w] = (
                model.globals.seq_type_log_proportions[r]
                + model.priors.warp_log_proportions[w]
                + event.summed_logZ[r, w]
                + gauss_info_logZ(
                    event.summed_potentials[r, w],
                    event.summed_precisions[r, w])
            )
        end
    end

    event.seq_type_posterior .-= logsumexp(event.seq_type_posterior)
end

"""
Create a singleton cluster / sequence event with spike `s`
and return new assignment index `k`.
"""
function add_event!(model::SeqModel, s::Spike)
    # Mark event k as non-empty.
    k = add_event!(model.sequence_events)
    # Add spike s to event k.
    return add_datapoint!(model, s, k)
end


text:
## Company Statement
Base.@kwdef mutable struct CompanyResults
    EBITDA::Float64 #Ebitda
    ghg_emissions::Float64 # GHG emissions in milions tonnes of CO2eq
    climate_costs::Union{Nothing,Float64} = nothing # Climate-related damages to society from the company, in Billion USD
end

company_a = CompanyResults(EBITDA = 25.2, ghg_emissions =  374)

## GHG emissions valuation
using Plots

function climate_valuation!(data::CompanyResults, scc::Float64)::CompanyResults
    data.climate_costs = - data.ghg_emissions * scc * 10^(-3) # to express in billion USD
    legend_axis = ["EBITDA", "GHG",
                "Air Pollution",
                "Water Consumption",
                "Water Pollution",
                "Land Use",
                "Waste","Incidents"]

    air_pollution = data.climate_costs * 0.2
    water_consumption = data.climate_costs * 0.1
    water_pollution = data.climate_costs * 0.08
    land_use = data.climate_costs * 0.05
    waste = data.climate_costs * 0.04
    incidents = data.climate_costs * 0.02

    list_data = [data.EBITDA,
                data.climate_costs,
                air_pollution,
                water_consumption,
                water_pollution,
                land_use,
                waste,
                incidents]
    list_data = reduce(vcat, list_data)
    println(list_data)
    display(bar(legend_axis, list_data, orientation = :horizontal, bar_width = 1, legend = nothing, xlabel = "Billion USD"))
    return data
end

climate_valuation!(company_a, 51.)

company_a.climate_costs / company_a.EBITDA * -1

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance


text:
function MatrixProfile.matrix_profile(T, m::Int, dist::DTW{<:Any, N}; showprogress=true) where N
    n   = lastlength(T)
    l   = n-m+1
    r   = dist.radius
    # n > 2m+1 || throw(ArgumentError("Window length too long, maximum length is $((n+1)÷2)"))
    P    = Vector{floattype(T)}(undef, l)
    I    = Vector{Int}(undef, l)
    prog = Progress((l - 1) ÷ 5, dt=1, desc="Matrix profile", barglyphs = BarGlyphs("[=> ]"), color=:blue)
    bsf = typemax(floattype(T))
    @inbounds for i = 1:l
        Ti = getwindow(T,m,i)
        res = dtwnn(Ti, T, dist.dist, r, N, transportcost=dist.transportcost, avoid=i-r:i+r)
        I[i] = res.loc
        P[i] = res.cost
        showprogress && i % 5 == 0 && next!(prog)
    end
    MatrixProfile.Profile(T, P, I, m, nothing)
end


text:
# abstract types
abstract type PointProcess end;

# specific types
mutable struct DeterminantalPointProcess <: PointProcess
    L::Symmetric
    Lfact::Eigen
    size::Int
    rng::AbstractRNG

    function DeterminantalPointProcess(L::Symmetric; seed::Int=42)
        Lfact = eigen(L)
        new(L, Lfact, length(Lfact.values), MersenneTwister(seed))
    end

    function DeterminantalPointProcess(Lfact::Eigen; seed::Int=42)
        L = Symmetric((Lfact.vectors .* Lfact.values') * Lfact.vectors')
        new(L, Lfact, length(Lfact.values), MersenneTwister(seed))
    end
end

mutable struct KroneckerDeterminantalPointProcess <: PointProcess
    # TODO
end

# aliases
const DPP = DeterminantalPointProcess
const KDPP = KroneckerDeterminantalPointProcess
const MCMCState = Tuple{BitArray{1}, Array{Float64, 2}}


text:
# SimpleMPSJulia
# Density matrix renormalization group (DMRG) in matrix product state (MPS)

# This file contains the defination of the Heisenberg model in MPO structure
# and solves its ground state
using LinearAlgebra

include("./mps.jl")
using .SimpleMPS

"""
Construct the building block of the Heisenberg XXZ model MPO.
The returned array has 4 dimensions and the physical dimensions
are in the middle.
"""
function heisenberg_mo(J, Jz, h)::Array
    # S^+
    Sp = [[0., 1.] [0., 0.]]
    # S^-
    Sm = [[0., 0.] [1., 0.]]

    Sz = [[0.5, 0.] [0., -0.5]]

    # zero matrix block
    S0 = zeros((2, 2))

    # identity matrix block
    S1 = Diagonal(ones(2))

    mpo_block1 = cat(S1, S0, S0, S0, S0, dims=3)
    mpo_block2 = cat(Sp, S0, S0, S0, S0, dims=3)
    mpo_block3 = cat(Sm, S0, S0, S0, S0, dims=3)
    mpo_block4 = cat(Sz, S0, S0, S0, S0, dims=3)
    mpo_block5 = cat(-h * Sz, J / 2 * Sm, J / 2 * Sp, Jz * Sz, S1, dims=3)

    # after the permutation, the first index corresponds to the vertical axis of the above matrix
    # and the last index corresponds to the horizontal axis of the above matrix
    single_mo = permutedims(cat(mpo_block1, mpo_block2, mpo_block3, mpo_block4, mpo_block5, dims=4), (4, 1, 2, 3))
end

const J = Jz = h = 1
mo = heisenberg_mo(J, Jz, h)

const site_num = 20
mpo = MatrixProductOperator(mo, site_num)

const bond_dimension = 20
mps, e = search_groundstate(mpo, bond_dimension)
println("Energies:")
println(e)


text:
###########Zero Rate Curve Types
#Abstract
abstract type AbstractZeroRateCurve end

#Scalar
struct ZeroRate{T2 <: Number} <: AbstractZeroRateCurve
    r::T2
    function ZeroRate(r::T2) where {T2 <: Number}
        return new{T2}(r)
    end
end

#Curve
struct ZeroRateCurve{num1 <: Number, num2 <: Number} <: AbstractZeroRateCurve
    r::Curve{num1, num2}
    function ZeroRateCurve(r_::Array{num1}, T::num2) where {num1 <: Number, num2 <: Number}
        new{num2, num1}(Curve(r_, T))
    end

    function ZeroRateCurve(r_::Curve{num1, num2}) where {num1 <: Number, num2 <: Number}
        new{num2, num1}(r_)
    end
end

#Common Constructors
function ZeroRate(r_::Array{num1}, T::num2) where {num1 <: Number, num2 <: Number}
    return ZeroRateCurve(r_, T)
end
function ImpliedZeroRate(r_::Array{num1}, T::num2) where {num1 <: Number, num2 <: Number}
    return ZeroRateCurve(ImpliedCurve(r_, T))
end
function ImpliedZeroRate(r_::Array{num1}, T::Array{num2}) where {num1 <: Number, num2 <: Number}
    return ZeroRateCurve(ImpliedCurve(r_, T))
end

export ZeroRateCurve;
export ZeroRate;


text:
@testset "Optimising" begin
    @testset "Parsing FlatZinc output format" begin
        @testset "one_solution.fzn" begin
            # Most simple output.
            out_string = "x = 10;\r\n\r\n----------\r\n==========\r\n"
            @test CP.FlatZinc._parse_to_assignments(out_string) == [
                Dict("x" => [10])
            ]
        end

        @testset "Infeasible" begin
            # Output for an infeasible model.
            out_string = "=====UNSATISFIABLE=====\r\n"
            @test CP.FlatZinc._parse_to_assignments(out_string) == []
        end

        @testset "basic.fzn" begin
            # Marker for the end of search: '=' ^ 10
            out_string = "x = 3;\r\n\r\n----------\r\n"
            @test CP.FlatZinc._parse_to_assignments(out_string) == [
                Dict("x" => [3])
            ]
        end

        @testset "basic.fzn with a float value" begin
            # Marker for the end of search: '=' ^ 10
            out_string = "x = 3.0;\r\n\r\n----------\r\n"
            @test CP.FlatZinc._parse_to_assignments(out_string) == [
                Dict("x" => [3.0])
            ]
        end

        @testset "several_solutions.fzn" begin
            # Several solutions (with CLI parameter -a)
            out_string = "xs = array1d(1..2, [2, 3]);\r\n\r\n----------\r\nxs = array1d(1..2, [1, 3]);\r\n\r\n----------\r\nxs = array1d(1..2, [1, 2]);\r\n\r\n----------\r\n==========\r\n"
            @test CP.FlatZinc._parse_to_assignments(out_string) == [
                Dict("xs" => [2, 3]),
                Dict("xs" => [1, 3]),
                Dict("xs" => [1, 2]),
            ]
        end

        @testset "puzzle.fzn" begin
            # 2D array
            out_string = "x = array2d(1..4, 1..4, [5, 1, 8, 8, 9, 3, 8, 6, 9, 7, 7, 8, 1, 7, 8, 9]);\r\n\r\n----------\r\n"
            @test CP.FlatZinc._parse_to_assignments(out_string) == [
                Dict("x" => [5, 1, 8, 8, 9, 3, 8, 6, 9, 7, 7, 8, 1, 7, 8, 9])
            ]
        end

        @testset "einstein.fzn" begin
            # Multiple variables
            out_string = "a = array1d(1..5, [5, 4, 3, 1, 2]);\r\nc = array1d(1..5, [3, 4, 5, 1, 2]);\r\nd = array1d(1..5, [2, 4, 3, 5, 1]);\r\nk = array1d(1..5, [3, 1, 2, 5, 4]);\r\ns = array1d(1..5, [3, 5, 2, 1, 4]);\r\n\r\n----------\r\n"
            @test CP.FlatZinc._parse_to_assignments(out_string) == [
                Dict(
                    "a" => [5, 4, 3, 1, 2],
                    "c" => [3, 4, 5, 1, 2],
                    "d" => [2, 4, 3, 5, 1],
                    "k" => [3, 1, 2, 5, 4],
                    "s" => [3, 5, 2, 1, 4],
                )
            ]
        end
    end
end


text:
let a = NDSparse([12,21,32], [52,41,34], [11,53,150]), b = NDSparse([12,23,32], [52,43,34], [56,13,10])
    p = vec(collect(IndexedTables.product(a, b)))
    @test p == [(11,56), (11,13), (11,10), (53,56), (53,13), (53,10), (150,56), (150,13), (150,10)]

    p = vec(collect(IndexedTables.product(a, b, a)))
    @test p == [(11,56,11),(11,56,53),(11,56,150),(11,13,11),(11,13,53),(11,13,150),(11,10,11),(11,10,53),(11,10,150),
                (53,56,11),(53,56,53),(53,56,150),(53,13,11),(53,13,53),(53,13,150),(53,10,11),(53,10,53),(53,10,150),
                (150,56,11),(150,56,53),(150,56,150),(150,13,11),(150,13,53),(150,13,150),(150,10,11),(150,10,53),(150,10,150)]
end

let a = [1:10;]
    @test IndexedTables._sizehint!(1:10, 20) == 1:10
    @test IndexedTables._sizehint!(a, 20) === a
end

@test Columns(([1,2], [3,4])) == Columns(([1,2], [3.0,4.0]))
@test Columns(([1,2], [3,4])) != Columns(([1,2], [3.0,4.1]))
@test Columns(([1,2], [3,4])) != Columns((a=[1,2], b=[3,4]))

function roundtrips(x)
    b = IOBuffer()
    serialize(b, x)
    seekstart(b)
    return deserialize(b) == x
end

@test roundtrips(Columns((rand(5), rand(5))))
@test roundtrips(Columns(c1 = rand(5), c2 = rand(5)))
@test roundtrips(convert(NDSparse, rand(3,3)))
@test roundtrips(NDSparse(Columns(y=rand(3), x=rand(3)), rand(3)))

let x = rand(3), y = rand(3), v = rand(3), w = rand(3)
    @test vcat((Columns((x,y))), Columns((v,w))) == Columns((vcat(x,v), vcat(y,w)))
    @test vcat(Columns(x=x,y=y), Columns(x=v,y=w)) == Columns(x=vcat(x,v), y=vcat(y,w))
end

@testset "compact_mem" begin
    v = Columns(a=rand(10), b = fill("string", 10))
    v_pooled = IndexedTables.replace_storage(v) do c
        isbitstype(eltype(c)) ? c : convert(PooledArray, c)
    end
    @test eltype(v) == eltype(v_pooled)
    @test all(v.a .== v_pooled.a)
    @test all(v.b .== v_pooled.b)
    @test !isa(v_pooled.a, PooledArray)
    @test isa(v_pooled.b, PooledArray)
    @test v_pooled == IndexedTables.compact_mem(v)
    s = WeakRefStrings.StringArray(["a", "b", "c"])
    @test IndexedTables.compact_mem(s) isa PooledArray{String}
    @test IndexedTables.compact_mem(s)[1] == "a"
    @test IndexedTables.compact_mem(s)[2] == "b"
    @test IndexedTables.compact_mem(s)[3] == "c"
    @test IndexedTables.compact_mem(IndexedTables.compact_mem(s)) == IndexedTables.compact_mem(s)
    s = WeakRefStrings.StringArray(["a", missing, "c"])
    @test IndexedTables.compact_mem(s) isa PooledArray{Union{String, Missing}}
    @test all(isequal.(s, IndexedTables.compact_mem(s)))
end

@testset "refs" begin
    a = WeakRefStrings.StringVector(["a", "b", "c"])
    b = PooledArrays.PooledArray(["1", "2", "3"])
    c = [:a, :b, :c]
    s = Columns(a=a, b=b, c=c)
    ref = IndexedTables.refs(s)
    @test ref[1].a isa WeakRefStrings.WeakRefString{UInt8}
    @test all(isequal.(ref.a, a))
    @test ref[1].b isa UInt8
    @test isequal(ref.b, UInt8.([1, 2, 3]))
    Base.permute!!(ref, sortperm(s))
    @test issorted(s)

    a = WeakRefStrings.StringVector(["a", missing, "c"])
    refa = IndexedTables.refs(a)
    @test refa isa WeakRefStrings.StringArray{Union{Missing, WeakRefStrings.WeakRefString{UInt8}}}
    @test all(isequal.(a, refa))
end


text:
using GeometryTypes, Visualize, Colors, GLAbstraction
using Visualize: Resolution, Open, Visible
using Visualize.GLRasterization: glwindow, renderloop
using Visualize: Model, Color, Position, LineVertex, LineSegments, Camera, BasicCamera

resolution = (1024, 1024)

window = glwindow(
    Resolution => resolution,
    Camera => BasicCamera(),
    Visualize.Debugging => true
);

scale = 10f0
res = Vec2f0(resolution)
radius = (min(resolution...) / 2f0) - scale

N = 32
ls = LineSegments((
    Position => [rand(Vec2f0) .* res for x in 1:N],
    Color => Vec4f0(1, 0,0,1)
))
ls
x = Visualize.visualize(window, ls)

renderloop(window)
window[Camera]
map(x-> window[Camera][ProjectionView] * Vec4f0(x..., 0, 1.0),  [rand(Vec2f0) .* res for x in 1:N])
# Visualize.destroy!(window)

using Visualize, GeometryTypes, Colors, FileIO, NIfTI

window = Visualize.GLRasterization.glwindow(
    Resolution => (1024, 1024),
    Visualize.Color => Vec4f0(0.0, 0.0, 0.0, 0.0),
)

niivol = Array(niread("/media/s/94BC3648BC362560/Users/sdani/juliastuff/demo/t2.nii"))
vol = niivol ./ maximum(niivol)

colorfun2(color, uniforms, intensities) = color

function accum2(accum, pos, stepdir, intensities, uniforms)
    alpha = intensities[pos][1] * 0.1f0
    sample = Vec4f0(0.9f0, 0f0, 0.4f0, alpha)
    alpha = 1f0 - ((1f0 - alpha) ^ 1.2f0)
    sample = Vec4f0(
        sample[1] * alpha,
        sample[2] * alpha,
        sample[3] * alpha,
        alpha
    )
    sample = (1f0 - accum[4]) * sample
    tmp = accum + sample
    tmp, tmp[4] > 0.95f0
end

volume = Visualize.Volume((
    Data => vol,
    AccumulationFunction => accum2,
    StartValue => Vec4f0(0),
    ColorFunction => colorfun2,
))

wf = wireframe(
    AABB(Vec3f0(-0.01), Vec3f0(1.02)), ()
)

visualize(window, wf)
visualize(window, volume)
@async renderloop(window)


text:
#---------------------------------------------------------------------------------------------------
# iTEBD
#---------------------------------------------------------------------------------------------------
# gate
function gate!(
    mat::AbstractMatrix,
    ΓA::AbstractArray{<:Number, 3},
    ΓB::AbstractArray{<:Number, 3},
    λl::AbstractVector{<:Real},
    d1::Integer,
    d2::Integer,
    bound::Integer,
    cutoff::Real
)
    α, β = size(ΓA, 1), size(ΓB, 3)
    tlmul!(λl, ΓA)
    Γ = tcmul(mat, tgroup(ΓA, ΓB))
    ΓA2, λc, ΓB2 = tsvd(Γ, α, d1, d2, β, bound=bound, cutoff=cutoff, renormalize=true)
    trmul!(ΓA2, λc)
    tlmul!(1.0 ./ λl, ΓA2)
    ΓA2, λc, ΓB2
end
#---------------------------------------------------------------------------------------------------
# 2-site iTEBD
export itebd2!
function itebd2!(
    mat::AbstractMatrix,
    ΓA::AbstractArray{<:Number, 3},
    ΓB::AbstractArray{<:Number, 3},
    λl::AbstractVector{<:Real};
    bound::Integer=100,
    cutoff::Real=1e-7
)
    d1, d2 = size(ΓA, 2), size(ΓB, 2)
    ΓA2, λ1, ΓB2 = gate!(mat,  ΓA,  ΓB, λl, d1, d2, bound, cutoff)
    ΓB3, λ2, ΓA3 = gate!(mat, ΓB2, ΓA2, λ1, d2, d1, bound, cutoff)
    ΓA3, ΓB3, λ2
end

#---------------------------------------------------------------------------------------------------
# 3-site iTEBD
export itebd3!
function itebd3!(
    mat::AbstractMatrix,
    ΓAB1::AbstractArray{<:Number, 3},
    ΓC1::AbstractArray{<:Number, 3},
    λl::AbstractVector{<:Real};
    bound::Integer=100,
    cutoff::Real=1e-7
)
    d1, d2, d3 = size(ΓA, 2), size(ΓB, 2),  size(ΓC, 2)
    ΓA2, λ1, ΓBC2 = gate!(mat, ΓAB1, ΓC1, λl, d1, d2 * d3, bound, cutoff)
    ΓB3, λ2, ΓCA3 = gate!(mat, ΓBC2, ΓA2, λ1, d2, d3 * d1, bound, cutoff)
    ΓC4, λ3, ΓAB4 = gate!(mat, ΓBC2, ΓA2, λ2, d3, d1 * d2, bound, cutoff)
    ΓAB4, ΓC4, λ3
end


text:
# TODO: Try to make as many types as possible immutable
# TODO: Generalize to arbitrary metric space

abstract AbstractBallTreeNode

type EmptyBallTreeNode <: AbstractBallTreeNode
end

type BallTreeNode <: AbstractBallTreeNode
	ball::AbstractBall
	index::Int # Index of ball's center point in source data
	parent::AbstractBallTreeNode
	left::AbstractBallTreeNode
	right::AbstractBallTreeNode
end

function BallTreeNode()
	BallTreeNode(EmptyBall(),
		         -1,
	             EmptyBallTreeNode(),
	             EmptyBallTreeNode(),
	             EmptyBallTreeNode())
end

immutable BallTree
	root::AbstractBallTreeNode
	n::Int # Number of leaves
	# Other metadata
end

function isleaf(n::AbstractBallTreeNode)
	return isa(n.left, EmptyBallTreeNode) &&
	       isa(n.right, EmptyBallTreeNode)
end

function Base.show(io::IO, n::EmptyBallTreeNode)
	print(io, EmptyBallTreeNode)
end

function Base.show(io::IO, n::BallTreeNode)
	print(io, n.ball)
end

function walk(t::BallTree)
	println("BallTree")
	walk(t.root)
end

function walk(n::BallTreeNode, depth::Integer = 0)
	for i in 1:depth
		print(' ')
	end
	println(n.ball)
	walk(n.left, depth + 1)
	walk(n.right, depth + 1)
	return
end

function walk(n::EmptyBallTreeNode, depth::Integer = 0)
	return
end

check(t::BallTree) = check(t.root)

function check(n::BallTreeNode)
	if n.ball.r == 0.0
		return true
	else
		return contains(n.ball, n.left.ball) &&
		       contains(n.ball, n.right.ball) &&
			   check(n.left) &&
			   check(n.right)
	end
end

check(n::EmptyBallTreeNode) = true

function BallTree{T <: Real}(X::Matrix{T})
	p, n = size(X)
	balls = Array(AbstractBallTreeNode, n)
	for i in 1:n
		balls[i] = BallTreeNode(Ball(X[:, i], 0.0),
			                    i,
			                    EmptyBallTreeNode(),
			                    EmptyBallTreeNode(),
			                    EmptyBallTreeNode())
	end
	return BallTree(wrap(balls, 1, n), n)
end

function coordselect!(balls::Vector{AbstractBallTreeNode},
	                  c::Integer,
	                  k::Integer,
	                  l::Integer,
	                  u::Integer)
	# Semi-sort balls[l:u] around k-th element on the c-th dimension
	#  * Items below have >=
	#  * Items above have <=
	while l < u
		r = rand(l:u)
		balls[r], balls[l] = balls[l], balls[r]
		m = l
		i = l + 1
		while i <= u
			if balls[i].ball.center[c] < balls[l].ball.center[c]
				m += 1
				balls[m], balls[i] = balls[i], balls[m]
			end
			i += 1
		end
		balls[l], balls[m] = balls[m], balls[l]
		if m <= k
			l = m + 1
		end
		if m >= k
			u = m - 1
		end
	end
	return
end

# Doing this exactly is too slow when data set is large
# Use a randomized algorithm instead
function findmax_spreadcoord(balls::Vector{AbstractBallTreeNode},
	                         l::Integer,
	                         u::Integer)
	# max_spread = -Inf
	# max_coord = -1
	p = length(balls[1].ball.center)
	return rand(1:p)
	# n = length(balls)
	# for itr in 1:3
	# 	c = rand(1:p)
	# 	min_c, max_c = Inf, -Inf
	# 	for i in 1:n
	# 		if balls[i].ball.center[c] < min_c
	# 			min_c = balls[i].ball.center[c]
	# 		end
	# 		if balls[i].ball.center[c] > max_c
	# 			max_c = balls[i].ball.center[c]
	# 		end
	# 	end
	# 	spread = max_c - min_c
	# 	if spread > max_spread
	# 		max_spread, max_coord = spread, c
	# 	end
	# end
	# return max_coord
end

# Construct a BallTreeNode that wraps the l-th through u-th elements of balls
function wrap(balls::Vector{AbstractBallTreeNode},
	          l::Integer,
	          u::Integer)
	if u == l
		return balls[u]
	else
		# Split left and right
		c = findmax_spreadcoord(balls, l, u)
		m = fld(l + u, 2) # May be off-by-1
		coordselect!(balls, c, m, l, u)

		res = BallTreeNode()

		res.left = wrap(balls, l, m)
		res.left.parent = res

		res.right = wrap(balls, m + 1, u)
		res.right.parent = res

		ball = Ball(Array(Float64, length(balls[1].ball.center)),
			        NaN)
		centroid!(ball.center,
			      res.left.ball.center,
			      res.right.ball.center)
		ball.r = max(euclidean(ball.center,
			                   res.left.ball.center) + res.left.ball.r,
			         euclidean(ball.center,
			         	       res.right.ball.center) + res.right.ball.r)
		res.ball = ball

		return res
	end
end

function k_nearest{T <: Real}(v::Vector{T},
	                          t::BallTree,
	                          k::Integer = 1,
	                          exclude::Integer = -1)
	if k > t.n
		error("k must be smaller than the size of the full data set")
	end

	# Should PriorityQueue's allow duplicate keys?
	pq = PriorityQueue{Int, Float64}(Base.Order.Reverse)
	for i in 1:k
		enqueue!(pq, -i, Inf)
	end

	# Don't search inside a sub-tree whose ball doesn't intersect search_ball
	search_ball = Ball(v, euclidean(v, t.root.ball.center) + t.root.ball.r)

	k_nearest_search!(v, t.root, pq, search_ball, exclude)

	return collect(keys(pq)), collect(values(pq))
end

function k_nearest_search!{T <: Real}(v::Vector{T},
                                      n::AbstractBallTreeNode,
                                      pq::PriorityQueue,
                                      search_ball::Ball,
                                      exclude::Integer)
	if isleaf(n)
		if n.index == exclude
			return
		end
		d = euclidean(search_ball.center, n.ball.center)
		if d < search_ball.r
			dequeue!(pq)
			enqueue!(pq, n.index, d)
			i_max, d_max = peek(pq)
			search_ball.r = d_max
		end
	else
		# Closest point in ball to center of search_ball
		ld = euclidean(n.left.ball.center,
			           search_ball.center) - n.left.ball.r
		rd = euclidean(n.right.ball.center,
			           search_ball.center) - n.right.ball.r
		if ld <= search_ball.r || rd <= search_ball.r
			if ld <= rd
				k_nearest_search!(v, n.left, pq, search_ball, exclude)
				if rd < search_ball.r
					k_nearest_search!(v, n.right, pq, search_ball, exclude)
				end
			else
				k_nearest_search!(v, n.right, pq, search_ball, exclude)
				if ld < search_ball.r
					k_nearest_search!(v, n.left, pq, search_ball, exclude)
				end
			end
		end
	end
	return
end

function inball{T <: Real}(v::Vector{T},
	                       t::BallTree,
	                       r::Real,
	                       exclude::Integer = -1)
	error("Not yet implemented")
end


text: # Implicit active contours


text:
mutable struct NormalSwitches{T1} <: AbstractDistributionSwitches{T1}
    pvec :: ParamVector
    # logNormal :: Bool
end

mutable struct Normal{T1} <: AbstractDistributionLH{T1}
    objId :: ObjectId
    switches :: NormalSwitches{T1}
    mean :: T1
    std :: T1
end

## ------------  Generic

Lazy.@forward Normal.switches (
    ModelParams.get_pvector
);

# function exp_draws!(drawV)
#     for (j, x) in enumerate(drawV)
#         drawV[j] = exp(x);
#     end
# end

function draw(u :: Normal{T1}, nDims, rng :: AbstractRNG) where T1
    drawV = u.mean .+ u.std .* randn(rng, T1, nDims...);
    # is_lognormal(u)  &&  exp_draws!(drawV);
    return drawV
end

function draw(u :: Normal{T1}, rng :: AbstractRNG)  where T1
    drawV = u.mean + u.std * rand(rng, T1);
    # is_lognormal(u)  &&  (drawV = exp(drawV));
    return drawV
end

function quantiles(u :: Normal{T1}, pctV)  where T1
    normalDistr = Distributions.Normal(u.mean, u.std);
    qV = Distributions.quantile.(normalDistr, pctV);
    # is_lognormal(u)  &&  exp_draws!(qV);
    return qV
end

# is_lognormal(switches :: NormalSwitches{T1}) where T1 = switches.logNormal;

Base.eltype(::Normal{T1}) where T1 = T1;

isbounded(u :: Normal{T1}) where T1 = false;

Base.show(io :: IO, um :: Normal{T1}) where T1 =
    print(io, "Normal{$T1}($(um.mean), $(um.std))");

# Validate draws
function check_draws(u :: Normal{T1}, drawM) where T1
    isValid = true;
    # if is_lognormal(u)
    #     isValid = isValid  &&  all(x -> x > zero(T1), drawM);
    # end
    return isValid
end

## ------------  Construction

ModelParams.get_pvector(switches :: NormalSwitches{T1}) where T1 = switches.pvec;
# ModelParams.get_pvector(nd :: Normal{T1}) where T1 = nd.switches.pvec;

init_distribution(objId :: ObjectId, switches :: NormalSwitches{T1}) where T1 =
    init_normal(switches);

function init_normal(switches :: NormalSwitches{T1}) where T1
    pvec = get_pvector(switches);
    return Normal{T1}(switches.pvec.objId, switches,
        param_default_value(pvec, :mean),
        param_default_value(pvec, :std))
end

function init_normal_switches(objId :: ObjectId;
    meanVal = 0.0, stdVal = 1.0, calMean = true, calStd = true)
    T1 = typeof(meanVal);
    pMean = init_mean(; meanVal, calMean);
    pStd = init_std(; stdVal, calStd);
    pvec = ParamVector(objId, [pMean, pStd]);
    return NormalSwitches{T1}(pvec)
end

function init_mean(; descr = "Mean", lsym = "μ", meanVal = 0.0,
    meanLb = -10.0, meanUb = 10.0, calMean = true)
    return Param(:mean, descr, lsym, meanVal, meanVal, meanLb, meanUb, calMean);
end

function init_std(; descr = "StdDev", lsym = "σ", stdVal = 1.0,
    stdLb = 0.0, stdUb = 10.0, calStd = true)
    return Param(:std, descr, lsym, stdVal, stdVal, stdLb, stdUb, calStd);
end

make_test_normal_switches() = init_normal_switches(ObjectId(:none);
    meanVal = 1.2, stdVal = 2.3);
make_test_normal() = init_normal(make_test_normal_switches());

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


text:
module Predictors

import ..Models: FMModel, predict_instance!

using FactorizationMachines: FMMatrix, FMFloat
using FactorizationMachines.Common: sigmoid
using FactorizationMachines.Tasks: PredictorTask, ClassificationTask, RegressionTask

type FMPredictor{T<:PredictorTask}
    task::T
    model::FMModel
end

"""Instance prediction specialized for classification"""
function predict_instance!(predictor::FMPredictor{ClassificationTask},
                           idx::StridedVector{Int64}, x::StridedVector{FMFloat},
                           f_sum::Vector{FMFloat}, sum_sqr::Vector{FMFloat})
    p = predict_instance!(predictor.model, idx, x, f_sum, sum_sqr)
    sigmoid(p)
end

"""Instance prediction specialized for regression"""
function predict_instance!(predictor::FMPredictor{RegressionTask},
                           idx::StridedVector{Int64}, x::StridedVector{FMFloat},
                           f_sum::Vector{FMFloat}, sum_sqr::Vector{FMFloat})
    p = predict_instance!(predictor.model, idx, x, f_sum, sum_sqr)
    max(min(p, predictor.task.target_max), predictor.task.target_min)
end

"""Predicts labels for each column of `X`"""
function predict(predictor::FMPredictor, X::FMMatrix)
    result = zeros(X.n)
    predict!(predictor, X, result)
    result
end

"""Predicts labels for each column of `X` and stores the results into `result`"""
function predict!(predictor::FMPredictor, X::FMMatrix, result::Vector{FMFloat})
    fill!(result, .0)
    f_sum = fill(.0, predictor.model.num_factors)
    sum_sqr = fill(.0, predictor.model.num_factors)
    for c in 1:X.n
        X_nzrange = nzrange(X, c)
        idx = sub(X.rowval, X_nzrange)
        x = sub(X.nzval, X_nzrange)
        result[c] = predict_instance!(predictor, idx, x, f_sum, sum_sqr)
    end
end

end


text:
module CovidNYT

using DocStringExtensions: SIGNATURES, FIELDS, TYPEDEF

const src_path = @__DIR__

"""

# rel_path_covidnyt

Relative path using the StatisticalRethinking src/ directory.

### Example to get access to the data subdirectory
```julia
rel_path_covid("..", "data")
```

"""
rel_path_covidnyt(parts...) = normpath(joinpath(src_path, parts...))

#include("refresh_dataframes.jl")

export
  rel_path_covidnyt

end # module


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

module GeoLearning

using Meshes
using GeoStatsBase

using Tables
using TableOperations
using MLJModelInterface
using Distributions

import GeoStatsBase: solve

# aliases
const TO = TableOperations
const MI = MLJModelInterface

include("pointwise.jl")
include("utils.jl")

export
  PointwiseLearn,
  learn, perform

end


text:
# Version 5.2
txt = "100"
for base = 2:21
    base10 = parse(Int, txt, base)
    println("String $txt in base $base is $base10 in base 10")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Used to allow custom dispatches to `_http_request`
abstract type AbstractBackend end

"""
    AWS.HTTPBackend <: AWS.AbstractBackend

An `HTTPBackend` can hold default `http_options::AbstractDict{Symbol,<:Any}`
to pass to HTTP.jl, which can be overwritten per-request by any `http_options`
supplied there.
"""
struct HTTPBackend <: AbstractBackend
    http_options::AbstractDict{Symbol,<:Any}
end

function HTTPBackend(; kwargs...)
    return if isempty(kwargs)
        HTTPBackend(LittleDict{Symbol,Any}())
    else
        HTTPBackend(LittleDict(kwargs))
    end
end
# populated in `__init__`
const DEFAULT_BACKEND = Ref{AbstractBackend}()

Base.@kwdef mutable struct Request
    service::String
    api_version::String
    request_method::String

    headers::AbstractDict{String,String} = LittleDict{String,String}()
    content::Union{String,Vector{UInt8}} = ""
    resource::String = ""
    url::String = ""

    # Note: User provided `IO` should support seeking in order to support API error handling
    response_stream::Union{IO,Nothing} = nothing

    http_options::AbstractDict{Symbol,<:Any} = LittleDict{Symbol,String}()
    backend::AbstractBackend = DEFAULT_BACKEND[]

    # Deprecated fields
    use_response_type::Bool = false
    return_stream::Union{Bool,Nothing} = nothing
    return_raw::Union{Bool,Nothing} = nothing
    response_dict_type::Union{Type{<:AbstractDict},Nothing} = nothing
end

"""
    submit_request(aws::AbstractAWSConfig, request::Request; return_headers::Bool=false)

Submit the request to AWS.

# Arguments
- `aws::AbstractAWSConfig`: AWSConfig containing credentials and other information for fulfilling the request, default value is the global configuration
- `request::Request`: All the information about making a request to AWS

# Keywords
- `return_headers::Bool=false`: Set to `true` if you want the headers from the response returned back. Only
  used if `request.use_response_type = false`.

# Returns
- `AWS.Response`: A struct containing the response details
"""
function submit_request(aws::AbstractAWSConfig, request::Request; return_headers=nothing)
    aws_response = nothing
    TOO_MANY_REQUESTS = 429
    EXPIRED_ERROR_CODES = ["ExpiredToken", "ExpiredTokenException", "RequestExpired"]
    REDIRECT_ERROR_CODES = [301, 302, 303, 304, 305, 307, 308]
    THROTTLING_ERROR_CODES = [
        "Throttling",
        "ThrottlingException",
        "ThrottledException",
        "RequestThrottledException",
        "TooManyRequestsException",
        "ProvisionedThroughputExceededException",
        "LimitExceededException",
        "RequestThrottled",
        "PriorRequestNotComplete",
    ]

    request.headers["User-Agent"] = user_agent[]
    request.headers["Host"] = HTTP.URI(request.url).host
    stream = @something request.response_stream IOBuffer()

    @repeat 3 try
        credentials(aws) === nothing || sign!(aws, request)

        aws_response = @mock _http_request(request.backend, request, stream)
        response = aws_response.response

        if response.status in REDIRECT_ERROR_CODES
            if HTTP.header(response, "Location") != ""
                request.url = HTTP.header(response, "Location")
            else
                e = HTTP.StatusError(response.status, response)
                throw(AWSException(e, stream))
            end
        end
    catch e
        if e isa HTTP.StatusError
            e = AWSException(e, stream)
        elseif !(e isa AWSException)
            rethrow(e)
        end

        @retry if occursin("Signature expired", e.message)
        end

        # Handle ExpiredToken...
        # https://github.com/aws/aws-sdk-go/blob/v1.31.5/aws/request/retryer.go#L98
        @retry if e isa AWSException && e.code in EXPIRED_ERROR_CODES
            check_credentials(credentials(aws); force_refresh=true)
        end

        # Throttle handling
        # https://github.com/boto/botocore/blob/1.16.17/botocore/data/_retry.json
        # https://docs.aws.amazon.com/general/latest/gr/api-retries.html
        @delay_retry if _http_status(e.cause) == TOO_MANY_REQUESTS ||
            e.code in THROTTLING_ERROR_CODES
        end

        # Handle BadDigest error and CRC32 check sum failure
        @retry if _header(e.cause, "crc32body") == "x-amz-crc32" ||
            e.code in ("BadDigest", "RequestTimeout", "RequestTimeoutException")
        end

        if occursin("Missing Authentication Token", e.message) &&
            aws.credentials === nothing
            return throw(
                NoCredentials(
                    "You're attempting to perform a request without credentials set."
                ),
            )
        end
    end

    if request.use_response_type
        return aws_response
    else
        return legacy_response(request, aws_response; return_headers=return_headers)
    end
end

function _http_request(http_backend::HTTPBackend, request::Request, response_stream::IO)
    http_options = merge(http_backend.http_options, request.http_options)

    # HTTP options such as `status_exception` need to be used when creating the stack
    http_stack = HTTP.stack(;
        redirect=false, retry=false, aws_authorization=false, http_options...
    )

    local buffer
    local response
    try
        @repeat 4 try
            # Use a sacrificial I/O stream so that we only write to the `response_stream`
            # once even with multiple attempted requests. Additionally this works around the
            # HTTP.jl issue (https://github.com/JuliaWeb/HTTP.jl/issues/543) where the
            # `response_stream` is closed automatically. Effectively, this works as if we're
            # not using streaming I/O at all, as we write all data at once, but only
            # returning data via I/O ensures we aren't relying on response's body being
            # populated.
            buffer = Base.BufferStream()

            response = @mock HTTP.request(
                http_stack,
                request.request_method,
                HTTP.URI(request.url),
                HTTP.mkheaders(request.headers),
                request.content;
                require_ssl_verification=false,
                response_stream=buffer,
                http_options...,
            )
        catch e
            # `Base.IOError` is needed because HTTP.jl can often have errors that aren't
            # caught and wrapped in an `HTTP.IOError`
            # https://github.com/JuliaWeb/HTTP.jl/issues/382
            @delay_retry if isa(e, Sockets.DNSError) ||
                isa(e, HTTP.ParseError) ||
                isa(e, HTTP.IOError) ||
                isa(e, Base.IOError) ||
                (isa(e, HTTP.StatusError) && _http_status(e) >= 500)
            end
        end
    finally
        # We're unable to read from the `Base.BufferStream` until it has been closed.
        # HTTP.jl will close passed in `response_stream` keyword. This ensures that it
        # is always closed (e.g. HTTP.jl 0.9.15)
        close(buffer)

        # Transfer the contents of the `BufferStream` into `response_stream` variable.
        write(response_stream, buffer)
    end

    return @mock Response(response, response_stream)
end

_http_status(e::HTTP.StatusError) = e.status
_header(e::HTTP.StatusError, k, d="") = HTTP.header(e.response, k, d)


text:
using ChainRulesTestUtils
using LinearAlgebra
using Zygote: ZygoteRuleConfig

# issue 897

test_rrule(ZygoteRuleConfig(), x->sum(sin, Diagonal(x)), ones(2); rrule_f=rrule_via_ad, check_inferred=false)
test_rrule(ZygoteRuleConfig(), x->sum(sin, Diagonal(x)), rand(3); rrule_f=rrule_via_ad, check_inferred=false)


text:
reponame = "arnavs/InstantiationTest"
tag1 = "master"
tag2 = "v0.2.0"
sha = "1daba6c" # 7 chars
pdir = pwd()

# Clean environment.
isdir(".projects") == false || rm(".projects", recursive=true)

prefix, path = activate_github(reponame)
@test prefix == "master"
@test path == joinpath(pdir, ".projects", "InstantiationTest-master")
@test Base.active_project() == joinpath(path, "Project.toml")

@warn "Shouldn't take too long on new precompilations, etc."
prefix, path = activate_github(reponame, tag=tag1)
@test prefix == "master"
@test path == joinpath(pdir, ".projects", "InstantiationTest-master")
@test Base.active_project() == joinpath(path, "Project.toml")

prefix, path = activate_github(reponame, tag=tag2)
@test prefix == "v0.2.0"
@test path == joinpath(pdir, ".projects", "InstantiationTest-v0.2.0")
@test Base.active_project() == joinpath(path, "Project.toml")

InstantiateFromURL.copy_env(reponame, "v0.2.0", "v0.2.0-edit")
prefix, path = activate_github(reponame, tag="v0.2.0-edit")
@test prefix == "v0.2.0-edit"
@test path == joinpath(pdir, ".projects", "InstantiationTest-v0.2.0-edit")
@test Base.active_project() == joinpath(path, "Project.toml")

prefix, path = activate_github(reponame, sha=sha)
@test prefix == sha
@test path == joinpath(pdir, ".projects", "InstantiationTest-$sha")
@test Base.active_project() == joinpath(path, "Project.toml")


text:
"""
    JacobiP(x, alpha, beta, N)

Evaluate Jacobi polynomial P_n^(α,β)(x)
Note that the order of arguments and return values are different from Jacobi.jl
"""
function JacobiP(x::T, alpha, beta, N) where {T<:Real}
    xp = x
    PL = zeros(N + 1)

    # P₀(x) and P₁(x)
    gamma0 =
        2^(alpha + beta + 1) / (alpha + beta + 1) * gamma(alpha + 1) * gamma(beta + 1) /
        gamma(alpha + beta + 1)
    PL[1] = 1.0 / sqrt(gamma0)
    if N == 0
        return PL[1]
    end
    gamma1 = (alpha + 1) * (beta + 1) / (alpha + beta + 3) * gamma0
    PL[2] = ((alpha + beta + 2) * xp / 2 + (alpha - beta) / 2) / sqrt(gamma1)
    if N == 1
        P = PL[2]
        return P
    end

    # forward recurrence using the symmetry of the recurrence
    aold = 2 / (2 + alpha + beta) * sqrt((alpha + 1) * (beta + 1) / (alpha + beta + 3))
    for i = 1:N-1
        h1 = 2 * i + alpha + beta
        anew =
            2 / (h1 + 2) * sqrt(
                (i + 1) * (i + 1 + alpha + beta) * (i + 1 + alpha) * (i + 1 + beta) /
                (h1 + 1) / (h1 + 3),
            )
        bnew = -(alpha^2 - beta^2) / h1 / (h1 + 2)
        PL[i+2] = 1 / anew * (-aold * PL[i] + (xp - bnew) * PL[i+1])
        aold = anew
    end

    P = PL[N+1]

    return P
end

function JacobiP(x::AbstractArray{T}, alpha, beta, N) where {T<:Real}
    xp = copy(x)
    PL = zeros(N + 1, length(xp))

    # P₀(x) and P₁(x)
    gamma0 =
        2^(alpha + beta + 1) / (alpha + beta + 1) * gamma(alpha + 1) * gamma(beta + 1) /
        gamma(alpha + beta + 1)
    PL[1, :] .= 1.0 / sqrt(gamma0)
    if N == 0
        P = PL[:]
        return P
    end
    gamma1 = (alpha + 1) * (beta + 1) / (alpha + beta + 3) * gamma0
    @. PL[2, :] = ((alpha + beta + 2) * xp / 2 + (alpha - beta) / 2) / sqrt(gamma1)
    if N == 1
        P = PL[N+1, :]
        return P
    end

    # forward recurrence using the symmetry of the recurrence
    aold = 2 / (2 + alpha + beta) * sqrt((alpha + 1) * (beta + 1) / (alpha + beta + 3))
    for i = 1:N-1
        h1 = 2 * i + alpha + beta
        anew =
            2 / (h1 + 2) * sqrt(
                (i + 1) * (i + 1 + alpha + beta) * (i + 1 + alpha) * (i + 1 + beta) /
                (h1 + 1) / (h1 + 3),
            )
        bnew = -(alpha^2 - beta^2) / h1 / (h1 + 2)
        @. PL[i+2, :] = 1 / anew * (-aold * PL[i, :] + (xp - bnew) * PL[i+1, :])
        aold = anew
    end

    P = PL[N+1, :]

    return P
end

function ∂JacobiP(r::T, alpha, beta, N) where {T<:Real}
    dP = 0.0
    if N != 0
        dP = sqrt(N * (N + alpha + beta + 1)) * JacobiP(r, alpha + 1, beta + 1, N - 1)
    end

    return dP
end

function ∂JacobiP(r::AbstractArray{T}, alpha, beta, N) where {T<:Real}
    dP = zero(r)
    if N != 0
        dP .= sqrt(N * (N + alpha + beta + 1)) .* JacobiP(r, alpha + 1, beta + 1, N - 1)
    end

    return dP
end