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

There are two situations under which we add binary constraints to the model in order to create a segmented cost curve
for a technology:
    1. When a technology has tax or investment incentives with upper capacity limits < tech.max_kw
        - first segment(s) have lower slope than last segment
    2. When a technology has multiple cost/size pairs (not implemented yet, used for CHP in v1)
        - we interpolate the slope between the cost/size points, typically with economies of scale pricing
We used to use cost curve segments for when a technology has a non-zero existing_kw by setting the first segment to a
zero cost (and slope) from zero kw to the existing_kw. Instead, we now have dvPurchaseSize >= dvSize - existing_kw.

To avoid unnecessary decision variables (and constraints) we use anonymous variables and register them in the model
manually. With this method the decision variable containers do not have to be rectangular.
For example, (as implemented in v1 of the API) the decision variable dvSegmentSystemSize is indexed on
[p.techs.all, p.Subdivision, p.Seg], where:
    - p.Subdivision is hardcoded to ["CapCost"] - and so it is unnecessary,
    - p.Seg is the maximum number of segments in any cost curve
    - and p.techs.all includes all technologies
We can instead construct dvSegmentSystemSize as follows:
```julia
    for t in p.SegmentedTechs
        dv = "dvSegmentSystemSize" * t
        m[Symbol(dv)] = @variable(m, [1:p.SegSize[t]], base_name=dv)
    end
```
Then, for example, the dvSegmentSystemSize for Wind at segment 2 is accessed via:
```julia
    m[:dvSegmentSystemSizeWind][2]
```
In this fashion only technologies with more than one cost curve segment, and only the necessary amount of segments, are
modeled using the following binary variables and constraints.
"""
function add_cost_curve_vars_and_constraints(m, p; _n="")

    for t in p.techs.segmented
        dv = "dvSegmentSystemSize" * t
        m[Symbol(dv)] = @variable(m, [1:p.n_segs_by_tech[t]], base_name=dv, lower_bound=0)

        dv = "binSegment" * t
        m[Symbol(dv)] = @variable(m, [1:p.n_segs_by_tech[t]], base_name=dv, binary=true)
    end

    ##Constraint (7f)-1: Minimum segment size
        @constraint(m, SegmentSizeMinCon[t in p.techs.segmented, s in 1:p.n_segs_by_tech[t]],
        m[Symbol("dvSegmentSystemSize"*t)][s] >= p.seg_min_size[t][s] * m[Symbol("binSegment"*t)][s]
    )

    ##Constraint (7f)-2: Maximum segment size
    @constraint(m, SegmentSizeMaxCon[t in p.techs.segmented, s in 1:p.n_segs_by_tech[t]],
        m[Symbol("dvSegmentSystemSize"*t)][s] <= p.seg_max_size[t][s] * m[Symbol("binSegment"*t)][s]
    )

    ##Constraint (7g):  Segments add up to system size
    @constraint(m, SegmentSizeAddCon[t in p.techs.segmented],
        sum(m[Symbol("dvSegmentSystemSize"*t)][s] for s in 1:p.n_segs_by_tech[t]) == m[:dvPurchaseSize][t]
    )

    ##Constraint (7h): At most one segment allowed
    @constraint(m, SegmentSelectCon[t in p.techs.segmented],
        sum(m[Symbol("binSegment"*t)][s] for s in 1:p.n_segs_by_tech[t]) <= 1
    )
end

max_stars_repo_name: NREL/REoptLite

max_issues_repo_name: NREL/REoptLite

max_forks_repo_name: NREL/REoptLite


text:
##########################################################################
# Precoding types
type PrecodingResults <: Results
    results::Dict{ASCIIString, Any}
end
PrecodingResults() =
    PrecodingResults(Dict{ASCIIString, Any}())
Base.getindex(p::PrecodingResults, k::ASCIIString) =
    getindex(p.results, k)
Base.setindex!(p::PrecodingResults, v, k::ASCIIString) =
    setindex!(p.results, v, k)

# AuxPrecodingParams is defined in network.jl due to cyclic dependencies
# between files.

##########################################################################
# Standard functions to calculate rates from optimal MSE weights.

# logdet rates assume an ML decoder with perfect CSI-R
calculate_logdet_rates(state) =
    calculate_weighted_logdet_rates(state, ones(length(state.W)))

function calculate_weighted_logdet_rates(state, alphas)
    K = length(state.W)
    ds = Int[ size(state.W[k], 2) for k = 1:K ]; max_d = maximum(ds)

    utilities = zeros(Float64, K, max_d)
    for k = 1:K; if ds[k] > 0
        # W is p.d., so we should only get abs eigenvalues. Numerically we may
        # get some imaginary noise however. Also, numerically the eigenvalues
        # may be less than 1, so we need to handle that to not get negative
        # rates.
        r = alphas[k]*log2(max(1, abs(eigvals(state.W[k]))))

        if ds[k] < max_d
            utilities[k,:] = cat(1, r, zeros(Float64, max_d - ds[k]))
        else
            utilities[k,:] = r
        end
    end; end

    return utilities
end

# MMSE rates assume an MMSE decoder with perfect CSI-R
calculate_MMSE_rates(state) =
    calculate_weighted_MMSE_rates(state, ones(length(state.W)))

function calculate_weighted_MMSE_rates(state, alphas)
    K = length(state.W)
    ds = Int[ size(state.W[k], 1) for k = 1:K ]; max_d = maximum(ds)

    MMSE_rates = zeros(Float64, K, max_d)
    for k = 1:K; if ds[k] > 0
        # Invert W to get MMSE matrix and pick out diagonal elements, to obtain
        # MSE performance of MMSE receiver. Then take log2 of the reciprocal
        # to get the rates for the streams.
        # N.B. it is a little wasteful to take the inverse here, but it makes
        # for cleaner code in the algorithms. Also, these matrices are typically
        # small (like 2-by-2 or 3-by-3), so the complexity isn't too bad.
        E = state.W[k]\eye(state.W[k])
        r = alphas[k]*log2(max(1, abs(1./diag(E))))

        if ds[k] < max_d
            MMSE_rates[k,:] = cat(1, r, zeros(Float64, max_d - ds[k]))
        else
            MMSE_rates[k,:] = r
        end
    end; end

    return MMSE_rates
end

function calculate_allocated_power(state)
    K = length(state.V)
    ds = Int[ size(state.V[k], 2) for k = 1:K ]; max_d = maximum(ds)

    allocated_power = zeros(Float64, K, max_d)

    for k = 1:K
        p = [ vecnorm(state.V[k][:,n])^2 for n = 1:ds[k] ]

        if ds[k] < max_d
            allocated_power[k,:] = cat(1, p, zeros(Float64, max_d - ds[k]))
        else
            allocated_power[k,:] = p
        end
    end

    return allocated_power
end

##########################################################################
# Initialization functions
function initial_receivers(channel::SinglecarrierChannel, Ps, sigma2s, ds,
    assignment, aux_params)

    U = Array(Matrix{Complex128}, channel.K)

    if haskey(aux_params, "initial_receivers")
        if aux_params["initial_receivers"] == "eigendirection"
            for i = 1:channel.I; for k in served_MS_ids(i, assignment)
                Utmp, _, _ = svd(channel.H[k,i])
                U[k] = Utmp[:, 1:ds[k]]
            end; end

            return U
        elseif aux_params["initial_receivers"] == "white"
            for i = 1:channel.I; for k in served_MS_ids(i, assignment)
                U[k] = eye(Complex128, channel.Ns[k], ds[k])
            end; end

            return U
        elseif aux_params["initial_receivers"] == "zeros"
            for i = 1:channel.I; for k in served_MS_ids(i, assignment)
                U[k] = zeros(Complex128, channel.Ns[k], ds[k])
            end; end

            return U
        end
    end

    # Default: decoders taken as columns of DFT matrix
    for i = 1:channel.I; for k in served_MS_ids(i, assignment)
        U[k] = fft(eye(channel.Ns[k], ds[k]), 1)
    end; end

    return U
end

function initial_MSE_weights(channel::SinglecarrierChannel, Ps, sigma2s, ds,
    assignment, aux_params)

    # So far only the trivial initial point available
    return [ Hermitian(eye(Complex128, ds[k])) for k = 1:channel.K ]
end

function initial_precoders(channel::SinglecarrierChannel, Ps, sigma2s, ds,
    assignment, aux_params)

    V = Array(Matrix{Complex128}, channel.K)

    if haskey(aux_params, "initial_precoders")
        if aux_params["initial_precoders"] == "eigendirection"
            for i = 1:channel.I
                served = served_MS_ids(i, assignment)
                Kc = length(served)

                for k in served
                    _, _, Vtmp = svd(channel.H[k,i])
                    V[k] = sqrt(Ps[i]/Kc)*Vtmp[:, 1:ds[k]]/vecnorm(Vtmp[:, 1:ds[k]])
                end
            end

            return V
        elseif aux_params["initial_precoders"] == "white"
            for i = 1:channel.I
                served = served_MS_ids(i, assignment)
                Kc = length(served)

                for k in served
                    V[k] = sqrt(Ps[i]/(ds[k]*Kc))*eye(Complex128, channel.Ms[i], ds[k])
                end
            end

            return V
        elseif aux_params["initial_precoders"] == "zeros"
            for i = 1:channel.I
                served = served_MS_ids(i, assignment)
                Kc = length(served)

                for k in served
                    V[k] = zeros(Complex128, channel.Ms[i], ds[k])
                end
            end

            return V
        end
    end

    # Default: precoders taken as columns of DFT matrix
    for i = 1:channel.I
        served = served_MS_ids(i, assignment)
        Kc = length(served)

        for k in served
            V[k] = sqrt(Ps[i]/(channel.Ms[i]*ds[k]*Kc))*fft(eye(channel.Ms[i], ds[k]), 1)
        end
    end

    return V
end

ext: jl


text:
############################################################################################
# Constants
"RNG for sampling based solutions"
const _rng = Random.GLOBAL_RNG   # shorthand
Random.seed!(_rng, 1)

"Tolerance for stochastic solutions"
const _TOL = 1.0e-6

"Number of samples"
N = 10^3

############################################################################################
"Normalize Vector inplace and return normalizing constant"
@inline function normalize!(vec::AbstractVector)
    norm = sum(vec)
    vec ./= norm
    return norm
end
"Faster max function for Viterbi algorithm in HMM. Not exported."
function vec_maximum(vec::AbstractVector)
    m = vec[1]
    @inbounds for i = Base.OneTo(length(vec))
        if vec[i] > m
            m = vec[i]
        end
    end
    return m
end

############################################################################################
"Compute Stationary distribution for given transition matrix"
function get_stationary!(Transition::AbstractMatrix{T}) where T<:Real
    # From: https://github.com/QuantEcon/QuantEcon.jl/blob/f454d4dfbaf52f550ddd52eff52471e4b8fddb9d/src/markov/mc_tools.jl
    # Grassmann-Taksar-Heyman (GTH) algorithm (Grassmann, Taksar, and Heyman 1985)
    n = size(Transition, 1)
    x = zeros(T, n)

    @inbounds for k in 1:n-1
        scale = sum(Transition[k, k+1:n])
        if scale <= zero(T)
            # There is one (and only one) recurrent class contained in
            # {1, ..., k};
            # compute the solution associated with that recurrent class.
            n = k
            break
        end
        Transition[k+1:n, k] /= scale

        for j in k+1:n, i in k+1:n
            Transition[i, j] += Transition[i, k] * Transition[k, j]
        end
    end

    # backsubstitution
    x[n] = 1
    @inbounds for k in n-1:-1:1, i in k+1:n
        x[k] += x[i] * Transition[i, k]
    end

    # normalisation
    x /= sum(x)

    return x
end
get_stationary(Transition::AbstractMatrix{T}) where {T<:Real} = get_stationary!(copy(Transition))

############################################################################################
"Check if ancestors are all in correct order - Not that first column can have different order as ancestors not resampled at 0."
function check_ancestors(ancestors::AbstractMatrix{I}) where {I<:Integer}
    for Ncol in 2:size(ancestors, 2)
        for Nrow in Base.OneTo(size(ancestors, 1))
            if ancestors[Nrow, Ncol] != Nrow
                return false
            end
        end
    end
    return true
end

############################################################################################
include("hmm.jl")
include("hmmMV.jl")
include("hsmm.jl")
include("HigherOrderMarkov.jl")

############################################################################################
resamplemethods = [Systematic(), Residual(), Stratified(), BaytesFilters.Multinomial()]
referencemethods = [Conditional(), Ancestral(), Marginal()]
objectives = [deepcopy(markov_objective), deepcopy(markov_MV_objective), deepcopy(semimarkov_objective)]
generating = [UpdateFalse(), UpdateTrue()]


text:
using ManifoldsBase
using Test
using ManifoldsBase:
    combine_allocation_promotion_functions, allocation_promotion_function, ℝ

struct AllocManifold{𝔽} <: Manifold{𝔽} end
AllocManifold() = AllocManifold{ℝ}()

function ManifoldsBase.exp!(::AllocManifold, v, x, y)
    v[1] .= x[1] .+ y[1]
    v[2] .= x[2] .+ y[2]
    v[3] .= x[3] .+ y[3]
    return v
end

@testset "Allocation" begin
    a = [[1.0], [2.0], [3.0]]
    b = [[2.0], [3.0], [-3.0]]
    M = AllocManifold()
    v = exp(M, a, b)
    @test v ≈ [[3.0], [5.0], [0.0]]
    @test allocate(([1.0], [2.0])) isa Tuple{Vector{Float64},Vector{Float64}}
    @test allocate(([1.0], [2.0]), Int) isa Tuple{Vector{Int},Vector{Int}}
    @test allocate([[1.0], [2.0]]) isa Vector{Vector{Float64}}
    @test allocate([[1.0], [2.0]], Int) isa Vector{Vector{Int}}

    a1 = allocate([1], 2, 3)
    @test a1 isa Matrix{Int}
    @test size(a1) == (2, 3)
    a2 = allocate([1], (2, 3))
    @test a2 isa Matrix{Int}
    @test size(a2) == (2, 3)
    a3 = allocate([1], Float64, 2, 3)
    @test a3 isa Matrix{Float64}
    @test size(a3) == (2, 3)
    a4 = allocate([1], Float64, (2, 3))
    @test a4 isa Matrix{Float64}
    @test size(a4) == (2, 3)

    @test allocation_promotion_function(M, exp, (a, b)) === identity
    @test combine_allocation_promotion_functions(identity, identity) === identity
    @test combine_allocation_promotion_functions(identity, complex) === complex
    @test combine_allocation_promotion_functions(complex, identity) === complex
    @test combine_allocation_promotion_functions(complex, complex) === complex

    @test number_eltype([2.0]) === Float64
    @test number_eltype([[2.0], [3]]) === Float64
    @test number_eltype([[2], [3.0]]) === Float64
    @test number_eltype([[2], [3]]) === Int
    @test number_eltype(([2.0], [3])) === Float64
    @test number_eltype(([2], [3.0])) === Float64
    @test number_eltype(([2], [3])) === Int
    @test number_eltype(Any[[2.0], [3.0]]) === Float64
end


text:
#=
VecE2: a 2d euclidean vector
=#

immutable VecE2 <: VecE
    x :: Float64
    y :: Float64

    VecE2() = new(0.0,0.0)
    VecE2(x::Real, y::Real) = new(x,y)
end

polar(r::Real, θ::Real) = VecE2(r*cos(θ), r*sin(θ))

Base.length(::VecE2) = 2
Base.copy(a::VecE2) = VecE2(a.x, a.y)
Base.convert(::Type{Vector{Float64}}, a::VecE2) = [a.x, a.y]
function Base.convert{R<:Real}(::Type{VecE2}, a::AbstractArray{R})
    @assert(length(a) == 2)
    VecE2(a[1], a[2])
end
Base.show(io::IO, a::VecE2) = @printf(io, "VecE2(%.3f, %.3f)", a.x, a.y)

@compat Base.:+(a::VecE2, b::Real)  = VecE2(a.x+b, a.y+b)
@compat Base.:+(b::Real,  a::VecE2) = VecE2(a.x+b, a.y+b)
@compat Base.:+(a::VecE2, b::VecE2) = VecE2(a.x+b.x, a.y+b.y)

@compat Base.:-(a::VecE2)  = VecE2(-a.x, -a.y)
@compat Base.:-(a::VecE2, b::Real)  = VecE2(a.x-b, a.y-b)
@compat Base.:-(a::VecE2, b::VecE2) = VecE2(a.x-b.x, a.y-b.y)

@compat Base.:*(a::VecE2, b::Real) = VecE2(a.x*b, a.y*b)
@compat Base.:*(b::Real, a::VecE2) = VecE2(a.x*b, a.y*b)

@compat Base.:/(a::VecE2, b::Real) = VecE2(a.x/b, a.y/b)

@compat Base.:^(a::VecE2, b::Integer) = VecE2(a.x^b, a.y^b)
@compat Base.:^(a::VecE2, b::AbstractFloat) = VecE2(a.x^b, a.y^b)

@compat Base.:%(a::VecE2, b::Real) = VecE2(a.x%b, a.y%b)

@compat Base.:(==)(a::VecE2, b::VecE2) = isequal(a.x, b.x) && isequal(a.y, b.y)
Base.isequal(a::VecE2, b::VecE2) = isequal(a.x, b.x) && isequal(a.y, b.y)

Base.isfinite(a::VecE2) = isfinite(a.x) && isfinite(a.y)
Base.isinf(a::VecE2) = isinf(a.x) || isinf(a.y)
Base.isnan(a::VecE2) = isnan(a.x) || isnan(a.y)

Base.round(a::VecE2) = VecE2(round(a.x), round(a.y))
Base.floor(a::VecE2) = VecE2(floor(a.x), floor(a.y))
Base.ceil(a::VecE2) = VecE2(ceil(a.x), ceil(a.y))
Base.trunc(a::VecE2) = VecE2(trunc(a.x), trunc(a.y))
Base.clamp(a::VecE2, lo::Real, hi::Real) = VecE2(clamp(a.x, lo, hi), clamp(a.y, lo, hi))

Base.abs(a::VecE2) = hypot(a.x, a.y)
Base.hypot(a::VecE2) = hypot(a.x, a.y)
Base.abs2(a::VecE2) = a.x*a.x + a.y*a.y
function Base.norm(a::VecE2)
    m = abs(a)
    VecE2(a.x/m, a.y/m)
end

Base.atan2(a::VecE2) = atan2(a.y, a.x)

dist(a::VecE2, b::VecE2) = hypot(a.x-b.x, a.y-b.y)
function dist2(a::VecE2, b::VecE2)
    Δx = a.x-b.x
    Δy = a.y-b.y
    Δx*Δx + Δy*Δy
end

Base.dot(a::VecE2, b::VecE2) = a.x*b.x + a.y*b.y
Base.cross(a::VecE2, b::VecE2) = a.x*b.y - a.y*b.x
proj(a::VecE2, b::VecE2, ::Type{Float64}) = (a.x*b.x + a.y*b.y) / hypot(b.x, b.y) # dot(a,b) / |b|
function proj(a::VecE2, b::VecE2, ::Type{VecE2})
    # dot(a,b) / dot(b,b) ⋅ b
    s = (a.x*b.x + a.y*b.y) / (b.x*b.x + b.y*b.y)
    VecE2(s*b.x, s*b.y)
end

lerp(a::VecE2, b::VecE2, t::Real) = VecE2(a.x + (b.x-a.x)*t, a.y + (b.y-a.y)*t)

Base.rot180(a::VecE2) = VecE2(-a.x, -a.y)
Base.rotl90(a::VecE2) = VecE2(-a.y,  a.x)
Base.rotr90(a::VecE2) = VecE2( a.y, -a.x)
function rot(a::VecE2, θ::Float64)

    #=
    Rotates counter-clockwise about origin
    =#

    c = cos(θ)
    s = sin(θ)

    VecE2(a.x*c - a.y*s, a.x*s+a.y*c)
end


text:
"""
   defineBold(df, signal; method = :Max, peak_to = :all, nmax = 4, window = 4:10, nanmath = true)
