text:
export Chebyshev

"""
`Chebyshev()` is the space spanned by the Chebyshev polynomials
```
    T_0(x),T_1(x),T_2(x),…
```
where `T_k(x) = cos(k*acos(x))`.  This is the default space
as there exists a fast transform and general smooth functions on `[-1,1]`
can be easily resolved.
"""
struct Chebyshev{D<:Domain,R} <: PolynomialSpace{D,R}
    domain::D
    function Chebyshev{D,R}(d) where {D,R}
        isempty(d) && throw(ArgumentError("Domain cannot be empty"))
        new(d)
    end
    Chebyshev{D,R}() where {D,R} = new(convert(D, ChebyshevInterval()))
end

Chebyshev(d::Domain) = Chebyshev{typeof(d),real(prectype(d))}(d)
Chebyshev() = Chebyshev(ChebyshevInterval())
Chebyshev(d) = Chebyshev(Domain(d))

Space(d::SegmentDomain) = Chebyshev(d)
function Space(d::AbstractInterval)
    a,b = endpoints(d)
    if isinf(norm(a)) && isinf(norm(b))
        Chebyshev(Line(d))
    elseif isinf(norm(a)) || isinf(norm(b))
        Chebyshev(Ray(d))
    else
        Chebyshev(d)
    end
end

setdomain(S::Chebyshev, d::Domain) = Chebyshev(d)

ones(::Type{T}, S::Chebyshev) where {T<:Number} = Fun(S,fill(one(T),1))
ones(S::Chebyshev) = Fun(S,fill(1.0,1))

function Base.first(f::Fun{<:Chebyshev})
    n = ncoefficients(f)
    n == 0 && return zero(cfstype(f))
    n == 1 && return f.coefficients[1]
    foldr(-,coefficients(f))
end

Base.last(f::Fun{<:Chebyshev}) = reduce(+,coefficients(f))

spacescompatible(a::Chebyshev,b::Chebyshev) = domainscompatible(a,b)
hasfasttransform(::Chebyshev) = true

## Transform

transform(::Chebyshev,vals::AbstractVector,plan) = plan*vals
itransform(::Chebyshev,cfs::AbstractVector,plan) = plan*cfs
plan_transform(::Chebyshev,vals::AbstractVector) = plan_chebyshevtransform(vals)
plan_itransform(::Chebyshev,cfs::AbstractVector) = plan_ichebyshevtransform(cfs)

## Evaluation

clenshaw(sp::Chebyshev, c::AbstractVector, x::AbstractArray) =
    clenshaw(c,x,ClenshawPlan(promote_type(eltype(c),eltype(x)),sp,length(c),length(x)))

function clenshaw(::Chebyshev,c::AbstractVector,x)
    N,T = length(c),promote_type(eltype(c),typeof(x))
    if N == 0
        return zero(x)
    elseif N == 1 # avoid issues with NaN x
        return first(c)*one(x)
    end

    x = 2x
    bk1,bk2 = zero(T),zero(T)
    @inbounds for k = N:-1:2
        bk2, bk1 = bk1, muladd(x,bk1,c[k]-bk2)
    end

    muladd(x/2,bk1,c[1]-bk2)
end

clenshaw(c::AbstractVector,x::AbstractVector,plan::ClenshawPlan{S,V}) where {S<:Chebyshev,V}=
    clenshaw(c,collect(x),plan)

#TODO: This modifies x, which is not threadsafe
function clenshaw(c::AbstractVector,x::Vector,plan::ClenshawPlan{S,V}) where {S<:Chebyshev,V}
    N,n = length(c),length(x)
    if isempty(c)
        return zeros(V,n)
    end

    bk=plan.bk
    bk1=plan.bk1
    bk2=plan.bk2

    @inbounds for i = 1:n
        x[i] = 2x[i]
        bk1[i] = zero(V)
        bk2[i] = zero(V)
    end

    @inbounds for k = N:-1:2
        ck = c[k]
        for i = 1:n
            bk[i] = muladd(x[i],bk1[i],ck-bk2[i])
        end
        bk2, bk1, bk = bk1, bk, bk2
    end

    ck = c[1]
    @inbounds for i = 1:n
        x[i] = x[i]/2
        bk[i] = muladd(x[i],bk1[i],ck-bk2[i])
    end

    bk
end

function clenshaw(c::AbstractMatrix{T},x::T,plan::ClenshawPlan{S,T}) where {S<:Chebyshev,T<:Number}
    bk=plan.bk
    bk1=plan.bk1
    bk2=plan.bk2

    m,n=size(c) # m is # of coefficients, n is # of funs

    @inbounds for i = 1:n
        bk1[i] = zero(T)
        bk2[i] = zero(T)
    end
    x = 2x

    @inbounds for k=m:-1:2
        for j=1:n
            ck = c[k,j]
            bk[j] = muladd(x,bk1[j],ck - bk2[j])
        end
        bk2, bk1, bk = bk1, bk, bk2
    end

    x = x/2
    @inbounds for i = 1:n
        ce = c[1,i]
        bk[i] = muladd(x,bk1[i],ce - bk2[i])
    end

    bk
end

function clenshaw(c::AbstractMatrix{T},x::AbstractVector{T},plan::ClenshawPlan{S,T}) where {S<:Chebyshev,T<:Number}
    bk=plan.bk
    bk1=plan.bk1
    bk2=plan.bk2

    m,n=size(c) # m is # of coefficients, n is # of funs

    @inbounds for i = 1:n
        x[i] = 2x[i]
        bk1[i] = zero(T)
        bk2[i] = zero(T)
    end

    @inbounds for k=m:-1:2
        for j=1:n
            ck = c[k,j]
            bk[j] = muladd(x[j],bk1[j],ck - bk2[j])
        end
        bk2, bk1, bk = bk1, bk, bk2
    end

    @inbounds for i = 1:n
        x[i] = x[i]/2
        ce = c[1,i]
        bk[i] = muladd(x[i],bk1[i],ce - bk2[i])
    end

    bk
end

# overwrite x

function clenshaw!(c::AbstractVector,x::AbstractVector,plan::ClenshawPlan{S,V}) where {S<:Chebyshev,V}
    N,n = length(c),length(x)

    if isempty(c)
        fill!(x,0)
        return x
    end

    bk=plan.bk
    bk1=plan.bk1
    bk2=plan.bk2

    @inbounds for i = 1:n
        x[i] = 2x[i]
        bk1[i] = zero(V)
        bk2[i] = zero(V)
    end

    @inbounds for k = N:-1:2
        ck = c[k]
        for i = 1:n
            bk[i] = muladd(x[i],bk1[i],ck - bk2[i])
        end
        bk2, bk1, bk = bk1, bk, bk2
    end

    ce = c[1]
    @inbounds for i = 1:n
        x[i] = x[i]/2
        x[i] = muladd(x[i],bk1[i],ce-bk2[i])
    end

    x
end

## Calculus

# diff T -> U, then convert U -> T
integrate(f::Fun{Chebyshev{D,R}}) where {D<:IntervalOrSegment,R} =
    Fun(f.space,fromcanonicalD(f,0)*ultraint!(ultraconversion(f.coefficients)))
differentiate(f::Fun{Chebyshev{D,R}}) where {D<:IntervalOrSegment,R} =
    Fun(f.space,1/fromcanonicalD(f,0)*ultraiconversion(ultradiff(f.coefficients)))

## Multivariate

# determine correct parameter to have at least
# N point
_padua_length(N) = Int(cld(-3+sqrt(1+8N),2))

function squarepoints(::Type{T}, N) where T
    pts = paduapoints(T, _padua_length(N))
    n = size(pts,1)
    ret = Array{Vec{2,T}}(undef, n)
    @inbounds for k=1:n
        ret[k] = Vec{2,T}(pts[k,1],pts[k,2])
    end
    ret
end

points(S::TensorSpace{<:Tuple{<:Chebyshev{<:ChebyshevInterval},<:Chebyshev{<:ChebyshevInterval}}}, N) =
    squarepoints(real(prectype(S)), N)

function points(S::TensorSpace{<:Tuple{<:Chebyshev,<:Chebyshev}},N)
    T = real(prectype(S))
    pts = squarepoints(T, N)
    pts .= fromcanonical.(Ref(domain(S)), pts)
    pts
end

plan_transform(S::TensorSpace{<:Tuple{<:Chebyshev,<:Chebyshev}},v::AbstractVector) =
    plan_paduatransform!(v,Val{false})

transform(S::TensorSpace{<:Tuple{<:Chebyshev,<:Chebyshev}},v::AbstractVector,
        plan=plan_transform(S,v)) = plan*copy(v)

plan_itransform(S::TensorSpace{<:Tuple{<:Chebyshev,<:Chebyshev}},v::AbstractVector) =
     plan_ipaduatransform!(eltype(v),sum(1:nblocks(Fun(S,v))),Val{false})

itransform(S::TensorSpace{<:Tuple{<:Chebyshev,<:Chebyshev}},v::AbstractVector,
         plan=plan_itransform(S,v)) = plan*pad(v,sum(1:nblocks(Fun(S,v))))

#TODO: adaptive
for op in (:(Base.sin),:(Base.cos))
    @eval ($op)(f::ProductFun{<:Chebyshev,<:Chebyshev}) =
        ProductFun(chebyshevtransform($op.(values(f))),space(f))
end

reverseorientation(f::Fun{<:Chebyshev}) =
    Fun(Chebyshev(reverseorientation(domain(f))),alternatesign!(copy(f.coefficients)))

include("ChebyshevOperators.jl")


text:
module FocalPointAnalysis

using Fermi
using Fermi.Output
using Fermi.Integrals: IntegralHelper
using Fermi.HartreeFock: RHF
using Fermi.Geometry: Molecule
using Fermi.CoupledCluster: RCCSD, RCCSDpT
using PrettyTables

abstract type AbstractFPA end

struct FPA{T} <: AbstractFPA
    basis::Array{String,1}
    methods::Array{String,1}
    data::Array{T,2}
    extrapolated_indexes::Array{Tuple{Int64,Int64},1}
end

struct ecFPA{T} <: AbstractFPA
    basis::Array{String,1}
    methods::Array{String,1}
    data::Array{T,2}
    extrapolated_indexes::Array{Tuple{Int64,Int64},1}
end

function print_fpa(fpa::T; incremental=false) where T <: AbstractFPA

    data = copy(fpa.data)
    if incremental
        for k in size(fpa.data,2):-1:2
           data[:,k] -= data[:,k-1]
        end
    end

    format = ft_printf("%15.10f")
    put_braket = (v,i,j) -> (i,j) in fpa.extrapolated_indexes ? "["*v*"]" : v

    out = pretty_table(String, data, fpa.methods, alignment=:c,
    formatters = (format, put_braket), row_names=fpa.basis,
    row_name_column_title="Basis Set")

    @output "\n{}\n" out
end

function Base.:*(N::T, F::FPA) where T <: Number
    return FPA(F.basis, F.methods, N.*F.data, F.extrapolated_indexes)
end

function Base.:+(A::FPA, B::FPA)
    @assert A.methods == B.methods
    @assert A.basis == B.basis
    @assert A.extrapolated_indexes == B.extrapolated_indexes
    return FPA(A.basis, B.methods, A.data+B.data, A.extrapolated_indexes)
end

function Base.:-(A::FPA, B::FPA)
    @assert A.methods == B.methods
    @assert A.basis == B.basis
    @assert A.extrapolated_indexes == B.extrapolated_indexes
    return FPA(A.basis, B.methods, A.data-B.data, A.extrapolated_indexes)
end

include("Extrapolation.jl")
include("FPA.jl")
include("ecFPA.jl")

end #module

ext: jl


text:
using DataFrames
using DataFramesMeta
using Dates
using DelimitedFiles
using CSV
using Query
using Statistics
using StatsBase

const FILE_ACCESS_OVER = "w+";
const FILE_ACCESS_APPEND = "a+";

DF_Model = CSV.read(".//outputs//postprocess//postproc_model_costs_01-01-2019_12-25-2019.csv", header=true, DataFrame);
DF_Inputs = CSV.read(".//outputs//postprocess//postproc_gens_cost_emiss_01-01-2019_12-25-2019.csv", header=true, DataFrame);

Start_date = "01/01/2019 07:00:00"
Final_date = "12/25/2019 06:00:00"

#Select a subset of values that fall into the months that the users specifies.
ini_date = DateTime(Start_date, "mm/dd/yyyy HH:MM:SS");
end_date = DateTime(Final_date, "mm/dd/yyyy HH:MM:SS");

rename!(DF_Model,:DATE => :BEGIN_DATE);

df_inputs_drange = @from i in DF_Inputs begin
        @where DateTime(i.BEGIN_DATE) >= ini_date && DateTime(i.BEGIN_DATE) <= end_date
        @select {i.BEGIN_DATE, i.PORTFOLIO_NAME, i.UNIT_NAME, i.UNIT_ID, i.UNIT_TYPE, i.TYPE_MAIN_FUEL,
                i.MW, i.GenCost, i.EmissCO2_ton, i.EmissNOx_kg, i.EmissSO2_kg}
        @collect DataFrame
   end;

df_model_drange = @from i in DF_Model begin
   @where DateTime(i.BEGIN_DATE) >= ini_date && DateTime(i.BEGIN_DATE) <= end_date
   @select {i.BEGIN_DATE, i.UNIT_ID, i.Output, i.GenCost, i.EmissCO2_ton,
            i.EmissNOx_kg, i.EmissSO2_kg}
   @collect DataFrame
end;

rename!(df_model_drange,:Output => :MW_Model, :GenCost => :GenCost_Model,
         :EmissCO2_ton => :EmissCO2_ton_Model, :EmissNOx_kg => :EmissNOx_kg_Model,
         :EmissSO2_kg => :EmissSO2_kg_Model);

DF_Merge = outerjoin(df_inputs_drange, df_model_drange,
            on = [:BEGIN_DATE => :BEGIN_DATE, :UNIT_ID => :UNIT_ID],
            validate=(true, true), makeunique=true);

##Estimate Error Metrics

# Mean Absolute Error (MAE) of power generation (MW)
insertcols!(DF_Merge, size(DF_Merge,2)+1, :Abs_Diff_MW => abs.(DF_Merge.MW_Model .- DF_Merge.MW));

insertcols!(DF_Merge, size(DF_Merge,2)+1, :MW_NoMiss => DF_Merge.MW);
insertcols!(DF_Merge, size(DF_Merge,2)+1, :MW_Model_NoMiss => DF_Merge.MW_Model);

DF_Merge.MW_NoMiss = collect(Missings.replace(DF_Merge.MW_NoMiss, 0))
DF_Merge.MW_Model_NoMiss = collect(Missings.replace(DF_Merge.MW_Model_NoMiss, 0))

DF_Merge.MW_NoMiss = convert.(Float64, DF_Merge.MW_NoMiss);
DF_Merge.MW_Model_NoMiss = convert.(Float64, DF_Merge.MW_Model_NoMiss);
insertcols!(DF_Merge, size(DF_Merge,2)+1, :Abs_Diff_MW_NoMiss => abs.(DF_Merge.MW_NoMiss .- DF_Merge.MW_Model_NoMiss));

