text:
###
### Structs that define the solver settings
###
abstract type SparseSolver end

# Sequential Least Squares methods
struct slstHard <: SparseSolver
    num_iter::Number
    hard_threshold::Number
    iterate_through_rows::Bool
end
slstHard(num_iter, hard_threshold) = slstHard(num_iter, hard_threshold, true)
slstHard(hard_threshold) = slstHard(1, hard_threshold)

struct slstQuantile <: SparseSolver
    num_iter::Number
    quantile_threshold::Number
    iterate_through_rows::Bool
end
slstQuantile(num_iter, quantile_threshold) =
    slstQuantile(num_iter, quantile_threshold, true)
slstQuantile(quantile_threshold) = slstQuantile(1, quantile_threshold)

struct slstNumber <: SparseSolver
    num_iter::Number
    num_terms::Tuple
    iterate_through_rows::Bool
end
slstNumber(num_iter, num_terms) = slstNumber(num_iter, num_terms, true)
slstNumber(num_terms) = slstNumber(1, num_terms)

struct denseSolver <: SparseSolver end

copy_optimizer(s::SparseSolver, options...) = typeof(s)(options...)

###
### Actually solving
###
"""
Uses SparseSolver objects and implements Sequential Least Squares Thresholding
    in various ways

This is the outer user-facing function that performs checks, namely whether
    recursion over rows is needed
"""
function sparse_regression(alg::SparseSolver, X::Matrix, y::Matrix)
                           # num_iter=2,
                           # quantile_threshold=0.1,
                           # hard_threshold=nothing,
                           # num_terms=nothing)
    # Initial check: should we do recursive?
    if size(y,1) > 1 && alg.iterate_through_rows
        A = zeros(size(y,1), size(X,1))
        # For the 2d predictor case, just loop over rows
        for i in 1:size(y,1)
            # if num_terms !== nothing
            #     # Constrain the number of terms row by row
            #     if length(num_terms)>1
            #         n = num_terms[i]
            #     else
            #         n = num_terms
            #     end
            # else
            #     n = nothing
            # end
            # Call inner function
            this_y = y[i:i,:]
            # println(size(this_y))
            A[i, :] = private_sparse_regression(alg, X, this_y, i)#,
                            # num_iter=num_iter,
                            # quantile_threshold=quantile_threshold,
                            # hard_threshold=hard_threshold,
                            # num_terms=n)
            # @show A
        end
    else
        A =  private_sparse_regression(alg, X, y, 0, nothing)
    end
    return A
end

###
### Private functions
###
"""
Fake "sparse" solver; just L2 solution
"""
function private_sparse_regression(alg::denseSolver, X::Matrix, y::Matrix,
                            A=nothing)
    return y/X
end

"""
Implements a hard threshold

see also: threshold_signal!
"""
function private_sparse_regression(alg::slstHard, X::Matrix, y::Matrix,
                            which_data_row::Number, A=nothing)
    if A == nothing
        A = y/X
    end
    # Redo, removing the rows in X that were zeroed out
    for i in 1:alg.num_iter
        threshold_signal!(A, hard_threshold=alg.hard_threshold)
        which_rows = vec(abs.(A) .> 0)
        A_subset = y/(X[which_rows, :])
        A[which_rows] = A_subset
        A[.!which_rows] .= 0.0 # Retain the previous zeros
    end

    return A
end

"""
Implements a quantile threshold

see also: sparsify_signal!
"""
function private_sparse_regression(alg::slstQuantile, X::Matrix, y::Matrix,
                            which_data_row::Number, A=nothing)
    if A == nothing
        A = y/X
    end
    # println(A)
    # Redo, removing the rows in X that were zeroed out
    for i in 1:alg.num_iter
        sparsify_signal!(A, quantile_threshold=alg.quantile_threshold)
        which_rows = vec(abs.(A) .> 0) # Find zeros from sparsification
        A_subset = y/(X[which_rows, :]) # Redo best fit
        A[which_rows] = A_subset # Update non-zero weights
        A[.!which_rows] .= 0.0 # Retain the previous zeros
    end

    return A
end

"""
Implements a hard cap on the number of terms allowed

see also: keep_n_terms
"""
function private_sparse_regression(alg::slstNumber, X::Matrix, y::Matrix,
                            which_data_row::Number, A=nothing)
    if A == nothing
        A = y/X
    end
    # Redo, removing the rows in X that were zeroed out
    for i in 1:alg.num_iter
        A = keep_n_terms(A, alg.num_terms[which_data_row])
        which_rows = vec(abs.(A) .> 0)
        A_subset = y/(X[which_rows, :])
        A[which_rows] = A_subset
        A[.!which_rows] .= 0.0 # Retain the previous zeros
    end

    return A
end

export slstHard, slstQuantile, slstNumber, denseSolver, SparseSolver,
    sparse_regression, copy_optimizer


text:
export FourierOperator, FourierOperatorMetaData
export filter_convolve, convolve, box_filter

import Base: *, ^, +, inv, -
# Derivative Struct and Operations
struct FourierOperator{T, S}
    op::T
    metadata::S
end

struct FourierOperatorMetaData{S}
    name::S
end

FourierOperator(a) = FourierOperator(a, nothing)

function *(∂x::FourierOperator, ϕ::AbstractArray)
    return ∂x.op .* ϕ
end

function (p::FourierOperator)(ϕ::AbstractArray)
    return *(p, ϕ)
end

## Operator Algebra
function ^(p::FourierOperator, α::Number)
    return FourierOperator(p.op .^(α))
end

function *(p::FourierOperator, α::Number)
    return FourierOperator(p.op.*(α))
end
*(q::Number, p::FourierOperator) = *(p, q)

function +(p::FourierOperator, q::FourierOperator)
    return FourierOperator(p.op .+ q.op)
end

function *(p::FourierOperator, q::FourierOperator)
    return FourierOperator(p.op .* q.op)
end

function +(p::FourierOperator, q::Number)
    return FourierOperator(p.op .+ q)
end

function -(p::FourierOperator)
    return FourierOperator( -(p.op))
end

-(p::FourierOperator, q::FourierOperator) = p + (-q)
+(q::Number, p::FourierOperator) = +(p, q)

## with names
function ^(p::FourierOperator{𝒮, 𝒯}, α::Number) where
    {𝒮, 𝒯 <: FourierOperatorMetaData}
    name = p.metadata.name * "^" * string(α)
    fomd = FourierOperatorMetaData(name)
    return FourierOperator(p.op .^(α), fomd)
end

function *(p::FourierOperator{𝒮, 𝒯}, α::Number) where
    {𝒮, 𝒯 <: FourierOperatorMetaData}
    name = string(α) * "*" * p.metadata.name
    fomd = FourierOperatorMetaData(name)
    return FourierOperator(p.op .* α, fomd)
end

function +(p::FourierOperator{𝒮, 𝒯}, q::FourierOperator{𝒮, 𝒯}) where
    {𝒮, 𝒯 <: FourierOperatorMetaData}
    name = "(" * p.metadata.name * "+" * q.metadata.name * ")"
    fomd = FourierOperatorMetaData(name)
    return FourierOperator(p.op .+ q.op, fomd)
end

function *(p::FourierOperator{𝒮, 𝒯}, q::FourierOperator{𝒮, 𝒯}) where
    {𝒮, 𝒯 <: FourierOperatorMetaData}
    name = "(" * p.metadata.name * "*" * q.metadata.name * ")"
    fomd = FourierOperatorMetaData(name)
    return FourierOperator(p.op .* q.op, fomd)
end

function +(p::FourierOperator{𝒮, 𝒯}, q::Number) where
    {𝒮, 𝒯 <: FourierOperatorMetaData}
    name = "(" * p.metadata.name * "+" * string(q) * ")"
    fomd = FourierOperatorMetaData(name)
    return FourierOperator(p.op .+ q, fomd)
end

function -(p::FourierOperator{𝒮, 𝒯}) where
    {𝒮, 𝒯 <: FourierOperatorMetaData}
    name = "-(" * p.metadata.name * ")"
    fomd = FourierOperatorMetaData(name)
    return FourierOperator( -(p.op), fomd)
end

"""
function inv(a::FourierOperator)

# Description
- Takes the inverse of a FourierOperator. If the inverse is zero,
the it returns zero by defualt

# Argument
- `L`: A FourierOperator Object

# Return
- `L⁻¹`: A FourierOperator that represents the inverse of the original
"""
function inv(a::FourierOperator)
    inv_op = 1 ./ a.op
    mask1 = abs.(inv_op) .!= Inf
    mask2 = .!(isnan.(norm.(inv_op)))
    inv_op = mask1 .* mask2 .* inv_op # false * nan or false * inf is false
    return FourierOperator(inv_op)
end

function inv(a::FourierOperator{𝒮, 𝒯}) where
    {𝒮, 𝒯 <: FourierOperatorMetaData}
    name = a.metadata.name * "⁻¹"
    fomd = FourierOperatorMetaData(name)
    inv_op = 1 ./ a.op
    mask1 = abs.(inv_op) .!= Inf
    mask2 = .!(isnan.(norm.(inv_op)))
    inv_op = mask1 .* mask2 .* inv_op # false * nan or false * inf is false
    return FourierOperator(inv_op, fomd)
end

# Filters
function box_filter(û; pm = 6)
    n = length(û)
    mid = div(n, 2) + 1
    pm  = div(n, pm)
    u = copy(û)
    u[mid - pm : mid + pm] .= 0.0
    return u
end

# Convolution operations
function convolve(û, v̂, P, iP)
    u = iP * û
    v = iP * v̂
    return P * (u .* v)
end

function filter_convolve(û, v̂, P, iP; pm = 6)
    u = box_filter(û, pm = pm)
    v = box_filter(v̂, pm = pm)
    u = iP * u
    v = iP * v
    return P * (u .* v)
end

function convolve!(u, v, w, ŵ, û, v̂, P, iP)
    mul!(u , iP , û)
    mul!(v , iP , v̂)
    @inbounds for i in eachindex(w)
        w[i] = u[i] * v[i]
    end
    mul!(ŵ, P, w)
    return nothing
end

##
function Base.show(io::IO, O::FourierOperator{𝒮, 𝒯}) where
    {𝒮, 𝒯 <: FourierOperatorMetaData}
    printstyled(io, O.metadata.name, color = 159)
end

ext: jl

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
# Default visualization functions

function cacheVariablePointEst!(dummyvis,
                                params::Dict{Symbol, Any},
                                fgl::FactorGraph;
                                meanmax=:max  )::Nothing
    #
    cachevars = params[:cachevars]

    # get all variables
    xx, ll = IIF.ls(fgl)
    vars = union(xx, ll)

    # update the variable point-estimate cache
    for vsym in vars

        # get vertex and estimate from the factor graph object
        vert = getVert(fgl, vsym)
        X = getKDE(vert)

        # TODO should be defined through params and not a keyword
        xmx = meanmax == :max ? getKDEMax(X) : getKDEMean(X)

        # get the variable type
        typesym = Caesar.getData(vert).softtype |> typeof |> Symbol

        # cache variable type and estimated value (slightly memory intensive)
        cachevars[vsym] = (typesym, [false;], xmx)
    end

    return nothing
end

"""
    $(SIGNATURES)

Draw variables (triads and points) assing the `cachevars` dictionary is populated with the latest information.  Function intended for internal module use.

`cachevars` === (softtype, [inviewertree=false;], [x,y,theta])
"""
function visualizeVariableCache!(vis::Visualizer,
                                 params::Dict{Symbol, Any},
                                 rose_fgl  )::Nothing
    #

    sessionId="Session"
    if isa(rose_fgl, FactorGraph) && (rose_fgl.sessionname != "" || rose_fgl.sessionname != "NA")
        sessionId = rose_fgl.sessionname
    elseif isa(rose_fgl, FactorGraph)
        sessionId = rose_fgl[2]
    else
        sessionId = params[:sessionId]
    end

    # draw all that the cache requires
    for (vsym, valpair) in params[:cachevars]
        # TODO -- consider upgrading to MultipleDispatch with actual softtypes
        if valpair[1] == :Point2
            visPoint2!(vis, sessionId, vsym, valpair[3],  updateonly=valpair[2][1])
        elseif valpair[1] == :Pose2
            visPose2!(vis, sessionId, vsym, valpair[3],  updateonly=valpair[2][1])
        elseif valpair[1] == :Point3
            visPoint3!(vis, sessionId, vsym, valpair[3],  updateonly=valpair[2][1])
        elseif valpair[1] == :Pose3
            visPose3!(vis, sessionId, vsym, valpair[3],  updateonly=valpair[2][1])
        else
            error("Unknown softtype symbol to visualize from cache.")
        end
        valpair[2][1] = true
    end

    return nothing