It adjust the Bold signal per genotype, mouse, stim, area and concentration
according to different algorythms define by method option.
:Peak = take all values after the peak until the index define by peak_to
    if peak_to == :all takes all values available (goes over every presentation
    if the presentation where not averaged before)
:Max = take the highest n values defined by nmax, default = 4
:Window = averages over the time points specified by window, default = 4:10
nanmath: if true (default) uses safemean and safe sem to calculate mean and sem,
else it uses regular mean and sem

"""
function define_bold(df0, signal; av_pres = false, method = :Peak, peak_to = :all, nmax = 4, window = 4:10, nanmath = true)
    fmean, fsem = nanmath ? (safemean, safesem) : (mean, sem)
    if method == :Peak
        if av_pres
            #average response over presentation
            df1 = average_presentation(df0,signal; nanmath = nanmath)
            signal = :Bold
            gd1 = groupby(df1,[:Genotype, :Mouse, :Stim, :Area, :Concentration])
        else
            #change groupby to find peak for each presentation
            gd1 = groupby(df0,[:Genotype, :Mouse, :Stim, :Area, :Concentration,:Repetition_nr])
        end
        #find peak in the bold signal and take all value from there until peak_to
        df2 = peak_bold(gd1, signal; peak_to = peak_to)
        df2.Bold = df2[:,signal]
    elseif method == :Max
        #average response over all presentation
        df1 = average_presentation(df0,signal; nanmath = nanmath)
        gd1 = groupby(df1,[:Genotype, :Mouse, :Stim, :Area, :Concentration])
        signal = :Bold
        #find the first n max values in the bold signal
        df2 = combine(gd1, signal => (s -> max_bold(s; max = nmax, nanmath = nanmath)) => :Bold)
    elseif method == :Window
        #average all bold values per mouse and concentration in a time window
        df1 = average_presentation(df0,signal; nanmath = nanmath)
        gd1 = groupby(df1,[:Genotype, :Mouse, :Stim, :Area, :Concentration])
        signal = :Bold
        df2 = combine(gd1, [signal, :Time] =>((s,t) -> window_bold(s,t; window = window, nanmath = nanmath)) => :Bold)
    elseif method == :None
        df2 = filter(r->r.Time in window, df0)
        df2.Bold = df2[:,signal]
    else
        error("unrecognised bold processing method:/n options are :Peak, :Max, :Timewindow")
    end
    return df2
end

function average_presentation(df0,signal; nanmath = true)
    fmean, fsem = nanmath ? (safemean, safesem) : (mean, sem)
    df0[:,signal] = [ismissing(x) ? NaN : x for x in df0[:,signal]]
    df1 = combine(groupby(df0, [:Time,:Concentration,:Area,:Stim,:Mouse,:Genotype]), signal => fmean => :Bold)
    return df1
end

function peak_bold(gd, signal; peak_to = :all)
    combine(gd) do dd
        #any(ismissing,dd[:,signal]) && println("got NaN ", dd.Mouse[1], " ", dd.Area[1], " ", dd.Stim[1], " ", dd.Concentration[1])
        if all(ismissing.(dd[:,signal]))
            return filter(r -> r.Area == "NoData",dd)
        else
            i = findmax(abs.(skipmissing(dd[:,signal])))[2]
            return peak_to == :all ? dd[i:end,:] : dd[i:peak_to,:]
        end
    end
end

function max_bold(v; max = 4, nanmath = true)
    fmean = nanmath ? safemean : mean
    v1 = [ismissing(x) ? NaN : x for x in v]
    v2 = sort(v1,rev = true)
    if length(v1) <= max
        return fmean(v2)
    else
        return fmean(v2[1:max])
    end
end

function window_bold(b,t; window = 4:10, nanmath = true)
    fmean = nanmath ? safemean : mean
    index = [x in window for x in t]
    b1 = [ismissing(x) ? NaN : x for x in b]
    fmean(b1[index])
end

"""
    plot_bold(df0,x; nanmath = true)
Function to check the average response over concentration per area.
df0 has to be a dataframe preprocessed with define_bold,
containing bold signal in a column named :Bold.
returns a vector of plots of bold over concentration per each area,
a dataframe with the mean and sem values of bold per area
and a filtered list of the areas that contain sufficient data for further analysis
"""

function plot_bold(df0; nanmath = true, fconc = true, lim = (-0.4,0.7))
    fmean = nanmath ? safemean : mean
    fsem = nanmath ? safesem : sem
    #calculate mean repsonse over concentrations per area
    df = combine(groupby(df0,[:Area, :Stim, :Concentration, :Genotype]),  :Bold .=> [fmean, fsem] .=>[:mean, :sem])
    # identify regions with missing concentration info
    res = combine(groupby(df,[:Area]), :mean => (m-> !(any(isnan.(m)))) => :to_keep)
    to_keep = filter(r -> r.to_keep, res)[:,:Area]
    # prepare a grid plot with the results for each regions
    df.Color = [x ? :red : :black for x in df.Stim]
    posdict = Dict(v=>p for (v,p) in zip(union(df.Concentration), collect(1:length(union(df.Concentration)))))
    df.Pos = [get(posdict, x,0) for x in df.Concentration]
    sort!(df,[:Stim,:Pos])
    fconc && filter!(r -> r.Concentration != 0.0001 ,df)
    df1 = filter(r-> r.Genotype == "sert",df)
    plt = []
        combine(groupby(df1, :Area)) do dd
            pp = @df dd scatter(:Pos, :mean, yerror = :sem, group = :Stim, color = :Color,
                title = :Area[1], legend = false, ylims = lim)
            @df dd plot!(:Pos, :mean, color = :Color, group = :Stim)
            Plots.abline!(0,0,color = :black, linestyle = :dash)
            push!(plt,pp)
        end
    return plt, df, to_keep
end

"""
    odour_response(df0, to_keep)