MAE_MW = mean(skipmissing(DF_Merge.Abs_Diff_MW));

#MAE_MW_NoMiss = mean(DF_Merge.Abs_Diff_MW_NoMiss);
MAE_MW_NoMiss = meanad(DF_Merge.MW_NoMiss, DF_Merge.MW_Model_NoMiss);

# Mean Percentage Error (MAPE) of power generation (MW)
insertcols!(DF_Merge, size(DF_Merge,2)+1, :Percent_Diff_MW => DF_Merge.Abs_Diff_MW ./ DF_Merge.MW);
insertcols!(DF_Merge, size(DF_Merge,2)+1, :Percent_Diff_MW_NoMiss => DF_Merge.Abs_Diff_MW_NoMiss ./ DF_Merge.MW_NoMiss);

#TODO: Check an alterantive to MAPE
MAPE_MW = (sum(skipmissing(DF_Merge.Percent_Diff_MW)) / size(DF_Merge,1)) * 100;
MAPE_MW_NoMiss = (sum(skipmissing(DF_Merge.Percent_Diff_MW_NoMiss)) / size(DF_Merge,1)) * 100;

#Root Mean Square Error (RMSE) of power generation
insertcols!(DF_Merge, size(DF_Merge,2)+1, :Sqr_Diff_MW => DF_Merge.Abs_Diff_MW.^2);
insertcols!(DF_Merge, size(DF_Merge,2)+1, :Sqr_Diff_MW_NoMiss => DF_Merge.Abs_Diff_MW_NoMiss.^2);

MSE_MW = mean(skipmissing(DF_Merge.Sqr_Diff_MW));
MSE_MW_NoMiss = msd(DF_Merge.MW_NoMiss, DF_Merge.MW_Model_NoMiss);

RMSE_MW = sqrt(mean(skipmissing(DF_Merge.Sqr_Diff_MW)));
RMSE_MW_NoMiss = rmsd(DF_Merge.MW_NoMiss, DF_Merge.MW_Model_NoMiss);
#RMSE_MW_NoMiss = sqrt(mean(DF_Merge.Sqr_Diff_MW_NoMiss));

# Mean Absolute Error (MAE) of generation cost ($)
insertcols!(DF_Merge, size(DF_Merge,2)+1, :Abs_Diff_GenCost => abs.(DF_Merge.GenCost_Model .- DF_Merge.GenCost));

insertcols!(DF_Merge, size(DF_Merge,2)+1, :GenCost_NoMiss => DF_Merge.GenCost);
insertcols!(DF_Merge, size(DF_Merge,2)+1, :GenCost_Model_NoMiss => DF_Merge.GenCost_Model);

DF_Merge.GenCost_NoMiss = collect(Missings.replace(DF_Merge.GenCost_NoMiss, 0))
DF_Merge.GenCost_Model_NoMiss = collect(Missings.replace(DF_Merge.GenCost_Model_NoMiss, 0))

DF_Merge.GenCost_NoMiss = convert.(Float64, DF_Merge.GenCost_NoMiss);
DF_Merge.GenCost_Model_NoMiss = convert.(Float64, DF_Merge.GenCost_Model_NoMiss);

insertcols!(DF_Merge, size(DF_Merge,2)+1, :Abs_Diff_GenCost_NoMiss => abs.(DF_Merge.GenCost_NoMiss .- DF_Merge.GenCost_Model_NoMiss));

MAE_GenCost = mean(skipmissing(DF_Merge.Abs_Diff_GenCost));
MAE_GenCost_NoMiss = mean(DF_Merge.Abs_Diff_GenCost_NoMiss);

# Mean Percentage Error (MAPE) of power generation (GenCost)
insertcols!(DF_Merge, size(DF_Merge,2)+1, :Percent_Diff_GenCost => DF_Merge.Abs_Diff_GenCost ./  DF_Merge.GenCost);
insertcols!(DF_Merge, size(DF_Merge,2)+1, :Percent_Diff_GenCost_NoMiss => (DF_Merge.Abs_Diff_GenCost_NoMiss ./  DF_Merge.GenCost_NoMiss));

MAPE_GenCost = (sum(skipmissing(DF_Merge.Percent_Diff_GenCost)) / size(DF_Merge,1)) * 100;
MAPE_GenCost_NoMiss = (sum(skipmissing(DF_Merge.Percent_Diff_GenCost_NoMiss)) / size(DF_Merge,1)) * 100;

#Root Mean Square Error (RMSE) of generation cost
insertcols!(DF_Merge, size(DF_Merge,2)+1, :Sqr_Diff_GenCost => DF_Merge.Abs_Diff_GenCost.^2);
insertcols!(DF_Merge, size(DF_Merge,2)+1, :Sqr_Diff_GenCost_NoMiss => DF_Merge.Abs_Diff_GenCost_NoMiss.^2);

MSE_GenCost = mean(skipmissing(DF_Merge.Sqr_Diff_GenCost));
MSE_GenCost_NoMiss = msd(DF_Merge.GenCost_NoMiss, DF_Merge.GenCost_Model_NoMiss);

RMSE_GenCost = sqrt(mean(skipmissing(DF_Merge.Sqr_Diff_GenCost)));
RMSE_GenCost_NoMiss = rmsd(DF_Merge.GenCost_NoMiss, DF_Merge.GenCost_Model_NoMiss);

DF_Merge_Byname = groupby(DF_Merge, [:UNIT_NAME, :UNIT_TYPE, :TYPE_MAIN_FUEL]);
DF_Merge_ByID = groupby(DF_Merge, [:UNIT_ID, :UNIT_NAME, :UNIT_TYPE, :TYPE_MAIN_FUEL]);
DF_Merge_Byfuel = groupby(DF_Merge, :TYPE_MAIN_FUEL);

DF_Summary_Byname = @combine(DF_Merge_Byname,
    MW = sum(:MW),
    GenCost = sum(:GenCost),
    EmissCO2_ton = sum(:EmissCO2_ton),
    EmissSO2_kg = sum(:EmissSO2_kg),
    EmissNOx_kg = sum(:EmissNOx_kg),
    MW_Model = sum(:MW_Model),
    GenCost_Model = sum(:GenCost_Model),
    EmissCO2_ton_Model = sum(:EmissCO2_ton_Model),
    EmissSO2_kg_Model = sum(:EmissSO2_kg_Model),
    EmissNOx_kg_Model = sum(:EmissNOx_kg_Model),
    #Sum_Abs_Diff_MW = sum(:Abs_Diff_MW),
    #Percent_Diff_over_MW = mean(:Abs_Diff_MW),
    #Sum_Abs_Diff_GenCost = sum(:Abs_Diff_GenCost),
    #Percent_Diff_over_GenCost = sum(:Abs_Diff_GenCost) / sum(:GenCost),
    #MAE_MW = mean(:Abs_Diff_MW_NoMiss),
    MAE_MW = meanad(:MW_NoMiss, :MW_Model_NoMiss),
    MSE_MW = msd(:MW_NoMiss, :MW_Model_NoMiss),
    #RMSE_MW = sqrt(mean(:Sqr_Diff_MW_NoMiss)),
    RMSE_MW = rmsd(:MW_NoMiss, :MW_Model_NoMiss),
    RMSE_MW_Normal = rmsd(:MW_NoMiss, :MW_Model_NoMiss; normalize=true),
    MAE_GenCost = meanad(:GenCost_NoMiss, :GenCost_Model_NoMiss),
    MSE_GenCost = msd(:GenCost_NoMiss, :GenCost_Model_NoMiss),
    RMSE_GenCost = rmsd(:GenCost_NoMiss, :GenCost_Model_NoMiss),
    RMSE_GenCost_Normal = rmsd(:GenCost_NoMiss, :GenCost_Model_NoMiss; normalize=true),
    N = size(:MW_NoMiss, 1));

DF_Summary_ByID = @combine(DF_Merge_ByID,
    MW = sum(:MW),
    GenCost = sum(:GenCost),
    EmissCO2_ton = sum(:EmissCO2_ton),
    EmissSO2_kg = sum(:EmissSO2_kg),
    EmissNOx_kg = sum(:EmissNOx_kg),
    MW_Model = sum(:MW_Model),
    GenCost_Model = sum(:GenCost_Model),
    EmissCO2_ton_Model = sum(:EmissCO2_ton_Model),
    EmissSO2_kg_Model = sum(:EmissSO2_kg_Model),
    EmissNOx_kg_Model = sum(:EmissNOx_kg_Model),
    MAE_MW = meanad(:MW_NoMiss, :MW_Model_NoMiss),
    MSE_MW = msd(:MW_NoMiss, :MW_Model_NoMiss),
    RMSE_MW = rmsd(:MW_NoMiss, :MW_Model_NoMiss),
    RMSE_MW_Normal = rmsd(:MW_NoMiss, :MW_Model_NoMiss; normalize=true),
    MAE_GenCost = meanad(:GenCost_NoMiss, :GenCost_Model_NoMiss),
    MSE_GenCost = msd(:GenCost_NoMiss, :GenCost_Model_NoMiss),
    RMSE_GenCost = rmsd(:GenCost_NoMiss, :GenCost_Model_NoMiss),
    RMSE_GenCost_Normal = rmsd(:GenCost_NoMiss, :GenCost_Model_NoMiss; normalize=true),
    N = size(:MW_NoMiss, 1));

DF_Summary_Byfuel = @combine(DF_Merge_Byfuel,
    MW = sum(:MW),
    GenCost = sum(:GenCost),
    EmissCO2_ton = sum(:EmissCO2_ton),
    EmissSO2_kg = sum(:EmissSO2_kg),
    EmissNOx_kg = sum(:EmissNOx_kg),
    MW_Model = sum(:MW_Model),
    GenCost_Model = sum(:GenCost_Model),
    EmissCO2_ton_Model = sum(:EmissCO2_ton_Model),
    EmissSO2_kg_Model = sum(:EmissSO2_kg_Model),
    EmissNOx_kg_Model = sum(:EmissNOx_kg_Model),
    MAE_MW = meanad(:MW_NoMiss, :MW_Model_NoMiss),
    MSE_MW = msd(:MW_NoMiss, :MW_Model_NoMiss),
    RMSE_MW = rmsd(:MW_NoMiss, :MW_Model_NoMiss),
    RMSE_MW_Normal = rmsd(:MW_NoMiss, :MW_Model_NoMiss; normalize=true),
    MAE_GenCost = meanad(:GenCost_NoMiss, :GenCost_Model_NoMiss),
    MSE_GenCost = msd(:GenCost_NoMiss, :GenCost_Model_NoMiss),
    RMSE_GenCost = rmsd(:GenCost_NoMiss, :GenCost_Model_NoMiss),
    RMSE_GenCost_Normal = rmsd(:GenCost_NoMiss, :GenCost_Model_NoMiss; normalize=true),
    N = size(:MW_NoMiss, 1) );

#Some of the operations in the @combine function add a total row. Delete this.
delete!(DF_Summary_Byname, size(DF_Summary_Byname,1))
delete!(DF_Summary_ByID, size(DF_Summary_ByID,1))
delete!(DF_Summary_Byfuel, size(DF_Summary_Byfuel,1))

#Store results
Start_date = replace(Start_date, "/" => "-");
Final_date = replace(Final_date, "/" => "-");

Start_date = Start_date[1:findfirst(isequal(' '), Start_date)-1]
Final_date = Final_date[1:findfirst(isequal(' '), Final_date)-1]

CSV.write(".//outputs//postprocess//gen_cost_emiss_inputANDmodel_$(Start_date)_$(Final_date).csv", DF_Merge)

CSV.write(".//outputs//postprocess//perform_error_byName_$(Start_date)_$(Final_date).csv", DF_Summary_Byname)
CSV.write(".//outputs//postprocess//perform_error_byID_$(Start_date)_$(Final_date).csv", DF_Summary_ByID)
CSV.write(".//outputs//postprocess//perform_error_byFuel_$(Start_date)_$(Final_date).csv", DF_Summary_Byfuel)