end


text:
#
# Generate an initial point in a grid, with a minimum distance defined
#

import Random

# Function that introduces a perturbation in the grid point
perturbation(step,σ) = -σ + 2*σ*rand()

function initial_point(N,side,tol;seed=123)

  # Use a different seed if desired
  if seed > 0
    Random.seed!(seed)
  else
    Random.seed!()
  end
  p = [ zeros(2) for i in 1:N ]

  # Properties of the regular grid
  l = ceil(Int,sqrt(N))
  step = (side-tol)/l
  if tol > step
    error("tol must be smaller than (side-tol)/√N = $step")
  end

  # The grid has l^2 points, but this is generally greater
  # than N. Thus, we generate list of points to be skiped
  skip = Int[]
  while length(skip) < l^2 - N
    i = rand(1:l^2)
    if findfirst(isequal(i), skip) == nothing
      push!(skip,i)
    end
  end

  # Maximum perturbation
  σ = (step-tol)/2

  # Generate points
  ip = 0
  igrid = 0
  x = [ tol/2, tol/2 ]
  for i in 1:l
    for j in 1:l
      igrid += 1
      if findfirst(isequal(igrid),skip) != nothing
        continue
      end
      # update x, with a perturbation of ±0.1*step
      ip += 1
      p[ip] .= x .+ perturbation(step,σ)
      x[2] = x[2] + step
    end
    x[2] = tol/2
    x[1] = x[1] + step
  end
  return p
end

# From data
initial_point(data::Data;tol=0.9,seed=123) = initial_point(data.N,data.side,tol,seed=seed)


text:
#this simple script gets IRS pdf data

function downloadIRSPDF(url::String, fileName::String, outPath::String,
  outName::String = fileName;
  refreshData::Bool=true, maxAttempts::Int = 12)::Nothing

  if !refreshData && isfile("$outPath\\$(outName).pdf")# if the file exists
    return nothing
  end

  attempt::Int = 0
  retry::Bool = true
  while attempt < maxAttempts && retry #loops and tries to download all the files
    attempt+=1
    retry = false
    #download the file
    oStream::IOStream = open("$outPath\\$outName.pdf", "w")
    #weird error if using the transcode codec

    #acquire the data
    try
      HTTP.open(io::IO->
        while !eof(io)
          write(oStream, readavailable(io))
        end,
        "GET", "$url/$fileName", retries=4, connection_limit=16,
        pipeline_limit=32, readtimeout=120)
      close(oStream)
    catch err
      @warn "DL Error with $outName. $(
        attempt<maxAttempts ? "Retrying ($attempt/$maxAttempts)" : "Max attempts Reached for file $outName.pdf.")
        \nErrorMessage:\n$err"
      retry = true
    end

    println("Downloaded IRS file $outName.")
  end

  return nothing
end

#vectorized version (useful for parallelization)
function downloadIRSPDFs(urls::Vector{String}, fileNames::Vector{String},
  outPath::String, outNames::Vector{String}, mode::Symbol = :serial,
  ;delay::Float64 = 0.0, refreshData::Bool=true)::Nothing

  #short-circuit for parallel operation
  if mode == :parallel
    return downloadIRSPDFs(urls, fileNames, outPath, outNames,
    delay=delay, refreshData=refreshData)
  end

  for i::Int ∈ 1:length(urls)

    if delay ≥ .001 #allows for a delay to avoid slamming the site
      sleep(delay)
    end

    downloadIRSPDF(urls[i], fileNames[i], outPath, outNames[i],
      refreshData=refreshData)
  end

  return nothing
end

#useful for parallelization
function downloadIRSPDFsPar(urls::Vector{String}, fileNames::Vector{String},
  outPath::String, outNames::Vector{String};
  delay::Float64 = 0.0, refreshData::Bool=true)::Nothing

  numFiles::Int = length(urls)
  pids = workers()
  np::Int = length(pids)

  #divide up the work
  assignments::Vector{Vector{Int}} = ((i::Int)->Vector{Int}()).(1:np)
  for i::Int ∈ 1:numFiles
    push!(assignments[i%np+1], i) #assigns a file to a pid
  end

  #set up a vector for pushing and receiving assignments
  futures::Vector{Future} = Vector{Future}()
  sizehint!(futures, np)

  for i::Int ∈ 1:np
    if length(assignments[i]) > 0
      push!(futures,
        (@spawnat pids[i] downloadIRSPDFs(urls[assignments[i]],
          fileNames[assignments[i]], outPath,
          outNames[assignments[i]], delay=delay,
          refreshData=refreshData)))
    end
  end

  #get the work
  (fetch).(futures)

  return nothing
end

function readIRSPDFsIndex(;inPath::String = DOCUMENTATION_PATH,
  filesIndexName::String = IRS_INDEX_NAME)::DataFrame
  return CSV.read("$inPath\\$filesIndexName.csv",# strings=:raw,
    categorical=false)
end

function downloadIRSPDFs(;outPath::String = IRS_PATH, mode::Symbol=:parallel,
  delay::Float64=0.0, refreshData::Bool=true)

  #read in the metadata
  fileDF::DataFrame = readIRSPDFsIndex()
  numFiles::Int = size(fileDF,1)

  #build the file names and urls
  fileNames::Vector{String} =  fileDF[:filePrefix]
  fileNames .*= ((string).(fileDF[:urlYears]))
  fileNames .*= ".pdf"
  urls::Vector{String} = fileDF[:urlStems]

  #now build the outnames
  outNames::Vector{String} = fileDF[:outPrefix]
  outNames .*= "_" .* ((string).(fileDF[:urlYears]))
  outNames .*= "_" .* fileDF[:outSuffix]

  return downloadIRSPDFs(urls, fileNames, outPath, outNames, mode, delay=delay,
    refreshData=refreshData)

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*--------------------------------------------------------------------
#* Models to simulate a power plant.
#*--------------------------------------------------------------------
#* Author: Argimiro R. Secchi
#* $Id: power_plant.mso 195 2007-03-07 20:30:12Z arge $
#*-------------------------------------------------------------------
# Declaracao de tipos
type CorrenteVap
	CorrenteVap()=begin
		PP=outers.PP
		new(
			CorrenteZ(),
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties"
			)),
			[
				:(H = PP.VapourEnthalpy(T, P, z)),
				:(S = PP.VapourEntropy(T, P, z)),
			],
			[
				"Vapour Enthalpy","Vapour Entropy",
			],
			[:PP,],
		)
	end
	_base_1::CorrenteZ
	PP::DanaPlugin
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export CorrenteVap
function setEquationFlow(in::CorrenteVap)
	addEquation(1)
	addEquation(2)