Use control data to identify which areas show a dose dependent response to
odour concentration with a mixed model.
df0 has to be a dataframe preprocessed with define_bold,
    containing bold signal in a column named :Bold.
filt is a vector of string specifying which area to keep and which to filter
    out from the analysis, default is :none which takes all areas
"""

function odour_response(df0; farea = :none, fconc = true)
    # use only control data to determine odour concentration dose dependency
    df1 = filter(r -> !r.Stim, df0)
    #remove lowest concentration
    fconc && filter!(r -> r.Concentration != 0.0001, df1)
    # filters area with insufficient response, identified from plot_bold function
    farea != :none && filter!(r -> r.Area in farea, df1)
    # identify if the model requires to incorporate time
    if iscolumn(df0,:Time)
        f0 = @formula(Bold ~ 1 + Time + Concentration + (1+Time|Mouse))
        f1 = @formula(Bold ~ 1 + Time + Concentration*Area + (1+Time|Mouse))
    else
        f0 = @formula(Bold ~ 1 + Concentration + (1|Mouse))
        f1 = @formula(Bold ~ 1 + Concentration*Area + (1|Mouse))
    end
    #replace NaN values used for plotting to missing for model
    # any(isnan.(df1.Bold)) && (df1.Bold = [isnan(x) ? missing : x for x in df1.Bold])
    # model fit and test
    m0 = fit(MixedModel, f0, df1)
    m1 = fit(MixedModel, f1, df1)
    mp = MixedModels.likelihoodratiotest(m0,m1)
    # identifies area that significantly respond to odour concentration
    res = DataFrame(Area = coefnames(m1), P = m1.pvalues, Z = m1.beta ./ m1.stderror)
    res2 = filter!(r -> occursin("Concentration & Area: ", r.Area), res)
    res2.Area = replace.(res2.Area, "Concentration & Area: " => "")
    to_keep = filter(r -> r.P < 0.05, res2)[:,:Area]

    return m1, mp, to_keep
end

"""
    stim_interaction(df0, farea; fconc = :none)
Use  data to identify which areas show an interaction between dose dependent
response to odour concentration and stimulation with a mixed model.
df0 has to be a dataframe preprocessed with define_bold,
    containing bold signal in a column named :Bold.
farea is a vector of string specifying which areas showed a significant dose
dependent response to odour concetration, found using odour_response function
fconc allows to remove lowest concetration from the dataset
"""
function stim_interaction(df0, farea; fconc = true, gamma = false)
    # filters area lacking concentration response identified by odour_response function
    df1 = filter(r -> r.Genotype == "sert" && r.Area in farea, df0)
    fconc && filter!(r -> r.Concentration != 0.0001, df0)
    # identify if the model requires to incorporate time
    if iscolumn(df0,:Time)
        f1 = @formula(Bold ~ 1 + Time + Concentration*Area + Stim + (1+Time|Mouse) + (1+Time|Area))
        f2 = @formula(Bold ~ 1 + Time + Concentration + Area + Stim + Concentration & Area + Stim & Concentration + Concentration & Area & Stim + (1+Time|Mouse) + (1+Time|Area))
        f3 = @formula(Bold ~ 1 + Time + Concentration*Area*Stim + (1+Time|Mouse) + (1+Time|Area))
        # f1 = @formula(Bold ~ 1 + Time + Concentration*Area + (1+Time|Mouse) + (1+Time|Area))
        # f2 = @formula(Bold ~ 1 + Time + Concentration + Area + Concentration & Area + Area & Stim + Concentration & Area & Stim + (1+Time|Mouse) + (1+Time|Area))
        # f3 = @formula(Bold ~ 1 + Time + Concentration*Area*Stim + (1+Time|Mouse) + (1+Time|Area))
    else
        f1 = @formula(Bold ~ 1 + Concentration*Area + Stim + (1|Mouse) + (1|Area))
        f2 = @formula(Bold ~ 1 + Concentration + Area + Stim + Concentration & Area + Stim & Concentration + Concentration & Area & Stim + (1|Mouse) + (1|Area))#(1 + Concentration + Area + Stim|Mouse))
        f3 = @formula(Bold ~ 1 + Concentration*Area*Stim + (1|Mouse) + (1|Area))
    end
    #replace NaN values used for plotting to missing for model
    # df1.Bold = [isnan(x) ? missing : x for x in df1.Bold]
    if gamma
        m1 = GeneralizedLinearMixedModel(f1, df1, Gamma())
        m2 = GeneralizedLinearMixedModel(f2, df1, Gamma())
        m3 = GeneralizedLinearMixedModel(f3, df1, Gamma())
    else
        m1 = fit(MixedModel, f1,df1)
        m2 = fit(MixedModel, f2,df1)
        m3 = fit(MixedModel, f3,df1)
    end
    p_simple = MixedModels.likelihoodratiotest(m1,m2)
    p_complex = MixedModels.likelihoodratiotest(m2,m3)
    res = model_res(m2)
    filter!(r -> occursin("Stim", r.Area), res)
    sort!(res,:P)
    return df1, m2, m3, p_simple, p_complex, res
end

function analysis_bold(df0, signal; fconc = true, av_pres = true, method = :Peak,
        peak_to = :all, nmax = 6, window = 4:12,
        nanmath = true, gamma = false, lim = (-0.4,0.7))
    #adjust bold signal
    bold_df = define_bold(df0,signal; method = method, av_pres = av_pres, peak_to = peak_to, nmax = nmax, window = window, nanmath = nanmath)
    if method == :None
        nanmath ? (fmean = safemean) : (fmean = mean)
        dfplt = combine(groupby(bold_df,[:Area, :Stim, :Concentration, :Genotype, :Mouse]),  :Bold => (b-> fmean(skipmissing(b))) => :Bold)
    else
        dfplt = bold_df
    end
    plt, plt_df, to_keep = plot_bold(dfplt; nanmath = true, fconc = fconc, lim = lim)
    display(plot(plt..., size = (900*1.5,1200*1.5)))#layout = (6,3))
    ## find areas that respond to concentration
    m1, mp, to_keep2 = odour_response(bold_df; farea = to_keep, fconc = fconc)
    # assess stim concentration interaction
    stim_df, m2, m3, p_simple, p_complex, res = stim_interaction(bold_df, to_keep2; fconc = fconc, gamma = gamma)
    show(m2)
    show(p_simple)
    show(m3)
    show(p_complex)
    return bold_df, plt, plt_df, to_keep, m1, mp, to_keep2, stim_df, m2, m3, p_simple, p_complex, res
end


text:
"""
    to_hue(B::SSHBinner[, s=1.0, v=1.0, a=1.0])

Maps each bin filling to a hue relative to the maximum bin filling.
"""
function to_hue(B::SSHBinner, s=1.0, v=1.0, a=1.0)
    # hue is cyclic, so the maximum should be < 1.0 (360)
    bins = 0.8 * B.bins ./ maximum(B.bins)
    map(bins) do k
        HSVA(360k, s, v, a)
    end
end

"""
    to_hue2(B::SSHBinner[, s=1.0, v=1.0, a=1.0])

Maps each bin filling to a hue relative to the maximum and minimum bin filling.
"""
function to_hue2(B::SSHBinner, s=1.0, v=1.0, a=1.0)
    # hue is cyclic, so the maximum should be < 1.0 (360)
    bins = 0.8 * (B.bins .- minimum(B.bins)) ./ (maximum(B.bins) - minimum(B.bins))
    map(bins) do k
        HSVA(360k, s, v, a)
    end
end

function to_alpha_hue(B::SSHBinner, s=1.0, v=1.0)
    # hue is cyclic, so the maximum should be < 1.0 (360)
    bins = 0.8 * (B.bins .- minimum(B.bins)) ./ (maximum(B.bins) - minimum(B.bins))
    map(bins) do k
        HSVA(360k, s, v, min(1.0, 2k))
    end
end

"""
    to_alpha(B::SSHBinner[, r=0.2, g=0.4, b=0.8])

Maps each bin filling to the alpha level of a color (given through r, g, b).
"""
function to_alpha(B::SSHBinner, r=0.2, g=0.4, b = 0.8)
    bins = B.bins ./ maximum(B.bins)
    map(bins) do k
        RGBA(r, g, b, k)
    end
end


text:
export fista

mutable struct FISTA{rT <: Real, vecT <: Union{AbstractVector{rT}, AbstractVector{Complex{rT}}}, matA, matAHA} <: AbstractLinearSolver where {rT}
  A::matA
  AᴴA::matAHA
  reg::Regularization
  x::vecT
  x₀::vecT
  xᵒˡᵈ::vecT
  res::vecT
  ρ::rT
  t::rT
  tᵒˡᵈ::rT
  iterations::Int64
  relTol::rT
  normalizeReg::Bool
  regFac::rT
  norm_x₀::rT
  rel_res_norm::rT
  verbose::Bool
end

"""
    FISTA(A, x::vecT=zeros(eltype(A),size(A,2))
          ; reg=nothing, regName=["L1"], λ=[0.0], kargs...)

creates a `FISTA` object for the system matrix `A`.

# Arguments
* `A`                       - system matrix
* `x::vecT`                 - array with the same type and size as the solution
* (`reg=nothing`)           - regularization object
* (`regName=["L1"]`)        - name of the Regularization to use (if reg==nothing)
* (`AᴴA=A'*A`)              - specialized normal operator, default is `A'*A`
* (`λ=0`)                   - regularization paramter
* (`ρ=1`)                   - step size for gradient step
* (`normalize_ρ=false`)     - normalize step size by the maximum eigenvalue of `AᴴA`
* (`t=1.0`)                 - parameter for predictor-corrector step
* (`relTol::Float64=1.e-5`) - tolerance for stopping criterion
* (`iterations::Int64=50`)  - maximum number of iterations
"""
function FISTA(A, x::AbstractVector{T}=Vector{eltype(A)}(undef,size(A,2)); reg=nothing, regName=["L1"]
              , AᴴA=A'*A
              , λ=0
              , ρ=0.95
              , normalize_ρ=true
              , t=1
              , relTol=eps(real(T))
              , iterations=50
              , normalizeReg=false
              , verbose = false
              , kargs...) where {rT <: Real, T <: Union{rT, Complex{rT}}}

  if reg == nothing
    reg = Regularization(regName, λ, kargs...)
  end

  x₀   = similar(x)
  xᵒˡᵈ = similar(x)
  res  = similar(x)
  res[1] = Inf # avoid spurious convergence in first iterations

  if normalize_ρ
    ρ /= abs(power_iterations(AᴴA))
  end

  return FISTA(A, AᴴA, vec(reg)[1], x, x₀, xᵒˡᵈ, res, rT(ρ),rT(t),rT(t),iterations,rT(relTol),normalizeReg,one(rT),one(rT),rT(Inf),verbose)
end

"""
    init!(it::FISTA, b::vecT
              ; A=solver.A
              , x::vecT=similar(b,0)
              , t::Float64=1.0) where T

(re-) initializes the FISTA iterator
"""
function init!(solver::FISTA{rT,vecT,matA,matAHA}, b::vecT
              ; x::vecT=similar(b,0)
              , t=1
              ) where {rT,vecT,matA,matAHA}

  solver.x₀ .= adjoint(solver.A) * b
  solver.norm_x₀ = norm(solver.x₀)

  if isempty(x)
    solver.x .= solver.ρ .* solver.x₀
  else
    solver.x .= x
  end

  solver.t = t
  solver.tᵒˡᵈ = t
  # normalization of regularization parameters
  if solver.normalizeReg
    solver.regFac = norm(solver.x₀,1)/length(solver.x₀)
  else
    solver.regFac = 1
  end
end

"""
    solve(solver::FISTA, b::Vector)

solves an inverse problem using FISTA.

# Arguments
* `solver::FISTA`                     - the solver containing both system matrix and regularizer
* `b::vecT`                           - data vector
* `A=solver.A`                        - operator for the data-term of the problem
* (`startVector::vecT=similar(b,0)`)  - initial guess for the solution
* (`solverInfo=nothing`)              - solverInfo object

when a `SolverInfo` objects is passed, the residuals are stored in `solverInfo.convMeas`.
"""
function solve(solver::FISTA, b; A=solver.A, startVector=similar(b,0), solverInfo=nothing, kargs...)
  # initialize solver parameters
  init!(solver, b; x=startVector)

  # log solver information
  solverInfo != nothing && storeInfo(solverInfo,solver.x,solver.res_norm)

  # perform FISTA iterations
  for (iteration, item) = enumerate(solver)
    solverInfo != nothing && storeInfo(solverInfo,solver.x,solver.res_norm)
  end

  return solver.x
end

"""
  iterate(it::FISTA, iteration::Int=0)