# Save performance error indicators for whole data set
performance_header    = ["Outcome", "Mean Absolute Error (MAE)", "Mean Square Error (MSE)", "Root Mean Square Error (RMSE)"]
open(".//outputs//postprocess//perform_error_TOTAL_$(Start_date)_$(Final_date).csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(performance_header), ',')
    writedlm(io, hcat("Power Output (MW)", MAE_MW, MSE_MW, RMSE_MW), ',')
    writedlm(io, hcat("Power Output (MW) [Missing replaced]", MAE_MW_NoMiss, MSE_MW_NoMiss, RMSE_MW_NoMiss), ',')
    writedlm(io, hcat("Generation Cost", MAE_GenCost, MSE_GenCost, RMSE_GenCost), ',')
    writedlm(io, hcat("Generation Cost [Missing replaced]", MAE_GenCost_NoMiss, MSE_GenCost_NoMiss, RMSE_GenCost_NoMiss), ',')
end; # closes file

#include(".//postprocessing//get_perform_err_metrics.jl")

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
using ReservoirComputing, Random

const res_size = 20
const ts = 0.:0.1:50.0
const data = sin.(ts)
const train_len = 400
const input_data = reduce(hcat, data[1:train_len-1])
const target_data = reduce(hcat, data[2:train_len])
const predict_len = 100
const test_data = reduce(hcat, data[train_len+1:train_len+predict_len])
const training_method = StandardRidge(10e-6)

states_types = [StandardStates, ExtendedStates, PaddedStates, PaddedExtendedStates]

for t in states_types
    Random.seed!(77)
    esn = ESN(input_data;
        reservoir=RandSparseReservoir(res_size, 1.2, 0.1))
    output_layer = train(esn, target_data, training_method)
    output = esn(Generative(predict_len), output_layer)
    @test maximum(abs.(test_data .- output)) ./ maximum(abs.(test_data)) < 0.1
end

max_stars_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_issues_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_forks_repo_name: MartinuzziFrancesco/EchoStateNetwork


text:
# to_expr
function to_expr(f::SAF)
    f = MOIU.canonical(f)
    if isempty(f.terms)
        return f.constant
    else
        linear_term_exprs = map(f.terms) do term
            :($(term.coefficient) * x[$(term.variable_index.value)])
        end
        expr = :(+($(linear_term_exprs...)))
        if !iszero(f.constant)
            push!(expr.args, f.constant)
        end
        return expr
    end
end

function to_expr(f::SQF)
    f = MOIU.canonical(f)
    linear_term_exprs = map(f.affine_terms) do term
        i = term.variable_index.value
        :($(term.coefficient) * x[$i])
    end
    quadratic_term_exprs = map(f.quadratic_terms) do term
        i = term.variable_index_1.value
        j = term.variable_index_2.value
        if i == j
            :($(term.coefficient / 2) * x[$i] * x[$j])
        else
            :($(term.coefficient) * x[$i] * x[$j])
        end
    end
    expr = :(+($(linear_term_exprs...), $(quadratic_term_exprs...)))
    if !iszero(f.constant)
        push!(expr.args, f.constant)
    end
    return expr
end

# check_variable_indices
function check_variable_indices(model::Optimizer, index::VI)
    @assert 1 <= index.value <= length(model.inner.variable_info)
end

function check_variable_indices(model::Optimizer, f::SAF)
    for term in f.terms
        check_variable_indices(model, term.variable_index)
    end
end

function check_variable_indices(model::Optimizer, f::SQF)
    for term in f.affine_terms
        check_variable_indices(model, term.variable_index)
    end
    for term in f.quadratic_terms
        check_variable_indices(model, term.variable_index_1)
        check_variable_indices(model, term.variable_index_2)
    end
end

function find_variable_info(model::Optimizer, vi::VI)
    check_variable_indices(model, vi)
    model.inner.variable_info[vi.value]
end


text:
using Test
using DistributedFactorGraphs
using Dates

#TODO implement with tests on
# TestFunctorInferenceType1
# TestCCW1

## Generated compare functions
# VariableNodeData
vnd1 = VariableNodeData(TestVariableType1())
vnd2 = deepcopy(vnd1)
vnd3 = VariableNodeData(TestVariableType2())

@test vnd1 == vnd2
push!(vnd1.val, [1.0;])
push!(vnd2.val, [1.0;])
@test vnd1 == vnd2
vnd2.val[1] = [0.1;]
@test !(vnd1 == vnd2)
@test !(vnd1 == vnd3)

# MeanMaxPPE
ppe1 = MeanMaxPPE(:default, [1.], [2.], [3.])
ppe2 = deepcopy(ppe1)
ppe3 = MeanMaxPPE(:default, [2.], [3.], [4.])

@test ppe1 == ppe2
ppe2.max[1] = 0.1
@test !(ppe1 == ppe2)
@test !(ppe1 == ppe3)

# DFGVariable
v1 = DFGVariable(:x1, TestVariableType1())
v2 = deepcopy(v1)
v3 = DFGVariable(:x2, TestVariableType2())

@test v1 == v2
v2.solvable = 0
@test !(v1 == v2)
@test !(v1 == v3)
@test !(DFGVariable(:x1, TestVariableType1()) == DFGVariable(:x1, TestVariableType2()))

# GenericFunctionNodeData
gfnd1 = GenericFunctionNodeData(true, true, [1,2], TestFunctorInferenceType1())
gfnd2 = deepcopy(gfnd1)
gfnd3 = GenericFunctionNodeData(true, true, [1,2], TestFunctorInferenceType2())

@test gfnd1 == gfnd2
@test !(gfnd1 == gfnd3)

# DFGFactor
f1 = DFGFactor(:f1, [:a, :b], gfnd1)
f2 = deepcopy(f1)
f3 = DFGFactor(:f1, [:b, :a], gfnd1)

@test f1 == f2
@test !(f1 == f3)

## Compare functions
vnd1 = VariableNodeData(TestVariableType1())
vnd2 = deepcopy(vnd1)
vnd3 = VariableNodeData(TestVariableType2())
@test compare(vnd1, vnd2)
@test !compare(vnd1, vnd3)

@test compare(vnd1, vnd2)
push!(vnd1.val, [1.0;])
push!(vnd2.val, [1.0;])
@test compare(vnd1, vnd2)
vnd2.val[1][1] = 0.1
@test !compare(vnd1, vnd2)
@test !compare(vnd1, vnd3)

gfnd1 = GenericFunctionNodeData(true, true, [1,2], TestFunctorInferenceType1())
gfnd2 = deepcopy(gfnd1)
gfnd3 = GenericFunctionNodeData(true, true, [1,2], PackedTestFunctorInferenceType1())

@test compare(gfnd1, gfnd2)
@test_broken !(compare(gfnd1, gfnd3))


text:
module OpSel

using Printf
using DataFrames
using JuMP
using Ipopt
using DelimitedFiles
using SparseArrays
using LinearAlgebra

# function testEqual(sp_csv, N, theta)
function testEqual(Z, N = 50)
    sp_csv = loadDataBySize(Z)
    (x_result, info_result) = steepestAscent(sp_csv, N)
    return (x_result, info_result, sp_csv)
end

thetaTable = spzeros(300000, 30000)
thetaTable[ 50,   200] = 3.3429999999999994e-02/2
thetaTable[ 50,  1050] = 6.2748061171875091e-02/2
thetaTable[ 50,  2045] = 7.1083439999999956e-02/2
thetaTable[ 50,  5050] = 1.0806060937499347e-01/2
thetaTable[ 50,  5255] = 2.4439005000000000e-02/2
thetaTable[ 50, 10100] = 7.0094325234373603e-02/2
thetaTable[ 50, 15100] = 6.7641368203122818e-02/2
thetaTable[ 50, 15222] = 3.8808090000000316e-02/2
thetaTable[100,   200] = 2.5817499999999993e-02/2
thetaTable[100,  1050] = 5.3925029121093773e-02/2
thetaTable[100,  2045] = 6.2823779999999996e-02/2
thetaTable[100,  5050] = 9.9463449550774755e-02/2
thetaTable[100,  5255] = 1.5510273749999984e-02/2
thetaTable[100, 10100] = 6.0998935781248584e-02/2
thetaTable[100, 15100] = 5.8678192968747774e-02/2
thetaTable[100, 15222] = 3.0043923750000281e-02/2

function steepestAscent(sp_csv, N=50)
    Z = size(sp_csv,1)
    try
        theta = thetaTable[N, Z]
    catch error1
        error(@sprintf("theta for N = %d, Z = %d is not available in opsel.thetaTable\n", N, Z))
    end

    info_result = Dict()
    info_result["build_time"] =  @elapsed begin
        Z = size(sp_csv,1)
        Ainv = makeAinv(sp_csv)
        g = sp_csv[:,4]
        l = zeros(Z)
        u = min.(sp_csv[:,5], 1/N)
        # A = inv(full(Ainv))
        A = inv(Matrix(Ainv))
        A = (A+A')/2
        if theta < 0
            println("*** Negative theta is given. Setting up new theta")
            theta = newTheta(Z, A, u, N)
        end

        println("Z = $Z, N = $N, theta = $theta")

        # JuMP
        model1 = Model(with_optimizer(Ipopt.Optimizer, acceptable_tol=1e-5))

        @variable(model1, y[1:Z])
        @constraint(model1, Ainv*y .<=u)
        @constraint(model1, l.<=Ainv*y)
        @constraint(model1, sum(Ainv*y) == 1.0)
        (I1,J1,V1) = findnz(Ainv)
        @constraint(model1, sum(V1[k]*y[I1[k]]*y[J1[k]] for k=1:length(V1)) <=2*theta)

        @objective(model1, Max, (Ainv*g)'*y)

    end # build _time

    info_result["solver_time"] = @elapsed optimize!(model1)

    y_value = value.(y)
    x_value = Ainv*y_value
    x_value = slimData(x_value,1.0e-5)

    @printf("=== SOCP Result Summary ===\n")
    @printf("JuMP status = %s\n", termination_status(model1))
    info_result["JuMPstatus"] = @sprintf("%s", termination_status(model1))
    info_result["socp_time"] = info_result["build_time"] + info_result["solver_time"]
    @printf("sum(x_socp) = %lf, gx = %lf, xAx = %lf, 2theta = %lf\n",
            sum(x_value), dot(g, x_value),
            dot(y_value, Ainv*y_value), 2*theta)

    info_result["steep_time"] = @elapsed begin
        # compute Lambda0 by Meuwissen's lagrange multiplier
        lambda0 = computeLambda0(Z, Ainv, g, theta);
        lambda0 = 2*lambda0; # mulitply with 2 is import for Z = 1050
        (x_sa, iter_sa, beforecobj_socp, aftercobj_socp) = localSwap(x_value, sp_csv, A, Ainv, N, theta, lambda0)
        info_result["steep_iter"] = iter_sa
    end

    @printf("sum(x_sa) = %lf, gx = %lf, xAx = %lf, 2theta = %lf\n",
            sum(x_sa), dot(g,x_sa),
            dot(x_sa, Ainv\x_sa), 2*theta)

    gx_value = dot(g, x_value)
    gx_sa    = dot(g, x_sa)

    @printf("SOCP = %f, STEEP = %f, gap = %f%%\n",
        gx_value, gx_sa, abs(gx_value - gx_sa)/max(gx_value, gx_sa)*100)

    info_result["total_time"] = info_result["socp_time"] + info_result["steep_time"]

    @printf("time(s): build = %.3f, solver = %.3f, steep = %.3f, total = %.3f\n",
            info_result["build_time"],
            info_result["solver_time"],
            info_result["steep_time"],
            info_result["total_time"])
    return (x_sa, info_result)
end

function testUnequal(Z)
    sp_csv = loadDataBySize(Z)
    if Z == 200
        println("Z = 200 is not solvable for N_s = 14 and N = 2800")
        return (0, 0, 0)
    end
    (x_result, info_result) = compactSOCP(sp_csv, 14, 2800)
    return (x_result, info_result, sp_csv)
end

function compactSOCP(sp_csv, N_s = 14, N = 2800)

    # N_s is the effective number
    # In the data, we fixed N_s = 14

    info_result = Dict()
    info_result["build_time"] =  @elapsed begin
        Z = size(sp_csv,1)
        theta = 0.5/N_s
        println("Z = $Z, N_s = $N_s, N = $N, theta = $theta")

        p = sp_csv[:,2]
        q = sp_csv[:,3]
        g = sp_csv[:,4]
        l = zeros(Z)
        u = sp_csv[:,5]/N
        inbreeding = sp_csv[:,6]
        Ainv = makeAinv(sp_csv)

        # JuMP
        model1 = Model(with_optimizer(Ipopt.Optimizer, acceptable_tol=1e-5))

        @variable(model1, y[1:Z])
        @constraint(model1, Ainv*y .<=u)
        @constraint(model1, l.<=Ainv*y)
        @constraint(model1, sum(Ainv*y) == 1.0)
        (I1,J1,V1) = findnz(Ainv)
        @constraint(model1, sum(V1[k]*y[I1[k]]*y[J1[k]] for k=1:length(V1)) <=2*theta)
        @objective(model1, Max, dot(Ainv*g, y))
        # @objective(model1, Max, dot(Ainv*g, y))
    end # build _time

    info_result["solver_time"] = @elapsed optimize!(model1)
    info_result["total_time"] = info_result["build_time"] + info_result["solver_time"]

    y_value = value.(y)
    x_value = slimData(Ainv*y_value,0)

    @printf("=== Result Summary ===\n")
    @printf("JuMP status = %s\n", termination_status(model1))
    info_result["JuMPstatus"] = @sprintf("%s", termination_status(model1))

    println("Z = $Z, N_s = $N_s, N = $N")
    @printf("gx = %lf, xAx = %lf, 2theta = %lf\n", dot(g, x_value),
            dot(y_value, Ainv*y_value), 2*theta)
    @printf("time(s): build = %.3f, solver = %.3f, total = %.3f\n",
            info_result["build_time"],
            info_result["solver_time"],
            info_result["total_time"])
    return (x_value, info_result)
end

function loadDataBySize(Z)
    filename = @sprintf("%s/../data/sorted%06d.csv",@__DIR__, Z)
    sp_csv = loadFile(filename)
    return sp_csv
end

function loadFile(filename)
    sp_csv = readdlm(filename,',')
    sp_csv = convert(Array{Float64,2}, sp_csv[2:end,:])
    # ignore the column names in the first line
    return sp_csv
end

function slimData(A, epsilon)
    if length(size(A)) == 1
        return A.*(abs.(A).>=epsilon)
    end
    (m,n) = size(A)
    (I1,J1,V1) = findnz(A)
    V2 = V1.*(abs(V1).>=epsilon)
    A2 = sparse(I1,J1,V2,m,n)
    return A2
end

# function [Ainv] = makeAinv2( sp_csv )
#   compute Ainverse by Henderson's formula

function makeAinv(sp_csv)
    Z = size(sp_csv,1)
    inbreeding = sp_csv[:,6]

    eleNo = 0
    for i = 1:Z
        p = convert(Int, sp_csv[i,2])
        q = convert(Int, sp_csv[i,3])
        if p != 0 && q != 0
            eleNo = eleNo + 9
            continue
        end
        if (p != 0 && q == 0) || (p == 0 && q != 0 )
            eleNo = eleNo + 4
            continue
        end
        if p == 0 && q == 0
            eleNo = eleNo + 1
        end
    end

    I = zeros(Int, eleNo)
    J = zeros(Int, eleNo)
    V = zeros(eleNo)

    eleNo = 0
    for i = 1:Z
        p = convert(Int, sp_csv[i,2])
        q = convert(Int, sp_csv[i,3])
        if p != 0 && q != 0
            b = 4.0/( (1-inbreeding[p]) + (1-inbreeding[q]) )
            I[eleNo + 1] = i; J[eleNo + 1] = i; V[eleNo + 1] =  1.0*b;
            I[eleNo + 2] = p; J[eleNo + 2] = i; V[eleNo + 2] = -0.5*b;
            I[eleNo + 3] = q; J[eleNo + 3] = i; V[eleNo + 3] = -0.5*b;
            I[eleNo + 4] = i; J[eleNo + 4] = p; V[eleNo + 4] = -0.5*b;
            I[eleNo + 5] = i; J[eleNo + 5] = q; V[eleNo + 5] = -0.5*b;
            I[eleNo + 6] = p; J[eleNo + 6] = p; V[eleNo + 6] = +0.25*b;
            I[eleNo + 7] = p; J[eleNo + 7] = q; V[eleNo + 7] = +0.25*b;
            I[eleNo + 8] = q; J[eleNo + 8] = p; V[eleNo + 8] = +0.25*b;
            I[eleNo + 9] = q; J[eleNo + 9] = q; V[eleNo + 9] = +0.25*b;
            eleNo = eleNo + 9
            continue
        end
        if (p != 0 && q == 0) || (p == 0 && q != 0 )
            p = max(p,q)
            b = 4.0/( 1.0*(1-inbreeding[p]) + 2.0*(1-0) )
            I[eleNo + 1] = i; J[eleNo + 1] = i; V[eleNo + 1] =  1.0*b;
            I[eleNo + 2] = p; J[eleNo + 2] = i; V[eleNo + 2] = -0.5*b;
            I[eleNo + 3] = i; J[eleNo + 3] = p; V[eleNo + 3] = -0.5*b;
            I[eleNo + 4] = p; J[eleNo + 4] = p; V[eleNo + 4] = +0.25*b;
            eleNo = eleNo + 4
            continue
        end
        if p == 0 && q == 0
            b = 4.0/( 2.0*(1-0) + 2.0*(1-0) )
            I[eleNo + 1] = i; J[eleNo + 1] = i; V[eleNo + 1] =  1.0*b;
            eleNo = eleNo + 1
        end
    end

    Ainv = sparse(I,J,V,Z,Z)
    # Ainv = Ainv + Ainv' - diag(diag(Ainv));
    return Ainv
end # end of makeAinv

function makeB(sp_csv, f)

    Z = size(sp_csv,1)
    p = convert(Array{Int,1}, sp_csv[:,2])
    q = convert(Array{Int,1}, sp_csv[:,3])

    indI = zeros(Int,3*Z)
    indJ = zeros(Int,3*Z)
    indV = zeros(3*Z)
    ind0 = 1

    for i=1:Z
        if  p[i] == 0 && q[i] == 0
            b = 4.0/(2*(1-0) + 2*(1-0))
            indI[ind0] = i
            indJ[ind0] = i
            indV[ind0] = 1*sqrt(b)
            ind0 = ind0+1
        elseif q[i] == 0
            b = 4.0/(1.0*(1-f[p[i]]) + 2.0*(1-0))
            indI[ind0] = i
            indJ[ind0] = i
            indV[ind0] = 1*sqrt(b)
            ind0 = ind0+1
            indI[ind0] = p[i]
            indJ[ind0] = i
            indV[ind0] = -0.5*sqrt(b)
            ind0 = ind0+1
        else
            b = 4.0/((1-f[p[i]])+(1-f[q[i]]));
            indI[ind0] = i
            indJ[ind0] = i
            indV[ind0] = 1.0*sqrt(b)
            ind0 = ind0+1
            indI[ind0] = p[i]
            indJ[ind0] = i
            indV[ind0] = -0.5*sqrt(b)
            ind0 = ind0+1
            indI[ind0] = q[i]
            indJ[ind0] = i
            indV[ind0] = -0.5*sqrt(b)
            ind0 = ind0+1
        end
    end

    B = sparse(indJ[1:ind0-1],indI[1:ind0-1],indV[1:ind0-1],Z,Z)
    return B
end # end of function makeB

function localSwap(x_value, sp_csv, A, Ainv, N, theta_relax, lambda0)

    Z = size(sp_csv,1)
    l = zeros(Z,1)
    g = sp_csv[:,4]
    u = min.(sp_csv[:,5],1/N)
    u_orig = u

    sort_index = sortperm(x_value, rev = true)
    sort_x = x_value[sort_index]
    xf = zeros(Z,1)
    F = convert(Array{Int,1},sort_index[1:N])
    xf[F] .= 1/N
    V = setdiff(findall(x -> x>=1/N, u_orig), F);

    take_max = 0; # usually zero
    # take_max = -inf; # usually zero
    c_obj = dot(g,xf) - lambda0*max(take_max, dot(xf,A*xf)/2 - theta_relax)
    before_cobj = c_obj
    @printf("======= Before swap c_obj = %f\n", c_obj);

    iter2 = 0
    xAx_prev = dot(xf, A*xf)/2
    diagA = diag(A)
    while true
        changed = 0
        iter2 = iter2+1
        A_FF = sum(A[F,F])
        gF = sum(g[F])
        best_f = 0
        best_v = 0
        for index_f = 1:length(F)
            f = F[index_f]
            xAx_vec =  (A_FF + diagA[f] - 2*sum(A[F,f])) .+
              (diagA[V] -2*A[V,f]+2*sum(A[V,F],dims=2))
            c_obj3 = 1/N*((gF - g[f]) .+ g[V]) -
                lambda0 * max.(take_max, 1/N^2*xAx_vec/2 .- theta_relax)
            (c_obj2, max_v) = findmax(c_obj3)
            if c_obj2 > c_obj
                changed = 1
                best_f = f
                best_v = V[max_v]
                c_obj = c_obj2
            end
        end
        if changed == 0
            @printf("No swap c_obj = %f, gx = %f, xAx = %f\n",
                    c_obj, dot(g,xf), dot(xf,A*xf));
            break
        end
        xf[best_f] = 0
        xf[best_v] = 1/N
        F = union(setdiff(F,best_f),best_v)
        V = union(setdiff(V,best_v),best_f)

        @printf("iter2 = %d, c_obj = %f, gx = %f, xAx = %f, best_f = %d, best_v = %d\n",
                iter2, c_obj, dot(g,xf), dot(xf, Ainv\xf), best_f, best_v)
    end
    after_cobj = c_obj
    @printf("Finished ... iter2 = %d, before_cobj = %f, after_cobj = %f\n",
            iter2, before_cobj, after_cobj)

    return  (xf, iter2, before_cobj, after_cobj)

end # end of localSswap

function computeLambda0(Z, Ainv, g, theta)
#  compute lambda0 based on Meuwissen's lagrange multiplier

    Q = ones(Z,1)
    r = ones(1,1)
    theta_min = dot(r,(Q'*Ainv*Q)\r)/2
    if theta_min > theta
        @printf("theta (%f) is too small. It must be >= %f\n", theta, theta_min)
        error("theta error")
    end
    Ainv_g = Ainv*g
    B = inv(Q'*Ainv*Q)
    l_bunshi = dot(g,Ainv_g) - dot(Q'*Ainv_g, B*(Q'*Ainv_g))
    l_bunbo = 8*theta - 4*dot(r, B*r)
    lambda0 = 2*sqrt(l_bunshi/l_bunbo)*2 # mulitply with 2 is import for Z = 1050

    return lambda0
end # end of computeLambda0

end # end of module


text:
using .Plots

"""
    pdata = pplot(x, w, y, yhat, a, t, pdata; kwargs...)
    pdata = pplot(pf, y, pdata; kwargs...)

To be called inside a particle filter, plots either particle density (`density=true`) or individual particles (`density=false`) \n
Will plot all the real states in `xindices` as well as the expected vs real measurements of `yindices`.
# Arguments:
- `x`: `Vector{Vector}(N)`. The states for each particle where `N` number of Particles
- `w`: `Vector(N)`. weight of each particle
- `y`: `Vector{Vector}(T)`. All true outputs. `T` is total number of time steps (will only use index `t`)
- `yhat`: `Vector{Vector}(N)` The expected output per particle. `N` number of Particles
- `a`, `Vector(N)`, reorderng of particles (e.g. `1:N`)
- `t`, Current time step
- `xreal`: `Vector{Vector}(T)`. All true states if available. `T` is total number of time steps (will only use index `t`)
- `xprev`: Same as `x`, but for previous time step, only used when `!density` to show states origins
- `pdata`: Persistant data for plotting. Set to `nothing` in first call and pdata on remaining \n
- `density = true` To only plot the particle trajectories, set (`leftonly=false`)\n
- `leftonly = true`: only plot the left column\n
- `xindices = 1:n_states`\n
- `yindices = 1:n_measurements`\n
Returns: `pdata`
"""
function pplot(pf, u, y, args...; kwargs...)
    s = state(pf)
    t = s.t[]
    pplot(s.x, s.we, u, y, LowLevelParticleFilters.measurement(pf).(s.x, Ref(u[t]), t), s.j, t, args...; xprev=s.xprev, kwargs...)
end

function pplot(x, w, u, y, yhat, a, t, pdata; xreal=nothing, xprev=nothing,  density = true, leftonly = true, xindices = 1:length(x[1]), yindices = 1:length(y[1]), lowpass=0.9)

    (t == 1 || t % 35 == 0) &&  @printf "Time     Surviving    Effective nbr of particles\n--------------------------------------------------------------\n"

    T = length(y)
    N = length(x)
    D = length(x[1])
    x = reduce(hcat, x)
    y = reduce(hcat, y)
    yhat = reduce(hcat, yhat)

    cols = leftonly ? 1 : 2
    grd = (r,c) -> (r-1)*cols+c
    @printf("t: %5d %7.3f %9.1f\n", t, (N-length(setdiff(Set(1:N),Set(a))))/N, effective_particles(w))
    plotvals = [x;yhat]
    realVals = xreal === nothing ? [fill(Inf,D,T);y] : [reduce(hcat,xreal);y] # Inf turns plotting off
    if !density
        plotvalsOld = xprev === nothing ? [Inf*x;yhat] : [reduce(hcat, xprev);yhat]
    end

    plotindices = [xindices; size(x,1) .+ yindices]
    if pdata === nothing
        pdata = plot(layout=(length(plotindices),cols)), zeros(length(plotindices),2)
    end
    p, minmax = pdata
    # dataMin = minimum(plotvals[plotindices,:])
    # dataMax = maximum(plotvals[plotindices,:])
    # minmax = [min.(minmax[:,1], dataMin)*lowpass .+ (1-lowpass)*dataMin max.(minmax[:,2], dataMax)*lowpass .+ (1-lowpass)*dataMax]

    for (i, pind) in enumerate(plotindices)
        #Plot the heatmap on the left plot
        if density
            heatboxplot!(plotvals[pind,:], t, w, subplot=grd(i,1), reuse=false) #ylims=minmax[i,:],
            if !leftonly
                densityplot!( plotvals[pind,:], w , c=:blue, subplot=grd(i,2), reuse=true, legend=false) #ylims = tuple(minmax[i,:]...)
            end
        else
            #Plot the line on the left plot
            plot!(repeat([t-0.5 t+0.5], N)', [plotvalsOld[pind,:] plotvalsOld[pind,:]]', legend=false, subplot=grd(i,1), l=(:black, 0.1))
            # plot!(repeat([t-1 t-0.5], N)', [plotvalsOld[pind,a] plotvals[pind,:]]', legend=false, subplot=grd(i,1), l=(:black, 0.1))
        end
        #Plot Real State Here
        scatter!( [t+0.5], [realVals[pind,t]], subplot=grd(i,1), legend=false, m=(:cyan,))
    end
    # gui(p)
    (p, minmax)
end

function commandplot(f)
    res = f(nothing)
    display(res[1])
    while true
        print("Waiting for command. q to Quit, s NN to skip NN steps:\n")
        line = readline()
        skip = 1
        if line[1] == 'q'
            return
        elseif occursin("s", line)
            if length(line) <= 2
                skip = 1
            else
                ss = split(strip(line,'\n'))
                skip = parse(Int,ss[2])
            end
        end
        foreach(1:skip) do i
            res = f(res)
        end
        display(res[1])
    end
end

"""
    commandplot(pf, u, y; kwargs...)

Produce a helpful plot. For customization options (`kwargs...`), see `?pplot`.
After each time step, a command from the user is requested.
- q: quit
- s n: step `n` steps
"""
function commandplot(pf, u, y; kwargs...)
    # pdata = nothing
    reset!(pf)
    pfp = pf isa AuxiliaryParticleFilter ? pf.pf : pf
    commandplot() do pdata
        pdata = pplot(pfp, u, y, pdata; kwargs...)
        t = index(pf)
        LowLevelParticleFilters.update!(pf,u[t],y[t])
        pdata
    end
end

"""
    debugplot(pf, u, y; runall=false, kwargs...)

Produce a helpful plot. For customization options (`kwargs...`), see `?pplot`.
- ` runall=false:` if true, runs all time steps befor displaying (faster), if false, displays the plot after each time step.

The generated plot becomes quite heavy. Initially, try limiting your input to 100 time steps to verify that it doesn't crash.
"""
function debugplot(pf, u, y; runall=false, kwargs...)
    pdata = nothing
    reset!(pf)
    pfp = pf isa AuxiliaryParticleFilter ? pf.pf : pf
    for i = 1:length(y)
        pdata = pplot(pfp, u, y, pdata; kwargs...)
        t = index(pf)
        LowLevelParticleFilters.update!(pf,u[t],y[t])
        runall || display(pdata[1])
    end
    display(pdata[1])
end


text:
# api.jl

module API

include("./spec.jl")

import ..AppliAR: Domain
using .Domain

using Dates

import AppliSales: Order # Order is not exported but is refered to in the next function
import AppliGeneralLedger: create_journal_entry

export create, conv2entry

# start invoice numbering
n = 1000

create(order::Order, invoice_id::String)::UnpaidInvoice = begin
    meta = MetaInvoice(order.id, order.training.id)
    header_invoice = Header(
		    invoice_id, order.org.name, order.org.address, order.org.zip, order.org.city, order.org.country, order.order_ref, order.contact_name, order.contact_email)
    body_invoice = OpentrainingItem(order.training.name, order.training.date, order.training.price, order.students)
	return UnpaidInvoice(invoice_id, meta, header_invoice, body_invoice)
end

create(invoice::UnpaidInvoice, bs::BankStatement)::PaidInvoice = begin
	id_inv = id(invoice)
	meta_inv = meta(invoice)
	header_inv = header(invoice)
	body_inv = body(invoice)
	s = bs
	return PaidInvoice(id_inv, meta_inv, header_inv, body_inv, s)
end

conv2entry(inv::UnpaidInvoice, from::Int, to::Int) = begin
    id = string(Date(now())) * "-" * string(global n += 1)
    customer_id = name(header(inv))
	inv_nbr = invoice_nbr(header(inv))
	b = body(inv)
    debit = price_per_student(b) * length(students(b))
    credit = 0.0
    vat = debit * vat_perc(b)
    descr = name_training(b)
    return create_journal_entry(id, customer_id, inv_nbr, from, to, debit, credit, vat, descr)
end

conv2entry(inv::PaidInvoice, from::Int, to::Int) = begin
    id = string(Date(now())) * "-" * string(global n += 1)
    customer_id = name(header(inv))
    inv_nbr = invoice_nbr(header(inv))
	b = body(inv)
    debit = amount(stm(inv))
    credit = 0.0
    vat = 0.0
    descr = name_training(b)
    return create_journal_entry(id, customer_id, inv_nbr, from, to, debit, credit, vat, descr)
end

end # module


text:
module TestConstructors

using Base.Test
using DataArrays
using DataFrames
using Dates

using TimeData

println("Running constructor tests")

## init test values
vals = rand(30, 4)
dats = Date[Date(2013, 7, ii) for ii=1:30]
tims = DateTime[DateTime(2010,01,ii,00,00,00)
                                 for ii=1:30]
idxs = [1:30]
## nams = ["A", "B", "C", "D"]
nams = [:A, :B, :C, :D]
valsDf = convert(DataFrame, vals)
rename!(valsDf, names(valsDf), nams)

allTypes = (:Timedata, :Timenum, :Timematr)

########################
## inner constructors ##
########################

for t in allTypes
    eval(quote
        td = $t(valsDf, dats)
        td = $t(valsDf, tims)
        td = $t(valsDf, idxs)
    end)
end

#######################################
## shortcuts from arrays without NAs ##
#######################################

for t in allTypes
    eval(quote
        ## from numeric array only
        td = $(t)(vals)

        ## from arrays without names
        td = $(t)(vals, dats)
        td = $(t)(vals, tims)
        td = $(t)(vals, idxs)

        ## from arrays without idx
        td = $(t)(vals, nams)

        ## from three arrays with names and idx
        td = $(t)(vals, nams, dats)
        td = $(t)(vals, nams, tims)
        td = $(t)(vals, nams, idxs)
    end)
end

###################
## NAs in values ##
###################

for t in allTypes
    eval(quote
        ## from DataFrame without idx
        td = $(t)(valsDf)

        ## from DataFrame with idx as array
        td = $(t)(valsDf, dats)
        td = $(t)(valsDf, tims)
        td = $(t)(valsDf, idxs)

        ## testing index values as DataArray
        td = $(t)(valsDf, convert(DataArray, idxs))
    end)
end

######################################
## throwing errors for wrong inputs ##
######################################

## idx initialized with wrong type
vals = DataFrame(a = [2, 3, 4])
invalidIdx = @data([1., NA, 3])
for t in allTypes
    eval(quote
        @test_throws Exception $(t)(vals, invalidIdx)
    end)
end

## idx and vals sizes not matching
idxVals = [Date(2013, 7, ii) for ii=1:30]
valsArr = rand(20, 4)
for t in allTypes
    eval(quote
        @test_throws Exception $(t)(valsArr, idxVals)
    end)
end

######################################
## check single column constructors ##
######################################

for t in allTypes
    eval(quote
        td = $t(rand(4))
    end)
end

end


text:
using Manopt, ManifoldsBase, Test

@testset "Gradient Plan" begin
    io = IOBuffer()
    M = ManifoldsBase.DefaultManifold(2)
    x = [4.0, 2.0]
    o = GradientDescentOptions(
        x; stopping_criterion=StopAfterIteration(20), stepsize=ConstantStepsize(1.0)
    )
    o.gradient = [1.0, 0.0]
    f(M, y) = distance(M, y, x) .^ 2
    gradf(M, y) = -2 * log(M, y, x)
    p = GradientProblem(M, f, gradf)
    @test get_initial_stepsize(p, o) == 1.0
    @test get_stepsize(p, o, 1) == 1.0
    @test get_last_stepsize(p, o, 1) == 1.0
    # Check Fallbacks of Problen
    @test get_cost(p, o.x) == 0.0
    @test get_gradient(p, o.x) == zero_vector(M, x)
    @test_throws MethodError get_proximal_map(p, 1.0, o.x, 1)
    @test_throws MethodError get_subgradient(p, o.x)
    # Additional Specific Debugs
    a1 = DebugGradient(false, io)
    a1(p, o, 1)
    @test String(take!(io)) == "gradF(x):[1.0, 0.0]"
    a1a = DebugGradient("s:", io)
    a1a(p, o, 1)
    @test String(take!(io)) == "s:[1.0, 0.0]"
    a2 = DebugGradientNorm(false, io)
    a2(p, o, 1)
    @test String(take!(io)) == "|gradF(x)|:1.0"
    a2a = DebugGradientNorm("s:", io)
    a2a(p, o, 1)
    @test String(take!(io)) == "s:1.0"
    a3 = DebugStepsize(false, io)
    a3(p, o, 1)
    @test String(take!(io)) == "s:1.0"
    a3a = DebugStepsize("S:", io)
    a3a(p, o, 1)
    @test String(take!(io)) == "S:1.0"
    # Additional Specific Records
    b1 = RecordGradient(o.gradient)
    b1(p, o, 1)
    @test b1.recorded_values == [o.gradient]
    b2 = RecordGradientNorm()
    b2(p, o, 1)
    @test b2.recorded_values == [1.0]
    b3 = RecordStepsize()
    b3(p, o, 1)
    b3(p, o, 2)
    b3(p, o, 3)
    @test b3.recorded_values == [1.0, 1.0, 1.0]
end


text:
"""
`calibForce(POSES, F, m0=0.3; offset=true)`
If result is bad, check if you send data in correct form;
POSES is always the position of the tool frame from the robot FK
cell array with 4x4 transformation matrices
F is a Nx3 vector of forces (accepts Nx6 matrix with torques also)
usage `Rf*force[i,1:3]' + forceoffs = POSES[1:3,1:3,i]'*[0, 0, mf*-9.82]`
Bagge
"""
function calibForce(POSES,F,m0=0.3; offset=true, useCVX = false)

    N = size(POSES,3)

    forceoffs = zeros(3)
    g = 9.82
    mg = m0*g # Initial guess is m0, the method accepts 5 orders of magnitude error at least

    I = eye(3)
    B = Array(Float64,3N)
    B2 = Array(Float64,3N)
    if offset # Include offset
        A = Array(Float64,3N,12)
        for i = 1:N
            RA = POSES[1:3,1:3,i]'
            At = [F[i,1]*I   F[i,2]*I    F[i,3]*I   I]
            b =  -RA[1:3,3]
            A[3(i-1)+1:3i,:] = At
            B[3(i-1)+1:3i] = b
        end
        info("Condition number of Gram matrix: ", cond(A'A))
        w = A\B.*mg
        if useCVX
            w = Variable(12)
            problem = minimize(norm(A*w-B.*mg,1))
            solve!(problem)
        end
        Rf = reshape(w[1:9],3,3)
        info("Determinant of Rf before orthonormalization: ", det(Rf)," (Should be close to 1, but don't be afraid if it's not, should however be positive!)")
        det(Rf) < 0 && warn("det(Rf) < 0, left handed coordinate system?")
        toOrthoNormal!(Rf)
        for ii = 1:2
            #        forceoffs = w(10:12)
            A2 = Array(Float64,3N,4)
            for i = 1:N
                RA = POSES[1:3,1:3,i]'
                At = [RA[:,3] I]
                b =  -Rf*F[i,1:3]'
                A2[3(i-1)+1:3i,:] = At
                B2[3(i-1)+1:3i] = b
            end
            w = A2\B2
            if useCVX
                w = Variable(4)
                problem = minimize(norm(A2*w-B2,1))
                solve!(problem)
            end
            mg = w[1]
            forceoffs = w[2:4]
        end
    else # Do not include offset
        A = Array(Float64,3N,9)
        for i = 1:N
            RA = POSES[1:3,1:3,i]'
            At = [F[i,1]*I, F[i,2]*I, F[i,3]*I]
            b =  -RA[1:3,3]
            A[3(i-1)+1:3i,:] = At
            B[3(i-1)+1:3i] = b
        end
        info("Condition number of Gram matrix: ", cond(A'A))
        w = A\B.*mg
        for ii = 1:4
            Rf = reshape(w[1:9],3,3)
            info("Determinant of Rf before orthonormalization: ", det(Rf))
            det(Rf) < 0 && warn("det(Rf) < 0, left handed coordinate system?")
            toOrthoNormal!(Rf)
            A2 = Array(Float64,3N,1)
            B2 = Array(Float64,3N)
            for i = 1:N
                RA = POSES[1:3,1:3,i]'
                At = RA[:,3]
                b =  -Rf*F[i,1:3]'
                A2[3(i-1)+1:3i,:] = At
                B2[3(i-1)+1:3i] = b
            end
            mg = A2\B2

        end
    end
    m = mg/g
    if m < 0
        warn("Estimated mass is negative, left handed coordinate system for sensor?")
    end
    if offset
        return Rf,m,forceoffs
    else
        return Rf,m
    end
end

# function calibForceCVX(T,f,m=0.3)
#
#     Rfv = Variable(3,3)
#     offsetv = Variable(3)
#     J = 0.0
#     m = 0.2205
#     for i = 1:10:N
#         J += sum(abs(Rfv*f[i,1:3]' + offsetv - T[1:3,1:3,i]'*[0, 0, m*-9.82]))
#     end
#     problem = minimize(J)
#     solve!(problem, SCSSolver(verbose=true,max_iters=10000))
#     Rf = Rfv.value
#     offset = Rf\offsetv.value
#     Rf = toOrthoNormal(Rf)
#     mv = Variable(1)
#     J = 0.0
#     for i = 1:10:N
#         J += sum(abs(Rf*(f[i,1:3]' + offset) - T[1:3,1:3,i]'*[0, 0, mv*-9.82]))
#     end
#     problem = minimize(J)
#     solve!(problem, SCSSolver(verbose=true,max_iters=10000))
#     m = mv.value
#
#     model = JuMP.Model(solver=IpoptSolver())
#     # model = JuMP.Model(solver=IpoptSolver(constr_viol_tol = 10., dual_inf_tol = 10., compl_inf_tol=10.))
#     JuMP.@defVar(model, g[1:3])
#     JuMP.@defVar(model, offset2[1:3])
#     JuMP.@defVar(model, Rf2[1:3,1:3])
#     # JuMP.setValue(Rf2,Rf)
#     JuMP.setValue(g,[0, 0, m*-9.82])
#     JuMP.setValue(offset2,offsetv.value)
#     J = 0.0
#
#     function absv(v::JuMP.Variable)
#         JuMP.@defVar(v.m, aux >= 0)
#         JuMP.@addConstraint(v.m, aux >= v)
#         JuMP.@addConstraint(v.m, aux >= -v)
#         return aux
#     end
#
#     function sumabsv(m,v::VecOrMat)
#         N = length(v)
#         JuMP.@defVar(m, aux[1:N] >= 0)
#         for i = 1:N
#             JuMP.@addConstraint(m, aux[i] >= v[i])
#             JuMP.@addConstraint(m, aux[i] >= -v[i])
#         end
#         return sum(aux)
#     end
#
#     for i = 1:10:N
#         X = Rf2*f[i,1:3]' + offset - T[1:3,1:3,i]'*g
#         J += sumabsv(model,X)
#     end
#
#     for i = 1:3
#         for j = 1:3
#             if i == j
#                 JuMP.@addConstraint(model, Rf2[1:3,i]'*Rf2[1:3,i] .== 1)
#                 # JuMP.@addNLConstraint(model, Rf2'Rf2 == eye(3))
#                 # JuMP.@addNLConstraint(model, sum{Rf2[i,k]*Rf2[k,j], k=1:3} == 1)
#             else
#                 JuMP.@addConstraint(model, Rf2[1:3,i]'*Rf2[1:3,j] .== 0)
#                 # JuMP.@addNLConstraint(model, sum{Rf2[i,k]*Rf2[k,j], k=1:3} == 0)
#             end
#
#         end
#     end
#     JuMP.@setObjective(model, Min, J)
#     JuMP.solve(model)
#     Rf2 = JuMP.getValue(Rf2)
#     g = JuMP.getValue(g)
#
#     err = cat(2,[Rf*(f[i,1:3]' + offset) - T[1:3,1:3,i]'*[0, 0, m*-9.82] for i = 1:N]...)'
#     err2 = cat(2,[Rf2*(f[i,1:3]' + offset) - T[1:3,1:3,i]'*g for i = 1:N]...)'
#     plot(f[:,1:3],lab="Force")
#     plot!(err,l=:dash,lab="Error")
#     plot!(err,l=:dash,lab="Error2")
#     println("Error:  ", round(rms(err),4))
#     println("Error2: ", round(rms(err2),4))
#
# end


text:
"""
$(TYPEDEF)

`ULELiouvillian` defines the Liouvillian operator corresponding the universal Lindblad equation.

# Fields

$(FIELDS)
"""
struct ULELiouvillian <: AbstractLiouvillian
    """Lindblad kernels"""
    kernels::Any
    """close system unitary"""
    unitary::Any
    """absolute error tolerance for integration"""
    atol::Float64
    """relative error tolerance for integration"""
    rtol::Float64
    """cache matrix for inplace unitary"""
    Ut::Union{Matrix,MMatrix}
    """cache matrix for inplace unitary"""
    Uτ::Union{Matrix,MMatrix}
    """cache matrix for integration"""
    Λ::Union{Matrix,MMatrix}
    """cache matrix for Lindblad operator"""
    LO::Union{Matrix,MMatrix}
    """tf minus coarse grain time scale"""
    Ta::Real
end

function ULELiouvillian(kernels, U, Ta, atol, rtol)
    m_size = size(kernels[1][2])
    Λ = m_size[1] <= 10 ? zeros(MMatrix{m_size[1],m_size[2],ComplexF64}) :
        zeros(ComplexF64, m_size[1], m_size[2])

    unitary = isinplace(U) ? U.func : (cache, t) -> cache .= U(t)
    ULELiouvillian(kernels, unitary, atol, rtol, similar(Λ),
        similar(Λ), similar(Λ), Λ, Ta)
end

function (L::ULELiouvillian)(du, u, p, t)
    s = p(t)
    LO = fill!(L.LO, 0.0)
    for (inds, coupling, cfun) in L.kernels
        for (i, j) in inds
            function integrand(cache, x)
                L.unitary(L.Ut, t)
                L.unitary(L.Uτ, x)
                L.Ut .= L.Ut * L.Uτ'
                mul!(L.Uτ, coupling[j](s), L.Ut')
                # The 5 arguments mul! will to produce NaN when it
                # should not. May switch back to it when this is fixed.
                # mul!(cache, L.Ut, L.Uτ, cfun[i, j](t, x), 0)
                mul!(cache, L.Ut, L.Uτ)
                lmul!(cfun[i, j](t, x), cache)
            end
            quadgk!(
                integrand,
                L.Λ,
                max(0.0, t - L.Ta),
                min(t + L.Ta, p.tf),
                rtol=L.rtol,
                atol=L.atol,
            )
            axpy!(1.0, L.Λ, LO)
        end
    end
    du .+= LO * u * LO' - 0.5 * (LO' * LO * u + u * LO' * LO)
end

"""
$(TYPEDEF)

The Liouvillian operator in Lindblad form.
"""
struct LindbladLiouvillian <: AbstractLiouvillian
    "1-d array of Lindblad rates"
    γ::Vector
    "1-d array of Lindblad operataors"
    L::Vector
    "size"
    size::Tuple
end

Base.length(lind::LindbladLiouvillian) = length(lind.γ)

function LindbladLiouvillian(L::Vector{Lindblad})
    if any((x) -> size(x) != size(L[1]), L)
        throw(ArgumentError("All Lindblad operators should have the same size."))
    end
    LindbladLiouvillian([lind.γ for lind in L], [lind.L for lind in L], size(L[1]))
end

function (Lind::LindbladLiouvillian)(du, u, p, t)
    s = p(t)
    for (γfun, Lfun) in zip(Lind.γ, Lind.L)
        L = Lfun(s)
        γ = γfun(s)
        du .+= γ * (L * u * L' - 0.5 * ( L' * L * u + u * L' * L))
    end
end

function update_cache!(cache, lind::LindbladLiouvillian, p, t::Real)
    s = p(t)
    for (γfun, Lfun) in zip(lind.γ, lind.L)
        L = Lfun(s)
        γ = γfun(s)
        cache .-= 0.5 * γ * L' * L
    end
end


text:
SRC_DIR = joinpath(Base.source_dir(), "..", "..", "..", "src")
OUT_DIR = joinpath(Base.source_dir(), "..", "..", "..", "output11")
mkpath(OUT_DIR)

#include(joinpath(SRC_DIR, "BridgeSDEInference.jl"))
#using Main.BridgeSDEInference
include(joinpath(SRC_DIR, "BridgeSDEInference_for_tests.jl"))

using StaticArrays, LinearAlgebra, GaussianDistributions
using Distributions
using Random
using Bridge
using CSV, DataFrames
DIR = "auxiliary"
include(joinpath(SRC_DIR, DIR, "data_simulation_fns.jl"))
include(joinpath(SRC_DIR, DIR, "utility_functions.jl"))
include(joinpath(SRC_DIR, DIR, "plotting_fns.jl"))

#=============================================================================
                Routines for setting up the MCMC sampler
=============================================================================#
function _prepare_setup(updt_order, ρ=0.96, num_mcmc_steps=4*10^3,
                        thin_path=10^0, save_path_every=1*10^3)
    model_setup = DiffusionSetup(P˟, P̃, PartObs())
    set_observations!(model_setup, Ls, Σs, obs_vals,
                      obs_time)
    set_imputation_grid!(model_setup, 1/100)
    set_x0_prior!(model_setup,
                  GsnStartingPt(x0, SArray{Tuple{4,4}, Float64}(I)), # prior over starting point
                  x0)
    set_auxiliary!(model_setup; skip_for_save=thin_path,
                   adaptive_prop=Adaptation(x0, fill(0.0, 50), fill(200, 50), 1))
    initialise!(eltype(x0), model_setup, Vern7(), false, NoChangePt())

    mcmc_setup = MCMCSetup(
        Imputation(NoBlocking(), ρ, Vern7()),
        ParamUpdate(MetropolisHastingsUpdt(), [1], fill(0.0, 8),
                    UniformRandomWalk(0.5, true), ExpUnif(-7.0, 2.0),
                    UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [1]))),
        ParamUpdate(MetropolisHastingsUpdt(), [2], fill(0.0, 8),
                    UniformRandomWalk(0.5, true), ExpUnif(-7.0, 2.0),
                    UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [2]))),
        ParamUpdate(MetropolisHastingsUpdt(), [3], fill(0.0, 8),
                    UniformRandomWalk(0.5, true), ExpUnif(-7.0, 2.0),
                    UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [3]))),
        ParamUpdate(MetropolisHastingsUpdt(), [4], fill(0.0, 8),
                    UniformRandomWalk(0.5, true), ExpUnif(-7.0, 2.0),
                    UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [4]))),
        ParamUpdate(MetropolisHastingsUpdt(), [5], fill(0.0, 8),
                    UniformRandomWalk(0.5, true), ExpUnif(-7.0, 2.0),
                    UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [5]))),
        #ParamUpdate(MetropolisHastingsUpdt(), [6], fill(0.0, 8),
        #            UniformRandomWalk(0.5, true), ExpUnif(-7.0, 2.0),
        #            UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [6]))),
        ParamUpdate(MetropolisHastingsUpdt(), [7], fill(0.0, 8),
                    UniformRandomWalk(0.5, true), ExpUnif(-7.0, 2.0),
                    UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [7]))),
        ParamUpdate(MetropolisHastingsUpdt(), [8], fill(0.0, 8),
                    UniformRandomWalk(0.5, true), ExpUnif(-7.0, 2.0),
                    UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [8]))),

    )
    schedule = MCMCSchedule(num_mcmc_steps, updt_order,
                            (save=save_path_every, verbose=10^2, warm_up=0,
                             readjust=(x->x%100==0), fuse=(x->x==2000)) )
    mcmc_setup, schedule, model_setup