end
function atributes(in::CorrenteVap,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Vapour Material Stream"
	fields[:Info]="Model for vapour material streams.
	This model should be used only when the phase of the stream
	is known ''a priori''."
	drive!(fields,_)
	return fields
end
CorrenteVap(_::Dict{Symbol,Any})=begin
	newModel=CorrenteVap()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using VulkanGen

const destfile = joinpath(dirname(dirname(@__DIR__)), "generated", "vulkan_wrapper.jl")
const docfile = joinpath(dirname(dirname(@__DIR__)), "generated", "vulkan_docs.jl")

vw = VulkanWrapper()
@info "Vulkan successfully wrapped."
write(vw, destfile, docfile)


text:
# maps.jl

"""
A map is any transformation of the form `y = f(x)`, where `x` has type `S` and
`y` has type `T`.
"""
abstract type AbstractMap{S,T}
end

domaintype(map::AbstractMap) = domaintype(typeof(map))
domaintype(::Type{AbstractMap{S,T}}) where {S,T} = S
domaintype(::Type{M}) where {M <: AbstractMap} = domaintype(supertype(M))

codomaintype(map::AbstractMap) = codomaintype(typeof(map))
codomaintype(::Type{AbstractMap{S,T}}) where {S,T} = T
codomaintype(::Type{M}) where {M <: AbstractMap} = codomaintype(supertype(M))

domain_space(map::AbstractMap) = GeometricSpace{domaintype(map)}()
range_space(map::AbstractMap) = GeometricSpace{codomaintype(map)}()

(*)(map::AbstractMap, x) = applymap(map, x)

(\)(map::AbstractMap, y) = inv(map) * y

"""
Apply the inverse of the given map. By default, this constructs the inverse of
`m`. However, some maps know its own inverse and may provide a more efficient
implementation of `apply_inverse`.
"""
apply_inverse(m::AbstractMap, y) = applymap(inv(m), y)

apply_left_inverse(m::AbstractMap, y) = applymap(left_inverse(m), y)
apply_right_inverse(m::AbstractMap, y) = applymap(right_inverse(m), y)

# Note that there is a difference between a map being invertible, and the map
# knowing its inverse explicitly and implementing `inv`.
inv(m::AbstractMap) = error("Map ", m, " does not have a known inverse.")

"""
Return a left inverse of the given map. This left inverse `mli` is not unique,
but in any case it is such that `(mli ∘ m) * x = x` for each `x` in the domain
of `m`.
"""# By default we invoke inv
left_inverse(m::AbstractMap) = inv(m)

"""
Return a right inverse of the given map. This right inverse `mri` is not unique,
but in any case it is such that `(m ∘ mri) * y = y` for each `y` in the range
of `m`.
"""# By default we invoke inv
right_inverse(m::AbstractMap) = inv(m)

isreal(m::AbstractMap) = isreal(domaintype(m)) && isreal(codomaintype(m))
isreal(::Type{Complex{T}}) where {T<:Real} = false
isreal(::Type{T}) where {T<:Real} = true
isreal(::Type{SVector{N,T}}) where {N,T} = isreal(T)
isreal(::Type{NTuple{N,T}}) where {N,T} = isreal(T)

"""
`return_type(map, U)` is a generic function that computes the return type when
the given map is applied to a variable of type `U`.

For any `AbstractMap{S,T}`, we have that `return_type(map, S) = T`. We also return
`T` if `U` can be promoted to `T` using Julia's promotion system.
"""# - If the map maps S to T and x has type S, then the result is T
return_type(map::AbstractMap{S,T}, ::Type{S}) where {S,T} = T
# - If x has a type different from S, perhaps it can be promoted to S?
return_type(map::AbstractMap{S,T}, ::Type{U}) where {S,T,U} = _return_type(map, T, promotes_to(U,S))
# -- Yes, it can, let's return T.
_return_type(map, ::Type{T}, ::Type{True}) where {T} = T
# -- No, it can not. We return Any.
_return_type(map, ::Type{T}, ::Type{False}) where {T} = Any

# The expected type of the Jacobian of a map for some combination of types
jac_type(::Type{S},::Type{T}) where {S<:Number,T<:Number} = T
jac_type(::Type{S},::Type{SVector{M,T}}) where {S<:Number,M,T} = SMatrix{M,1,T,M}
jac_type(::Type{SVector{N,S}},::Type{T}) where {N,S,T<:Number} = SVector{N,S}
jac_type(::Type{SVector{N,S}},::Type{SVector{M,T}}) where {N,S,M,T} = SMatrix{M,N,T,M*N}


text:
# ---
# title: 946. Validate Stack Sequences
# id: problem946
# author: zhwang
# date: 2022-03-11
# difficulty: Medium
# categories: Stack
# link: <https://leetcode.com/problems/validate-stack-sequences/description/>
# hidden: true
# ---
#
# Given two sequences `pushed` and `popped`  **with distinct values** , return
# `true` if and only if this could have been the result of a sequence of push
# and pop operations on an initially empty stack.
#
#
#
# **Example 1:**
#
#
#
#     Input: pushed = [1,2,3,4,5], popped = [4,5,3,2,1]
#     Output: true
#     Explanation: We might do the following sequence:
#     push(1), push(2), push(3), push(4), pop() -> 4,
#     push(5), pop() -> 5, pop() -> 3, pop() -> 2, pop() -> 1
#
#
# **Example 2:**
#
#
#
#     Input: pushed = [1,2,3,4,5], popped = [4,3,5,1,2]
#     Output: false
#     Explanation: 1 cannot be popped before 2.
#
#
#
#
# **Constraints:**
#
#   * `0 <= pushed.length == popped.length <= 1000`
#   * `0 <= pushed[i], popped[i] < 1000`
#   * `pushed` is a permutation of `popped`.
#   * `pushed` and `popped` have distinct values.
#
#
## @lc code=start
using LeetCode
function validate_stack_sequences(pushed::Vector{Int}, popped::Vector{Int})
    return validate_stack_sequences!(copy(pushed), popped)
end
function validate_stack_sequences!(pushed::Vector{Int}, popped::Vector{Int})::Bool
    stack = Int[]
    for num in popped
        while isempty(stack) || last(stack) != num
            isempty(pushed) && return false
            push!(stack, popfirst!(pushed))
        end
        pop!(stack)
    end
    return true
end

## @lc code=end


text:
struct CuthillMcKeePermuter
    reverse::Bool
    assumesymmetry::Bool
    checksymmetry::Bool
end

function symrcm(A::SparseMatrixCSC{Tv, Ti}) where {Tv, Ti}
    @inbounds begin
        m, n = size(A)

        deg = undefs(Ti, max(m + 2, n))
        bkt = undefs(Ti, n)
        ord = undefs(Ti, n)
        deg[1] = 1
        for i = 2:m + 2
            deg[i] = 0
        end
        for j = 1:n
            deg[A.colptr[j + 1] - A.colptr[j] + 2] += 1
        end
        for i = 1:m + 1
            deg[i + 1] += deg[i]
        end
        for j = 1:n
            d = A.colptr[j + 1] - A.colptr[j] + 1
            k = deg[d]
            ord[k] = j
            deg[d] = k + 1
        end

        k_start = 1

        prm = undefs(Ti, n)
        vst = falses(n)
        k_current = 0
        k_frontier = k_current
        while k_frontier < n
            k_current += 1
            if k_current > k_frontier
                while vst[ord[k_start]]
                    k_start += 1
                end
                j = ord[k_start]
                prm[k_current] = j
                k_frontier = k_current
                vst[j] = true
            else
                j = prm[k_current]
            end
            k_frontier′ = k_frontier
            for q = A.colptr[j]:A.colptr[j + 1] - 1
                i = A.rowval[q]
                if !vst[i]
                    k_frontier′ += 1
                    prm[k_frontier′] = i
                    vst[i] = true
                end
            end

            if k_frontier′ > k_frontier
                let i = prm[k_frontier′]
                    deg[k_frontier′] = A.colptr[i + 1] - A.colptr[i]
                    bkt[k_frontier′] = k_frontier′
                end
                for k = k_frontier′-1:-1:k_frontier + 1
                    i = prm[k]
                    d = A.colptr[i + 1] - A.colptr[i]
                    while k != k_frontier′ && d > deg[k + 1]
                        deg[k] = deg[bkt[k + 1]]
                        prm[k] = prm[bkt[k + 1]]
                        bkt[k] = bkt[k + 1] - (d != deg[k + 1])
                        k = bkt[k + 1]
                    end
                    prm[k] = i
                    deg[k] = d
                    bkt[k] = k
                end

                k_frontier = k_frontier′
            end
        end
        return prm
    end
end

function rctrcm(A::SparseMatrixCSC{Tv, Ti}, adj_A::SparseMatrixCSC{Tv, Ti}) where {Tv, Ti}
    @inbounds begin
        m, n = size(A)

        deg = undefs(Ti, max(m + 2, n + 2, m + n))
        bkt = undefs(Ti, m + n)
        ord = undefs(Ti, m + n)
        deg[1] = 1
        for v = 2:max(n + 2, m + 2)
            deg[v] = 0
        end
        for i = 1:m
            deg[adj_A.colptr[i + 1] - adj_A.colptr[i] + 2] += 1
        end
        for j = 1:n
            deg[A.colptr[j + 1] - A.colptr[j] + 2] += 1
        end
        for v = 1:max(m + 1, n + 1)
            deg[v + 1] += deg[v]
        end
        for i = 1:m
            d = adj_A.colptr[i + 1] - adj_A.colptr[i] + 1
            k = deg[d]
            ord[k] = i
            deg[d] = k + 1
        end
        for j = 1:n
            d = A.colptr[j + 1] - A.colptr[j] + 1
            k = deg[d]
            ord[k] = m + j
            deg[d] = k + 1
        end

        k_start = 1

        prm = undefs(Ti, m + n)
        vst = falses(m + n)
        k_current = 0
        k_frontier = k_current
        while k_frontier < m + n
            k_current += 1
            if k_current > k_frontier
                while vst[ord[k_start]]
                    k_start += 1
                end
                v = ord[k_start]
                prm[k_current] = v
                k_frontier = k_current
                vst[v] = true
            else
                v = prm[k_current]
            end
            k_frontier′ = k_frontier
            if v <= m
                i = v
                for q = adj_A.colptr[i]:adj_A.colptr[i + 1] - 1
                    j = adj_A.rowval[q]
                    if !vst[m + j]
                        k_frontier′ += 1
                        prm[k_frontier′] = m + j
                        vst[m + j] = true
                    end
                end
            else
                j = v - m
                for q = A.colptr[j]:A.colptr[j + 1] - 1
                    i = A.rowval[q]
                    if !vst[i]
                        k_frontier′ += 1
                        prm[k_frontier′] = i
                        vst[i] = true
                    end
                end
            end

            if k_frontier′ > k_frontier
                let v = prm[k_frontier′]
                    if v <= m
                        i = v
                        deg[k_frontier′] = adj_A.colptr[i + 1] - adj_A.colptr[i]
                    else
                        j = v - m
                        deg[k_frontier′] = A.colptr[j + 1] - A.colptr[j]
                    end
                    bkt[k_frontier′] = k_frontier′
                end
                for k = k_frontier′-1:-1:k_frontier + 1
                    v = prm[k]
                    if v <= m
                        i = v
                        d = adj_A.colptr[i + 1] - adj_A.colptr[i]
                    else
                        j = v - m
                        d = A.colptr[j + 1] - A.colptr[j]
                    end
                    while k != k_frontier′ && d > deg[k + 1]
                        deg[k] = deg[bkt[k + 1]]
                        prm[k] = prm[bkt[k + 1]]
                        bkt[k] = bkt[k + 1] - (d != deg[k + 1])
                        k = bkt[k + 1]
                    end
                    prm[k] = v
                    deg[k] = d
                    bkt[k] = k
                end

                k_frontier = k_frontier′
            end
        end
        return prm
    end
end

CuthillMcKeePermuter(; reverse=false, assumesymmetry=false, checksymmetry=true) = CuthillMcKeePermuter(reverse, assumesymmetry, checksymmetry)

function permute_stripe(A::Symmetric{Tv, SparseMatrixCSC{Tv, Ti}}, method::CuthillMcKeePermuter) where {Tv, Ti}
    prm = symrcm(SparseMatrixCSC(A))
    if method.reverse
        reverse!(prm)
    end
    return DomainPermutation(prm)
end

function permute_stripe(A::SparseMatrixCSC, method::CuthillMcKeePermuter)
    return permute_plaid(A, method)[end]
end

function permute_plaid(A::Symmetric{Tv, SparseMatrixCSC{Tv, Ti}}, method::CuthillMcKeePermuter) where {Tv, Ti}
    prm = symrcm(SparseMatrixCSC(A))
    if method.reverse
        reverse!(prm)
    end
    return (DomainPermutation(prm), DomainPermutation(prm))
end

function permute_plaid(A::SparseMatrixCSC{Tv, Ti}, method::CuthillMcKeePermuter; adj_A = nothing) where {Tv, Ti}
    @inbounds begin
        if method.assumesymmetry || (method.checksymmetry && issymmetric(pattern(A)))
            prm = symrcm(A)
            if method.reverse
                reverse!(prm)
            end
            return (DomainPermutation(prm), DomainPermutation(prm))
        else
            m, n = size(A)
            if adj_A === nothing
                adj_A = adjointpattern(A)
            end
            prm = rctrcm(A, adj_A)
            iprm = undefs(Ti, m)
            jprm = undefs(Ti, n)
            if method.reverse
                i = m
                j = n
                for v in prm
                    if v <= m
                        iprm[i] = v
                        i -= 1
                    else
                        jprm[j] = v - m
                        j -= 1
                    end
                end
            else
                i = 1
                j = 1
                for v in prm
                    if v <= m
                        iprm[i] = v
                        i += 1
                    else
                        jprm[j] = v - m
                        j += 1
                    end
                end
            end
            return (DomainPermutation(iprm), DomainPermutation(jprm))
        end
    end
end


text:
using GeoStats
using ImageQuilting
using GeoStatsImages
using Plots; gr()
using Base.Test
using VisualRegressionTests

# setup GR backend for Travis CI
ENV["GKSwstype"] = "100"
ENV["PLOTS_TEST"] = "true"

# list of maintainers
maintainers = ["juliohm"]

# environment settings
ismaintainer = "USER" ∈ keys(ENV) && ENV["USER"] ∈ maintainers
istravislinux = "TRAVIS" ∈ keys(ENV) && ENV["TRAVIS_OS_NAME"] == "linux"
datadir = joinpath(@__DIR__,"data")

@testset "ImageQuilting.jl" begin
  @testset "Basic checks" begin
    # the output of a homogeneous image is also homogeneous
    TI = ones(20,20,20)
    reals = iqsim(TI, 10, 10, 10, size(TI)...)
    @test reals[1] == TI

    # categories are obtained from training image only
    ncateg = 3; TI = rand(RandomDevice(), 0:ncateg, 20, 20, 20)
    reals = iqsim(TI, 10, 10, 10, size(TI)..., simplex=true)
    @test Set(reals[1]) ⊆ Set(TI)
  end

  @testset "Soft data" begin
    # trends with soft data
    TI = [zeros(10,20,1); ones(10,20,1)]
    trend = [zeros(20,10,1) ones(20,10,1)]
    reals = iqsim(TI, 10, 10, 1, size(TI)..., soft=[(trend,TI)], tol=1)
    @test mean(reals[1][:,1:10,:]) ≤ mean(reals[1][:,11:20,:])

    # no side effects with soft data
    TI = ones(20,20,20)
    TI[:,5,:] = NaN
    aux = ones(TI)
    iqsim(TI, 10, 10, 10, size(TI)..., soft=[(aux,aux)])
    @test aux == ones(TI)

    # auxiliary variable with integer type
    TI = ones(20,20,20)
    aux = [i for i in 1:20, j in 1:20, k in 1:20]
    iqsim(TI, 10, 10, 10, size(TI)..., soft=[(aux,aux)])
    @test aux == [i for i in 1:20, j in 1:20, k in 1:20]
  end

  @testset "Hard data" begin
    # hard data is honored everywhere
    TI = ones(20,20,20)
    obs = rand(size(TI))
    data = HardData((i,j,k)=>obs[i,j,k] for i=1:20, j=1:20, k=1:20)
    reals = iqsim(TI, 10, 10, 10, size(TI)..., hard=data)
    @test reals[1] == obs

    # multiple realizations with hard data
    TI = ones(20,20,20)
    data = HardData((20,20,20)=>10)
    reals = iqsim(TI, 10, 10, 10, size(TI)..., hard=data, nreal=3)
    for real in reals
      @test real[20,20,20] == 10
    end
  end

  @testset "Masked grids" begin
    # masked simulation domain
    TI = ones(20,20,20)
    shape = HardData()
    active = trues(TI)
    for i=1:20, j=1:20, k=1:20
      if (i-10)^2 + (j-10)^2 + (k-10)^2 < 25
        push!(shape, (i,j,k)=>NaN)
        active[i,j,k] = false
      end
    end
    reals = iqsim(TI, 10, 10, 10, size(TI)..., hard=shape)
    @test all(isnan.(reals[1][.!active]))
    @test all(.!isnan.(reals[1][active]))

    # masked training image
    TI = ones(20,20,20)
    TI[:,5,:] = NaN
    reals = iqsim(TI, 10, 10, 10, size(TI)...)
    @test reals[1] == ones(TI)
    TI[1,5] = 0
    reals = iqsim(TI, 10, 10, 10, size(TI)..., simplex=true)
    @test reals[1] == ones(TI)

    # masked domain and masked training image
    TI = ones(20,20,20)
    TI[:,5,:] = NaN
    aux = ones(TI)
    shape = HardData((i,j,k)=>NaN for i=1:20, j=5, k=1:20)
    reals = iqsim(TI, 10, 10, 10, size(TI)..., hard=shape)
    @test all(isnan.(reals[1][:,5,:]))
    @test all(reals[1][:,1:4,:] .== 1)
    @test all(reals[1][:,6:20,:] .== 1)
    reals = iqsim(TI, 10, 10, 10, size(TI)..., hard=shape, soft=[(aux,aux)])
    @test all(isnan.(reals[1][:,5,:]))
    @test all(reals[1][:,1:4,:] .== 1)
    @test all(reals[1][:,6:20,:] .== 1)
  end

  @testset "Minimum error cut" begin
    # 3D cut
    TI = ones(20,20,20)
    for cut in [:dijkstra,:boykov]
      _, _, voxs = iqsim(TI, 10, 10, 10, size(TI)..., overlapx=1/3, overlapy=1/3, overlapz=1/3, cut=cut, debug=true)
      @test 0 ≤ voxs[1] ≤ 1
    end
  end

  @testset "Simulation paths" begin
    # different simulation paths
    for kind in [:rasterup,:rasterdown,:dilation,:random]
      path = ImageQuilting.genpath((10,10,10), kind)
      @test length(path) == 1000
    end

    # datum is visited first if present
    path = ImageQuilting.genpath((10,10,10), :datum, [(1,1,1),(10,10,10)])
    @test path[1:2] == [1,1000] || path[1:2] == [1000,1]
  end

  @testset "Voxel reuse" begin
    # mean voxel reuse is in range [0,1]
    TI = rand(20,20,20)
    μ, σ = voxelreuse(TI, 10, 10, 10, nreal=1)
    @test 0 ≤ μ ≤ 1
  end

  # visual regression tests are very hard to get
  # right given the frequent changes in the Julia
  # plotting packages, therefore we only run them
  # on maintainers machines and on Travis CI
  if ismaintainer || istravislinux
    @testset "Visual tests" begin
      for TIname in ["StoneWall","WalkerLake"]
        function plot_voxel_reuse(fname)
          srand(2017)
          TI = training_image(TIname)[1:20,1:20,:]
          voxelreuseplot(TI)
          png(fname)
        end
        refimg = joinpath(datadir, "Voxel"*TIname*".png")

        @test test_images(VisualTest(plot_voxel_reuse, refimg), popup=!istravislinux) |> success
      end

      for TIname in ["Strebelle","StoneWall"]
        function plot_reals(fname)
          srand(2017)
          TI = training_image(TIname)[1:50,1:50,:]
          reals = iqsim(TI, 30, 30, 1, size(TI)..., nreal=4)
          ps = []
          for real in reals
            push!(ps, heatmap(real[:,:,1]))
          end
          plot(ps...)
          png(fname)
        end
        refimg = joinpath(datadir, "Reals"*TIname*".png")

        @test test_images(VisualTest(plot_reals, refimg), popup=!istravislinux) |> success
      end
    end
  end

  @testset "GeoStats.jl API" begin
    geodata = GeoDataFrame(DataFrames.DataFrame(x=[25.,50.,75.], y=[25.,75.,50.], variable=[1.,0.,1.]), [:x,:y])
    grid = RegularGrid{Float64}(100,100)
    problem = SimulationProblem(geodata, grid, :variable, 3)

    TI = training_image("Strebelle")
    inactive = [(i,j,1) for i in 1:30 for j in 1:30]
    solver = ImgQuilt(:variable => @NT(TI=TI, template=(30,30,1), inactive=inactive))

    srand(2017)
    solution = solve(problem, solver)
    @test keys(solution.realizations) ⊆ [:variable]

    incomplete_solver = ImgQuilt()
    @test_throws AssertionError solve(problem, incomplete_solver)

    if ismaintainer || istravislinux
      function plot_solution(fname)
        plot(solution, size=(1000,300))
        png(fname)
      end
      refimg = joinpath(datadir, "GeoStatsAPI.png")

      @test test_images(VisualTest(plot_solution, refimg), popup=!istravislinux) |> success
    end
  end

  if ImageQuilting.cl ≠ nothing && ImageQuilting.clfft ≠ nothing
    @testset "GPU support" begin
      # CPU and GPU give same results
      TI = ones(20,20,20)
      TI[10:end,:,:] = 2
      srand(0); realscpu = iqsim(TI, 10, 10, 10, size(TI)..., gpu=false)
      srand(0); realsgpu = iqsim(TI, 10, 10, 10, size(TI)..., gpu=true)
      @test realscpu[1] == realsgpu[1]
    end
  end
end


text:
infile = "./src/Day4/input.txt"
a = split(join(readlines(infile, keep = true)), "\n\n")

validcount = 0
invalidcount = 0
passport = Set{String}()

# part 1
requiredfields = ["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"] # cid
for row in a
    fields = split(row, ['\n', ' '])
    fieldnames = first.(split.(fields, ":"))
    fieldset = Set(fieldnames)
    if all(in(fieldset).(requiredfields))
        validcount += 1
    else
        invalidcount += 1
    end
end
validcount

# part 2

validcount = 0
invalidcount = 0
# passport = Dict{String, String}()
for row in a
    fields = split(row, ['\n', ' '])
    fieldnames = first.(split.(fields, ":"))
    fieldvals =  last.(split.(fields, ":"))
    fieldset = Set(fieldnames)
    if all(in(fieldset).(requiredfields))
        # then check values
        valid = true
        passport = Dict(fieldnames .=> fieldvals)
        # byr (Birth Year) - four digits; at least 1920 and at most 2002.
        if length(passport["byr"]) != 4
            valid = false
        end
        byr = parse(Int, passport["byr"])
        if (byr < 1920) | (byr > 2020)
            valid = false
        end
        # iyr (Issue Year) - four digits; at least 2010 and at most 2020.
        if length(passport["iyr"]) != 4
            valid = false
        end
        iyr = parse(Int, passport["iyr"])
        if (iyr < 2010) | (iyr > 2020)
            valid = false
        end

        # eyr (Expiration Year) - four digits; at least 2020 and at most 2030.
        if length(passport["eyr"]) != 4
            valid = false
        end
        eyr = parse(Int, passport["eyr"])
        if (eyr < 2020) | (eyr > 2030)
            valid = false
        end
        # hgt (Height) - a number followed by either cm or in:
        # If cm, the number must be at least 150 and at most 193.
        # If in, the number must be at least 59 and at most 76.
        # println(passport["hgt"])
        hgt_string = passport["hgt"]
        if length(hgt_string) > 2
            hgt_units = hgt_string[(end-1):end]
            hgt_value = parse(Int, hgt_string[1:(end-2)])
            if hgt_units == "cm"
                if (hgt_value < 150) | (hgt_value > 193)
                    valid = false
                end
            elseif hgt_units == "in"
                if (hgt_value < 59) | (hgt_value > 76)
                    valid = false
                end
            end
        else
            valid = false
        end
        # hcl (Hair Color) - a # followed by exactly six characters 0-9 or a-f.
        if !occursin(r"^#[0-9a-f]{6}$", passport["hcl"])
            valid = false
        end

        # ecl (Eye Color) - exactly one of: amb blu brn gry grn hzl oth.
        if !(passport["ecl"] in ["amb","blu","brn","gry","grn","hzl","oth"])
            valid = false
        end

        # pid (Passport ID) - a nine-digit number, including leading zeroes.
        if !occursin(r"^\d{9}$", passport["pid"])
            valid = false
        end
        # cid (Country ID) - ignored, missing or not.

        if valid
            validcount += 1
        else
            invalidcount += 1
        end
    else
        invalidcount += 1
    end
end
validcount
invalidcount


text:
using LinearAlgebra;
using Printf;

function linspace(start::Number, stop::Number, n::Int64)
	linx = zeros(Float64, n);
	delta = (stop - start) / (n - 1);
	for i = 1:n
		linx[i] = start + delta * (i - 1)
	end
	return linx;
end

lx = 1.0;
t = 0.0
nx = 11;
dx = lx / (nx - 1);
dt = 1e-2;
tlim = 0.1;
kappa = 1.0
beta = dt * kappa / 2.0dx^2;
t_save = linspace(0.0, tlim, nx);
save_count = 1;

function solve_linear(n::Int64, a::Float64, b::Float64, c::Float64, r::Array)
	A = zeros(n, n);
	for i = 2:n+1
		for j = 2:n+1
			if(i-1 == j-1-1)
				A[i-1, j-1] = c;
			elseif (i-1==j-1)
				A[i-1,j-1] = b;
			elseif (i-1-1 == j-1)
				A[i-1,j-1] = a;
			end
		end
	end
	return A \ r;
end

function crank_nicolson(u::Array)
	a = Float64(0.0);
	b = Float64(0.0);
	c = Float64(0.0);
	r = zeros(Float64, nx - 2);
	for i = 2:nx-1
		a = -beta;
		b = 1.0 + 2.0beta;
		c = -beta;
		r[i-1] = u[i] + beta*(u[i+1] - 2.0u[i] + u[i-1]);
	end
	r[begin] = r[begin] - a*u[begin]
	r[end] = r[end] - c * u[end];
	u[2:end-1] = solve_linear(nx - 2, a, b, c, r)
end

function output(x::Array, u::Array, t::Float64)
	filename = @sprintf("cn_t=%4.3f.dat", t);
	open(filename, "w") do io
		for i = 1:nx
			write(io, @sprintf("%.5f\t%.5f\n", x[i], u[i]));
		end
	end;
end

function main()
	@printf("beta = %.5f, dt = %.5f, dx = %.5f\n", beta, dt, dx);
	u = ones(Float64, nx);
	u[begin] = u[end] = 0.0;
	x = linspace(0, lx, nx);
	while t < tlim
		output(x, u, t);
		crank_nicolson(u);
		global t = t + dt;
	end
end

main()


text:
# mini-max fit to a generic function on the interval [-1,1].
# Ignores the domain outside of this interval.

using Profile

using FFTW

using Distributed
using SharedArrays

@everywhere import DelimitedFiles

@everywhere import Random
@everywhere import LinearAlgebra
@everywhere import PyPlot
@everywhere import Printf

@everywhere import Calculus
@everywhere import ForwardDiff

@everywhere include("../src/misc/declarations.jl")
@everywhere include("../src/minimax/Chebyshev.jl")
#@everywhere include("../src/minimax/APF.jl")
@everywhere include("../src/minimax/interpolators.jl")
@everywhere include("../src/minimax/eval.jl")
@everywhere include("../src/minimax/filter.jl")

@everywhere include("../src/misc/utilities.jl")
@everywhere include("../src/misc/template_signals.jl")
@everywhere include("../src/misc/bounds.jl")
@everywhere include("../src/minimax/exchange.jl")
@everywhere include("../src/minimax/exchange_strategy.jl")
@everywhere include("../src/misc/IO.jl")

@everywhere include("../src/minimax/engine.jl")
@everywhere include("../src/minimax/extrema_search.jl")

@everywhere include("../src/FIR/design.jl")
@everywhere include("../src/FIR/type_conversion.jl")

@everywhere include("../src/misc/frontend.jl")
@everywhere include("../test/visualize.jl")

fig_num = 1
PyPlot.close("all")

Random.seed!(25)

fig_num = 1

L = 25
f = xx->sin(xx*π+3*xx)+xx

𝓧_new = frontendoptimization(f, L)

xq = convert(Vector{BigFloat},collect(LinRange(-2,1,10*L)))
f_itp_xq, f_itp_𝓧, 𝑒_xq, 𝑒_𝓧, fig_num = visualizeoptimzationsolution(xq, 𝓧_new, f, fig_num)

@assert 1==232

L = 25 # 55 # set to odd.
# L_interpolant = L
# if isodd(L)
#     L_interpolant += 1
# end
L_interpolant = 4

max_iters = 200
filter_type_num = 2
tol_converged = convert(BigFloat,1e-3)
tol_𝓧_spacing = convert(BigFloat,1e-6/(L+2))
tol_derivative_zero = convert(BigFloat,0.2)
tol_no_update = convert(BigFloat,1e-6)
tol_h = convert(BigFloat,1e-17)

N_samples_interval = 10
N_test_positions = 200

verbose_flag = true
plot_flag = false #true

savefig_flag = true
save_delay = 0.5

config = MiniMaxConfigType( tol_converged,
                            tol_𝓧_spacing,
                            tol_derivative_zero,
                            tol_no_update,
                            tol_h,
                            N_samples_interval,
                            N_test_positions,
                            max_iters,
                            L_interpolant,
                            verbose_flag,
                            plot_flag)

candidate_multiple = 20

## set up functions.
wfunc = getweightfunc(filter_type_num)

zfunc = xx->sin(xx*π+3*xx)+xx
dfunc = xx->zfunc(acos(xx))
#dfunc = xx->sin(xx*π+3*xx)+xx

#dfunc = xx->rectfunc(xx, BigFloat("0.5"))
#dfunc = xx->sinc(xx)

# β = BigFloat("1")
# 𝑇 = BigFloat("2")
# dfunc = xx->raisedcosinefreqrsp(xx,β,𝑇)

𝓧_selection_string = "Chebyshev"

@time 𝓧_new, h_history,
        iters_ran, status_msg, fig_num = designfilter(  config,
                                                        L,
                                                        candidate_multiple,
                                                        wfunc,
                                                        dfunc,
                                                        𝓧_selection_string,
                                                        fig_num)
#

println("Status: ", status_msg)

include("visualize.jl")

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
using Surrogates

#1D
lb = 0.0
ub = 5.0
n = 20
x = sample(n,lb,ub,SobolSample())
f = x->2*x+x^2
y = f.(x)
my_ear1d = EarthSurrogate(x,y,lb,ub)
val = my_ear1d(3.0)
add_point!(my_ear1d,6.0,48.0)

#ND
lb = [0.0,0.0]
ub = [10.0,10.0]
n = 30
x = sample(n,lb,ub,SobolSample())
f = x -> x[1]*x[2] + x[1]
y = f.(x)
my_earnd = EarthSurrogate(x,y,lb,ub)
val = my_earnd((2.0,2.0))
add_point!(my_earnd,(2.0,2.0),6.0)


text:
using DRIPs
using Test
using LinearAlgebra
@testset "Endogenous Feedback" begin
    ρ   = 0.6;        #persistence of money growth
    σ_u = 0.1;        #std. deviation of shocks to money growth
    α   = 0.8;        #degree of strategic complementarity
    L   = 40;         #length of truncation
    Hq  = ρ.^(0:L-1); #state-space rep. of Δq

    ## specifying the primitives of the drip
    ω   = 0.2;
    β   = 0.99;
    A   = [1 zeros(1,L-2) 0; Matrix(I,L-1,L-1) zeros(L-1,1)];
    Q   = [σ_u; zeros(L-1,1)];

    ## Function to iterate over ge
    function ge_drip(ω,β,A,Q,          #primitives of drip except for H because H is endogenous
                     α,                #strategic complementarity
                     Hq,               #state space rep. of Δq
                     L;                #length of truncation
                     H0       = Hq,    #optional: initial guess for H (Hq is the true solution when α=0)
                     maxit    = 200,   #optional: max number of iterations for GE code
                     tol      = 1e-4)  #optional: tolerance for iterations
        # set primitives
        err   = 1;
        iter  = 0;
        M     = [zeros(1,L-1) 0; Matrix(I,L-1,L-1) zeros(L-1,1)];
        # iterate on GE
        while (err > tol) & (iter < maxit)
                if iter == 0
                    global ge  = Drip(ω,β,A,Q,H0, w = 0.9);
                else
                    global ge  = Drip(ω,β,A,Q,H0;Ω0 = ge.ss.Ω ,Σ0 = ge.ss.Σ_1);
                end

                XFUN(jj) = ((I-ge.ss.K*ge.ss.Y')*ge.A)^jj * (ge.ss.K*ge.ss.Y') * (M')^jj
                X = DRIPs.infinitesum(XFUN; maxit=L, start = 0);  #E[x⃗]=X×x⃗

                XpFUN(jj) = α^jj * X^(jj)
                Xp = DRIPs.infinitesum(XpFUN; maxit=L, start = 0);

                H1 = (1-α)*Xp'*Hq;
                err= 0.5*norm(H1-H0,2)/norm(H0)+0.5*err;
                H0 = H1;
                iter += 1;
        end
        return(ge,err)
    end;
    println("Solving RI with Endogenous Feedback ...")
    ge,err = ge_drip(ω,β,A,Q,α,Hq,L) # remove suppress to see convergence log
    @test ge.ss.err < 1e-4
    @test err < 1e-4
    ge,err = ge_drip(0.001,β,A,Q,α,Hq,L) # remove suppress to see convergence log
    @test ge.ss.err < 1e-4
    @test err < 1e-4
    ge,err = ge_drip(0.1,β,A,Q,α,Hq,L) # remove suppress to see convergence log
    @test ge.ss.err < 1e-4
    @test err < 1e-4
    ge,err = ge_drip(1,β,A,Q,α,Hq,L) # remove suppress to see convergence log
    @test ge.ss.err < 1e-4
    @test err < 1e-4
end


text:
import Revise
import GLRenderer as GL
import ThreeDP3 as T
import Images as I
import MiniGSG as S
import Rotations as R
import PoseComposition: Pose
using Gen
import MeshCatViz

import MeshCatViz
MeshCatViz.setup_visualizer()

YCB_DIR = "/home/nishadg/mcs/ThreeDVision.jl/data/ycbv2"
world_scaling_factor = 100.0
id_to_cloud, id_to_shift, id_to_box  = T.load_ycbv_models_adjusted(YCB_DIR, world_scaling_factor);
all_ids = sort(collect(keys(id_to_cloud)));

IDX =400
@show T.get_ycb_scene_frame_id_from_idx(YCB_DIR,IDX)
gt_poses, ids, gt_rgb_image, gt_depth_image, cam_pose, camera = T.load_ycbv_scene_adjusted(
    YCB_DIR, IDX, world_scaling_factor, id_to_shift
);
img = I.colorview(I.Gray, gt_depth_image ./ maximum(gt_depth_image))
rgb_image = I.colorview(I.RGB, permutedims(Float64.(gt_rgb_image)./255.0, (3,1,2)))

nominal_colors = [I.colorant"red",I.colorant"blue",I.colorant"green",I.colorant"yellow",I.colorant"white",
I.colorant"black"];
diffs = cat([I.colordiff.(rgb_image, c) for c in nominal_colors]...,dims=3);

argmaxs = map(x->x[3],argmin(diffs, dims=3))[:,:,1];

mod_image = copy(rgb_image)
for (idx,c) in enumerate(nominal_colors)
    mod_image[argmaxs .== idx] .= c
end
mod_image

import Plots as P

P.heatmap

# +
resolution = 0.5

renderer = GL.setup_renderer(camera, GL.DepthMode())
for id in all_ids
    cloud = id_to_cloud[id]
    v,n,f = GL.mesh_from_voxelized_cloud(GL.voxelize(cloud, resolution), resolution)
    GL.load_object!(renderer, v, f)
end
depth_image = GL.gl_render(
    renderer, ids, gt_poses, T.IDENTITY_POSE)
depth_image[depth_image .== 50000.0] .= 200.0
img = I.colorview(I.Gray, depth_image ./ maximum(depth_image))
# -

renderer_color = GL.setup_renderer(camera, GL.RGBMode())
for id in all_ids
    cloud = id_to_cloud[id]
    v,n,f = GL.mesh_from_voxelized_cloud(GL.voxelize(cloud, resolution), resolution)
    GL.load_object!(renderer_color, v, n, f)
end
colors = I.distinguishable_colors(length(ids), I.colorant"green")
# colors = [I.colorant"red", I.colorant"green", I.colorant"cyan"]
colors  = [
    I.colorant"yellow", I.colorant"cyan", I.colorant"lightgreen",
    I.colorant"red", I.colorant"purple", I.colorant"orange"
]
rgb_image, depth_image = GL.gl_render(
    renderer_color, ids, gt_poses, colors[1:length(ids)], T.IDENTITY_POSE)
depth_image[depth_image .== 50000.0] .= 200.0
img = I.colorview(I.Gray, depth_image ./ maximum(depth_image))
I.colorview(I.RGBA, permutedims(rgb_image,(3,1,2)))

gt_poses

renderer_color = GL.setup_renderer(camera, GL.RGBMode())
for id in all_ids
    cloud = id_to_cloud[id]
    v,n,f = GL.mesh_from_voxelized_cloud(GL.voxelize(cloud, resolution), resolution)
    GL.load_object!(renderer_color, v, n, f)
end
colors = I.distinguishable_colors(length(ids), I.colorant"green")
# colors = [I.colorant"red", I.colorant"green", I.colorant"cyan"]
colors  = [
    I.colorant"yellow", I.colorant"cyan", I.colorant"lightgreen",
    I.colorant"red", I.colorant"purple", I.colorant"orange"
]
rand_poses = [T.uniformPose(-20.0,20.0,-20.0,20.0,97.0,120.0) for _ in 1:length(gt_poses)]
rgb_image, depth_image = GL.gl_render(
    renderer_color, ids, rand_poses, colors[1:length(ids)], T.IDENTITY_POSE)
depth_image[depth_image .== 50000.0] .= 200.0
img = I.colorview(I.Gray, depth_image ./ maximum(depth_image))
I.colorview(I.RGBA, permutedims(rgb_image,(3,1,2)))

# +
renderer_texture = GL.setup_renderer(camera, GL.TextureMode())
obj_paths = T.load_ycb_model_obj_file_paths(YCB_DIR)
texture_paths = T.load_ycb_model_texture_file_paths(YCB_DIR)
for id in all_ids
    v,n,f,t = renderer_texture.gl_instance.load_obj_parameters(
        obj_paths[id]
    )
    v = v * world_scaling_factor
    v .-= id_to_shift[id]'

    GL.load_object!(renderer_texture, v, n, f, t,
        texture_paths[id]
    )
end

rgb_image, depth_image = GL.gl_render(
    renderer_texture, ids, gt_poses, T.IDENTITY_POSE)
depth_image[depth_image .== 50000.0] .= 200.0
img = I.colorview(I.Gray, depth_image ./ maximum(depth_image))
I.colorview(I.RGBA, permutedims(rgb_image,(3,1,2)))
# -
focus_idx = 6
rgb_image, depth_image = GL.gl_render(
    renderer_texture, [ids[focus_idx]], [gt_poses[focus_idx]], T.IDENTITY_POSE)
depth_image[depth_image .== 50000.0] .= 200.0
img = I.colorview(I.Gray, depth_image ./ maximum(depth_image))
I.colorview(I.RGBA, permutedims(rgb_image,(3,1,2)))

idxs = [4,6]
rgb_image, depth_image = GL.gl_render(
    renderer_texture, ids[idxs], gt_poses[idxs], T.IDENTITY_POSE)
depth_image[depth_image .== 50000.0] .= 200.0
img = I.colorview(I.Gray, depth_image ./ maximum(depth_image))
I.colorview(I.RGBA, permutedims(rgb_image,(3,1,2)))

focus_idx = 5
rgb_image, depth_image = GL.gl_render(
    renderer_texture, [ids[focus_idx]], [Pose([0.0, 50.0, 80.0], R.RotXYZ(pi/2, pi/2, 0.0))], Pose(zeros(3),R.RotX(-0.5)))
depth_image[depth_image .== 50000.0] .= 200.0
img = I.colorview(I.Gray, depth_image ./ maximum(depth_image))
I.colorview(I.RGBA, permutedims(rgb_image,(3,1,2)))


text:
#=
    This file contains all of the functions that a solver needs to implement
    in order to use LQOI.

       min/max: c'x + x'Qx
    subject to:
                # Variable bounds
                abᵢ <=  xᵢ              <= bbᵢ, i=1,2,...Nx
                # Linear Constraints
                llᵢ <= alᵢ'x             <= ulᵢ, i=1,2,...Nl
                # Quadratic Constraints
                lqᵢ <= aqᵢ' x + x'Qqᵢ' x <= uqᵢ, i=1,2,...Nq
                # SOS1, SOS2 constraints
                # Binary Constraints
                xᵢ ∈ {0, 1}
                # Integer Constraints
                xᵢ ∈ Z
=#

"""
    backend_type(m, ::MOI.AbstractSet)::CChar

An overloadable type for dispatching the appropriate types to the backends.

For example, GLPK.jl uses `'E'` for `a'x=b` constraints, where as Gurobi.jl uses
`==`.

Three are special cases:
 - `Val{:Continuous}`: for the type of a continuous variable
 - `Val{:Upperbound}`: for the upper bound of a variable
 - `Val{:Lowerbound}`: for the lower bound of a variable

### Defaults

    MOI.GreaterThan  - 'G'
    MOI.LessThan     - 'L'
    MOI.EqualTo      - 'E'

    MOI.Zeros        - 'E'
    MOI.Nonpositives - 'L'
    MOI.Nonnegatives - 'G'

    MOI.ZeroOne        - 'B'
    MOI.Integer        - 'I'
    MOI.Semicontinuous - 'S'
    MOI.Semiinteger    - 'N'

    MOI.SOS1 - :SOS1  # '1'
    MOI.SOS2 - :SOS2  # '2'

    Val{:Continuous} - 'C'
    Val{:Upperbound} - 'U'
    Val{:Lowerbound} - 'L'
"""
backend_type(m::LinQuadOptimizer, ::MOI.GreaterThan{T}) where T = Cchar('G')
backend_type(m::LinQuadOptimizer, ::MOI.LessThan{T}) where T    = Cchar('L')
backend_type(m::LinQuadOptimizer, ::MOI.EqualTo{T}) where T     = Cchar('E')
# Implemented separately
# backend_type(m::LinQuadOptimizer, ::MOI.Interval{T}) where T    = Cchar('R')

backend_type(m::LinQuadOptimizer, ::MOI.Zeros)        = Cchar('E')
backend_type(m::LinQuadOptimizer, ::MOI.Nonpositives) = Cchar('L')
backend_type(m::LinQuadOptimizer, ::MOI.Nonnegatives) = Cchar('G')

backend_type(m::LinQuadOptimizer, ::MOI.ZeroOne) = Cchar('B')
backend_type(m::LinQuadOptimizer, ::MOI.Integer) = Cchar('I')

backend_type(m::LinQuadOptimizer, ::MOI.SOS1{T}) where T = :SOS1  # Cchar('1')
backend_type(m::LinQuadOptimizer, ::MOI.SOS2{T}) where T = :SOS2  # Cchar('2')

backend_type(m::LinQuadOptimizer, ::MOI.Semicontinuous{T}) where T = Cchar('S')
backend_type(m::LinQuadOptimizer, ::MOI.Semiinteger{T}) where T    = Cchar('N')

backend_type(m::LinQuadOptimizer, ::Val{:Continuous}) = Cchar('C')
backend_type(m::LinQuadOptimizer, ::Val{:Upperbound}) = Cchar('U')
backend_type(m::LinQuadOptimizer, ::Val{:Lowerbound}) = Cchar('L')

"""
    LinearQuadraticModel(M::Type{<:LinQuadOptimizer}, env)

Initializes a model given a model type `M` and an `env` that might be a `nothing`
for some solvers.
"""
function LinearQuadraticModel end
@deprecate LinQuadModel LinearQuadraticModel

"""
    supported_constraints(m)::Vector{
        Tuple{MOI.AbstractFunction, MOI.AbstractSet}
    }

Get a list of supported constraint types in the model `m`.

For example, `[(LQOI.Linear, LQOI.EQ)]`
"""
function supported_constraints end
@deprecate lqs_supported_constraints supported_constraints

"""
    supported_objectives(m)::Vector{MOI.AbstractScalarFunction}

Get a list of supported objective types in the model `m`.

For example, `[LQOI.Linear, LQOI.Quad]`
"""
function supported_objectives end
@deprecate lqs_supported_objectives supported_objectives

# Constraints

"""
    change_variable_bounds!(m, cols::Vector{Int}, values::Vector{Float64}, senses::Vector)

Change the bounds of the variable. The sense of the upperbound
is given by `backend_type(m, Val{:Upperbound}())`. The sense
of the lowerbound is given by `backend_type(m, Val{:Lowerbound}())`
"""
function change_variable_bounds! end
@deprecate lqs_chgbds! change_variable_bounds!

"""
    get_variable_lowerbound(m, col::Int)::Float64

Get the lower bound of the variable in 1-indexed column `col` of the model `m`.
"""
function get_variable_lowerbound end
@deprecate lqs_getlb get_variable_lowerbound

"""
    get_variable_upperbound(m, col::Int)::Float64

Get the upper bound of the variable in 1-indexed column `col` of the model `m`.
"""
function get_variable_upperbound end
@deprecate lqs_getub get_variable_upperbound

"""
    get_number_linear_constraints(m)::Int

Get the number of linear constraints in the model `m`.
"""
function get_number_linear_constraints end
@deprecate lqs_getnumrows get_number_linear_constraints

"""
    add_linear_constraints!(m, rows::Vector{Int}, cols::Vector{Int},
        coefs::Vector{Float64},
        sense::Vector{Cchar}, rhs::Vector{Float64})::Void

Adds linear constraints of the form `Ax (sense) rhs` to the model `m`.

The A matrix is given in triplet form `A[rows[i], cols[i]] = coef[i]` for all
`i`, and `length(rows) == length(cols) == length(coefs)`.

The `sense` is given by `backend_type(m, set)`.

Ranged constraints (`set=MOI.Interval`) should be added via `add_ranged_constraint!`
instead.
"""
function add_linear_constraints! end
@deprecate lqs_addrows! add_linear_constraints!

"""
    add_ranged_constraints!(m, rows::Vector{Int}, cols::Vector{Int},
        coefs::Vector{Float64}, lowerbound::Vector{Float64}, upperbound::Vector{Float64})

Adds linear constraints of the form `lowerbound <= Ax <= upperbound` to the
model `m`.

The A matrix is given in triplet form `A[rows[i], cols[i]] = coef[i]` for all
`i`,

This is a special case compared to standard `add_linear_constraints!` since it
is often implemented via multiple API calls.
"""
function add_ranged_constraints! end

"""
    modify_ranged_constraint!(m, rows::Vector{Int}, lowerbound::Vector{Float64}, upperbound::Vector{Float64})

Modify the lower and upperbounds of a ranged constraint in the model `m`.

This is a special case compared to standard the `change_coefficient!` since it
is often implemented via multiple API calls.
"""
function modify_ranged_constraints! end

"""
    get_rhs(m, row::Int)::Float64

Get the right-hand side of the linear constraint in the 1-indexed row `row` in
the model `m`.
"""
function get_rhs end
@deprecate lqs_getrhs get_rhs

"""
    get_linear_constraint(m, row::Int)::Tuple{Vector{Int}, Vector{Float64}}

Get the linear component of the constraint in the 1-indexed row `row` in
the model `m`. Returns a tuple of `(cols, vals)`.
"""
function get_linear_constraint end
@deprecate lqs_getrows get_linear_constraint

"""
    change_coefficient!(m, row, col, coef)

Set the linear coefficient of the variable in column `col`, constraint `row` to
`coef`.
"""
function change_coefficient! end
@deprecate lqs_chgcoef! change_coefficient!

"""
    delete_linear_constraints!(m, start_row::Int, end_row::Int)::Void

Delete the linear constraints `start_row`, `start_row+1`, ..., `end_row` from
the model `m`.
"""
function delete_linear_constraints! end
@deprecate lqs_delrows! delete_linear_constraints!

"""
    delete_quadratic_constraints!(m, start_row::Int, end_row::Int)::Void

Delete the quadratic constraints `start_row`, `start_row+1`, ..., `end_row` from
the model `m`.
"""
function delete_quadratic_constraints! end

"""
    lqs_chgctype(m, cols::Vector{Int}, types):Void

Change the variable types. Type is the output of one of:
 - `backend_type(m, ::ZeroOne)`, for binary variables;
 - `backend_type(m, ::Integer)`, for integer variables; and
 - `backend_type(m, Val{:Continuous}())`, for continuous variables.
"""
function change_variable_types! end
@deprecate lqs_chgctype! change_variable_types!

"""
    change_linear_constraint_sense!(m, rows::Vector{Int}, sense::Vector{Symbol})::Void

Change the sense of the linear constraints in `rows` to `sense`.

`sense` is the output of `backend_type(m, set)`, where `set`
is the corresponding set for the row `rows[i]`.

`Interval` constraints require a call to `change_range_value!`.
"""
function change_linear_constraint_sense! end
@deprecate lqs_chgsense! change_linear_constraint_sense!

"""
    make_problem_type_integer(m)::Void

If an explicit call is needed to change the problem type integer (e.g., CPLEX).
"""
function make_problem_type_integer(m::LinQuadOptimizer)
    nothing  # default
end
@deprecate lqs_make_problem_type_integer make_problem_type_integer

"""
    make_problem_type_continuous(m)::Void

If an explicit call is needed to change the problem type continuous (e.g., CPLEX).
"""
function make_problem_type_continuous(m::LinQuadOptimizer)
    nothing  # default
end
@deprecate lqs_make_problem_type_continuous make_problem_type_continuous

"""
    add_sos_constraint!(m, cols::Vector{Int}, vals::Vector{Float64}, typ::Symbol)::Void

Add the SOS constraint to the model `m`. `typ` is either `:SOS1` or `:SOS2`.
"""
function add_sos_constraint! end
@deprecate lqs_addsos! add_sos_constraint!

"""
    delete_sos!(m, start_idx::Int, end_idx::Int)::Void

Delete the SOS constraints `start_idx`, `start_idx+1`, ..., `end_idx` from
the model `m`.
"""
function delete_sos! end
@deprecate lqs_delsos! delete_sos!

"""
    get_sos_constraint(m, idx::Int)::Tuple{Vector{Int}, Vector{Float64}, Symbol}

Get the SOS constraint `idx` from the model `m`. Returns the triplet
    `(cols, vals, typ)`.
"""
function get_sos_constraint end

"""
    get_number_quadratic_constraints(m)::Int

Get the number of quadratic constraints in the model `m`.
"""
function get_number_quadratic_constraints end
@deprecate lqs_getnumqcosntrs get_number_quadratic_constraints

"""
    add_quadratic_constraint!(m, cols::Vector{Int}, coefs::Vector{Float64}, rhs::Float64,
        sense, I::Vector{Int}, J::Vector{Int}, V::Vector{Float64})::Void

Add a quadratic constraint `a'x + 0.5 x' Q x`.
See `add_linear_constraints!` for information of linear component.
Arguments `(I,J,V)` given in triplet form for the Q matrix in `0.5 x' Q x`.
"""
function add_quadratic_constraint! end

"""
    set_quadratic_objective!(m, I::Vector{Int}, J::Vector{Int}, V::Vector{Float64})::Void

Set the quadratic component of the objective. Arguments given in triplet form
for the Q matrix in `0.5 x' Q x`.
"""
function set_quadratic_objective! end
@deprecate lqs_copyquad! set_quadratic_objective!

"""
    set_linear_objective!(m, cols::Vector{Int}, coefs::Vector{Float64})::Void

Set the linear component of the objective.
"""
function set_linear_objective! end
@deprecate lqs_chgobj! set_linear_objective!

"""
    change_objective_sense!(m, sense::Symbol)::Void

Change the optimization sense of the model `m` to `sense`. `sense` must be
`:min` or `:max`.
"""
function change_objective_sense! end
@deprecate lqs_chgobjsen! change_objective_sense!

"""
    get_linear_objective!(m, x::Vector{Float64})

Change the linear coefficients of the objective and store
in `x`.
"""
function get_linear_objective! end
@deprecate lqs_getobj get_linear_objective!

"""
    get_objectivesense(m)::MOI.OptimizationSense

Get the optimization sense of the model `m`.
"""
function get_objectivesense end
@deprecate lqs_getobjsen get_objectivesense

"""
    solve_mip_problem!(m)::Void

Solve a mixed-integer model `m`.
"""
function solve_mip_problem! end
@deprecate lqs_mipopt! solve_mip_problem!

"""
    solve_quadratic_problem!(m)::Void

Solve a model `m` with quadratic components.
"""
function solve_quadratic_problem! end
@deprecate lqs_qpopt! solve_quadratic_problem!

"""
    solve_linear_problem!(m)::Void

Solve a linear program `m`.
"""
function solve_linear_problem! end
@deprecate lqs_lpopt! solve_linear_problem!

"""
    get_variable_primal_solution!(m, x::Vector{Float64})

Get the primal solution for the variables in the model `m`, and
store in `x`. `x`must have one element for each variable.
"""
function get_variable_primal_solution! end
@deprecate lqs_getx! get_variable_primal_solution!

"""
    get_linear_primal_solution!(m, x::Vector{Float64})

Given a set of linear constraints `l <= a'x <= b` in the model `m`, get the
constraint primal `a'x` for each constraint, and store in `x`.
`x` must have one element for each linear constraint.
"""
function get_linear_primal_solution! end
@deprecate lqs_getax! get_linear_primal_solution!

"""
    get_quadratic_primal_solution!(m, x::Vector{Float64})

Given a set of quadratic constraints `l <= a'x + x'Qx <= b` in the model `m`,
get the constraint primal `a'x + x'Qx` for each constraint, and store in `x`.
`x` must have one element for each quadratic constraint.
"""
function get_quadratic_primal_solution! end
@deprecate lqs_getqcax! get_quadratic_primal_solution!

"""
    get_variable_dual_solution!(m, x::Vector{Float64})

Get the dual solution (reduced-costs) for the variables in the model `m`, and
store in `x`. `x`must have one element for each variable.
"""
function get_variable_dual_solution! end
@deprecate lqs_getdj! get_variable_dual_solution!

"""
    get_linear_dual_solution!(m, x::Vector{Float64})

Get the dual solution for the linear constraints in the model `m`, and
store in `x`. `x`must have one element for each linear constraint.
"""
function get_linear_dual_solution! end
@deprecate lqs_getpi! get_linear_dual_solution!

"""
    get_quadratic_dual_solution!(m, x::Vector{Float64})

Get the dual solution for the quadratic constraints in the model `m`, and
store in `x`. `x`must have one element for each quadratic constraint.
"""
function get_quadratic_dual_solution! end
@deprecate lqs_getqcpi! get_quadratic_dual_solution!

"""
    get_objective_value(m)

Get the objective value of the solved model `m`.
"""
function get_objective_value end
@deprecate lqs_getobjval get_objective_value

"""
    get_objective_bound(m)

Get the objective bound of the model `m`.
"""
function get_objective_bound end
@deprecate lqs_getbestobjval get_objective_bound

"""
    get_relative_mip_gap(m)

Get the relative MIP gap of the solved model `m`.
"""
function get_relative_mip_gap end
@deprecate lqs_getmiprelgap get_relative_mip_gap

"""
    get_iteration_count(m)

Get the number of simplex iterations performed during the most recent
optimization of the model `m`.
"""
function get_iteration_count end
@deprecate lqs_getitcnt get_iteration_count

"""
    get_barrier_iterations(m)

Get the number of barrier iterations performed during the most recent
optimization of the model `m`.
"""
function get_barrier_iterations end
@deprecate lqs_getbaritcnt get_barrier_iterations

"""
    get_node_count(m)

Get the number of branch-and-cut nodes expolored during the most recent
optimization of the model `m`.
"""
function get_node_count end
@deprecate lqs_getnodecnt get_node_count

"""
    get_farkas_dual!(m, x::Vector{Float64})

Get the farkas dual (certificate of primal infeasiblility) for the linear
constraints in the model `m`, and store in `x`. `x`must have one element for
each linear constraint.
"""
function get_farkas_dual! end
@deprecate lqs_dualfarkas! get_farkas_dual!

"""
    get_unbounded_ray!(m, x::Vector{Float64})

Get the unbounded ray (certificate of dual infeasiblility) for the linear
constraints in the model `m`, and store in `x`. `x`must have one element for
each variable.
"""
function get_unbounded_ray! end
@deprecate lqs_getray! get_unbounded_ray!

"""
    get_termination_status(m)

Get the termination status of the model `m`.
"""
function get_termination_status end
@deprecate lqs_terminationstatus get_termination_status

"""
    get_primal_status(m)

Get the primal status of the model `m`.
"""
function get_primal_status end
@deprecate lqs_primalstatus get_primal_status

"""
    get_dual_status(m)

Get the dual status of the model `m`.
"""
function get_dual_status end
@deprecate lqs_dualstatus get_dual_status

"""
    get_number_variables(m)::Int

Get the number of variables in the model `m`.
"""
function get_number_variables end
@deprecate lqs_getnumcols get_number_variables

"""
    add_variables!(m, n::Int)::Void

Add `n` new variables to the model `m`.
"""
function add_variables! end
@deprecate lqs_newcols! add_variables!

"""
    delete_variables!(m, start_col::Int, end_col::Int)::Void

Delete the columns `start_col`, `start_col+1`, ..., `end_col` from the model `m`.
"""
function delete_variables! end
@deprecate lqs_delcols! delete_variables!

"""
    add_mip_starts!(m, cols::Vector{Int}, x::Vector{Float64})::Void

Add the MIP start `x` for the variables in the columns `cols` of the model `m`.
"""
function add_mip_starts! end
@deprecate lqs_addmipstarts! add_mip_starts!


text:
module TransportLayerSecurity

# implementation packages must adhere to the following:

# implement types
# type SSLConfig
#     # hold configuration options
#     # cert/key values, etc.
# end

# constructor that takes cert and key files as inputs and returns a "default" configuration
# function SSLConfig#(cert_file, key_file)
#
# end

# type SSLContext
#     # a TLS-enabled socket for reading requests and sending responses
# end

# to "initialize" a TLS socket w/ a configuration
function setup!#(tls::SSLContext, config::SSLConfig)

end

# "upgrade" a regular TCPSocket to use the TLS SSLContext
function associate!#(tls::SSLContext, tcp::TCPSocket)

end

# perform the https handshake on a configured, initialized TLS SSLContext
function handshake!#(tls::SSLContext)

end

end # module


text:
f = open("big_input.txt","w+")
write(f,fill("1\n",4500000))
flush(f)
close(f)


text:
function calculate_ratios(pos::Array{T, 1})::Array{T, 2} where {T <: AbstractFloat}
    n::Int = length(pos)
    ratios::Array{T, 2} = zeros(Float64, n, n)
    @fastmath @inbounds for i in 1:n
        for j in i+1:n
            ratios[i,j] = pos[i] / pos[j]
        end
    end
    return ratios
end
function calculate_abs_diffs(pos::Array{T, 1})::Array{T, 2} where {T <: AbstractFloat}
    n::Int = length(pos)
    ratios::Array{T, 2} = zeros(Float64, n, n)
    @fastmath @inbounds for i in 1:n
        for j in i+1:n
            ratios[i, j] = pos[j] - pos[i]
        end
    end
    return ratios
end

function determine_calibration_constant_through_peak_ratios(fPositionX::Array{<:Real,1}, photon_lines::Array{<:Real,1}; α::Float64=0.01, rtol::Real = 5e-3)
    c::Float64 = 0
    photon_lines = Float64.(photon_lines)
    photon_lines = sort(photon_lines)
    fPositionX = sort(fPositionX)
    theo_ratios = calculate_ratios(photon_lines)
    exp_ratios = calculate_ratios(fPositionX)
    Δtheo = calculate_abs_diffs(photon_lines)
    n_peaks = length(fPositionX)
    n_lines = length(photon_lines)
    pcfs = Float64[] # pre calibration factors
    calibration_constant::Float64 = 0
    Δexp::Float64 = 0
    for i in 1:n_peaks
        for j in i + 1:n_peaks
            for it in 1:n_lines-1
                for jt in it + 1:n_lines
                    if 1 - α <= theo_ratios[it, jt] / exp_ratios[i, j] <= 1 + α
                        calibration_constant = (photon_lines[it] - photon_lines[jt]) / (fPositionX[i] - fPositionX[j])
                        Δexp = calibration_constant * (fPositionX[j] - fPositionX[i])
                        exp_peak_j = fPositionX[j] * calibration_constant
                        exp_peak_i = fPositionX[i] * calibration_constant
                        peak_i_rel_diff = abs(exp_peak_i - photon_lines[it]) / photon_lines[it]
                        peak_j_rel_diff = abs(exp_peak_j - photon_lines[jt]) / photon_lines[jt]
                        if peak_i_rel_diff < rtol && peak_j_rel_diff < rtol
                            push!(pcfs, calibration_constant )
                        end
                    end
                end
            end
        end
    end
    if isempty(pcfs)
        error("No peaks identified. Try different values for keyword parameters: `α`, `rtol`")
    end

    c_abs_min_diffs = zeros(Float64, length(pcfs))

    for (ic, c) in enumerate(pcfs)
        foundlines = c * fPositionX
        for (il, l) in enumerate(photon_lines)
            diffs = (foundlines .- l) ./ l
            abs_diffs = abs.(diffs)
            min_abs_diff, nearest_idx = findmin(abs_diffs)
            c_abs_min_diffs[ic] += min_abs_diff
        end
    end
    c = pcfs[findmin(c_abs_min_diffs)[2]]
    return c
end

function _filter_peaks_from_peakfinder(h::Histogram{<:Real, 1}, peakPositions::Vector{T}, σ::Real) where {T <: Real}
    pfits = []
    for p in peakPositions
        fitrange = (p - 5σ*step(h.edges[1]), p + 5σ*step(h.edges[1]))
        d = fit(NormalPeakUvD, subhist(h, fitrange))[1]
        push!(pfits, d)
    end
    d = map(i -> abs(pfits[i].UvNormal.μ - peakPositions[i]), eachindex(peakPositions));
    σs = map(f -> f.UvNormal.σ, pfits);
    hd = harmmean(d);
    hσ = harmmean(σs);
    accepted_peakPositions = T[]
    for (ipf, pf) in enumerate(pfits)
        if d[ipf] < 4 * hd && abs(pf.UvNormal.σ) < 2hσ
            push!(accepted_peakPositions, pf.UvNormal.μ)
        end
    end
    accepted_peakPositions
end

function determine_calibration_constant_through_peak_ratios(h::Histogram{<:Real, 1, E}, photon_lines::Array{<:Real, 1};
            min_n_peaks::Int = length(photon_lines), max_n_peaks::Int = 5 * length(photon_lines), threshold::Real = 10., α::Real = 0.01, σ::Real = 3.0, rtol::Real = 5e-3) where {E}
    h_deconv, peakPositions = peakfinder(h, threshold=threshold, σ = σ)
    fitted_peak_Positions = _filter_peaks_from_peakfinder(h, peakPositions, σ)
    while length(fitted_peak_Positions) > max_n_peaks
        threshold *= 1.1
        h_deconv, peakPositions = peakfinder(h, threshold=threshold, σ = σ)
        fitted_peak_Positions = _filter_peaks_from_peakfinder(h, peakPositions, σ)
    end
    while length(fitted_peak_Positions) < min_n_peaks
        threshold *= 0.75
        h_deconv, peakPositions = peakfinder(h, threshold=threshold, σ = σ)
        fitted_peak_Positions = _filter_peaks_from_peakfinder(h, peakPositions, σ)
    end
    c = determine_calibration_constant_through_peak_ratios(peakPositions, photon_lines, α = α, rtol = rtol)
    return c, h_deconv, peakPositions, threshold
end

function determine_calibration_constant_through_peak_fitting(h::Histogram{<:Real, 1, E}, photon_lines::Vector{<:Real}, c_pre::Real=1.0) where {E}
    c_pre = Float64(c_pre)
    peak_fits = []
    photon_lines = Float64.(photon_lines)
    processed_photon_lines = Float64[]
    for pl in photon_lines
        line::Float64 = pl
        fitrange = (line - 20, line + 20)
        fitrange = fitrange ./ c_pre

        first_bin, last_bin = StatsBase.binindex(h, fitrange[1]), StatsBase.binindex(h, fitrange[2])
        if last_bin > length(h.weights) || first_bin < 1
            @warn("$pl keV line not sufficiently contained in Histogram, skipping...");
            continue
        end
        d = fit(NormalPeakUvD, subhist(h, fitrange))[1]
        push!(peak_fits, d)
        push!(processed_photon_lines, pl)
    end
    photon_lines = intersect(photon_lines, processed_photon_lines)
    fitted_peak_positions = [ fr.UvNormal.μ for fr in peak_fits ]

    @fastmath function _linear_function_fixed_offset_at_zero(p::Vector{T}) where {T}
        return sqrt( sum(((fitted_peak_positions * p[1]) .- photon_lines).^2) )
    end

    opt_result = Optim.optimize(
        _linear_function_fixed_offset_at_zero,
        Float64[0], Float64[Inf], Float64[c_pre],
        Fminbox(BFGS());
        autodiff=:forward )

    c = opt_result.minimizer[1]

    return c, peak_fits, opt_result
end

"""
    calibrate_spectrum(h_uncal::Histogram, photon_lines::Array{Real, 1}; <keyword arguments>)

Returns the calibrated histogram, the deconvoluted spectrum, the found (uncalibrated) peak positions and the final `threshold` value.

# Keywords
- `σ::Real = 2.0`: The expected sigma of a peak in the spectrum. In units of bins.
- `threshold::Real = 10.0`: Threshold for being identified as a peak in the deconvoluted spectrum. A single bin is identified as a peak when its weight exceeds the `threshold` and the previous bin was not identified as an peak.
- `min_n_peaks::Int = 0`: If the number of found peaks is smaller than `min_n_peaks` the functions lowers the parameter `threshold` until enough peaks are found.
- `max_n_peaks::Int = 50`: Use only the first (strongest) `max_n_peaks` peaks for peak identification.
- `α::Real = 0.005`:  = 0.5%. Acceptance level in the comparison of the peak position ratios in the peak indentification step. When the difference between the ratio of two found peak positions and the ratio of two photon lines (`photon_lines`) is smaller than `α`, the found peaks are identified as the two photon lines.
- `rtol::Real = 5e-3`:  = 5e-3. Acceptance level for tolerance of the absolute difference between true and found line position.

Calibrate the spectrum `h_uncal`. This is done by:
1) finding peaks through devoncolution
2) identifying them through comparison of the ratios of their positions with the ratios of the known `lines`
3) fitting all identified peaks (with a gaussian plus first order polynomial) to get their position more precisely
4) performe a linear fit (offset forced to 0) of these positions vs the true positions (`lines`) to get the calibration constant
"""
function calibrate_spectrum(h_uncal::Histogram{<:Real, 1}, photon_lines::Vector{<:Real};
            σ::Real = 3.0, threshold::Real = 50.0, min_n_peaks::Int = length(photon_lines),
            max_n_peaks::Int = 4 * length(photon_lines), α::Real = 0.01, rtol::Real = 5e-3)
    c_precal, h_deconv, peakPositions, threshold = determine_calibration_constant_through_peak_ratios(h_uncal, photon_lines,
                                min_n_peaks = min_n_peaks, max_n_peaks = max_n_peaks, threshold=threshold, α=α, σ=σ, rtol=rtol)
    c = determine_calibration_constant_through_peak_fitting(h_uncal, photon_lines, c_precal)[1]
    h_cal = Histogram( h_uncal.edges[1] .* c, :left )
    h_cal.weights = h_uncal.weights
    return h_cal, h_deconv, peakPositions, threshold, c, c_precal
end


text:
"""
    BaselliModel()

A simple structure that stores a Baselli-model in form of transfer functions. It has a recipe defined to get some bode plots.

The default initialization is not proper, to get a working model use 'getModel(sim)', which gives one of the four configurations from the original paper.

# Examples
```julia
julia> model = getModel(1)
[ Info: Simulation 1: Respiration enters on S (LF are generated in the S, S loop)
BaselliModel...
```
"""
@with_kw mutable struct BaselliModel
    # Transfer Functions
    G_SI::TransferFunction = tf([0], 1)
    G_IS::TransferFunction = tf([0], 1)
    G_SS::TransferFunction = tf([0], 1)
    MS::TransferFunction = tf([0], 1)
    MI::TransferFunction = tf([0], 1)
    Mρ::TransferFunction = tf([0], 1)
    R_Sρ::TransferFunction = tf([0], 1)
    R_Iρ::TransferFunction = tf([0], 1)
    Gs::TransferFunction = tf([0], 1)
    Gt::TransferFunction = tf([0], 1)
    DS::TransferFunction = tf([0], 1)
    DI::TransferFunction = tf([0], 1)
    # Input GWN
    wS::Distribution = Normal(0,1)
    wI::Distribution = Normal(0,1)
    wρ::Distribution = Normal(0,1)

    # Data associsated with fitting, to calculate quality measures later on
    # Order
    n::Integer = 10
    # Data length (beats)
    N::Integer = 0
    # Witheness of residuals according to 5% Anderson TransferFunction
    pS_anderson::Real = -1
    pI_anderson::Real = -1
    # residuals for further tests
    res_S::Vector{<:Real} = Float64[]
    res_I::Vector{<:Real} = Float64[]
end

"""
    predict(model::BaselliModel, N::Int; burnIn::Int = 50)

Predicts 'N' values of the cardiovascular variables for the respective model. With 'burnIn' a certain number af values can be dropped in the beginning.

# Examples
```julia
 julia> S, I, ρ = predict(model, 100)
 (S = [...], I = [...], ρ = [...])
 ```
 """
function predict(model::BaselliModel, N::Int; burnIn::Int = 50)
    Random.seed!(1234)
    # Generate white noise input
    wρ = rand(model.wρ, N+burnIn)
    wS = rand(model.wS, N+burnIn)
    wI = rand(model.wI, N+burnIn)
    # time base for simulation
    T = 1:N+burnIn
    # respiration
    ρ = lsim(model.Mρ, wρ, T)[1][:]
    # w = [wS wI ρ]

    # a priori known
    rs = lsim(model.R_Sρ, ρ, T)[1][:]
    us = lsim(model.MS, wS, T)[1][:]
    cs = rs .+ us

    rt = lsim(model.R_Iρ, ρ, T)[1][:]
    ut = lsim(model.MI, wI, T)[1][:]
    ct = rt .+ ut

    # feedback dependend generate S and I
    Gs_cs = 1 / (1 - model.G_SS - (model.G_SI * model.G_IS))
    Gs_ct = model.G_SI / (1 - model.G_SS - (model.G_SI * model.G_IS))
    S = lsim([Gs_cs Gs_ct], [cs ct], T)[1][:]

    Gt_cs = model.G_IS / (1 - model.G_SS - (model.G_SI * model.G_IS))
    Gt_ct = 1 / (1 - ((model.G_IS * model.G_SI) / (1 - model.G_SS)))
    I = lsim([Gt_cs Gt_ct], [cs ct], T)[1][:]

    # return S, I and ρ
    return (S = S[burnIn+1:end], I = I[burnIn+1:end], ρ = ρ[burnIn+1:end])
end

# Simulations from original paper Baselli et al. 1988
"""
    getModel(sim::Int)

Gives a proper Baselli-model taken from the original paper, where 'sim' can take the value from 1:4, with the meaning:

* 1: Respiration enters on S (LF are generated in the S, S loop)
* 2: Respiration enters on I (LF are generated in the S, S loop)
* 3: The LF are generated by S regulating feedback in the S, S loop (respiration enters both on S and on I)
* 4: LF are generated externally from the closed-loop system (i.e., via a central command-respiration enters both on S and on I)

"""
function getModel(sim::Int)
    # Same for all simulations
    gSI = [0,-15]
    G_SI = impRes2tf(gSI)
    gIS = [range(0.005, stop = 0, length = 6);]
    G_IS = impRes2tf(gIS)
    gSS = [0, 0.5, 0, 0, -0.05, -0.1, -0.2, -0.1, -0.05]

    MS = tf([1], [1, 0, -0.25], 1)
    wS = Normal(0, sqrt(5))

    MI = tf([1], [1, -0.7], 1)
    wI = Normal(0, sqrt(0.00005))
    # AR(2) process
    wρ = Normal(0,1)
    # Mρ = zpk(Int[], [pol2imag(0.9, 0.44pi), pol2imag(0.9, -0.44pi)], 1, 1)
    Mρ = tf([1], [1, -0.337, 0.81], 1)

    gainFactor_S = 1/6
    gainFactor_I = 1/6

    # Unique properties
    if sim == 1
        @info "Simulation 1: Respiration enters on S (LF are generated in the S, S loop)"
        gainFactor_I = 0.1/6
    elseif sim == 2
        @info "Simulation 2: Respiration enters on I (LF are generated in the S, S loop)"
        gainFactor_S = 0.1/6
    elseif sim == 3
        @info "Simulation 3: The LF are generated by S regulating feedback in the S, S loop (respiration enters both on S and on I)"
    elseif sim == 4
        @info "Simulation 4: LF are generated externally from the closed-loop system (i.e., via a central command-respiration enters both on S and on I)"
        gSS = [0, 0.5]
        MS = zpk( Int[], [pol2imag(0.95, 0.14π), pol2imag(0.95, -0.14π)], 1, 1)
    else
        throw(DomainError(sim, "Only simulation 1:4 exist. Choose one of them."))
    end
    impulse_S = gainFactor_S .* [3:-1:1;]
    R_Sρ = impRes2tf([0; impulse_S])
    impulse_I = gainFactor_I .* [3:-1:1;]
    R_Iρ = impRes2tf([0; impulse_I])
    G_SS = impRes2tf(gSS)

    return BaselliModel(G_SI = G_SI, G_IS = G_IS, G_SS = G_SS, MS = MS, wS = wS, MI = MI, wI = wI, Mρ = Mρ, wρ = wρ, R_Sρ = R_Sρ, R_Iρ = R_Iρ)
end

# Helper to get the tf from the impulse response
function impRes2tf(h::Vector{<:Real})
    den = zeros(length(h))
    den[1] = 1
    return tf(h, den, 1) # discrete
    # return tf(h, den) # continous
end
#turns poles that are given in polar coordinates into a complex pole
function pol2imag(M::Real, ϕ::Real)
    return Complex(M * cos(ϕ), M * sin(ϕ))
end

@recipe function f(model::BaselliModel; plotphase = false)
    w = range(0, stop = π, length = 256)

    layout := (3,4)
    # layout := @layout [_ G_IS R_Iρ MI
     #G_SS G_SI R_Sρ MS
     #wS wI wρ Mρ]

    x = 1

    @series begin
        mag, phase, w = bode(model.G_IS, w)
        title := "G_IS"
        subplot := x+1
        (w./2π), plotphase ? vec(phase) : vec(mag)
    end

    @series begin
        mag, phase, w = bode(model.R_Iρ, w)
        title := "R_Iρ"
        subplot := x+2
        (w./2π), plotphase ? vec(phase) : vec(mag)
    end

    @series begin
        mag, phase, w = bode(model.MI, w)
        title := "MI"
        subplot := x+3
        (w./2π), plotphase ? vec(phase) : vec(mag)
    end

    @series begin
        mag, phase, w = bode(model.G_SS, w)
        title := "G_SS"
        subplot := x+4
        (w./2π), plotphase ? vec(phase) : vec(mag)
    end

    @series begin
        mag, phase, w = bode(model.G_SI, w)
        title := "G_SI"
        subplot := x+5
        (w./2π), plotphase ? vec(phase) : round.(vec(mag), digits = 2)
    end

    @series begin
        mag, phase, w = bode(model.R_Sρ, w)
        title := "R_Sρ"
        subplot := x+6
        (w./2π), plotphase ? vec(phase) : vec(mag)
    end

    @series begin
        mag, phase, w = bode(model.MS, w)
        title := "MS"
        subplot := x+7
        (w./2π), plotphase ? vec(phase) : vec(mag)
    end

    @series begin
        title := "wS"
        subplot := x+8
        model.wS
    end
    @series begin
        title := "wI"
        subplot := x+9
        model.wI
    end
    @series begin
        title := "wρ"
        subplot := x+10
        model.wρ
    end
    @series begin
        mag, phase, w = bode(model.Mρ, w)
        title := "Mρ"
        subplot := x+11
        (w./2π), plotphase ? vec(phase) : vec(mag)
    end
end