performs one fista iteration.
"""
function iterate(solver::FISTA, iteration::Int=0)
  if done(solver, iteration) return nothing end

  solver.xᵒˡᵈ .= solver.x

  # calculate residuum and do gradient step
  # solver.x .-= solver.ρ .* (solver.AᴴA * solver.x .- solver.x₀)
  mul!(solver.res, solver.AᴴA, solver.xᵒˡᵈ)
  solver.res .-= solver.x₀
  solver.x .-= solver.ρ .* solver.res

  solver.rel_res_norm = norm(solver.res) / solver.norm_x₀
  solver.verbose && println("Iteration $iteration; rel. residual = $(solver.rel_res_norm)")

  # the two lines below are equivalent to the ones above and non-allocating, but require the 5-argument mul! function to implemented for AᴴA, i.e. if AᴴA is LinearOperator, it requires LinearOperators.jl v2
  # mul!(solver.x, solver.AᴴA, solver.xᵒˡᵈ, -solver.ρ, 1)
  # solver.x .+= solver.ρ .* solver.x₀

  # proximal map
  solver.reg.prox!(solver.x, solver.regFac*solver.ρ*solver.reg.λ; solver.reg.params...)

  # predictor-corrector update
  solver.tᵒˡᵈ = solver.t
  solver.t = (1 + sqrt(1 + 4 * solver.tᵒˡᵈ^2)) / 2
  solver.x .+= ((solver.tᵒˡᵈ-1)/solver.t) .* (solver.x .- solver.xᵒˡᵈ)

  # return the residual-norm as item and iteration number as state
  return solver, iteration+1
end

@inline converged(solver::FISTA) = (solver.rel_res_norm < solver.relTol)

@inline done(solver::FISTA,iteration) = converged(solver) || iteration>=solver.iterations


text:
function local_update!(model::AshkinTeller, T::Real, Jsigma, Jtau, K)
    rng = model.rng
    nsites = numsites(model)
    nbonds = numbonds(model)
    mbeta = -1.0/T

    @inbounds for kind in 1:2
        Js = ifelse(kind==1, Jsigma, Jtau)
        dual_kind = 3-kind
        for site in 1:nsites
            center = model.spins[kind, site]
            dual_center = model.spins[dual_kind, site]
            de = 0.0
            for (n,b) in neighbors(model, site)
                de += 2center * model.spins[kind, n] * (Js[bondtype(model,b)]
                                                        + K[bondtype(model,b)]*model.spins[dual_kind, n]*dual_center)
            end
            if rand(rng) < exp(mbeta*de)
                model.spins[kind, site] *= -1
            end
        end
    end

    return nothing
end

function SW_update!(model::AshkinTeller, T::Real, Jsigma, Jtau, K)
    rng = model.rng
    m2b = -2.0/T
    spins = model.spins
    nsites = numsites(model)
    nbonds = numbonds(model)
    nbt = numbondtypes(model)
    ps = zeros(nbt,2)
    @inbounds for kind in 1:2
        uf = UnionFind(nsites)
        Js = ifelse(kind==1, Jsigma, Jtau)
        ps[:,1] .= -expm1.(m2b.*(Js .+ K))
        ps[:,2] .= -expm1.(m2b.*(Js .- K))
        dual_kind = 3-kind
        for bond in bonds(model)
            s1,s2 = source(bond), target(bond)
            bt = bondtype(bond)
            if spins[kind, s1] == spins[kind, s2] &&
                    rand(rng) < ps[bt, ifelse(spins[dual_kind, s1] == spins[dual_kind, s2], 1,2)]
                unify!(uf, s1,s2)
            end
        end
        nc = clusterize!(uf)
        clusterspin = rand(rng,[1,-1], nc)

        @inbounds for site in 1:nsites
            id = clusterid(uf, site)
            model.spins[kind, site] = clusterspin[id]
        end
    end
    return nothing
end

function Wolff_update!(model::AshkinTeller, T::Real, Jsigma, Jtau, K)
    rng = model.rng
    m2beta = -2.0/T
    nsites = numsites(model)
    ps = zeros(numbondtypes(model), 2)

    @inbounds for kind in 1:2
        Js = ifelse(kind==1, Jsigma, Jtau)
        ps[:,1] .= -expm1.(m2beta.*(Js.+K))
        ps[:,2] .= -expm1.(m2beta.*(Js.-K))
        dual_kind = 3-kind
        clustersize = 0
        st = Stack(Deque{Int}())
        center = rand(rng, 1:nsites)
        sp = model.spins[kind, center]
        model.spins[kind, center] *= -1
        push!(st, center)
        while !isempty(st)
            clustersize += 1
            s = pop!(st)
            for (n,b) in neighbors(model, s)
                bt = bondtype(model,b)
                if (model.spins[kind, n] == sp &&
                    rand(rng) < ps[bt, ifelse(model.spins[dual_kind,n] == model.spins[dual_kind, s], 1,2)])
                    model.spins[kind, n] *= -1
                    push!(st, n)
                end
            end
        end
    end

    return nothing
end

@gen_convert_parameter(AshkinTeller, ("T", 1, 1.0), ("Jsigma", numbondtypes, 1.0),
                       ("Jtau", numbondtypes, 1.0), ("K", numbondtypes, 0.0))


text:
include("QR.jl")
# using QR
using Lehmann
using StaticArrays, Printf

const Float = FQR.Float
const Double = FQR.Double
const DotF = FQR.Float
const Tiny = DotF(1e-5)

struct FreqGrid{D} <: FQR.Grid
    sector::Int                       # sector
    omega::SVector{D,Float}             # actual location of the grid point
    coord::SVector{D,Int}         # integer coordinate of the grid point on the fine meshes
end

Base.show(io::IO, grid::FreqGrid{2}) = print(io, "ω$(grid.sector) = ($(@sprintf("%12.4f", grid.omega[1])), $(@sprintf("%12.4f", grid.omega[2])))")

struct FreqFineMesh{D} <: FQR.FineMesh
    color::Int                            # D+1 sectors
    symmetry::Int                         # symmetrize colors and (omega1, omega2) <-> (omega2, omega1)
    candidates::Vector{FreqGrid{D}}       # vector of grid points
    selected::Vector{Bool}
    residual::Vector{Double}

    ## for frequency mesh only ###
    fineGrid::Vector{Float}         # fine grid for each dimension
    cache1::Vector{DotF}            # cache for exp(-x)
    cache2::Matrix{DotF}            # cache for exp(-x-y)

    function FreqFineMesh{D}(Λ, rtol; sym = 0) where {D}
        # initialize the residual on fineGrid with <g, g>
        _finegrid = Float.(fineGrid(Λ, rtol))
        separationTest(D, _finegrid)
        Nfine = length(_finegrid)

        _cache1 = zeros(DotF, Nfine)
        _cache2 = zeros(DotF, (Nfine, Nfine))
        for (xi, x) in enumerate(_finegrid)
            _cache1[xi] = exp(-DotF(x))
            for (yi, y) in enumerate(_finegrid)
                _cache2[xi, yi] = exp(-DotF(x) - DotF(y))
            end
        end

        color = D + 1
        # color = 1
        mesh = new{D}(color, sym, [], [], [], _finegrid, _cache1, _cache2)

        if D == 2
            for (xi, x) in enumerate(_finegrid)
                for (yi, y) in enumerate(_finegrid)
                    coord = (xi, yi)
                    for sector in 1:color
                        if irreducible(D, sector, coord, sym)  # if grid point is in the reducible zone, then skip residual initalization
                            g = FreqGrid{D}(sector, (x, y), coord)
                            push!(mesh.candidates, g)
                            push!(mesh.residual, FQR.dot(mesh, g, g))
                            push!(mesh.selected, false)
                        end
                    end
                end
            end
            # elseif D == 3
        else
            error("not implemented!")
        end
        println("fine mesh initialized.")
        return mesh
    end
end

"""
composite expoential grid
"""
function fineGrid(Λ, rtol)
    ############## use composite grid #############################################
    # degree = 8
    # ratio = Float(1.4)
    # N = Int(floor(log(Λ) / log(ratio) + 1))
    # panel = [Λ / ratio^(N - i) for i in 1:N]
    # grid = Vector{Float}(undef, 0)
    # for i in 1:length(panel)-1
    #     uniform = [panel[i] + (panel[i+1] - panel[i]) / degree * j for j in 0:degree-1]
    #     append!(grid, uniform)
    # end
    # append!(grid, Λ)
    # println(grid)
    # println("Composite expoential grid size: $(length(grid))")
    # return grid

    ############# DLR based fine grid ##########################################
    dlr = DLRGrid(Euv = Float64(Λ), beta = 1.0, rtol = Float64(rtol) / 100, isFermi = true, symmetry = :ph, rebuild = true)
    # println("fine basis number: $(dlr.size)\n", dlr.ω)
    degree = 4
    grid = Vector{Double}(undef, 0)
    panel = Double.(dlr.ω)
    for i in 1:length(panel)-1
        uniform = [panel[i] + (panel[i+1] - panel[i]) / degree * j for j in 0:degree-1]
        append!(grid, uniform)
    end

    println("fine grid size: $(length(grid)) within [$(grid[1]), $(grid[2])]")
    return grid
end

"""
Test the finegrids do not overlap
"""
function separationTest(D, finegrid)
    if D == 2
        epsilon = eps(DotF(1)) * 10
        for (i, f) in enumerate(finegrid)
            # either zero, or sufficiently large
            @assert abs(f) < epsilon || abs(f) > Tiny "$i: $f should either smaller than $epsilon or larger than $Tiny"
            for (j, g) in enumerate(finegrid)
                # two frequencies are either the same, or well separated
                @assert abs(f - g) < epsilon || abs(f - g) > Tiny "$i: $f and $j: $g should either closer than $epsilon or further than $Tiny"
                fg = f + g
                for (k, l) in enumerate(finegrid)
                    @assert abs(l - fg) < epsilon || abs(l - fg) > Tiny "$i: $f + $j: $g = $fg and $k: $l should either closer than $epsilon or further than $Tiny"
                end
            end
        end
    elseif D==1
        return
    else
        error("not implemented!")
    end
end

function coord2omega(mesh::FreqFineMesh{dim}, coord) where {dim}
    fineGrid = mesh.fineGrid
    if dim == 1
        return fineGrid[coord[1]]
    elseif dim == 2
        return (fineGrid[coord[1]], fineGrid[coord[2]])
    elseif dim == 3
        return (fineGrid[coord[1]], fineGrid[coord[2]], fineGrid[coord[3]])
    else
        error("not implemented!")
    end
end

function irreducible(D, sector, coord, symmetry)
    if symmetry == 0
        return true
    else
        if D == 2
            # return (coord[1] <= coord[2]) && (sector == 1)
            return (coord[1] <= coord[2])
        elseif D == 3
            # return (coord[1] <= coord[2] <= coord[3]) && (sector == 1)
            return (coord[1] <= coord[2] <= coord[3])
        else
            error("not implemented!")
        end
    end
end

function FQR.irreducible(grid::FreqGrid{D}) where {D}
    return irreducible(D, grid.sector, grid.coord, mesh.symmetry)
end

function FQR.mirror(mesh::FreqFineMesh{D}, idx) where {D}
    grid = mesh.candidates[idx]
    coord, s = grid.coord, grid.sector
    if mesh.symmetry == 0
        return []
    end
    if D == 2
        x, y = coord
        coords = unique([(x, y), (y, x),])
        # println(coords)
    elseif D == 3
        x, y, z = coord
        coords = unique([(x, y, z), (x, z, y), (y, x, z), (y, z, x), (z, x, y), (z, y, x)])
    else
        error("not implemented!")
    end
    newgrids = FreqGrid{D}[]
    # for s in 1:mesh.color
    for c in coords
        if s != grid.sector || c != Tuple(grid.coord)
            push!(newgrids, FreqGrid{D}(s, coord2omega(mesh, c), c))
        end
    end
    # end
    return newgrids
end

"""
F(x) = (1-exp(-y))/(x-y)
"""
# @inline function G2d(a::T, b::T, expa::T, expb::T) where {T}
#     if abs(a - b) > Tiny
#         return (expa - expb) / (b - a)
#     else
#         return (expa + expb) / 2
#     end
# end

"""
G(x, y) = (exp(-x)-exp(-y))/(x-y)
G(x, x) = -exp(-x)
"""
@inline function G2d(a::T, b::T, expa::T, expb::T) where {T}
    if abs(a - b) > Tiny
        return (expa - expb) / (b - a)
    else
        return (expa + expb) / 2
    end
end

"""
F(a, b, c) = (G(a, c)-G(a, c))/(a-b) where a != b, but a or b could be equal to c
"""
@inline function F2d(a::T, b::T, c::T, expa::T, expb::T, expc::T) where {T}
    @assert abs(a - b) > Tiny "$a - $c > $Tiny"
    return (G2d(a, c, expa, expc) - G2d(b, c, expb, expc)) / (b - a)
end

"""
F(any, any, 0)
"""
@inline function Fii2d(ω1::T, ω2::T, expω1::T, expω2::T) where {T}
    if ω1 < Tiny && ω2 < Tiny
        return T(1) / 2
    elseif ω1 < Tiny && ω2 > Tiny
        return (1 - ω2 - expω2) / ω2 / (ω1 - ω2)
    elseif ω1 > Tiny && ω2 < Tiny
        return (1 - ω1 - expω1) / ω1 / (ω2 - ω1)
    elseif abs(ω1 - ω2) < Tiny
        # @assert abs(ω1 - ω2) < eps(Float(1)) * 1000 "$ω1 - $ω2 = $(ω1-ω2)"
        ω = (ω1 + ω2) / 2
        expω = (expω1 + expω2) / 2
        return T((1 - expω * (1 + ω)) / ω / ω)
    else
        return T((ω1 - ω2 + expω1 * ω2 - expω2 * ω1) / (ω1 * ω2 * (ω1 - ω2)))
    end
end

"""
F(a,b,c)
"""
@inline function Fij2d(a::T, b::T, c::T, expa::T, expb::T, expc::T) where {T}
    if abs(a - b) > Tiny #a!=b
        return F2d(a, b, c, expa, expb, expc)
    else # a=b
        if abs(a - c) > Tiny # a=b != c
            return F2d(a, c, b, expa, expc, expb)
        else # a==b==c: exp(-a)/2
            return (expa + expb + expc) / 6
        end
    end
end

"""
basis dot for 2D
"""
function FQR.dot(mesh::FreqFineMesh{2}, g1::FreqGrid{2}, g2::FreqGrid{2})
    # println("dot: ", g1, ", ", g2)
    cache1 = mesh.cache1
    cache2 = mesh.cache2
    s1, s2 = g1.sector, g2.sector
    c1, c2 = g1.coord, g2.coord
    if s1 == s2  # F11, F22, F33
        ω1, ω2 = g1.omega[1] + g2.omega[1], g1.omega[2] + g2.omega[2]
        expω1 = cache2[c1[1], c2[1]]
        expω2 = cache2[c1[2], c2[2]]
        return Fii2d(ω1, ω2, expω1, expω2)
    elseif (s1 == 1 && s2 == 2) || (s1 == 2 && s2 == 3) || (s1 == 3 && s2 == 1) #F12, F23, F31
        a, b, c = g2.omega[2], g1.omega[1], g1.omega[2] + g2.omega[1]
        ea, eb, ec = cache1[c2[2]], cache1[c1[1]], cache2[c1[2], c2[1]]
        return Fij2d(a, b, c, ea, eb, ec)
    else  #F21, F32, F13
        a, b, c = g1.omega[2], g2.omega[1], g2.omega[2] + g1.omega[1]
        ea, eb, ec = cache1[c1[2]], cache1[c2[1]], cache2[c2[2], c1[1]]
        return Fij2d(a, b, c, ea, eb, ec)
    end
end

if abspath(PROGRAM_FILE) == @__FILE__

    D = 2

    lambda, rtol = 10, 1e-4
    mesh = FreqFineMesh{D}(lambda, rtol, sym = 0)

    # KK = zeros(3, 3)
    # n = (2, 2)
    # o = (mesh.fineGrid[n[1]], mesh.fineGrid[n[2]])
    # for i in 1:3
    #     g1 = FreqGrid{2}(i, o, n)
    #     for j in 1:3
    #         g2 = FreqGrid{2}(j, o, n)
    #         println(g1, ", ", g2)
    #         KK[i, j] = FQR.dot(mesh, g1, g2)
    #     end
    # end
    # display(KK)
    # println()

    basis = FQR.Basis{D,FreqGrid{D}}(lambda, rtol, mesh)
    FQR.qr!(basis, verbose = 1)

    lambda, rtol = 1000, 1e-8
    mesh = FreqFineMesh{D}(lambda, rtol, sym = 0)
    basis = FQR.Basis{D,FreqGrid{D}}(lambda, rtol, mesh)
    @time FQR.qr!(basis, verbose = 1)

    FQR.test(basis)

    mesh = basis.mesh
    grids = basis.grid
    open("basis.dat", "w") do io
        for (i, grid) in enumerate(grids)
            if grid.sector == 1
                println(io, grid.omega[1], "   ", grid.omega[2])
            end
        end
    end
    Nfine = length(mesh.fineGrid)
    open("finegrid.dat", "w") do io
        for i in 1:Nfine
            println(io, basis.mesh.fineGrid[i])
        end
    end
    open("residual.dat", "w") do io
        # println(mesh.symmetry)
        residual = zeros(Double, Nfine, Nfine)
        for i in 1:length(mesh.candidates)
            if mesh.candidates[i].sector == 1
                x, y = mesh.candidates[i].coord
                residual[x, y] = mesh.residual[i]
                # println(x, ", ", y, " -> ", length(mirror(mesh, i)))

                for grid in FQR.mirror(mesh, i)
                    if grid.sector == 1
                        xp, yp = grid.coord
                        residual[xp, yp] = residual[x, y]
                        # println(xp, ", ", yp)
                    end
                end
            end
        end

        for i in 1:Nfine
            for j in 1:Nfine
                println(io, residual[i, j])
            end
        end
    end
end


text:
#=
Lattice is implemented in a generic way, such that it should work with any
Nodes, Edges and Graphs as long as they fulfill a few requirements:

- They must inherit from AbstractNode, AbstractEdge or AbstractGraph
  respectively
- They must contain the fields StandardSite, StandardBond and StandardGraph contain
- They must have a constructor which only requires these fields.

StandardSite, StandardBond and StandardGraph are also generic in the sense that they
should work with any type inheriting from AbstractNode, AbstractEdge or
AbstractGraph.

E.g. the following works to add positions:

struct PosSite{D, T <: AbstractEdge} <: AbstractNode
    neighbors::Vector{Vector{T}}
    pos::Point{D}
    uvw::NTuple{D+1, Int64}
end
function PosSite(
        uc::UnitCell{D, T},
        neighbors::Vector{Vector{ET}},
        iuvw::NTuple
    ) where {D, T <: AbstractFloat, ET <: AbstractEdge}
    PosSite(neighbors, get_pos(B, iuvw[1], uvw[2:end]))
end
=#

abstract type AbstractNode end
abstract type AbstractEdge end
abstract type AbstractGraph end

struct StandardSite{N, T <: AbstractEdge} <: AbstractNode
    neighbors::Vector{Vector{T}}
    uvw::NTuple{N, Int64}
end

function StandardSite(
        uc::UnitCell,
        neighbors::Vector{Vector{T}},
        uvw::NTuple{N, Int64}
    ) where {N, T <: AbstractEdge}
    StandardSite(neighbors, uvw)
end

struct StandardBond{T <: AbstractNode} <: AbstractEdge
    node1::T
    node2::T
end

struct StandardGraph{
        Dimensions,
        T,
        NodeType <: AbstractNode,
        EdgeType <: AbstractEdge
    } <: AbstractGraph

    unitcell::UnitCell{Dimensions, T}
    nodes::Array{NodeType}
    edges::Vector{Vector{EdgeType}}
end

"""
    Lattice(uc[; kwargs...])