end

#==============================================================================
                        Auxiliary routines for saving
==============================================================================#

function save_paths(tt, paths, filename)
    d = length(paths[1][1])
    xx = copy(tt)
    for i in 1:length(paths)
        for j in 1:d
            xx = hcat(xx, [p[j] for p in paths[i]])
        end
    end
    CSV.write(joinpath(OUT_DIR, filename), DataFrame(xx))
end

function save_marginals(tt, paths, filename, indices)
    d = length(paths[1][1])
    xx = copy(tt[indices])
    N = length(paths)
    stride = div(N, 100)
    for i in 1:length(paths)
        if i % stride == 0
            print(div(i, 100), "% done...\n")
        end
        for j in 1:d
            xx = hcat(xx, [p[j] for p in paths[i][indices]])
        end
    end
    CSV.write(joinpath(OUT_DIR, filename), DataFrame(xx'))
end

function save_param_chain(chain, filename)
    d = length(chain[1])
    xx = [c[1] for c in chain]
    for i in 2:d
        xx = hcat(xx, [c[i] for c in chain])
    end
    CSV.write(joinpath(OUT_DIR, filename), DataFrame(xx))
end

function save_obs(tt, paths, filename)
    d = length(paths[1][1])
    xx = copy(tt)
    for i in 1:length(paths)
        for j in 1:d
            xx = hcat(xx, [p[j] for p in paths[i]])
        end
    end
    CSV.write(joinpath(OUT_DIR, filename), DataFrame(xx))
end

function save_history(to_save, filename)
    n, k = length(to_save), length(to_save[1])
    xx = zeros(n, k)
    for i in 1:n
        xx[i,1:end] .= to_save[i]
    end
    CSV.write(joinpath(OUT_DIR, filename), DataFrame(xx))
end

#================================================================================
                            Generate the process
==============================================================================#

_obs = open(joinpath(OUT_DIR, "../output/prokaryote_custom.dat")) do f
    [map(x->parse(Float64, x), split(l, ' '))[[2,3]] for (i,l) in enumerate(eachline(f)) if i != 1]
end

#_obs = open(joinpath(OUT_DIR, "prokaryote_custom.dat")) do f
#    [map(x->parse(Float64, x), split(l, ' '))[2] for (i,l) in enumerate(eachline(f)) if i != 1]
#end
#obs_vals = vcat([0.0],[o for o in _obs])

_obs_time = collect(range(1.0, length(_obs), step=1))

#obs_vals = vcat([0.0],[o[2] for o in _obs])
obs_vals = vcat([[0.0]],[(i % 8 == 0 ? o : [o[2]]) for (i,o) in enumerate(_obs)])
obs_time = vcat(0.0, _obs_time)

# Permutation of observations
Random.seed!(5)
obs_RNA = [o[1] for (i,o) in enumerate(obs_vals) if i % 8 == 1][2:end]
obs_RNA_permuted = obs_RNA[randperm(length(obs_RNA))]
obs_vals = [((i % 8 == 1 && i != 1) ? [obs_RNA_permuted[div(i,8)], o[2]] : o) for (i,o) in enumerate(obs_vals)]

# let's start from the true values that generated the data
θ_init = [0.05, 0.05, 0.05, 0.05, 0.05, 0.5, 0.05, 0.05]#[0.1, 0.7, 0.35, 0.2, 0.1, 0.9, 0.3, 0.1]#
K = 10.0
P˟ = Prokaryote(θ_init..., K)
x0 = ℝ{4}([8.0,8.0,8.0,5.0])

auxFlag = Val{:custom}()
start_v = @SVector[5.0, 5.0, 5.0, 5.0]
P̃ = [ProkaryoteAux(θ_init..., K, t₀, u, T, v, auxFlag, start_v) for (t₀, T, u, v)
     in zip(obs_time[1:end-1], obs_time[2:end], obs_vals[1:end-1], obs_vals[2:end])]

Σ = @SMatrix[4.0]
L = @SMatrix[0.0 1.0 2.0 0.0]
#Ls = [L for _ in P̃]
#Σs = [Σ for _ in P̃]

Σextra = @SMatrix[1.0 0.0; 0.0 4.0]
Lextra = @SMatrix[1.0 0.0 0.0 0.0; 0.0 1.0 2.0 0.0]
Ls = [(i % 8 == 0 ? Lextra : L) for i in 1:length(P̃)]
Σs = [(i % 8 == 0 ? Σextra : Σ) for i in 1:length(P̃)]

#=
model_setup = DiffusionSetup(P˟, P̃, PartObs())
set_observations!(model_setup, [L for _ in P̃], [Σ for _ in P̃], obs_vals,
                  obs_time)
set_imputation_grid!(model_setup, 1/100)
set_x0_prior!(model_setup,
              GsnStartingPt(x0, SArray{Tuple{4,4}, Float64}(I)), # prior over starting point
              x0)
set_auxiliary!(model_setup; skip_for_save=1,
               adaptive_prop=Adaptation(x0, fill(0.0, 10), fill(100, 10), 1))
initialise!(eltype(x0), model_setup, Vern7(), false, NoChangePt())
model_setup

mcmc_setup = MCMCSetup(
    Imputation(NoBlocking(), 0.95, Vern7()),
    ParamUpdate(MetropolisHastingsUpdt(), [1], fill(0.0, 8),
                UniformRandomWalk(0.5, true), ImproperPosPrior(),
                UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [1]))),
    ParamUpdate(MetropolisHastingsUpdt(), [2], fill(0.0, 8),
                UniformRandomWalk(0.5, true), ImproperPosPrior(),
                UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [2]))),
    ParamUpdate(MetropolisHastingsUpdt(), [3], fill(0.0, 8),
                UniformRandomWalk(0.5, true), ImproperPosPrior(),
                UpdtAuxiliary(Vern7(), check_if_rec->falseompute_ODEs(P̃, [3]))),
    ParamUpdate(MetropolisHastingsUpdt(), [7], fill(0.0, 8),
                UniformRandomWalk(0.5, true), ImproperPosPrior(),
                UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [7]))))
schedule = MCMCSchedule(1*10^4, [[1,2,3,4,5]],
                        (save=10^2, verbose=10^2, warm_up=0,
                         readjust=(x->x%100==0), fuse=(x->false)) )
out = mcmc(mcmc_setup, schedule, model_setup)
=#
setup = _prepare_setup([[1,2,3,4,5,6,7,8]], 0.9, 1*10^5, 1, 10^3)
Random.seed!(6)
out, elapsed = @timeit mcmc(setup...)
plot_chains(out[2]; truth=[0.1, 0.7, 0.35, 0.2, 0.1, 0.9, 0.3, 0.1])
plot_paths(out[1], out[2], setup[2])

save_paths(out[1].time, out[1].paths, "prokaryote_paths.csv")
save_param_chain(out[2].θ_chain, "prokaryote_chain.csv")
save_history(out[2].updates[1].accpt_history, "prokaryote_accpt_hist_1.csv")
save_history(out[2].updates[2].accpt_history, "prokaryote_accpt_hist_2.csv")
save_history(out[2].updates[3].accpt_history, "prokaryote_accpt_hist_3.csv")
save_history(out[2].updates[4].accpt_history, "prokaryote_accpt_hist_4.csv")
save_history(out[2].updates[1].ρs_history, "prokaryote_rho_hist_1.csv")