Generates a lattice of size L^d based on a UnitCell uc, which can be given
as a ::Function, a ::UnitCell.

Additional keyword arguments (kwargs) include
- `L::Int64 = 8`: The size of the lattice
- `Ls::NTuple{D, Int64} = (L, ..., L)`: The size of the lattice in D dimensions.
- `do_periodic::Bool = true`: Generate lattice with periodic bonds.
- `N_neighbors::Int64 = 1`: The number neighbor levels used.
- `nodetype::Type{N} = StandardSite`: Type of Node used to construct the lattice.
- `edgetype::Type{E} = StandardBond`: Type of Edge used to construct the lattice.
- `graphtype::Type{G} = StandardGraph`: Type of Graph used to construct the lattice.
"""
function Lattice(
        uc::UnitCell{D, T};
        L::Int64 = 8,
        Ls::NTuple{D, Int64} = ([L for _ in 1:D]...),
        do_periodic::Bool = true,
        N_neighbors::Int64 = 1,
        nodetype::Type{N} = StandardSite,
        edgetype::Type{E} = StandardBond,
        graphtype::Type{G} = StandardGraph
    ) where {D, T, N <: AbstractNode, E <: AbstractEdge, G <: AbstractGraph}

    dimsize = reduce(*, Ls)

    # Search for neighbors
    relative_offsets = get_neighbors(uc, N_neighbors)

    # Generate Nodes/Sites without edges
    # shape: (#number of Bravais lattices, L₁, ..., Lₙ) with n dimensions
    nodes = reshape([
        nodetype(
            uc,
            [edgetype[] for _ in relative_offsets[i]],
            (i, ind2sub(Ls, j)...)
        ) for i in eachindex(uc.basis), j in 1:dimsize
    ], length(uc.basis), Ls...)

    # Generate and connect edges
    edges = [edgetype[] for _ in 1:N_neighbors]

    for i in eachindex(nodes)
        from = ind2sub(nodes, i)

        for (lvl, lvl_offsets) in enumerate(relative_offsets[from[1]])
            for (j, offset) in enumerate(lvl_offsets)
                accept, to = move(from, offset, do_periodic, Ls)
                !accept && continue

                e = edgetype(nodes[from...], nodes[to...])
                if !(e in nodes[from...].neighbors[lvl])
                    push!(nodes[from...].neighbors[lvl], e)
                    push!(nodes[to...].neighbors[lvl], e)
                    push!(edges[lvl], e)
                end
            end
        end
    end

    graphtype(uc, nodes, edges)
end

# Extra methods
function Lattice(uc_func::Function, args...; kwargs...)
    Lattice(uc_func(), args...; kwargs...)
end

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

# Required for in
# Checks whether if e1 == e2, ignoring the direction of e1 and e2
function ==(e1::T, e2::T) where {T <: AbstractEdge}
    ((e1.node1 == e2.node1) && (e1.node2 == e2.node2)) ||
    ((e1.node1 == e2.node2) && (e1.node2 == e2.node1))
end

# Checks if edge is in edgelist
function in(edge::T, edgelist::Vector{T}) where {T <: AbstractEdge}
    for e in edgelist
        if edge == e
            return true
        end
    end
    false
end

# This calculates to = from + by with or without periodic boundaries.
# If the result is within the boundaries, it'll return true, to.
function move(
        from::NTuple{N, T},     # starting indices
        by::NTuple{N, T},       # offset indices
        do_periodic::Bool,      # use periodic boundaries
        Ls::NTuple{D, Int64}    # maximum index value
    ) where {D, N, T}

    to = [from[i] + by[i] for i in 1:N]

    if do_periodic
        for i in 2:length(to)
            (to[i] < 1) && (to[i] += Ls[i-1])
            (to[i] > Ls[i-1]) && (to[i] -= Ls[i-1])
        end
        return true, (to...)
    else
        return all(xL -> 1 <= xL[1] <= xL[2], zip(to[2:end], Ls)), (to...)
    end
end

# Pretty printing to avoid printing circular references
function show(io::IO, g::T) where {T <: AbstractGraph}
    D = dims(g.unitcell)
    print(io, "$D-dimensional lattice graph with ")
    print(io, length(g.nodes))
    print(io, " sites and ")
    print(io, mapreduce(length, +, g.edges))
    print(io, " bonds.")
end

# Pretty printing to avoid printing circular references
function show(io::IO, n::T) where {T <: AbstractNode}
    print(io, "Node $(n.uvw) with ")
    print(io, "$(map(length, n.neighbors)) edges.")
end
# Pretty printing to avoid printing circular references
function show(io::IO, e::T) where {T <: AbstractEdge}
    print(io, "Edge connecting Node $(e.node1.uvw) and $(e.node2.uvw)")
end


text:
cd("/home/astupidbear/.julia/v0.4/JuLearn/src/BLAS/examples/mlp/regression")
ENV["JuLearnBackend"]="BLAS"
using JuLearn#,Plots
@Eval lr=0.01 δ=0. momentum=0. batchsize=10 ▽max=0.01 λ=0.0
include("dataProvider.jl")
include("predict.jl")

inputBase,outputBase=dataProvider();
net=chain([:d,:h,:h,:l],[1 1;1 1000;1000 1;1 1],:regression);
@time err=backprop(net,Array(inputBase),Array(outputBase),nIter=1000)

# plot(err)
# plot(inputBase',Array(predict(net,Array(inputBase)))')
# plot!(inputBase',outputBase')


text:
using FastaIO;

DNA_ONE_HOT = Dict("A" => [1,0,0,0], "C" => [0,1,0,0], "G" => [0,0,1,0], "T" => [0,0,0,1]);

function parse_training(pos_seq_fp, total_seq_fp, balance)

    pos_mat = parse_seq(pos_seq_fp);

    neg_seqs = readfasta(total_seq_fp)[1:100];
    neg_mat = zeros(Int, 1, size(pos_mat, 2));

    for seq in neg_seqs
        seq_segs = create_17b_segments(seq[2]);
        for seg in seq_segs
            oh = dna_to_one_hot(seg);
            neg_mat = vcat(neg_mat, oh.');
        end
    end
    neg_mat = neg_mat[2:end, :];

    neg_mat = filter_mats(pos_mat, neg_mat)

    neg_mat = neg_mat[rand(1:size(neg_mat, 1), size(pos_mat, 1)*balance),:]

    labels = ones(size(pos_mat, 1))
    labels = vcat(labels, zeros(size(neg_mat, 1)))

    training_data = vcat(pos_mat, neg_mat)

    training_data, labels

end

function parse_testing(test_fp)

    test_mat = parse_seq(test_fp)

    test_mat

end

function create_17b_segments(seq)

    seqs = []
    for i in 1:17:length(seq)
        if i + 16 < length(seq)
            push!(seqs, seq[i:i+16])
        end
    end

    seqs

end

# Hash function to compare sequences
function seq_hash(seq)

    tmp = seq .* 2.^collect(0:(length(seq)-1))
    sum(tmp)

end

function filter_mats(pos_mat, neg_mat)

    pos_hash = [seq_hash(pos_mat[i,:]) for i in 1:size(pos_mat, 1)]
    neg_hash = [seq_hash(neg_mat[i,:]) for i in 1:size(neg_mat, 1)]

    to_remove = []
    for i in 1:size(neg_hash, 1)
        if contains(==, neg_hash[i], pos_hash)
            push!(to_remove, i)
        end

        if length(neg_hash[neg_hash .== neg_hash[i]]) > 1
            push!(to_remove, i)
        end

    end

    neg_mat = neg_mat[setdiff(1:end, to_remove), :];

    neg_mat

end

function parse_seq(fp)

    f = open(fp, "r");
    l = collect(readlines(f));

    mat = zeros(Int, length(l), length(l[1]) * 4)

    for i in 1:length(l)

        seq = l[i]
        mat[i,:] = dna_to_one_hot(seq);

    end

    close(f);

    mat

end

convert_dna(x) = DNA_ONE_HOT[x];

function dna_to_one_hot(seq)
    one_hot = zeros(Int, 4, length(seq))
    for b in 1:length(seq)
        one_hot[:,b] = convert_dna(string(seq[b]))
    end

    collect(Iterators.flatten(one_hot));
end


text:
using Genie
using Genie.Renderer.Html
using Genie.Requests

const content = "abcd"
const content2 = "efgh"

@testset "Basic rendering" begin
  r = Requests.HTTP.Response()

  @testset "Empty string" begin
    r = html("")

    @test String(r.body) == ""
  end;

  @testset "Empty string force parse" begin
    r = html("", forceparse = true)

    @test String(r.body) == ""
  end;

  @testset "String no spaces" begin
    r = html(content)

    @test String(r.body) == "$content"
  end;

  @testset "String no spaces force parse" begin
    r = html(content, forceparse = true)

    @test String(r.body) == "abcd"
  end;

  @testset "String with 2 spaces" begin
    r = html("  $content  ")

    @test String(r.body) == "  abcd  "
  end;

  @testset "String with 2 spaces force parse" begin
    r = html("  $content  ", forceparse = true)

    @test String(r.body) == "abcd  "
  end;

  @testset "String with  " begin
    r = html("      ")

    @test String(r.body) == "      "
  end;

  @testset "String with   force parse" begin
    r = html("      ", forceparse = true)

    @test String(r.body) == "    "
  end;

  @testset "String with 2   and 2 spaces" begin
    r = html("  $content  ")

    @test String(r.body) == "  abcd  "
  end;

  @testset "String with 2   and 2 spaces force parse" begin
    r = html("  $content  ", forceparse = true)

    @test String(r.body) == "  abcd  "
  end;

  @testset "String with newline" begin
    r = html("$content  \n    $content2")

    @test String(r.body) == "abcd  \n    efgh"
  end;

  @testset "String with newline force parse" begin
    r = html("$content  \n    $content2", forceparse = true)

    @test String(r.body) == "abcd  \n  efgh"
  end;

  @testset "String with quotes" begin
    r = html("He said \"wow!\"")

    @test String(r.body) == "He said \"wow!\""
  end;

  @testset "String with quotes force parse" begin
    r = html("He said \"wow!\"", forceparse = true)

    @test String(r.body) == "He said \"wow!\"\0"
  end;

  @testset "String with quotes" begin
    r = html(""" "" """)

    @test String(r.body) == " \"\" "
  end;

  @testset "String with quotes force parse" begin
    r = html(""" "" """, forceparse = true)

    @test String(r.body) == "\"\" "
  end;

  @testset "String with quotes" begin
    r = html("\"\"")

    @test String(r.body) == "\"\""
  end;

  @testset "String with quotes force parse" begin
    r = html("\"\"", forceparse = true)

    @test String(r.body) == "\"\"\0"
  end;

  @testset "String with interpolated vars" begin
    r = html("$(reverse(content))")

    @test String(r.body) == "dcba"
  end;

  @testset "String with interpolated vars force parse" begin
    r = html("$(reverse(content))", forceparse = true)

    @test String(r.body) == "dcba"
  end;

  @test r.status == 200
  @test r.headers[1]["Content-Type"] == "text/html; charset=utf-8"
end;


text:
## Partitions of Sample Space Ω
## ============================

"""Approximate partition of a set: both an `under` approximation and `rest`.
`under ∪ rest` is an overapproximation of set."""
type ApproxPartition{D<:Domain}
  under::Vector{D}
  rest::Vector{D}
end

"Lower and upper bounds for measure of partition"
function measure(p::ApproxPartition)
  lb = isempty(p.under) ? 0 : sum([measure(event) for event in p.under])
  rest = isempty(p.rest) ? 0 : sum([measure(event) for event in p.rest])
  ub = lb + rest
  Interval(lb,ub)
end

"`under ∪ rest` - over approximation"
collect{D}(p::ApproxPartition{D}) = vcat(p.under, p.rest)

## Sampling
## ========

"A partition which is efficient for drawing many samples"
type SampleablePartition{D<:Domain}
  over::Vector{D}
  last_under::Int
  cat::Categorical

  function SampleablePartition{D}(under::Vector, rest::Vector) where {D<:Domain}
    # Compute all measures and categorical weights only once
    over = vcat(under, rest)
    vols = Float64[measure(box) for box in over]
    pnormalize!(vols)
    cat = Categorical(vols, Distributions.NoArgCheck())
    new(over,length(over),cat)
  end
end

SampleablePartition{D}(p::ApproxPartition{D}) =
  SampleablePartition{D}(p.under, p.rest)

"Point sample from preimage - may be invalid point due to approximations"
abstract_sample(p::SampleablePartition) = p.over[rand(p.cat)]

"`n` Point samples from preimage - may be invalid point due to approximations"
function point_sample(p::SampleablePartition, n::Integer)
  rand_indices = rand(p.cat, n)
  events = [p.over[i] for i in rand_indices]
  [rand(event) for event in events]
end

point_sample(p::SampleablePartition) = point_sample(p,1)[1]

"Do refined rejection sampling from preimage"
function point_sample_exact(p::SampleablePartition, Y::RandVar{Bool};
                            maxtries = 1E7)
  for i = 1:maxtries
    sample = rand(p)
    if Y(sample) return sample end
  end
  error("Could not get sample in $maxtries tries")
end

"Generic preimage partition interface - partition of Y^-1{true}."
function pre_partition{T}(
    Y::SymbolicRandVar{Bool},
    ::Type{T};
    RandVarType::Type = default_randvar(),
    precision = default_precision(),
    args...)

  init_box = unit_box(LazyBox{Float64}, dims(Y))
  Y_conv = convert(RandVarType{Bool}, Y)
  set_precision!(Y_conv, precision)
  pre_partition(Y_conv, init_box, T; precision=precision, args...)
end


text:
# https://microsoft.github.io/language-server-protocol/specifications/specification-3-17/#textDocument_semanticTokens
const SemanticTokenKind = String
const SemanticTokenKinds = (
    # Namespace = "namespace",
    # Type = "type",
    # Class = "class",
    # Enum = "enum",
    # Interface = "interface",
	Struct = "struct",
	TypeParameter = "typeParameter",
	Parameter = "parameter",
	Variable = "variable",
	Property = "property",
	# EnumMember = "enumMember",
	# Event = "event",
	Function = "function",
	# Method = "method",
	Macro = "macro",
	Keyword = "keyword",
	# Modifier = "modifier",
	Comment = "comment",
	String = "string",
	Number = "number",
	Regexp = "regexp",
	Operator = "operator"
)

const SemanticTokenModifiersKind = String
const SemanticTokenModifiersKinds = (
    Declaration = "declaration",
	Definition = "definition",
	# Readonly = "readonly",
	# Static = "static",
	# Deprecated = "deprecated",
	# Abstract = "abstract",
	# Async = "async",
	Modification = "modification",
	Documentation = "documentation",
	DefaultLibrary = "defaultLibrary"

)

struct SemanticTokensLegend <: Outbound
	# /**
	#  * The token types a server uses.
	#  */
	tokenTypes::Vector{String}

	# /**
	#  * The token modifiers a server uses.
	#  */
	tokenModifiers::Vector{String}
end

const JuliaSemanticTokensLegend = SemanticTokensLegend(
		collect(values(SemanticTokenKinds)),
		collect(values(SemanticTokenModifiersKinds))
)

function semantic_token_encoding(token :: String) :: UInt32
	for (i, type) in enumerate(JuliaSemanticTokensLegend.tokenTypes)
		if token == type
			return i - 1 # -1 to shift to 0-based indexing
		end
	end
end
# function SemanticTokensLegend() :: SemanticTokensLegend
# 	SemanticTokensLegend(
# 	)
# end
# const Tok

@dict_readable struct SemanticTokensFullDelta <: Outbound
    delta::Union{Bool,Missing}
end

@dict_readable struct SemanticTokensClientCapabilitiesRequests <: Outbound
    range::Union{Bool,Missing}
    full::Union{Bool,Missing,SemanticTokensFullDelta}

end
@dict_readable struct SemanticTokensClientCapabilities <: Outbound
    dynamicRegistration::Union{Bool,Missing}
    tokenTypes::Vector{String}
    tokenModifiers::Vector{String}
    formats::Vector{String}
    overlappingTokenSupport::Union{Bool,Missing}
    multilineTokenSupport::Union{Bool,Missing}
end

struct SemanticTokensOptions <: Outbound
	legend::SemanticTokensLegend
	range::Union{Bool,Missing}
	full::Union{Bool,SemanticTokensFullDelta,Missing}
end

struct SemanticTokensRegistrationOptions <: Outbound
    documentSelector::Union{DocumentSelector,Nothing}
    # workDoneProgress::Union{Bool,Missing}
end

@dict_readable struct SemanticTokensParams <: Outbound
    textDocument::TextDocumentIdentifier
    # position::Position
    workDoneToken::Union{Int,String,Missing} # ProgressToken
    partialResultToken::Union{Int,String,Missing} # ProgressToken
end

struct SemanticTokens <: Outbound
    resultId::Union{String,Missing}
    data::Vector{UInt32}
end

SemanticTokens(data::Vector{UInt32}) = SemanticTokens(missing, data)

struct SemanticTokensPartialResult <: Outbound
    data::Vector{UInt32}
end

struct SemanticTokensDeltaParams <: Outbound
	workDoneToken::Union{Int,String,Missing}
    partialResultToken::Union{Int,String,Missing} # ProgressToken
    textDocument::TextDocumentIdentifier
	previousResultId::String
end
struct SemanticTokensEdit <: Outbound
	start::UInt32
	deleteCount::Int
	data::Union{Vector{Int},Missing}
end
struct SemanticTokensDelta <: Outbound
	resultId::Union{String,Missing}
	edits::Vector{SemanticTokensEdit}
end

struct SemanticTokensDeltaPartialResult <: Outbound
	edits::Vector{SemanticTokensEdit}
end

struct SemanticTokensRangeParams <: Outbound
	workDoneToken::Union{Int,String,Missing}
    partialResultToken::Union{Int,String,Missing} # ProgressToken
    textDocument::TextDocumentIdentifier
	range::Range
end

struct SemanticTokensWorkspaceClientCapabilities <: Outbound
	refreshSupport::Union{Bool,Missing}
end


text:
using AlmostBandedMatrices, FillArrays, Compat.Test, Compat.LinearAlgebra
using Compat

# Constructions:
@testset "construction" begin
    let n = 10, v = rand(n, 2), w = rand(n, 2), y = rand(n, 3), z = rand(n, 3), U = LowRankMatrix(v, w), L = LowRankMatrix(y, z)
        @test SemiSeparableMatrix(U, L, 1, 2) isa SemiSeparableMatrix{Float64}
        @test Matrix(SemiSeparableMatrix(U, L, 1, 2)) == triu(U, 2)+tril(L, -3)
        @test_throws AssertionError SemiSeparableMatrix(U, L, 10, 2)
        @test_throws AssertionError SemiSeparableMatrix(U, L, 1, 10)
    end

    let v = rand(n, 2), w = rand(n, 2), y = rand(n+1, 3), z = rand(n+1, 3), U = LowRankMatrix(v, w), L = LowRankMatrix(y, z)
        @test_throws AssertionError SemiSeparableMatrix(U, L, 1, 2)
    end
end


text:
using Documenter
using ClinicalTrialSampleSize

makedocs()

deploydocs(
    repo = "github.com/ContaTP/ClinicalTrialSampleSize.jl.git",
    julia  = "0.5",
    deps = Deps.pip("mkdocs", "python-markdown-math", "mkdocs-cinder"),
)


text:
InputFolder = './Images/PointClusters/';
OutputFolder = './Results/Images/Voronoi/';
Rescale = 1;

@iA = '*.tif';

@fxs_gVoronoi3D [iA] > [G];
params.ovs = 1;
/endf

/show iA >;
/show G >;
/keep G > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
using JSON
using Logging
using ArgParse

include("src/config.jl")
include("src/network.jl")
include("src/gym.jl")

argtable = ArgParseSettings()
@add_arg_table(
    argtable,
    "--seed", arg_type = Int, default = 0,
    "--log", arg_type = String, default = "gym.log",
    "--id", arg_type = String, default = "CartPole-v0",
    "--plot", arg_type = String, default = "",
    "--n_trials", arg_type = Int64, default = 1,
    "--connectivity", arg_type = Float64, default = 0.3,
    "--tinput", arg_type = Int64, default = 4,
    "--toutput", arg_type = Int64, default = 4,
    "--trest", arg_type = Int64, default = 7,
    "--ma_rate", arg_type = Float64, default = 0.8,
    "--fr", arg_type = Float64, default = 0.5,
    "--n", arg_type = Int64, default = 12,
    "--nratio", arg_type = Float64, default = 2.8,
    "--onpolicy", arg_type = Bool, default = true,
    "--dmethod", arg_type = Int64, default = 1
)

args, cfg = get_config(argtable)

Logging.configure(filename=args["log"], level=INFO)
Logging.info(@sprintf("Arguments: %s", JSON.json(args)))
Logging.info(@sprintf("Config: %s", JSON.json(cfg)))

env = gym.make(args["id"])
nin = env[:observation_space][:shape][1] * args["n"]
naction = env[:action_space][:n]
nout = naction * args["n"]
nn = Int64(round((nin + nout) * args["nratio"]))

srand(args["seed"])
n = Network(nn, nin, nout, args["connectivity"], cfg)

pcfg = Dict()
pcfg["n_actions"] = naction
for k in ["tinput", "toutput", "trest"]
    pcfg[k] = args[k] * 10
end
for k in ["ma_rate", "fr", "n_trials", "seed", "id", "onpolicy", "dmethod"]
    pcfg[k] = args[k]
end

fit = repeat_trials(n, env, pcfg)
Logging.info(@sprintf("E%0.6f", -fit))

if length(args["plot"]) > 0
    include("src/plot.jl")
    plot_spikes(n, args["plot"])
end


text:
using Pkg
Pkg.activate(".")

#using Revise
using Dash, DashHtmlComponents, DashCoreComponents, DashBootstrapComponents
using Statistics, DataFrames, CSVFiles
using Dates
import TOML

include("utils.jl")

# App begins
df = load_dataframe()

app = dash(external_stylesheets = ["assets/style.css", dbc_themes.COSMO])

app.title = "Julia Packages Requests"

title = dbc_row(
    dbc_col(
        dbc_alert(
            color = "primary",
            html_center(
                html_div(
                    style = Dict("vertical-alight" => "middle"),
                    html_h1("Julia Packages requests"),
                ),
            ),
            className = "bmi-jumbotron",
        ),
    ),
    className = "mt-3",
)

ttip =
    dbc_tooltip("Enter a Julia package name and press the search button", target = "hola")

graph = dcc_graph(id = "graph", figure = plot_default())

graph_content = html_div([
    dbc_row([dbc_col(graph, width = 12)]),
    html_br(),
    dbc_row([
        dbc_col(html_div(id = "info-content1", ["Info1"]), width = 6),
        dbc_col(html_div(id = "info-content2", ["Info2"]), width = 6),
    ]),
])

tabs = dbc_row([
    dbc_col(
        html_div([
            dbc_tabs(
                [
                    dbc_tab(label = "Users", tab_id = "user"),
                    dbc_tab(label = "CI runs", tab_id = "ci"),
                ],
                id = "tabs",
                active_tab = "user",
            ),
            html_div(id = "graph-content", graph_content),
        ]),
        width = 12,
    ),
],);

date_picker = dcc_datepickerrange(
    id = "dates",
    start_date = df."date"[begin],
    end_date = df."date"[end],
)

date_input = dbc_inputgroup([
    dbc_inputgrouptext(id = "date-range", "Date range"),
    date_picker,
    dbc_tooltip("Select a date range for the request's data", target = "date-range"),
])

search_input = dbc_inputgroup([
    dbc_inputgrouptext("Package", id = "package-label"),
    dbc_input(id = "search-input", placeholder = "name"),
    dbc_tooltip(
        "Enter a Julia package name and press the search button",
        target = "package-label",
    ),
    dbc_button("Search", color = "primary", id = "search-button"),
])

drop = dcc_dropdown(
    id = "drop",
    options = [(label = "PlotlyJS", value = "PlotlyJS")],
    style = Dict("width" => "300px"),
    multi = true,
    value = nothing,
)

drop2 = dcc_dropdown(
    id = "drop2",
    options = [(label = "Regular", value = "REG"), (label = "Cumulative", value = "CMT")],
    style = Dict("width" => "300px"),
    value = "REG",
)

mode2 = dbc_inputgroup([
    dbc_inputgrouptext("Mode", id = "mode2"),
    drop2,
    dbc_tooltip("Select plot type", target = "mode2"),
])

packages_input = dbc_row([
    dbc_col(search_input, width = 4),
    dbc_col(dbc_inputgroup([dbc_inputgrouptext("Select", id = "multi"), drop]), width = 8),
    dbc_tooltip(
        "Select one or more of the searched package(s) to render",
        target = "multi",
    ),
])

drop_input = dbc_row(dbc_col(date_input, width = 6))
opts = []

callback!(
    app,
    Output("drop", "options"),
    Output("drop", "value"),
    Output("search-input", "value"),
    Input("search-button", "n_clicks"),
    State("search-input", "value"),
) do clicks, input_value
    value = input_value
    if value !== nothing
        value = split(input_value, ".")[1] # To strip .jl part
    end
    package_df = package_dataframe(value)
    if package_df === nothing
        return opts, nothing, ""
    end
    append!(opts, [(label = value, value = value)])
    opts_set = Set(opts)
    opts_new = [o for o in opts_set]
    opts_new, value, ""
end

checks = dbc_checklist(
    id = "checks",
    options = [
        Dict("label" => "Regression Line", "value" => "RL"),
        Dict("label" => "Seven days moving averages", "value" => "MA"),
    ],
    value = [],
    labelStyle = Dict("display" => "inline-block"),
)

function package_global_stats(name, stats)
    total_requests = stats["total"]
    interval = stats["interval"]
    precision = 2
    daily_average = round(total_requests / interval, digits = precision)
    α = stats["α"]
    β = stats["β"]
    dbc_card(
        [
            dbc_cardbody([
                html_h4("$(name) global statistics", className = "card-title"),
                html_div(
                    [
                        html_p("Total requests: $total_requests"),
                        html_p("Time interval: $interval days"),
                        html_p("Daily average: $daily_average"),
                        html_p("Regression slope: $β"),
                    ],
                    className = "card-text",
                ),
            ]),
        ],
        style = Dict("width" => "12rem"),
    )
end

function packages_info_table(packages, stats)
    table_header = [
        html_thead(
            html_tr([
                html_th("Package"),
                html_th("Daily average"),
                html_th("Regression slope"),
            ]),
        ),
    ]
    rows = []
    for (i, name) in enumerate(packages)
        s = stats[i]
        total = s["total"]
        interval = s["interval"]
        precision = 2
        daily_average = round(total / interval, digits = precision)
        β = s["β"]
        row = html_tr([html_td(name), html_td("$daily_average"), html_td("$β")])
        append!(rows, [row])
    end

    table_body = [html_tbody(rows)]
    dbc_table([table_header; table_body], bordered = true)
end

callback!(
    app,
    Output("graph", "figure"),
    Output("info-content1", "children"),
    Output("info-content2", "children"),
    Input("tabs", "active_tab"),
    Input("drop", "value"),
    Input("dates", "start_date"),
    Input("dates", "end_date"),
    Input("checks", "value"),
    Input("drop2", "value"),
) do at, options, sd, ed, check_vals, drop2_val
    println(options)
    if (options === nothing) || (length(options) == 0)
        return plot_default(), [""], [""]
    end
    if typeof(options) === String
        options = [options]
    end
    figure, c_figure, stats_info =
        plot_graphs(options, check_vals, at, start_date = sd, end_date = ed)

    fout = figure
    if drop2_val == "CMT"
        fout = c_figure
    end
    info1 = [""]
    info2 = [""]
    if options !== nothing
        if length(options) == 1
            info1 = [package_global_stats(options[1], stats_info[1])]
        else
            vals = [s["total"] for s in stats_info]
            figure = plot(pie(values = vals, labels = options, marker_colors = colors))
            info1 = [dcc_graph(figure = figure)]
            info2 = [packages_info_table(options, stats_info)]
        end
    end
    return fout, info1, info2
end

app.layout = dbc_container(
    [
        title
        html_br()
        packages_input
        html_br()
        dbc_row([dbc_col(checks, width = 4), dbc_col(mode2, width = 8)])
        html_br()
        drop_input
        html_br()
        html_br()
        dbc_spinner(tabs)
        html_br()
        html_br()
    ],
)

run_server(app, "0.0.0.0", debug = false)


text:
# # Using Atomistic-Compatible Packages

# There are three main steps to using an implementation of the Atomistic API.
# This simple toy example will walk you through each step using [Molly.jl](https://github.com/JuliaMolSim/Molly.jl) for the dynamics.
# Some of the values from this example were taken from [this guide](https://ase.tufts.edu/chemistry/lin/images/FortranMD_TeachersGuide.pdf).

# ## Step 0: Load Dependencies

using Atomistic
using InteratomicPotentials
using Molly

# ## Step 1A: Configuring the System
# First, we must create an [AtomsBase.jl](https://github.com/JuliaMolSim/AtomsBase.jl)-style system.
# Here we will create a small Argon cluster with periodic boundary conditions.
# We are using a helper function to initialize the system automatically with a particular temperature.

N = 864
box_size = 3.47786u"nm"
reference_temp = 94.4u"K"

initial_system = generate_atoms_in_cubic_cell(N, :Ar, box_size, reference_temp)

# ## Step 1B: Configuring the Simulator
# Second, we initialize our simulator. In this case we are using an AndersenThermostat which is provided by Molly.

Δt = 1e-2u"ps"
steps = 2000
thermostat = Molly.AndersenThermostat(reference_temp, Δt * 10)
simulator = MollySimulator(Δt, steps, coupling=thermostat)

# ## Step 1C: Configuring the Potential
# Lastly, we specify the interatomic potential that we will use for the simulation, Lennard-Jones in this case.

potential = InteratomicPotentials.LennardJones(1.657e-21u"J", 0.34u"nm", 0.765u"nm", [:Ar])

# ## Step 2: Running the Simulation

result = simulate(initial_system, simulator, potential)

# ## Step 3: Analyzing the Results
# We can now analyze the simulation results.
# The Atomistic API exposes a variety of quantities from each timestep of the simulation (more details [here](@ref MolecularDynamicsResult_Specification)).
# In this example we will look at the temperature, energy, and radial distribution function.

# ### Temperature
plot_temperature(result, simulator.steps ÷ 200)

# ### Energy
plot_energy(result, simulator.steps ÷ 200)

# ### Radial Distribution Function (RDF)
plot_rdf(result, potential.σ, Int(0.95 * steps))


text:
## these are differently named than SymPy or missing or ...

Base.abs(x::SymbolicObject) = sympy.Abs(x)
Base.abs2(x::SymbolicObject) = abs(x)^2
Base.max(x::Sym, a) = sympy.Max(x, a)
Base.min(x::Sym, a) = sympy.Min(x, a)

Base.cbrt(x::Sym) = x^(1//3)
Base.ceil(x::Sym) = ceiling(x)

## Trig
Base.asech(z::Sym) = log(sqrt(1/z-1)*sqrt(1/z+1) + 1/z)
Base.acsch(z::Sym) = log(sqrt(1+1/z^2) + 1/z) ## http://mathworld.wolfram.com/InverseHyperbolicCosecant.html
Base.atan(y::Sym, x) = sympy.atan2(y,x)

Base.sinc(x::Sym) = sympy.sinc(PI*x)
cosc(x::Sym) = diff(sinc(x))

Base.sincos(x::Sym) = (sin(x), cos(x))
Base.sinpi(x::Sym) = sympy.sin(x*PI)
Base.cospi(x::Sym) = sympy.cos(x*PI)
degree_variants = (:sind, :cosd, :tand, :cotd, :secd, :cscd,
                   :asind, :acosd, :atand, :acotd, :asecd, :acscd)

for  methvar in degree_variants
    meth = Symbol(String(methvar)[1:end-1])
    @eval begin
        (Base.$methvar)(ex::SymbolicObject) = ($meth)(PI*ex)
    end
end

Base.rad2deg(x::Sym) = (x * 180) / PI
Base.deg2rad(x::Sym) = (x * PI) / 180

## exponential
Base.log(x::Sym) = sympy.log(x)
Base.log(b::Number, x::Sym) = sympy.log(x, b)
Base.log2(x::SymbolicObject) = log(2,x)
Base.log10(x::SymbolicObject) = log(10,x)

## calculus.
## use a pair for limit x=>0
limit(x::SymbolicObject, xc::Pair, args...;kwargs...) = limit(x, xc[1], xc[2], args...;kwargs...)
## allow a function
limit(f::Function, x::Sym, c;kwargs...) = limit(f(x), x, c;kwargs...)
limit(f::Function, c;kwargs...) = limit(f, symbols("x"), c;kwargs...)

function Base.diff(f::Function, n::Int=1)
    x = symbols("x")
    diff(f(x), x, n)
end

## integrate(ex,a,b)
function integrate(ex::SymbolicObject, a::Number, b::Number)
    fs = free_symbols(ex)
    if length(fs) !== 1
        @warn "Need exactly on free symbol. Use `integrate(ex, (x, a, b))` instead"
        return
    end
    integrate(ex, (fs[1], a, b))
end
function integrate(f::Function, a::Number, b::Number)
    x = symbols("x")
    integrate(f(x), (x, a, b))
end
function integrate(f::Function)
    x = symbols("x")
    integrate(f(x))
end

## Add interfaces for solve, nonlinsolve when vector of equations passed in
solve(V::Vector{T}, args...; kwargs...) where {T <: SymbolicObject} =
    sympy.solve(V, args...; kwargs...)

"""
    nonlinsolve

Note: if passing variables in use a tuple (e.g., `(x,y)`) and *not* a vector (e.g., `[x,y]`).
"""
nonlinsolve(V::AbstractArray{T,N}, args...; kwargs...) where {T <: SymbolicObject, N} =
    sympy.nonlinsolve(V, args...; kwargs...)

linsolve(V::AbstractArray{T,N}, args...; kwargs...) where {T <: SymbolicObject, N} =
    sympy.linsolve(V, args...; kwargs...)
linsolve(Ts::Tuple, args...; kwargs...) where {T <: SymbolicObject} =
    sympy.linsolve(Ts, args...; kwargs...)

## dsolve allowing initial condiation to be specified

"""
   dsolve(eqn, var, args..,; ics=nothing, kwargs...)

Call `sympy.dsolve` with possible difference for initial condition specification.

For problems with an initial condition, the `ics` argument may be specified. This is *different* from the `ics` argument of `sympy.dsolve`. (Call directly if that is preferred.)

Here `ics` allows the specification of a term like `f(x0) = y0` as a tuple `(f, x0, y0)`. Similarly, a term like `f''(x0)=y0` is specified through `(f'', x0, y0)`. If more than one initial condition is needed, a tuple of tuples is used, as in `((f,x0,y0), (f',x0,z0))`.