text:
module TestAutoMLPipeline
using Test

# suppress warnings
@info "suppressing PyCall warnings"
using PyCall
warnings = pyimport("warnings")
warnings.filterwarnings("ignore")

# test modules
include("test_skpreprocessing.jl")
include("test_sklearner.jl")
include("test_skcrossvalidator.jl")

end


text:
"""
    F = JopExp(spc, c)

where `F` is the exponential operator e^(c/x) with domain and range given by `spc::JetSpace`, and scalar value c.
"""
JopExp(spc::JetSpace{T}, c = 1) where {T} = JopNl(dom = spc, rng = spc, f! = JopExp_f!, df! = JopExp_df!, df′! = JopExp_df′!, s = (c=T(c),))
export JopExp

JopExp_f!(d::AbstractArray, m::AbstractArray; c) = d .= exp.(m ./ c)
JopExp_df!(δd::AbstractArray, δm::AbstractArray; mₒ, c) = δd .= (δm ./ c) .* exp.(mₒ ./ c)
JopExp_df′!(δm::AbstractArray, δd::AbstractArray; mₒ, c) = δm .= (δd ./ conj(c)) .* conj.(exp.(mₒ ./ c))


text:
using Oceananigans
using Oceananigans.Units
using GLMakie

grid = RectilinearGrid(size=(128, 128), x=(0, 128), z=(-64, 0), topology=(Periodic, Flat, Bounded))

b_bcs = FieldBoundaryConditions(top=FluxBoundaryCondition(1e-8))

@inline growth_func(x, y, z, t, p) = 1 / p.τ * exp(z / p.h)
growth = Forcing(growth_func, parameters=(τ=1hour, h=4.0))
sinking = AdvectiveForcing(WENO5(), w=-1)

model = NonhydrostaticModel(; grid,
                            tracers = (:b, :P),
                            buoyancy = BuoyancyTracer(),
                            boundary_conditions = (; b=b_bcs),
                            forcing = (; P = (growth, sinking)))

bᵢ(x, y, z) = 1e-5 * z + 1e-9 * rand()
set!(model, b=bᵢ)

Δz = grid.Δzᵃᵃᶜ
Δt = 0.1 * Δz # for a sinking velocity w=1
simulation = Simulation(model; Δt, stop_iteration = 1000)

fig = Figure()
ax = Axis(fig[1, 1])

P = model.tracers.P
hm = heatmap!(ax, interior(P, :, 1, :))

function update!(sim)
    hm.input_args[1][] = interior(P, :, 1, :)
    return nothing
end

simulation.callbacks[:plot] = Callback(update!, IterationInterval(100))

display(fig)

run!(simulation)


text:
module Calcul

export interest, rate

function interest(amount, rate)
    return amount * (1 + rate)
end

function rate(amount, interest)
    return interest / amount
end

ende

max_stars_repo_name: jlvvlj/Calcul

max_issues_repo_name: jlvvlj/Calcul

max_forks_repo_name: jlvvlj/Calcul


text:
export train_chainq, quantize_chainq, experiment_chainq

"Quantize using Vitebi algorithm implemented in c++"
function quantize_chainq_cpp!(
  CODES::Matrix{Int16},  # out. Where to save the result
  X::Matrix{T},                # in. d-by-n matrix to encode
  C::Vector{Matrix{T}},        # in. m-long vector with d-by-h codebooks
  binaries::Vector{Matrix{T}}, # in. Binary terms
  IDX::UnitRange{Int64}) where T <: AbstractFloat # in. Index to save the result

  # Get unaries
  unaries = get_unaries( X, C )

  h, n = size(unaries[1])
  m    = length(C)

  if h != 256
    error("The C++ implementation of chain quantization encoding only supports
    codebooks with 256 entries")
  end

  CODES2 = zeros(Cuchar, m, n)
  unaries2, binaries2 = vcat(unaries...), hcat(binaries...)

  ccall(("viterbi_encoding", encode_icm_so), Nothing,
    (Ptr{Int16}, Ptr{Cfloat}, Ptr{Cfloat}, Cint, Cint),
    CODES2, unaries2, binaries2, n, m)

  CODES2 = convert(Matrix{Int16}, CODES2) .+ one(Int16)
  CODES[:] = CODES2[:]
end

"Quantize using Vitebi algorithm implemented in julia"
function quantize_chainq!(
  CODES::SharedMatrix{Int16},  # out. Where to save the result
  X::Matrix{T},                # in. d-by-n matrix to encode
  C::Vector{Matrix{T}},        # in. m-long vector with d-by-h codebooks
  binaries::Vector{Matrix{T}}, # in. Binary terms
  IDX::UnitRange{Int64}) where T <: AbstractFloat # in. Index to save the result

  # Get unaries
  unaries = get_unaries( X, C )

  h, n = size( unaries[1] )
  m    = length( binaries ) + 1

  # We need a matrix to keep track of the min and argmin
  mincost = zeros(T, h)
  minidx  = zeros(Int32, h, m )

  # Allocate memory for brute-forcing each pair
  cost = zeros( T, h )

  U = zeros(T, h, m)

  minv = typemax(T)
  mini = 1

  # unaries = vcat(unaries...)
  # @show typeof(unaries)

  uidx = 1
  @inbounds for idx = IDX # Loop over the datapoints

    # Put all the unaries of this item together
    # U[:] = unaries[:,idx]
    for i = 1:m
      ui = unaries[i]
      @simd for j = 1:h
        U[j,i] = ui[j,uidx]
      end
    end

    # Forward pass
    for i = 1:(m-1) # Loop over states

      # If this is not the first iteration, add the precomputed costs
      if i > 1
        @simd for j = 1:h
          U[j,i] += mincost[j]
        end
      end

      bb = binaries[i]
      for j = 1:h # Loop over the cost of going to j
        @simd for k = 1:h # Loop over the cost of coming from k
          ucost   =  U[k, i] # Pay the unary of coming from k
          bcost   = bb[k, j] # Pay the binary of going from j-k
          cost[k] = ucost + bcost
        end

        # findmin -- julia's is too slow
        minv = cost[1]
        mini = 1
        for k = 2:h
          costi = cost[k]
          if costi < minv
            minv = costi
            mini = k
          end
        end

        mincost[j] = minv
         minidx[j, i] = mini
      end
    end

    # @show mincost, minidx

    @simd for j = 1:h
      U[j,m] += mincost[j]
    end

    _, mini = findmin( U[:,m] )

    # Backward trace
    backpath = [ mini ]
    for i = (m-1):-1:1
      push!(backpath, minidx[ backpath[end], i ])
    end

    # Save the inferred code
    CODES[:, idx] = reverse( backpath )
    uidx = uidx + 1;
  end # for idx = IDX
end

"Quantize using Vitebi algorithm implemented in julia. This algorithm is batched and thus may be faster, but in practice
it's about as fast as the sequential version."
function quantize_chainq_batched!(
  CODES::SharedMatrix{Int16},  # out. Where to save the result
  X::Matrix{T},                # in. d-by-n matrix to encode
  C::Vector{Matrix{T}},        # in. m-long vector with d-by-h codebooks
  binaries::Vector{Matrix{T}}, # in. Binary terms
  IDX::UnitRange{Int64}) where T <: AbstractFloat # in. Index to save the result

  # Get unaries
  # @Profile.profile begin
  unaries = get_unaries( X, C )

  h, n = size( unaries[1] )
  m    = length( binaries ) + 1

  # We need a matrix to keep track of the min and argmin
  mincost = zeros(T, h, n)
  minidx  = zeros(Int32, h, m, n)

  # Allocate memory for brute-forcing each pair
  cost = zeros(T, h, n)

  minv = typemax(T)
  mini = 1

  # Put all the unaries of this item together
  # unaries = vcat(unaries...)

  # Forward pass
  @inbounds for i = 1:(m-1) # Loop over states

    if i > 1; unaries[i] .+= mincost; end
    ucost = unaries[i]

    bb = binaries[i]
    for j = 1:h # Loop over the cost of going to j
      bcost = bb[:,j]
      # cost  = ucost .+ bcost

      for kk = IDX
        @simd for k=1:h
          cost[k,kk] = ucost[k,kk] + bcost[k]
          # ui[k,kk] += bcost[k]
        end
      end

      # Findmin
      minv, mini = findmin(cost,1)
      mincost[j,:]  = minv
      minidx[j,i,:] = rem.(mini-1,h) + 1
    end
  end

  unaries[m] .+= mincost
  _, mini = findmin(unaries[m],1)
  mini = rem.(mini-1,h) + 1

  @inbounds for idx = IDX # Loop over the datapoints
    # Backward trace
    backpath = [ mini[idx] ]
    for i = (m-1):-1:1
      push!(backpath, minidx[ backpath[end], i, idx ])
    end

    # Save the inferred code
    CODES[:, idx] = reverse( backpath )
  end # for idx = IDX
  # end
end

"Quantize using Vitebi algorithm implemented in CUDA"
function quantize_chainq_cuda!(
  CODES::Matrix{Int16},  # out. Where to save the result
  X::Matrix{T},                # in. d-by-n matrix to encode
  C::Vector{Matrix{T}},        # in. m-long vector with d-by-h codebooks
  binaries::Vector{Matrix{T}}, # in. Binary terms
  IDX::UnitRange{Int64}) where T <: AbstractFloat # in. Index to save the result

  d, n = size( X )
  _, h = size( C[1] )
  m    = length( binaries ) + 1

  # We need a matrix to keep track of the min and argmin
  mincost = zeros(T, h, n)
  # minidx  = zeros(Int32, h, m, n)
  minidx = convert(Array{Int32}, rand(1:h, h, m, n))

  # Allocate memory for brute-forcing each pair
  cost = zeros(T, h, n)

  minv = zeros(T, n)
  mini = zeros(Int32, n)

  # Setup GPU stuff
  gpuid = 0
  CudaUtilsModule.init(gpuid, cudautilsptx)

  # Make space for unaries
  d_X = CuArrays.CuArray(X)
  d_unaries   = Vector{CuArrays.CuArray{T}}(undef, m)
  d_codebooks = Vector{CuArrays.CuArray{T}}(undef, m)
  for j = 1:m
    d_codebooks[j] = CuArrays.CuArray(C[j])

    # -2 * C' * X
    # d_unaries[j] = -2.0f0 * d_codebooks[j]' * d_RX <-- thus runs out of memory real fast
    d_unaries[j] = CuArrays.BLAS.gemm('T', 'N', -2.0f0, d_codebooks[j], d_X)

    # Add self-codebook interactions || C_{i,i} ||^2
    CudaUtilsModule.vec_add(n, (1,h), d_unaries[j].buf, CuArrays.CuArray(sum(C[j].^2, dims=1)).buf, Cint(n), Cint(h))
  end

  d_binaries = Vector{CuArrays.CuArray{T}}(undef, length(binaries))
  for i = 1:length(binaries)
    d_binaries[i] = CuArrays.CuArray(binaries[i])
  end

  d_mincost = CuArrays.CuArray(mincost)
  d_mini = CuArrays.CuArray(mini)

  # Forward pass
  @inbounds for i = 1:(m-1) # Loop over states
    if i > 1; CuArrays.BLAS.axpy!(n * h, 1.0f0, d_mincost, 1, d_unaries[i], 1); end

    for j = 1:h # Loop over the cost of going to j
      CudaUtilsModule.viterbi_forward(n, (1, h), d_unaries[i].buf, d_binaries[i].buf, d_mincost.buf, d_mini.buf, Cint(n), Cint(j-1))

      Mem.download!(mini, d_mini.buf)
      minidx[j,i,:] .= mini .+ one(eltype(mini))
    end
  end

  CuArrays.BLAS.axpy!(n * h, 1.0f0, d_mincost, 1, d_unaries[m], 1)
  CudaUtilsModule.viterbi_forward(n, (1, h), d_unaries[m].buf, CuArrays.CuArray(zeros(Float32, h, h)).buf, d_mincost.buf, d_mini.buf, Cint(n), Cint(0))
  Mem.download!(mini, d_mini.buf)
  mini .+= one(eltype(mini))

  # Backward trace
  @inbounds for idx = IDX # Loop over the datapoints

    backpath = [ mini[idx] ]
    for i = (m-1):-1:1
      push!(backpath, minidx[ backpath[end], i, idx ])
    end

    # Save the inferred code
    CODES[:, idx] .= reverse!( backpath )
  end

  CudaUtilsModule.finit()

  nothing
end

"""
    quantize_chainq(X, C, use_cuda=false, use_cpp=false) -> B, ellapsed

Given data and chain codebooks, find codes using the Viterbi algorithm chain quantizer.

# Arguments
- `X::Matrix{T}`: `d`-by-`n` data to quantize
- `C::Vector{Matrix{T}}`: `m`-long vector with `d`-by-`h` matrices. Each matrix is a pretrained codebook of size approximately `d`-by-`h`.
- `use_cuda::Bool`: whether to use a CUDA implementation
- `use_cpp::Bool`: whether to use a c++ implementation

If both `use_cuda` and `use_cpp` are `true`, the CUDA implementation is used.

# Returns
- `B::Matrix{Int16}`: `m`-by-`n` matrix with the codes
- `ellapsed::Float64`: The time spent encoding
"""
function quantize_chainq(
  X::Matrix{Float32},         # d-by-n matrix. Data to encode
  C::Vector{Matrix{Float32}}, # m-long vector with d-by-h codebooks
  use_cuda::Bool=false,
  use_cpp::Bool=false)

  start_time = time_ns()
  d, n = size(X)
  m    = length(C)

  # Compute binary tables
  binaries = Vector{Matrix{Float32}}(undef, m-1)
  for i = 1:(m-1)
    binaries[i] = 2 * C[i]' * C[i+1]
  end

  CODES = SharedArray{Int16,2}(m, n)
  CODES .= zero(Int16)

  if use_cuda # CUDA version
    quantize_chainq_cuda!(sdata(CODES), X, C, binaries, 1:n)
    GC.gc()
  elseif nworkers() == 1
    if use_cpp # C++ parallel version
      quantize_chainq_cpp!(sdata(CODES), X, C, binaries, 1:n)
      # quantize_chainq_batched!(CODES, X, C, binaries, 1:n)
    else # Julia single-threaded version
      quantize_chainq!(CODES, X, C, binaries, 1:n)
    end
  else # Julia-paralell version
    paridx = splitarray(1:n, nworkers())
    @sync begin
      for (i,wpid) in enumerate(workers())
        @async begin
          Xw = X[:,paridx[i]]
          remotecall_wait(quantize_chainq!, wpid, CODES, Xw, C, binaries, paridx[i])
        end
      end
    end
  end

  B = sdata(CODES)
  B, (time_ns() - start_time)/1e9
end

"""
    train_chainq(X, m, h, R, B, C, niter, V=false) -> C, B, R, error

Train a chain quantizer.
This method is typically initialized by [Optimized product quantization (OPQ)](@ref)

# Arguments
- `X::Matrix{T}`: `d`-by-`n` data to quantize
- `m::Integer`: Number of codebooks
- `h::Integer`: Number of entries in each codebook (typically 256)
- `R::Matrix{T}`: `d`-by-`d` rotation matrix for initialization
- `B::Matrix{Int16}`: `m`-by-`n` matrix with pre-trained codes for initialization
- `C::Vector{Matrix{T}}`: `m`-long vector with `d`-by-`h` matrices. Each matrix is a pretrained codebook of size approximately `d`-by-`h`.
- `niter::Integer`: Number of iterations to use
- `V::Bool`: Whether to print progress

# Returns
- `C::Vector{Matrix{T}}`: `m`-long vector with `d`-by-`h` matrix entries. Each matrix is a codebook of size approximately `d`-by-`h`.
- `B::Matrix{Int16}`: `m`-by-`n` matrix with the codes
- `R::Matrix{T}`: `d`-by-`d` optimized rotation matrix
- `error::T`: The quantization error after training
"""
function train_chainq(
  X::Matrix{T},             # d-by-n matrix of data points to train on.
  m::Integer,               # number of codebooks
  h::Integer,               # number of entries per codebook
  R::Matrix{T},             # Init rotation matrix
  B::Matrix{Int16},         # Init codes
  C::Vector{Matrix{T}},     # Init codebooks
  niter::Integer,           # number of optimization iterations
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  if V; @printf("Training a chain quantizer\n"); end

  d, n = size(X)
  obj  = zeros(T, niter+1)

  # TODO expose these
  use_cuda = true
  use_cpp = false

  CB = zeros(T, size(X))
  RX = R' * X

  # Initialize C
  # C, Ctime = update_codebooks_chain(RX, B, h)
  C, Ctime = update_codebooks_chain_bin(RX, B, h)
  if V; @printf("%3d %e... %.2f secs updating C\n", -2, qerror(RX, B, C), Ctime); end

  # Initialize B
  B, Btime = quantize_chainq(RX, C, use_cuda, use_cpp)
  if V; @printf("%3d %e... %.2f secs updating B\n", -1, qerror(RX, B, C), Btime); end

  for iter = 0:niter
    if V; start_time = time_ns(); end # Take time if asked to

    obj[iter+1] = qerror(RX, B, C)
    if V; @printf("%3d %e... ", iter, obj[iter+1]); end

    # update CB
    CB .= zero(T)
    for i = 1:m; CB .+= C[i][:, vec(B[i,:])]; end

    # update R
    U, S, VV = svd(X * CB', full=false)
    R = U * VV'

    # update R*X
    RX = R' * X

    # Update the codebooks #
    # C, Ctime = update_codebooks_chain(RX, B, h)
    C, Ctime = update_codebooks_chain_bin(RX, B, h)

    # Update the codes with lattice search
    B, Btime = quantize_chainq(RX, C, use_cuda, use_cpp)
    if V; @printf("done in %.2f secs. %.2f secs updating B with %d workers. %.2f secs updating C\n", (time_ns()-start_time)/1e9, Btime, nworkers(), Ctime); end
  end

  return C, B, R, obj
end


text:
# default colors
orange = RGBA(1.0, 153.0 / 255.0, 51.0 / 255.0, 1.0)
cyan = RGBA(51.0 / 255.0, 1.0, 1.0, 1.0)
magenta = RGBA(255.0 / 255.0, 0.0, 255.0 / 255.0, 1.0)
gray = RGBA(0.7, 0.7, 0.7, 1.0)

# light colors
orange_light = RGBA(1.0, 153.0 / 255.0 + 0.2, 51.0 / 255.0 + 0.2, 1.0)
cyan_light = RGBA(51.0 / 255.0 + 0.6, 1.0, 1.0, 1.0)
magenta_light = RGBA(255.0 / 255.0, 0.75, 255.0 / 255.0, 1.0)
gray_light = RGBA(0.9, 0.9, 0.9, 1.0)


text:
function levendist(s::AbstractString, t::AbstractString)
    ls, lt = length.((s, t))
    ls == 0 && return lt
    lt == 0 && return ls

    s₁, t₁ = s[2:end], t[2:end]
    ld₁ = levendist(s₁, t₁)
    s[1] == t[1] ? ld₁ : 1 + min(ld₁, levendist(s, t₁), levendist(s₁, t))
end

@show levendist("kitten", "sitting") # 3
@show levendist("rosettacode", "raisethysword") # 8

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
struct DynamicsConstraints <: Constraints
    n::Int
    ineq
    w
end

function dynamics_constraints(model, T;
    w = [zeros(model.d) for t = 1:T-1])

    n = model.n * (T - 1)
    ineq = (1:0)

    return DynamicsConstraints(n, ineq, w)
end

function constraints!(c, Z, con::DynamicsConstraints, model, idx, h, T)
    n = model.n

    # note: x1 and xT constraints are handled as simple bound constraints
    #       e.g., x1 <= x <= x1, xT <= x <= xT

    for t = 1:T-1
        x = view(Z, idx.x[t])
        u = view(Z, idx.u[t])
        x⁺ = view(Z, idx.x[t + 1])

        c[(t-1) * n .+ (1:n)] = fd(model, x⁺, x, u, con.w[t], h, t)
    end

    return nothing
end

function constraints_jacobian!(∇c, Z, con::DynamicsConstraints, model, idx, h, T)
    n = model.n
    m = model.m

    # note: x1 and xT constraints are handled as simple bound constraints
    #       e.g., x1 <= x <= x1, xT <= x <= xT

    shift = 0

    for t = 1:T-1
        x = view(Z, idx.x[t])
        u = view(Z, idx.u[t])
        x⁺ = view(Z, idx.x[t + 1])

        dyn_x(z) = fd(model, x⁺, z, u, con.w[t], h, t)
        dyn_u(z) = fd(model, x⁺, x, z, con.w[t], h, t)
        dyn_x⁺(z) = fd(model, z, x, u, con.w[t], h, t)

        r_idx = (t-1) * n .+ (1:n)

        s = n * n
        ∇c[shift .+ (1:s)] = vec(ForwardDiff.jacobian(dyn_x, x))
        shift += s

        s = n * m
        ∇c[shift .+ (1:s)] = vec(ForwardDiff.jacobian(dyn_u, u))
        shift += s

        s = n * n
        ∇c[shift .+ (1:s)] = vec(ForwardDiff.jacobian(dyn_x⁺, x⁺))
        shift += s
    end

    return nothing
end

function constraints_sparsity(con::DynamicsConstraints, model, idx, T;
    shift_row = 0, shift_col = 0)
    n = model.n

    row = []
    col = []

    for t = 1:T-1
        r_idx = shift_row + (t-1) * n .+ (1:n)
        row_col!(row, col, r_idx, shift_col .+ idx.x[t])
        row_col!(row, col, r_idx, shift_col .+ idx.u[t])
        row_col!(row, col, r_idx, shift_col .+ idx.x[t + 1])
    end

    return collect(zip(row, col))
end


text:
# Julia wrapper for header: /usr/local/include/libavutil/channel_layout.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_get_channel_layout,
    av_get_channel_layout_string,
    av_get_channel_layout_nb_channels,
    av_get_default_channel_layout,
    av_get_channel_layout_channel_index,
    av_channel_layout_extract_channel,
    av_get_channel_name

function av_get_channel_layout(name)
    ccall((:av_get_channel_layout,libavutil),Uint64,(Ptr{Uint8},),name)
end

function av_get_channel_layout_string(buf,buf_size::Integer,nb_channels::Integer,channel_layout::Uint64)
    ccall((:av_get_channel_layout_string,libavutil),Void,(Ptr{Uint8},Cint,Cint,Uint64),buf,buf_size,nb_channels,channel_layout)
end

function av_get_channel_layout_nb_channels(channel_layout::Uint64)
    ccall((:av_get_channel_layout_nb_channels,libavutil),Cint,(Uint64,),channel_layout)
end

function av_get_default_channel_layout(nb_channels::Integer)
    ccall((:av_get_default_channel_layout,libavutil),Uint64,(Cint,),nb_channels)
end

function av_get_channel_layout_channel_index(channel_layout::Uint64,channel::Uint64)
    ccall((:av_get_channel_layout_channel_index,libavutil),Cint,(Uint64,Uint64),channel_layout,channel)
end

function av_channel_layout_extract_channel(channel_layout::Uint64,index::Integer)
    ccall((:av_channel_layout_extract_channel,libavutil),Uint64,(Uint64,Cint),channel_layout,index)
end

function av_get_channel_name(channel::Uint64)
    ccall((:av_get_channel_name,libavutil),Ptr{Uint8},(Uint64,),channel)
end


text:
################
# OuterProduct #
################
struct OuterProduct{P,N,S,K<:Ket,B<:Bra} <: DiracOp{P,N}
    scalar::S
    kt::K
    br::B
    # OuterProduct(scalar::S, kt::Ket{P,N}, br::Bra{P,N}) = new(scalar, kt, br)
end

OuterProduct(scalar::S, kt::Ket{P,N}, br::Bra{P,N}) where {P,N,S} = OuterProduct{P,N,S,typeof(kt),typeof(br)}(scalar, kt, br)

ptype(op::OuterProduct) = ptype(op.kt)

Base.copy(op::OuterProduct) = OuterProduct(copy(op.scalar), copy(op.kt), copy(op.br))

function Base.convert(::Type{OpSum}, op::OuterProduct{P,N}) where {P,N}
    result = OpSum(ptype(op), cons_outer!(OpDict{N,eltype(op)}(), op.kt, op.br))
    return scale!(op.scalar, result)
end

Base.promote_rule(::Type{G}, ::Type{O}) where {G<:OpSum, O<:OuterProduct} = OpSum

#######################
# Dict-Like Functions #
#######################
Base.eltype(op::OuterProduct) = promote_type(typeof(op.scalar), eltype(op.kt), eltype(op.br))

# these equality/hash functions are pretty inefficient...
Base.:(==)(a::OuterProduct{P,N}, b::OuterProduct{P,N}) where {P,N} = convert(OpSum, a) == convert(OpSum, b)
Base.hash(op::OuterProduct) = hash(convert(OpSum, op))
Base.hash(op::OuterProduct, h::UInt64) = hash(hash(op), h)

Base.length(op::OuterProduct) = length(op.kt)*length(op.br)

Base.getindex(op::OuterProduct, k, b) = op.scalar * op.kt[k] * op.br[b]
Base.getindex(op::OuterProduct, o::OpLabel) = op[klabel(o), blabel(o)]

# would be great if the below worked with normal indexing
# notation (e.g. op[k,:]) but the Colon is apparently
# special and doesn't dispatch directly to getindex...
# Base.getindex(op::OuterProduct, k, ::Colon) = (op.scalar * op.kt[k]) * op.br
# Base.getindex(op::OuterProduct, ::Colon, b) = (op.scalar * op.br[b]) * op.kt
# Base.getindex(op::OuterProduct, ::Colon, ::Colon) = convert(OpSum, op)

getbra(op::OuterProduct, k) = (op.scalar * get(op.kt,k)) * op.br
getket(op::OuterProduct, b) = (op.scalar * get(op.br,b)) * op.kt

Base.haskey(op::OuterProduct, k, b) = hasket(op, k) && hasbra(op, b)
Base.haskey(op::OuterProduct, o::OpLabel) = haskey(op, klabel(o), blabel(o))

Base.get(op::OuterProduct, k, b, default=predict_zero(eltype(op))) = haskey(op, k, b) ? op[k,b] : default
Base.get(op::OuterProduct, o::OpLabel, default=predict_zero(eltype(op))) = get(op, klabel(o), blabel(o), default)

Base.collect(op::OuterProduct{P,N}) where {P,N} = collect_pairs!(Array((OpLabel{N}, eltype(op)), length(op)), op)

function collect_pairs!(result, op::OuterProduct)
    i = 1
    for (k,kc) in dict(op.kt)
        for (b,bc) in dict(op.br)
            result[i] = (OpLabel(k, b), op.scalar * kc * bc')
            i += 1
        end
    end
    return result
end

##############
# ctranspose #
##############
Base.adjoint(op::OuterProduct) = OuterProduct(op.scalar', op.br', op.kt')

#########
# inner #
#########
inner(br::Bra, op::OuterProduct) = op.scalar * inner(br, op.kt) * op.br
inner(op::OuterProduct, kt::Ket) = op.scalar * op.kt * inner(op.br, kt)
inner(a::OuterProduct, b::OuterProduct) = OuterProduct(a.scalar * b.scalar * inner(a.br,b.kt), a.kt, b.br)
inner(a::OuterProduct, b::AbsOpSum) = a.scalar * a.kt * inner(a.br, b)
inner(a::AbsOpSum, b::OuterProduct) = inner(a, b.kt) * b.br * b.scalar

##########
# act_on #
##########
act_on(op::OuterProduct, kt::Ket, i) = op.scalar * op.kt * act_on(op.br, kt, i)
act_on(op::OuterProduct, br::Bra, i) = op.scalar * act_on(op.kt, br, i) * op.kt

##########
# tensor #
##########
tensor(kt::Ket, br::Bra) = OuterProduct(1, kt, br)
tensor(br::Bra, kt::Ket) = tensor(kt, br)
tensor(a::OuterProduct, b::OuterProduct) = OuterProduct(a.scalar * b.scalar, tensor(a.kt,b.kt), tensor(a.br, b.br))

###########
# Scaling #
###########
scale!(c::Number, op::OuterProduct) = (op.scalar = c*op.scalar; return op)
scale!(op::OuterProduct, c::Number) = (op.scalar = op.scalar*c; return op)

scale(c::Number, op::OuterProduct) = OuterProduct(c * op.scalar, copy(op.kt), copy(op.br))
scale(op::OuterProduct, c::Number) = OuterProduct(op.scalar * c, copy(op.kt), copy(op.br))

###########
# + and - #
###########
Base.:-(op::OuterProduct) = scale(-1, op)
Base.:+(a::OuterProduct, b::OuterProduct) = convert(OpSum, a) + convert(OpSum, b)

#################
# Normalization #
#################
function LinearAlgebra.norm(op::OuterProduct)
    result = predict_zero(eltype(op))
    for v in values(dict(op.kt)), c in values(dict(op.br))
        result += abs2(op.scalar * v * c')
    end
    return sqrt(result)
end

#########
# Trace #
#########
function LinearAlgebra.tr(op::OuterProduct)
    result = predict_zero(eltype(op))
    for (k,v) in dict(op.kt), (b,c) in dict(op.br)
        if b == k
            result += v * c'
        end
    end
    return op.scalar * result
end

#################
# Partial Trace #
#################
function ptrace_dict!(result, op::OuterProduct, over)
    for k in keys(dict(op.kt)), b in keys(dict(op.br))
        if k[over] == b[over]
            add_to_dict!(result, traceout(k, b, over), op[k,b])
        end
    end
    return result
end

function ptrans_dict!(result, op::OuterProduct, over)
    for k in keys(dict(op.kt)), b in keys(dict(op.br))
        add_to_dict!(result, ptranspose(k, b, over), op[k,b])
    end
    return result
end

########################
# Misc. Math Functions #
########################
nfactors(::OuterProduct{P,N}) where {P,N} = N
xsubspace(op::OuterProduct,x) = xsubspace(convert(OpSum, op), x)
filternz(op::OuterProduct) = filternz(convert(OpSum, op))
switch(op::OuterProduct, i, j) = switch(convert(OpSum, op), i, j)
permute(op::OuterProduct, perm::Vector) = permute(convert(OpSum, op), perm)

######################
# Printing Functions #
######################
labelrepr(op::OuterProduct, k, b, pad) = "$pad$(op[k,b]) $(ktstr(k))$(brstr(b))"

Base.summary(op::OuterProduct) = "OuterProduct with $(length(op)) operator(s); $(typeof(op.kt)) * $(typeof(op.br))"

function Base.show(io::IO, op::OuterProduct)
    print(io, summary(op)*":")
    pad = "  "
    maxlen = 10
    i = 1
    for k in keys(dict(op.kt)), b in keys(dict(op.br))
        if i > maxlen
            break
        else
            i += 1
            println(io)
            print(io, labelrepr(op, k, b, pad))
        end
    end
    if length(op) > maxlen
        println(io)
        print(io, "$pad$vdots")
    end
end

export getbra,
    getket,
    hasket,
    hasbra,
    ptrace,
    purity,
    filternz,
    xsubspace,
    nfactors


text:
using BinDeps
using CMakeWrapper
using Pkg

@BinDeps.setup

function cflags_validator(pkg_names...)
    return (name, handle) -> begin
        for pkg_name in pkg_names
            try
                run(`pkg-config --cflags $(pkg_name)`)
                return true
            catch ErrorException
            end
        end
        false
    end
end

# Explicitly disallow global LCM installations
validate_lcm(name, handle) = occursin(dirname(@__FILE__), name)

@static if Sys.islinux()
    python = library_dependency("python", aliases=["libpython2.7.so", "libpython3.2.so", "libpython3.3.so", "libpython3.4.so", "libpython3.5.so", "libpython3.6.so", "libpython3.7.so", "libpython3.8.so"], validate=cflags_validator("python", "python2", "python3"))
    glib = library_dependency("glib", aliases=["libglib-2.0-0", "libglib-2.0", "libglib-2.0.so.0"], depends=[python], validate=cflags_validator("glib-2.0"))
    lcm = library_dependency("lcm", aliases=["liblcm", "liblcm.1"], depends=[glib],
                             validate=validate_lcm)

    provides(AptGet, Dict("python-dev" => python, "libglib2.0-dev" => glib))
else
    glib = library_dependency("glib", aliases = ["libglib-2.0-0", "libglib-2.0", "libglib-2.0.so.0"])
    lcm = library_dependency("lcm", aliases=["liblcm", "liblcm.1"], depends=[glib],
                             validate=validate_lcm)
end

prefix = joinpath(BinDeps.depsdir(lcm), "usr")

lcm_cmake_arguments = [
    "-DCMAKE_BUILD_TYPE=Release",
    "-DLCM_ENABLE_TESTS:BOOL=OFF",
    "-DLCM_ENABLE_EXAMPLES:BOOL=OFF",
    "-DLCM_ENABLE_JAVA:BOOL=OFF",
    "-DLCM_ENABLE_LUA:BOOL=OFF"
]

@static if Sys.isapple()
    using Homebrew
    provides(Homebrew.HB, "glib", glib, os=:Darwin)
    push!(lcm_cmake_arguments,
        "-DCMAKE_LIBRARY_PATH=$(joinpath(dirname(pathof(Homebrew)), "..", "deps", "usr", "lib"))",
        "-DCMAKE_INCLUDE_PATH=$(joinpath(dirname(pathof(Homebrew)), "..", "deps", "usr", "include"))")
end

provides(Yum,
    Dict("glib" => glib))

lcm_sha = "0e6472188332bf5683d1ccdd29ba6855ac3cd593"
lcm_folder = "lcm-$(lcm_sha)"

provides(Sources,
    URI("https://github.com/lcm-proj/lcm/archive/$(lcm_sha).zip"),
    lcm,
    unpacked_dir=lcm_folder)

provides(BuildProcess, CMakeProcess(cmake_args=lcm_cmake_arguments),
         lcm,
         onload="""
const lcm_prefix = "$prefix"
""")

@BinDeps.install Dict(:lcm => :liblcm)


text:
using Distributions

function gibbs(N::Int64,thin::Int64)
  xa=Array(Float64,N)
  ya=Array(Float64,N)
  x=0.0
  y=0.0
  for i in 1:N
    for j in 1:thin
      x=rand(Gamma(3.0,1.0/(y*y+4.0)))
      y=rand(Normal(1.0/(x+1.0),1.0/sqrt(2.0*x+2.0)))
    end
    xa[i]=x
    ya[i]=y
  end
  return(xa,ya)
end

# Warm up
xaj, yaj = gibbs(1,1)

# Time
@elapsed xaj,yaj=gibbs(50000,1000)

max_stars_repo_name: ruqianl/distributions

max_issues_repo_name: ruqianl/distributions

max_forks_repo_name: ruqianl/distributions