Example:

```
x = Sym("x")
y = SymFunction("y")
eqn = y''(x) - y(x) - exp(x)
dsolve(eqn, y(x), ics=((y,0,1), (y, 1, 1//2)))
```
"""
function dsolve(eqn::Sym, args...; ics=nothing, kwargs...)
    if isa(ics, Nothing)
        sympy.dsolve(eqn, args...; kwargs...)
    else
        if isempty(args)
            var = first(free_symbols(eqn))
        else
            var = first(args)
        end
        # var might be f(x) or x, we want `x`
        if Introspection.classname(var) != "Symbol"
            var = first(var.args)
        end
        ## if we have one initial condition, can be passed in a (u,x0,y0) *or* ((u,x0,y0),)
        ## if more than one a tuple of tuples
        if eltype(ics) <: Tuple
            _dsolve(eqn, var, ics; kwargs...)
        else
            _dsolve(eqn, var, (ics,); kwargs...)
        end
    end
end

function _dsolve(eqn::Sym, var::Sym, ics; kwargs...)
    if length(ics) == 0
        throw(ArgumentError("""Some initial value specification is needed.
Specifying the function, as in `dsolve(ex, f(x))`, is deprecated.
Use `sympy.dsolve(ex, f(x); kwargs...)` directly for that underlying interface.
"""))
    end

    out = sympy.dsolve(eqn; kwargs...)
    ord = sympy.ode_order(eqn, var)

    ## `out` may be an array of solutions. If so we do each one.
    ## we want to use an array for output only if needed
    if !isa(out, Array)
        return _solve_ivp(out, var, ics,ord)
    else
        output = Sym[]
        for o in out
            a = _solve_ivp(o, var, ics,ord)
            a != nothing && push!(output, a)
        end
        return length(output) == 1 ? output[1] : output
    end
end

rhs(x::SymbolicObject) = pycall_hasproperty(x, :rhs) ? x.rhs : x
lhs(x::SymbolicObject) = pycall_hasproperty(x, :lhs) ? x.lhs : x

export dsolve, rhs, lhs

## Helper.
## out is an equation in var with constants. Args are intial conditions
## Return `nothing` if initial condition is not satisfied (found by `solve`)
function _solve_ivp(out, var, args, o)

    eqns = Sym[(diff(out.rhs(), var, f.n))(var=>x0) - y0 for (f, x0, y0) in args]
    sols = solve(eqns, Sym["C$i" for i in 1:o], dict=true)
    if length(sols) == 0
       return nothing
    end

    ## massage output
    ## Might have more than one solution, though unlikely. But if we substitute a variable
    ## for y0 we will get an array back from solve which may have length 1.
    if isa(sols, Array)
        if length(sols) == 1
            sols = sols[1]
        else
            return [out([Pair(k,v) for (k,v) in sol]...) for sol in sols]
        end
    end

    out([Pair(k,v) for (k,v) in sols]...)
end


text:
using Conception, Samsara

include("test/test_conception_for_cirularsys.jl")

max_stars_repo_name: leikanger/Samsara

max_issues_repo_name: leikanger/Samsara

max_forks_repo_name: leikanger/Samsara


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#    Tests of default printing.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Default" begin
    expected = """
\\begin{table}
  \\begin{tabular}{rrrr}
    \\hline\\hline
    \\textbf{Col. 1} & \\textbf{Col. 2} & \\textbf{Col. 3} & \\textbf{Col. 4} \\\\\\hline
    1 & false & 1.0 & 1 \\\\
    2 & true & 2.0 & 2 \\\\
    3 & false & 3.0 & 3 \\\\
    4 & true & 4.0 & 4 \\\\
    5 & false & 5.0 & 5 \\\\
    6 & true & 6.0 & 6 \\\\\\hline\\hline
  \\end{tabular}
\\end{table}
"""

    result = pretty_table(String, data, backend = :latex)
    @test result == expected
end

@testset "Pre-defined formats" begin

    # Default
    # ==========================================================================

    expected = """
\\begin{table}
  \\begin{tabular}{|r|r|r|r|}
    \\hline\\hline
    \\textbf{Col. 1} & \\textbf{Col. 2} & \\textbf{Col. 3} & \\textbf{Col. 4} \\\\\\hline
    1 & false & 1.0 & 1 \\\\\\hline
    2 & true & 2.0 & 2 \\\\\\hline
    3 & false & 3.0 & 3 \\\\\\hline
    4 & true & 4.0 & 4 \\\\\\hline
    5 & false & 5.0 & 5 \\\\\\hline
    6 & true & 6.0 & 6 \\\\\\hline\\hline
  \\end{tabular}
\\end{table}
"""

    result = pretty_table(String, data,
                          hlines = :all,
                          vlines = :all,
                          tf = tf_latex_default)

    @test result == expected

    # Simple
    # ==========================================================================

    expected = """
\\begin{table}
  \\begin{tabular}{|r|r|r|r|}
    \\hline
    \\textbf{Col. 1} & \\textbf{Col. 2} & \\textbf{Col. 3} & \\textbf{Col. 4} \\\\\\hline
    1 & false & 1.0 & 1 \\\\\\hline
    2 & true & 2.0 & 2 \\\\\\hline
    3 & false & 3.0 & 3 \\\\\\hline
    4 & true & 4.0 & 4 \\\\\\hline
    5 & false & 5.0 & 5 \\\\\\hline
    6 & true & 6.0 & 6 \\\\\\hline
  \\end{tabular}
\\end{table}
"""

    result = pretty_table(String, data,
                          hlines = :all,
                          vlines = :all,
                          tf = tf_latex_simple)

    @test result == expected

    # Modern
    # ==========================================================================

    expected = """
\\begin{table}
  \\begin{tabular}{!{\\vrule width 2pt}r!{\\vrule width 1pt}r!{\\vrule width 1pt}r!{\\vrule width 1pt}r!{\\vrule width 2pt}}
    \\noalign{\\hrule height 2pt}
    \\textbf{Col. 1} & \\textbf{Col. 2} & \\textbf{Col. 3} & \\textbf{Col. 4} \\\\\\noalign{\\hrule height 2pt}
    1 & false & 1.0 & 1 \\\\\\noalign{\\hrule height 1pt}
    2 & true & 2.0 & 2 \\\\\\noalign{\\hrule height 1pt}
    3 & false & 3.0 & 3 \\\\\\noalign{\\hrule height 1pt}
    4 & true & 4.0 & 4 \\\\\\noalign{\\hrule height 1pt}
    5 & false & 5.0 & 5 \\\\\\noalign{\\hrule height 1pt}
    6 & true & 6.0 & 6 \\\\\\noalign{\\hrule height 2pt}
  \\end{tabular}
\\end{table}
"""

    result = pretty_table(String, data,
                          hlines = :all,
                          vlines = :all,
                          tf = tf_latex_modern)

    @test result == expected

    # booktabs
    # ==========================================================================

    expected = """
\\begin{table}
  \\begin{tabular}{rrrr}
    \\toprule
    \\textbf{Col. 1} & \\textbf{Col. 2} & \\textbf{Col. 3} & \\textbf{Col. 4} \\\\\\midrule
    1 & false & 1.0 & 1 \\\\\\midrule
    2 & true & 2.0 & 2 \\\\\\midrule
    3 & false & 3.0 & 3 \\\\\\midrule
    4 & true & 4.0 & 4 \\\\\\midrule
    5 & false & 5.0 & 5 \\\\\\midrule
    6 & true & 6.0 & 6 \\\\\\bottomrule
  \\end{tabular}
\\end{table}
"""

    result = pretty_table(String, data,
                          hlines = :all,
                          vlines = :all,
                          tf = tf_latex_booktabs)

    @test result == expected
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct PinsApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_pinsAdd(_api::PinsApi, channel::String; timestamp=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/pins.add", ["slackAuth"])
    Swagger.set_param(_ctx.form, "channel", channel)  # type String
    Swagger.set_param(_ctx.form, "timestamp", timestamp)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Pins an item to a channel.
Param: channel::String (required)
Param: timestamp::String
Return: Dict{String, Any}
"""
function pinsAdd(_api::PinsApi, channel::String; timestamp=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_pinsAdd(_api, channel; timestamp=timestamp, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function pinsAdd(_api::PinsApi, response_stream::Channel, channel::String; timestamp=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_pinsAdd(_api, channel; timestamp=timestamp, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_pinsList(_api::PinsApi, channel::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Any, "/pins.list", ["slackAuth"])
    Swagger.set_param(_ctx.query, "channel", channel)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Lists items pinned to a channel.
Param: channel::String (required)
Return: Any
"""
function pinsList(_api::PinsApi, channel::String; _mediaType=nothing)
    _ctx = _swaggerinternal_pinsList(_api, channel; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function pinsList(_api::PinsApi, response_stream::Channel, channel::String; _mediaType=nothing)
    _ctx = _swaggerinternal_pinsList(_api, channel; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_pinsRemove(_api::PinsApi, channel::String; timestamp=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/pins.remove", ["slackAuth"])
    Swagger.set_param(_ctx.form, "channel", channel)  # type String
    Swagger.set_param(_ctx.form, "timestamp", timestamp)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Un-pins an item from a channel.
Param: channel::String (required)
Param: timestamp::String
Return: Dict{String, Any}
"""
function pinsRemove(_api::PinsApi, channel::String; timestamp=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_pinsRemove(_api, channel; timestamp=timestamp, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function pinsRemove(_api::PinsApi, response_stream::Channel, channel::String; timestamp=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_pinsRemove(_api, channel; timestamp=timestamp, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export pinsAdd, pinsList, pinsRemove