text:
using Dates
using Test

using TimeSeries
using MarketData
using DataFrames
using Tables
using CSV

@testset "Tables.jl integration" begin

@testset "interface" begin

    @testset "single column" begin
        @test Tables.istable(cl)
        @test Tables.istable(typeof(cl))
        @test Tables.rowaccess(cl)
        @test Tables.rowaccess(typeof(cl))
        @test Tables.columnaccess(cl)
        @test Tables.columnaccess(typeof(cl))

        sch = Tables.schema(cl)
        @test sch.names == (:timestamp, :Close)
        @test sch.types == (Date, Float64)
    end

    @testset "multiple column" begin
        @test Tables.istable(ohlc)
        @test Tables.istable(typeof(ohlc))
        @test Tables.rowaccess(ohlc)
        @test Tables.rowaccess(typeof(ohlc))
        @test Tables.columnaccess(ohlc)
        @test Tables.columnaccess(typeof(ohlc))

        sch = Tables.schema(ohlc)
        @test sch.names == (:timestamp, :Open, :High, :Low, :Close)
        @test sch.types == (Date, Float64, Float64, Float64, Float64)
    end
end

@testset "iterator" begin
    @testset "single column" begin
        @test Tables.columnnames(cl)           == [:timestamp; colnames(cl)]
        @test Tables.getcolumn(cl, 1)          == timestamp(cl)
        @test Tables.getcolumn(cl, :timestamp) == timestamp(cl)
        @test Tables.getcolumn(cl, 2)          == values(cl)
        @test Tables.getcolumn(cl, :Close)     == values(cl)

        # column iterator
        i = Tables.columns(cl)
        @test size(i)      == (length(cl), 2)
        @test length(i)    == 2  # timestamp and the close columns
        @test i[100, 1]    == timestamp(cl)[100]
        @test i[100, 2]    == values(cl)[100]
        @test i[end]       == values(cl)
        @test i[end, 1]    == timestamp(cl)[end]
        @test i[end, 2]    == values(cl)[end]
        @test timestamp(i) == timestamp(cl)
        @test colnames(i)  == colnames(cl)
        @test i.Close      == values(cl)
        @test Tables.getcolumn(i, 1)          == timestamp(cl)
        @test Tables.getcolumn(i, :timestamp) == timestamp(cl)
        @test Tables.getcolumn(i, 2)          == values(cl)
        @test Tables.getcolumn(i, :Close)     == values(cl)

        @test propertynames(i) == (:timestamp, :Close)

        # row iterator
        i = Tables.rows(cl)
        for r ∈ i
            @test r.timestamp == timestamp(cl)[1]
            @test r.Close     == values(cl)[1]
            break
        end
    end

    @testset "multi column" begin
        @test Tables.columnnames(ohlc)           == [:timestamp; colnames(ohlc)]
        @test Tables.getcolumn(ohlc, 1)          == timestamp(ohlc)
        @test Tables.getcolumn(ohlc, :timestamp) == timestamp(ohlc)
        @test Tables.getcolumn(ohlc, 3)          == values(ohlc[:High])
        @test Tables.getcolumn(ohlc, :High)      == values(ohlc[:High])

        # column iterator
        i = Tables.columns(ohlc)
        @test size(i)      == (length(ohlc), 5)
        @test length(i)    == 5
        @test i[100, 1]    == timestamp(ohlc)[100]
        @test i[100, 3]    == values(ohlc)[100, 2]
        @test i[end]       == values(ohlc[:Close])
        @test i[end, 1]    == timestamp(ohlc)[end]
        @test i[end, 3]    == values(ohlc)[end, 2]
        @test timestamp(i) == timestamp(ohlc)
        @test colnames(i)  == colnames(ohlc)
        @test i.Open       == values(ohlc.Open)
        @test Tables.getcolumn(i, 1)          == timestamp(ohlc)
        @test Tables.getcolumn(i, :timestamp) == timestamp(ohlc)
        @test Tables.getcolumn(i, 3)          == values(ohlc[:High])
        @test Tables.getcolumn(i, :High)      == values(ohlc[:High])

        @test propertynames(i) == (:timestamp, :Open, :High, :Low, :Close)

        # row iterator
        i = Tables.rows(ohlc)
        for r ∈ i
            @test r.timestamp == timestamp(ohlc)[1]
            @test r.Open      == values(ohlc)[1, 1]
            @test r.High      == values(ohlc)[1, 2]
            @test r.Low       == values(ohlc)[1, 3]
            @test r.Close     == values(ohlc)[1, 4]
            break
        end
    end
end  # @testset "iterator"

@testset "DataFrames.jl" begin
    @testset "single column" begin
        df = DataFrame(cl)
        @test names(df)    == [:timestamp; colnames(cl)]
        @test df.timestamp == timestamp(cl)
        @test df.Close     == values(cl.Close)
    end

    @testset "multi column" begin
        df = DataFrame(ohlc)
        @test names(df)    == [:timestamp; colnames(ohlc)]
        @test df.timestamp == timestamp(ohlc)
        @test df.Open      == values(ohlc.Open)
        @test df.High      == values(ohlc.High)
        @test df.Low       == values(ohlc.Low)
        @test df.Close     == values(ohlc.Close)
    end

    @testset "column name collision" begin
        ta = TimeArray(ohlc, colnames = [:Open, :High, :timestamp, :Close])
        df = DataFrame(ta)
        @test names(df)      == [:timestamp, :Open, :High, :timestamp_1, :Close]
        @test df.timestamp   == timestamp(ta)
        @test df.Open        == values(ta.Open)
        @test df.High        == values(ta.High)
        @test df.timestamp_1 == values(ta.timestamp)
        @test df.Close       == values(ta.Close)

        # no side effect on column renaming
        @test colnames(ta)  == [:Open, :High, :timestamp, :Close]
    end

    @testset "DataFrame to TimeArray" begin
        ts = Date(2018, 1, 1):Day(1):Date(2018, 1, 3)
        df = DataFrame(A  = [1., 2, 3],
                       B  = [4, 5, 6],
                       C  = [7, 8, 9],
                       ts = ts)
        ta = TimeArray(df; timestamp = :ts)

        @test timestamp(ta) == ts
        @test colnames(ta)  == [:A, :B, :C]
        @test meta(ta)       ≡ df
        @test values(ta.A)  == [1., 2, 3]
        @test values(ta.B)  == [4, 5., 6]
        @test values(ta.C)  == [7, 8, 9.]
    end
end  # @testset "DataFrames.jl"

@testset "CSV.jl" begin
    @testset "single column" begin
        ta = TimeArray(cl[1:5], values = [1.1, 2.2, 3.3, 4.4, 5.5])
        io = IOBuffer()
        CSV.write(io, ta)
        @test String(take!(io)) == "timestamp,Close\n" *
                                   "2000-01-03,1.1\n" *
                                   "2000-01-04,2.2\n" *
                                   "2000-01-05,3.3\n" *
                                   "2000-01-06,4.4\n" *
                                   "2000-01-07,5.5\n"
    end

    @testset "multi column" begin
        ta = TimeArray(ohlc[1:5], values = reshape(1.05:.1:2.95, 5, :))
        io = IOBuffer()
        CSV.write(io, ta)
        @test String(take!(io)) == "timestamp,Open,High,Low,Close\n" *
                                   "2000-01-03,1.05,1.55,2.05,2.55\n" *
                                   "2000-01-04,1.15,1.65,2.15,2.65\n" *
                                   "2000-01-05,1.25,1.75,2.25,2.75\n" *
                                   "2000-01-06,1.35,1.85,2.35,2.85\n" *
                                   "2000-01-07,1.45,1.95,2.45,2.95\n"
    end

    @testset "read csv into TimeArray, single column" begin
        file = "timestamp,Close\n" *
               "2000-01-03,111.94\n" *
               "2000-01-04,102.5\n" *
               "2000-01-05,104\n" *
               "2000-01-06,95\n" *
               "2000-01-07,99.5\n"
        io = IOBuffer(file)
        csv = CSV.File(io)
        ta = TimeArray(csv, timestamp = :timestamp)
        ans = cl[1:5]
        @test timestamp(ta)   == timestamp(ans)
        @test values(ta.Close) == values(ans.Close)
        @test meta(ta)        ≡ csv
    end

    @testset "read csv into TimeArray, multi column" begin
        file = "timestamp,Open,High,Low,Close\n" *
               "2000-01-03,104.88,112.5,101.69,111.94\n" *
               "2000-01-04,108.25,110.62,101.19,102.5\n" *
               "2000-01-05,103.75,110.56,103,104\n" *
               "2000-01-06,106.12,107,95,95\n" *
               "2000-01-07,96.5,101,95.5,99.5\n"
        io = IOBuffer(file)
        csv = CSV.File(io)
        ta = TimeArray(csv, timestamp = :timestamp)
        ans = ohlc[1:5]
        @test timestamp(ta)    == Date(2000, 1, 3):Day(1):Date(2000, 1, 7)
        @test values(ta.Open)  == values(ans.Open)
        @test values(ta.High)  == values(ans.High)
        @test values(ta.Low)   == values(ans.Low)
        @test values(ta.Close) == values(ans.Close)
        @test meta(ta)         ≡ csv
    end

    @testset "issue #442" begin
        file = "date,high,low,open,close,volume\n" *
               "1438992000,50.0,0.00262,50.0,0.00312499,1205.80332085\n" *
               "1439078400,0.0041,0.0024,0.00299999,0.00258069,898.12343401\n" *
               "1439164800,0.0029022,0.0022,0.00264996,0.00264498,718.36526568\n" *
               "1439251200,0.0044,0.002414,0.00264959,0.00395009,3007.27411094\n"
        io = IOBuffer(file)
        csv = CSV.File(io)
        ta = TimeArray(csv, timestamp = :date, timeparser = Date ∘ unix2datetime)
        @test timestamp(ta) == Date(2015, 8, 8):Day(1):Date(2015, 8, 11)
        @test colnames(ta) == [:high, :low, :open, :close, :volume]
    end
end  # @testset "CSV.jl"

end  # @testset "Tables.jl integration


text:
"""
    format_md(text::AbstractString; style::AbstractStyle = DefaultStyle(), kwargs...)

Normalizes the Markdown source and formats Julia code blocks.

See [`format_text`](@ref) for description of formatting options.
"""
function format_md(text::AbstractString; style::AbstractStyle = DefaultStyle(), kwargs...)
    return format_md(text, style; kwargs...)
end

function format_md(text::AbstractString, style::AbstractStyle; kwargs...)
    isempty(text) && return text
    opts = Options(; merge(options(style), kwargs)...)
    return format_md(text, style, opts)
end

function format_md(text::AbstractString, style::AbstractStyle, opts::Options)
    markdown(enable!(
        Parser(),
        [
            AdmonitionRule(),
            FootnoteRule(),
            MathRule(),
            TableRule(),
            FormatRule(style, opts),
        ],
    )(
        text,
    ))
end

ext: jl


text:
"""
    Frames

Stores the actual computed frames and the user input
which can be `:same` or `RFrames(10)`.
The `frames` are computed in [`render`](@ref).
"""
mutable struct Frames{T}
    frames::Union{Nothing,UnitRange}
    user::T
end

Base.convert(::Type{Frames}, x) = Frames(nothing, x)
Base.convert(::Type{Frames}, x::UnitRange) = Frames(x, x)

Base.copy(f::Frames) = Frames(f.frames, f.user)

"""
    set_frames!(parent, elem, last_frames::UnitRange; is_first=false)

Compute the frames based on a.frames and `last_frames`.
Save the result in `a.frames.frames` which can be accessed via [`get_frames`](@ref).

# Arguments
- `parent` is either nothing or the Object for the Action
- `elem` is the Object or Action
- `last_frames` holds the frames of the previous object or action.
- `is_first` defines whether this is the first child of the parent (for actions)
"""
function set_frames!(parent, elem, last_frames::UnitRange; is_first = false)
    frames = elem.frames.user
    elem.frames.frames = get_frames(parent, elem, frames, last_frames; is_first = is_first)
end

"""
    get_frames(elem)

Return `elem.frames.frames` which holds the computed frames for the AbstractObject or AbstractAction `a`.
"""
get_frames(elem) = elem.frames.frames

"""
    get_frames(parent, elem, frames::Symbol, last_frames::UnitRange; is_first=false)

Get the frames based on a symbol (currently only `same`) and the `last_frames`.
Throw `ArgumentError` if symbol is unknown
"""
function get_frames(parent, elem, frames::Symbol, last_frames::UnitRange; is_first = false)
    if frames === :same
        if elem isa AbstractAction && is_first
            return 1:length(last_frames)
        end
        return last_frames
    else
        throw(
            ArgumentError(
                "Currently the only symbol supported for defining frames is `:same`",
            ),
        )
    end
end

"""
    get_frames(parent, elem, relative::RFrames, last_frames::UnitRang; is_first=falsee)

Return the frames based on a relative frames [`RFrames`](@ref) object and the `last_frames`.
"""
function get_frames(
    parent,
    elem,
    relative::RFrames,
    last_frames::UnitRange;
    is_first = false,
)
    start_frame = last(last_frames) + first(relative.frames)
    last_frame = last(last_frames) + last(relative.frames)
    return start_frame:last_frame
end

"""
    get_frames(parent, elem, glob::GFrames, last_frames::UnitRange)

Return the frames based on a global frames [`GFrames`](@ref) object and the `last_frames`.
If `is_action` is false this is the same as defining the frames as just a unit range.
Inside an action it's now defined globally though.
"""
function get_frames(parent, elem, glob::GFrames, last_frames::UnitRange; is_first = false)
    if elem isa AbstractAction
        return glob.frames .- first(get_frames(parent)) .+ 1
    end
    return glob.frames
end


text:
struct Scale{T, N} <: Bijector{N}
    a::T
end

function Scale(a::Union{Real,AbstractArray}; dim::Val{D} = Val(ndims(a))) where D
    return Scale{typeof(a), D}(a)
end

(b::Scale)(x) = b.a .* x
(b::Scale{<:Real})(x::AbstractArray) = b.a .* x
(b::Scale{<:AbstractMatrix})(x::AbstractVecOrMat) = b.a * x
(b::Scale{<:AbstractVector{<:Real}, 2})(x::AbstractMatrix{<:Real}) = b.a .* x

(ib::Inversed{<:Scale})(y) = Scale(inv(ib.orig.a))(y)
(ib::Inversed{<:Scale{<:AbstractVector}})(y) = Scale(inv.(ib.orig.a))(y)
function (ib::Inversed{<:Scale{<:AbstractMatrix, 1}})(y::AbstractVecOrMat)
    return ib.orig.a \ y
end

# inv(b::Scale{T, D}) where {T, D} = Scale(inv(b.a); dim = Val(D))
# inv(b::Scale{<:AbstractVector, D}) where {D} = Scale(inv.(b.a); dim = Val(D))

# We're going to implement custom adjoint for this
logabsdetjac(b::Scale{T, N}, x) where {T, N} = _logabsdetjac_scale(b.a, x, Val(N))

_logabsdetjac_scale(a::Real, x::Real, ::Val{0}) = log(abs(a))
_logabsdetjac_scale(a::Real, x::AbstractVector, ::Val{0}) = fill(log(abs(a)), length(x))
_logabsdetjac_scale(a::Real, x::AbstractVector, ::Val{1}) = log(abs(a)) * length(x)
_logabsdetjac_scale(a::Real, x::AbstractMatrix, ::Val{1}) = fill(log(abs(a)) * size(x, 1), size(x, 2))
_logabsdetjac_scale(a::AbstractVector, x::AbstractVector, ::Val{1}) = sum(log.(abs.(a)))
_logabsdetjac_scale(a::AbstractVector, x::AbstractMatrix, ::Val{1}) = fill(sum(log.(abs.(a))), size(x, 2))
_logabsdetjac_scale(a::AbstractMatrix, x::AbstractVector, ::Val{1}) = log(abs(det(a)))
_logabsdetjac_scale(a::AbstractMatrix, x::AbstractMatrix{T}, ::Val{1}) where {T} = log(abs(det(a))) * ones(T, size(x, 2))


text:
# common directories and files for all tests
const fake_data = "../data/foo.conllu"
const real_tdata = "/state/partition1/tmp/okirnap/UD_English/en-ud-train.conllu"
const chmodel = "../data/english_chmodel.jld" # converted chmodelfile
const real_ddata = "/state/partition1/tmp/okirnap/UD_English/en-ud-dev.conllu"

const sc_tdata = "/ai/data/nlp/conll17/ud-treebanks-v2.0/UD_English/en-ud-train.conllu"
const sc_ddata = "/ai/data/nlp/conll17/ud-treebanks-v2.0/UD_English/en-ud-dev.conllu"
const sc_chmodel = "/Users/okirnap/charNL/trained_models/converged/english_chmodel.jld"

max_stars_repo_name: kirnap/KU-dependencyParser

max_issues_repo_name: kirnap/KU-dependencyParser

max_forks_repo_name: kirnap/KU-dependencyParser


text:
module TestOperators

using Timestamps
using Test

# clstamp = convert(TimestampArray, cl)

# facts("timestamp mathematical operators between two Timestamps") do

#     context("mathematical + works") do
#         @fact (clstamp[1] + clstamp[2]).value     --> 214.44
#         @fact (clstamp[1] + clstamp[2]).timestamp --> Day(1)
#     end

#     context("mathematical - works") do
#         @fact (clstamp[1] - clstamp[3]).value     --> roughly(7.94, atol=.01)
#         @fact (clstamp[1] - clstamp[3]).timestamp --> Day(2)
#     end

#     context("mathematical * works") do
#         @fact (clstamp[1] * clstamp[2]).value     --> 11473.85
#         @fact (clstamp[1] * clstamp[2]).timestamp --> Day(1)
#     end

#     context("mathematical / works") do
#         @fact (clstamp[1] / clstamp[2]).value     --> roughly(1.0920975609756098)
#         @fact (clstamp[1] / clstamp[2]).timestamp --> Day(1)
#     end
# end

# facts("timestamp mathematical operators between one Timestamps and a number") do

#     context("mathematical + works") do
#         @fact (clstamp[1] + 100).value     --> 211.94
#         @fact (clstamp[1] + 100).timestamp --> clstamp[1].timestamp
#     end

#     context("mathematical - works") do
#         @fact (clstamp[1] - 100).value     --> roughly(11.94)
#         @fact (clstamp[1] - 100).timestamp --> clstamp[1].timestamp
#     end

#     context("mathematical * works") do
#         @fact (clstamp[1] * 100).value     --> 11194.0
#         @fact (clstamp[1] * 100).timestamp --> clstamp[1].timestamp
#     end

#     context("mathematical / works") do
#         @fact (clstamp[1] / 100).value     --> 1.1194
#         @fact (clstamp[1] / 100).timestamp --> clstamp[1].timestamp
#     end
# end

# facts("timestamp comparison operators between two Timestamps") do

#     context("comparison > works") do
#         @fact (clstamp[1] > clstamp[2]).value  --> true
#     end

#     context("comparison < works") do
#         @fact (clstamp[1] < clstamp[3]).value  --> false
#     end

# #     context("comparison == works") do
# #         @fact clstamp[1] == clstamp[2].value --> false   # this does work
# #         @fact (clstamp[1] == clstamp[2].value) --> false # this should work
# #     end

#     context("comparison <= works") do
#         @fact (clstamp[1] <= clstamp[2]).value --> false
#     end

#     context("comparison >= works") do
#         @fact (clstamp[1] >= clstamp[2]).value --> true
#     end
# end

# facts("timestamp comparison operators between one Timestamp and a number") do

#     context("comparison > works") do
#         @fact (clstamp[1] > 100).value  --> true
#     end

#     context("comparison < works") do
#         @fact (clstamp[1] < 100).value  --> false
#     end

#     context("comparison == works") do
# #        @fact (clstamp[1] == 100).value --> false
#     end

#     context("comparison <= works") do
#         @fact (clstamp[1] <= 100).value --> false
#     end

#     context("comparison >= works") do
#         @fact (clstamp[1] >= 100).value --> true
#     end
# end

end # module


text:
using Adcv2020
using Test

@testset "Adcv2020.jl" begin
    # Write your tests here.
end


text:
using GaussianRandomFields
using DelimitedFiles
using LinearAlgebra
using DelimitedFiles
using Plots
using NPZ

#-------------------------------------------------
# 	PLOT A SLICE OF THE SAVED ARRAY

function dir__plotslice(dir)
	dims 				=	npzread("$(dir)/dimensions.npy")
	dims 				=	round.(Int64, dims)
	births 				=	npzread("$(dir)/pixel_births.npy")

	println(length(dims))

	if length(dims) == 2
		heatmap(births)
	elseif length(dims) ==3
		heatmap(births[:,:,1])
	end

end

#-------------------------------------------------
# 	FINAL RUN

function dir_pixperaxis_numaxes_covariancetype__save(;
														dir = "/Users/gh10/a/r/ac/re/pu/ar/202006_ffwg_fastforwardalgorithmforgenerators/computation/data/for_paper/cube",
														pixperaxis = 10,
														numaxes = 2,
														covariancetype = "exponential"
													)

	# 	GENERATE DATA
	# 	-------------

	dims 					=	fill(pixperaxis, numaxes)

	# 	DETERMINE COVARIANCE FUNCTION
	if covariancetype == "exponential"
		cov 				=	CovarianceFunction(numaxes, Exponential(.5))
	elseif covariancetype == "anisotropic"
		A 					=	fill( 800, (numaxes, numaxes) ) + 200*I
		cov = CovarianceFunction(numaxes, AnisotropicExponential(A))
	end

	# 	SET PIXEL RANGES
	pts 					=	range(0, stop=1, length=pixperaxis)

	# 	ASSEMBLE GENERATOR
	if numaxes == 2
		grf 				= 	GaussianRandomField(cov, CirculantEmbedding(), pts, pts, minpadding=numaxes*pixperaxis)
	elseif numaxes == 3
		grf 				= 	GaussianRandomField(cov, CirculantEmbedding(), pts, pts, pts, minpadding=numaxes*pixperaxis)
	end

	# 	SAMPLE
	field 					=	sample(grf)

	# 	SAVE DATA
	# 	---------

	# 	CREATE FOLDER
	foldername 				=	"grf_$(covariancetype)$(numaxes)d_$(pixperaxis)pixperaxis"
	dir_save				=	"$(dir)/$(foldername)"
	mkdir(dir_save)

	# 	WRITE CSV
	writedlm("$(dir_save)/dimensions.csv", dims)
	# writedlm("$(dir_save)/pixel_births.csv", arr__rowflattened(field))  # REMEMBER TO FLATTEN DATA

	# 	WRITE NUMPY
	npzwrite("$(dir_save)/pixel_births.npy", field)
	npzwrite("$(dir_save)/dimensions.npy", dims)

end

#-------------------------------------------------
# 	FLATTEN ARRAY IN ROW-MAJOR ORDER

function arr__rowflattened(arr)
	dims 					=	size(arr)
	flattened 				=	zeros(prod(dims))

	counter 				=	0

	if length(dims) == 2
		for i0 = 1:dims[1]
			for i1 = 1:dims[2]
				counter 	=	counter+1
				flattened[counter] 	=	arr[i0, i1]
			end
		end

	elseif length(dims) == 3
		for i0 = 1:dims[1]
			for i1 = 1:dims[2]
				for i2 = 1:dims[3]
					counter 	=	counter+1
					flattened[counter] 	=	arr[i0, i1, i2]
				end
			end
		end
	end

	return flattened
end

#-------------------------------------------------
# 	UNFLATTEN ARRAY IN ROW-MAJOR ORDER

function arr_dims__rowUNflattened(arr, dims)
	unflattened				=	zeros(dims...)

	counter 				=	0

	if length(dims) == 2
		for i0 = 1:dims[1]
			for i1 = 1:dims[2]
				counter 	=	counter+1
				unflattened[i0, i1] 	=	arr[counter]
			end
		end

	elseif length(dims) == 3
		for i0 = 1:dims[1]
			for i1 = 1:dims[2]
				for i2 = 1:dims[3]
					counter 	=	counter+1
					unflattened[i0, i1, i2]		=	arr[counter]
				end
			end
		end
	end

	return unflattened
end


text:
"""
DR <: ProjectionsMethod
Douglas-Rachford method

Reflect with respect one set, than to another and move there halfway: x → (x + R₂R₁x)/2
Stop after fixed number of iterations or if the desired accuracy is achieved
"""
abstract type DR <: ProjectionsMethod end

Base.@kwdef struct DRparam <: DR
    x⁰ = missing
    maxϵ::Union{Float64, Missing} = missing
    maxit::Union{Missing,Int64} = missing
    keephistory::Bool = false
    snapshots::Array{Int64} = Int64[]
end

# DRparam() = DRparam(missing,missing,missing, false, Int64[])

initial(alg::DRparam) = alg.x⁰
tolerance(alg::DRparam) = alg.maxϵ
maxit(alg::DRparam) = alg.maxit
keephistory(alg::DRparam) = alg.keephistory
snapshots(alg::DRparam) = alg.snapshots

function solve(p::TwoSetsFP, alg::DRparam, x⁰, maxϵ, maxit, keephistory::Bool, snapshots::Vector{Int64})
    A = p.A
    B = p.B

    # process the default parameters
    !ismissing(x⁰) || ( x⁰ = getelement(A) )
    !ismissing(maxϵ) || (maxϵ = 1e-15)
    !ismissing(maxit) || (maxit = 100)

    k = 0
    ϵ = Inf

    xᵏ = copy(x⁰)
    xᵏ⁺¹ = similar(xᵏ)
    x̃ᵏ⁺¹ = similar(xᵏ)
    ỹᵏ = similar(xᵏ)
    yᵏ = similar(xᵏ)

    #  preallocated error vector
    err = similar(xᵏ)

    if keephistory
        errhist = Vector{Float64}(undef, maxit)
        disthist = Vector{Float64}(undef, maxit)
    else
        errhist = Float64[]
        disthist = Float64[]
    end

    if length(snapshots) != 0
        xhist = [copy(x⁰) for i in 1:length(snapshots)]
    else
        xhist = typeof(x⁰)[]
    end
    j = 1

    while k < maxit && ϵ > maxϵ

        reflect!(yᵏ, xᵏ, B)
        reflect!(x̃ᵏ⁺¹, yᵏ, A)
        @. xᵏ⁺¹ = (x̃ᵏ⁺¹+ xᵏ) / 2

        err .= xᵏ⁺¹ .- xᵏ # This doesn't say much in infeasible case, but is OK in case of binary aperture
        # dist .= xᵏ⁺¹ .- yᵏ # this calculates true error but can stay large in case of infeasible case
        ϵ = LinearAlgebra.norm(err)
        xᵏ .= xᵏ⁺¹
        k += 1

    #         println(ϵ)
        if keephistory
            errhist[k] = ϵ
            disthist[k] = LinearAlgebra.norm(xᵏ⁺¹ .- yᵏ)
        end

        if k ∈ snapshots
            println("Saving snapshot # $j, iteration # $k")
            xhist[j] .= xᵏ
             j += 1
        end

    end

    println("Using $(supertype(typeof(alg))): to converge with $ϵ accuracy, it took me $k iterations")
    # if keephistory
    #     if length(snapshots) != 0
    #         return xᵏ, errhist, xhist
    #     else
    #         return xᵏ, errhist, Vector{T, (0,)}
    #     end
    # else
    #     return xᵏ, Vector{Float64, (0,)}, Vector{T, (0,)}
    # end
    return xᵏ, (lasty = yᵏ, errhist = errhist, xhist = xhist[1:j - 1], disthist = disthist, k= k)

end

export DR, DRparam


text:
# Apply kinetic operator to wave function in reciprocal space

function op_K( pw::PWGrid, psi::Array{ComplexF64,2} )
    #
    out = zeros(ComplexF64,size(psi))
    Ncol = size(psi)[2]

    Ngwx = pw.gvecw.Ngwx
    G2   = pw.gvec.G2[pw.gvecw.idx_gw2r]

    for is = 1:Ncol
        for ig = 1:Ngwx
            out[ig,is] = psi[ig,is]*G2[ig]
        end
    end

    return 0.5*out # two minus signs
end

function op_K( pw::PWGrid, psi::Array{ComplexF64,1} )
    out = zeros(ComplexF64,size(psi))
    Ngwx = pw.gvecw.Ngwx
    G2   = pw.gvec.G2[pw.gvecw.idx_gw2r]
    for ig = 1:Ngwx
        out[ig] = psi[ig]*G2[ig]
    end
    return 0.5*out # two minus signs
end


text:
# ---
# title: 1448. Count Good Nodes in Binary Tree
# id: problem1448
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/count-good-nodes-in-binary-tree/description/>
# hidden: true
# ---
#
# Given a binary tree `root`, a node _X_ in the tree is named  **good** if in
# the path from root to _X_ there are no nodes with a value _greater than_ X.
#
# Return the number of **good** nodes in the binary tree.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2020/04/02/test_sample_1.png)**
#
#
#
#     Input: root = [3,1,4,3,null,1,5]
#     Output: 4
#     Explanation: Nodes in blue are **good**.
#     Root Node (3) is always a good node.
#     Node 4 -> (3,4) is the maximum value in the path starting from the root.
#     Node 5 -> (3,4,5) is the maximum value in the path
#     Node 3 -> (3,1,3) is the maximum value in the path.
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/04/02/test_sample_2.png)**
#
#
#
#     Input: root = [3,3,null,4,2]
#     Output: 3
#     Explanation: Node 2 -> (3, 3, 2) is not good, because "3" is higher than it.
#
# **Example 3:**
#
#
#
#     Input: root = [1]
#     Output: 1
#     Explanation: Root is considered as **good**.
#
#
#
# **Constraints:**
#
#   * The number of nodes in the binary tree is in the range `[1, 10^5]`.
#   * Each node's value is between `[-10^4, 10^4]`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
struct Entity <: AbstractEntity
    id::Int
end

function Entity(m::AbstractLedger)
	if !isempty(free_entities(m))
		e = pop!(free_entities(m))
		entities(m)[e.id] = e
		return e
	end
	n = length(entities(m)) + 1
	e = Entity(n)
	push!(entities(m), e)
	return e
end

function Entity(m::AbstractLedger, datas::ComponentData...)
	e = Entity(m)
	for d in datas
		m[e] = d
	end
	return e
end

function Entity(c::AbstractComponent, i::Integer)
    return Entity(c.indices.packed[i])
end

function Entity(m::AbstractLedger, parent::Entity, datas::ComponentData...)
	e = Entity(m)
	for d in datas
		m[e] = d
	end
	for (T, component) in components(m)
		if component isa PooledComponent && in(parent, component) && !in(e, component)
			component[e] = parent
		end
	end
	return e
end

Base.iterate(e::Entity, state=1) = state > 1 ? nothing : (e, state+1)

const EMPTY_ENTITY = Entity(0)

Base.:(==)(e1::AbstractEntity, e2::AbstractEntity) = e1.id == e2.id


text:
import Laplacian
import ConjugateGradient

# Test delta
a = [1 2 3; 4 5 6; 7 8 9];
println("Testing `δ`")
@assert Laplacian.δ(i->i,j->j,3) * a == a "identity matrix operational fail"
@assert Laplacian.δ(i->i,j->j,10) == eye(10,10) "identity matrix equality fail"
@assert Laplacian.δ(i->i,j->mod1(j+1,5),5) * [1:5] == [5,1:4] "permute fail"
println("`δ` passed tests.")

# Test opsplit
println("Testing `opsplit`")
@assert Laplacian.opsplit(a, [0,1,0]) == ([1 0 3; 4 0 6; 7 0 9], [0 2 0; 0 5 0; 0 8 0]) "split fail"
println("`opsplit` passed tests")

# Test conjgrad
println("Testing `conjgrad`")
c = [1 2 3; 2 4 5; 3 5 6]
for i in 1:3
    x = rand(3)
    b = c*x
    xx = ConjugateGradient.conjgrad(c,b)
    @assert dot(xx-x, xx-x) < 1e-10 "inaccurate conjugate gradient"
end
println("`conjgrad` passed tests")


text:
@testset "CTMCTran1" begin
    Q = @expr [
        -x x 0;
        y -(y+z) z;
        0 1 -1
    ]
    x = @expr [1.0, 0, 0]
    r = @expr [1.0, 1, 0]
    ma = exrt(LinRange(0, 1, 10), Q, x, r)
    @bind begin
        x = 2.0
        y = 1.0
        z = 5.0
    end
    result = seval(ma)
    println(result)
    # ex = exp(Matrix(seval(Q, test, SymbolicCache())))' * [1,0,0]
    # @test isapprox(ex, result)
end

@testset "CTMCTran2" begin
    Q = @expr [
        -x x 0;
        y -(y+z) z;
        0 1 -1
    ]
    x = @expr [1.0, 0, 0]
    r = @expr [1.0, 1, 0]
    ma = exrt(LinRange(0, 1, 10), Q, x, r)

    @bind begin
        x = 2.0
        y = 1.0
        z = 5.0
    end
    result = seval(ma, :x)
    println(result)
    # ex = exp(Matrix(seval(Q, test, SymbolicCache())))' * [1,0,0]
    # @test isapprox(ex, result)
end

@testset "CTMCTran3" begin
    Q = @expr [
        -x x 0;
        y -(y+z) z;
        0 1 -1
    ]
    x = @expr [1.0, 0, 0]
    r = @expr [1.0, 1, 0]
    ma = exrt(LinRange(0, 1, 10), Q, x, r)

    @bind begin
        x = 2.0
        y = 1.0
        z = 5.0
    end
    result = seval(ma, (:x,:y))
    println(result)
    # ex = exp(Matrix(symboliceval(Q, test, SymbolicCache())))' * [1,0,0]
    # @test isapprox(ex, result)
end

@testset "CTMCTran4" begin
    Q = @expr [
        -x x 0;
        y -(y+z) z;
        0 1 -1
    ]
    x = @expr [1.0, 0, 0]
    r = @expr [1.0, 1, 0]
    ma = exrt(1.0, Q, x, r)

    @bind begin
        x = 2.0
        y = 1.0
        z = 5.0
    end
    result = seval(ma, (:x,:y))
    println(result)
    # ex = exp(Matrix(symboliceval(Q, test, SymbolicCache())))' * [1,0,0]
    # @test isapprox(ex, result)
end

@testset "CTMCTran4" begin
    Q = @expr [
        -x x 0;
        y -(y+z) z;
        0 1 -1
    ]
    x = @expr [1.0, 0, 0]
    r = @expr [1.0, 1, 0]
    ma = cexrt(LinRange(0, 1, 10), Q, x, r)

    @bind begin
        x = 2.0
        y = 1.0
        z = 5.0
    end
    result = seval(ma, (:x,:y))
    println(result)
    # ex = exp(Matrix(symboliceval(Q, test, SymbolicCache())))' * [1,0,0]
    # @test isapprox(ex, result)
end


text:
@testset "JDIS" begin
    mesh = buildSquare(0,1,0,1,2,2)

    sys = System{Float64}(mesh)

    dog = addVariable!(sys, "dog")
    cat = addVariable!(sys, "cat")

    initialize!(sys)

    jdis = JuliaDenseImplicitSolver(sys)

    # Make sure it's not iniitialized by default
    @test !jdis.initialized

    initialize!(jdis)

    # Now it should be
    @test jdis.initialized

    n_dofs = sys.n_dofs

    @test size(jdis.mat) == (n_dofs, n_dofs)
    @test length(jdis.rhs) == n_dofs
    @test length(jdis.solution) == n_dofs

    # Put something in the matrix so we can do a solve
    # Just make it an identity matrix
    # Note: we actually want to fill the matrix with identity values
    #   not replace it with an identity matrix
    for i in 1:n_dofs
        jdis.mat[i,i] = 1
    end

    for i in 1:n_dofs
        jdis.rhs[i] = i
    end

    # Now solve:
    solve!(jdis, assemble=false)

    # Should be the case!
    @test jdis.solution == -jdis.rhs
end

function testFullSolve(diffusion_system)
    u = addVariable!(diffusion_system, "u")

    diffusion_kernel = Diffusion(u)
    addKernel!(diffusion_system, diffusion_kernel)

    left_boundary = DirichletBC(u, [4], 0.0)
    addBC!(diffusion_system, left_boundary)

    right_boundary = DirichletBC(u, [2], 1.0)
    addBC!(diffusion_system, right_boundary)

    initialize!(diffusion_system)

    solver = JuliaDenseImplicitSolver(diffusion_system)

    solve!(solver)

    for i in 1:length(solver.solution)
        @test abs(solver.solution[i] - [0.0, 0.5, 0.5, 0.0, 1.0, 1.0, 0.5, 0.0, 1.0][i]) < 1e-9
    end
end

# Test a 'Full Solve' using JDIS
@testset "JDISFS" begin
    mesh = buildSquare(0, 1, 0, 1, 2, 2)

    diffusion_system = System{Float64}(mesh)

    testFullSolve(diffusion_system)
end

# Test a 'Full Solve' using JDIS and Automatic Differentiation
@testset "JDISFD" begin
    mesh = buildSquare(0, 1, 0, 1, 2, 2)

    diffusion_system = System{Dual{4, Float64}}(mesh)

    testFullSolve(diffusion_system)
end


text:
using CSV, Distributions, HypothesisTests

data = CSV.read("../data/machine1.csv", header=false)[:,1]
n, s, alpha = length(data), std(data), 0.1
ci = (  (n-1)*s^2/quantile(Chisq(n-1),1-alpha/2),
        (n-1)*s^2/quantile(Chisq(n-1),alpha/2)  )

println("Point estimate for the variance: ", s^2)
println("Confidence interval for the variance: ", ci)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Redis

function returnIfExists(client::RedisConnection, problem::Problem, opts::Opts, sd::StepData, expIndices::Vector{Int64},numOutputsFromOutputsFunction; outputLevel::Int64=0)

	longKey = problem.name*":"*problem.lossFunctionString*":"*string(problem.L2reg)*":"sd.stepString*":"*string(opts.stepSizePower)

	outputLevel>1 && println("longKey = $longKey")

	if ~exists(client, longKey*":gnum")
		outputLevel>0 && println("RUN: key does not exist")
		return false
	else
		existingListOfgunum=Int64[]
		existingListRes1=Float64[]
		existingListResEnd=Float64[]
		existingListOfk=Int64[]
		existingListOfOrigWant=Int64[]

		arrayOfStrings  = lrange(client, longKey*":gnum", 0, -1)
		for i in arrayOfStrings
			push!(existingListOfgunum, parse(Int,i))
		end
		arrayOfStrings  = lrange(client, longKey*":1", 0, -1)
		for i in arrayOfStrings
			push!(existingListRes1, parse(Float64,i))
		end

		arrayOfStrings  = lrange(client, longKey*":$numOutputsFromOutputsFunction", 0, -1)
		for i in arrayOfStrings
			push!(existingListResEnd, parse(Float64,i))
		end

		arrayOfStrings  = lrange(client, longKey*":k", 0, -1)
		for i in arrayOfStrings
			push!(existingListOfk, parse(Int,i))
		end
		arrayOfStrings  = lrange(client, longKey*":origWant", 0, -1)
		for i in arrayOfStrings
			push!(existingListOfOrigWant, parse(Int,i))
		end

		# CHECK IF HAVE AN OUTPUT FOR EACH INDEX IN EXPINDICES!!!!!
		i,j = 1,1
		while i<=length(expIndices)
			val = expIndices[i]
			while existingListOfgunum[j]<val
				outputLevel>1 && println("existingListOfgunum[j] = $(existingListOfgunum[j])")
				j+=1
				if j<=length(existingListOfgunum) && (isnan(existingListRes1[j]) || isnan(existingListResEnd[j]))
					outputLevel>0 && println("Found NaN at $j, when existingListOfgunum[$j] = $(existingListOfgunum[j]) and val = $val")
					i = -1
					break
				end
				if j>length(existingListOfgunum)
					outputLevel>0 && println("RE-RUN: could not find element >= $val in existingListOfgunum")
					return false
				end
			end
			if i==-1
				break
			end
			if existingListOfOrigWant[j]<=val
				i+=1
				outputLevel>1 && println("val = $val")
				outputLevel>1 && println("existingListOfOrigWant[j] = $(existingListOfOrigWant[j])")
			else
				outputLevel>0 && println("RE-RUN: for val = $val existingListOfOrigWant[j] is $(existingListOfOrigWant[j]), it is not small enough")
				return false
			end
		end

		i=1
		while Redis.exists(client, longKey*":$i")
			i+=1
		end
		numOutputs = i-1

		kk = zeros(length(arrayOfStrings),0)
		for outputNumber in 1:numOutputs
			arrayss = Float64[]
			fullRange = lrange(client, longKey*":$outputNumber",0,-1)
			for i in fullRange
				push!(arrayss, float(i))
			end
			kk=hcat(kk, arrayss)
		end

		outputLevel>0 && println("NO RE-RUN! Using existing")

		return (
		"Already exists!"
		,
		[parse(Int64,s) for s = lrange(client, longKey*":k", 0,-1)]
		,
		[parse(Int64,s) for s = lrange(client, longKey*":gnum",0, -1)]
		,
		kk
		)
	end
end

function writeFunction(client::RedisConnection, problem::Problem, opts::Opts, sd::StepData,  k::Int64, gnum::Int64, origWant::Int64, fromOutputsFunction; outputLevel::Int64=0)

	longKey = problem.name*":"*problem.lossFunctionString*":"*string(problem.L2reg)*":"sd.stepString*":"*string(opts.stepSizePower)

	if ~exists(client, longKey*":k")
		outputLevel>0 && println("Key doesn't exist, so pushing new")
		rpush(client, longKey*":k",k)
		rpush(client, longKey*":gnum",gnum)
		rpush(client, longKey*":origWant",origWant)
		for j in 1:length(fromOutputsFunction.resultLine)
			rpush(client, longKey*":$j",fromOutputsFunction.resultLine[j])
		end
	else
		arrayOfStrings  = lrange(client, longKey*":gnum", 0, -1)
		currentIndex=0
		for i in 1:length(arrayOfStrings)
			if gnum==parse(Int,arrayOfStrings[i])
				outputLevel>0 && println("Not writing because exists exact same with gnum = $gnum, at i = $i")
				ah = parse(Int,lindex(client,longKey*":origWant",i-1))
				if origWant < ah
					outputLevel>0 && println("origWant < ah, $origWant < $ah, so replacing!")
					lset(client, longKey*":origWant", i-1, origWant)
				end
				return
			end
			if parse(Int,arrayOfStrings[i])>gnum
				outputLevel>0 && println("Found where to insert: before $(i-1)")
				currentIndex=i-1
				break
			end
		end
		if currentIndex ==0
			outputLevel>0 && println("Didn't find where to insert, so attaching at the end")
			rpush(client, longKey*":k",k)
			rpush(client, longKey*":gnum",gnum)
			rpush(client, longKey*":origWant",origWant)
			for j in 1:length(fromOutputsFunction.resultLine)
				rpush(client, longKey*":$j",fromOutputsFunction.resultLine[j])
			end
		else
			ah = lindex(client,longKey*":k",currentIndex  )
			lset(client, longKey*":k", currentIndex, "oh well")
			linsert(client, longKey*":k" ,"after", "oh well" ,ah)
			linsert(client, longKey*":k" ,"after", "oh well" ,k)
			lrem(client, longKey*":k" ,0, "oh well")

			ah = lindex(client,longKey*":gnum",currentIndex  )
			lset(client, longKey*":gnum", currentIndex, "oh well")
			linsert(client, longKey*":gnum" ,"after", "oh well" ,ah)
			linsert(client, longKey*":gnum" ,"after", "oh well" ,gnum)
			lrem(client, longKey*":gnum" ,0, "oh well")

			ah = lindex(client,longKey*":origWant",currentIndex  )
			lset(client, longKey*":origWant", currentIndex, "oh well")
			linsert(client, longKey*":origWant" ,"after", "oh well" ,ah)
			linsert(client, longKey*":origWant" ,"after", "oh well" ,origWant)
			lrem(client, longKey*":origWant" ,0, "oh well")

			for j in 1:length(fromOutputsFunction.resultLine)
				ah = lindex(client,longKey*":$j",currentIndex  )
				lset(client, longKey*":$j", currentIndex, "oh well")
				linsert(client, longKey*":$j" ,"after", "oh well" ,ah)
				linsert(client, longKey*":$j" ,"after", "oh well" ,fromOutputsFunction.resultLine[j])
				lrem(client, longKey*":$j" ,0, "oh well")
			end
		end
	end
end


text:
import Base: download, nameof
import HTTP

abstract type AbstractBackend end

struct CURL <: AbstractBackend end
nameof(::CURL) = "cURL"

function download(::CURL, url, filename; verbose::Bool=false)
    curl = Sys.which("curl")
    curl === nothing && error("The `curl` executable was not found.")
    try
        if verbose
            run(`$curl -o $filename -L $url`)
        else
            run(`$curl -s -o $filename -L $url`)
        end
    catch err
        if (isdefined(Base, :ProcessFailedException) &&
            err isa ProcessFailedException) || err isa ErrorException
            throw(DownloadError(sprint(showerror, err)))
        else
            rethrow(err)
        end
    end
end

struct Wget <: AbstractBackend end
nameof(::Wget) = "wget"

function download(::Wget, url, filename; verbose::Bool=false)
    wget = Sys.which("wget")
    wget === nothing && error("The `wget` executable was not found.")
    try
        if verbose
            run(`$wget -O $filename $url`)
        else
            run(`$wget -q -O $filename $url`)
        end
    catch err
        if (isdefined(Base, :ProcessFailedException) &&
            err isa ProcessFailedException) || err isa ErrorException
            throw(DownloadError(sprint(showerror, err)))
        else
            rethrow(err)
        end
    end
end

struct Http <: AbstractBackend end
nameof(::Http) = "HTTP.jl"

const HEADERS = ["User-Agent" => "RemoteFiles.jl/0.3 (+https://github.com/helgee/RemoteFiles.jl)", "Accept" => "*/*"]

function download(::Http, url, filename; verbose::Bool=false)
    try
        if verbose
            HTTP.download(url, filename, HEADERS)
        else
            HTTP.download(url, filename, HEADERS; update_period=typemax(Int))
        end
    catch err
        throw(DownloadError((sprint(showerror, err))))
    end
end


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

@userplot DistPlot2D

@recipe function f(dp::DistPlot2D; quantiles=[0.25,0.50,0.75],
                                   bandwidthx=100, bandwidthy=100)
  # retrieve inputs
  sdata = dp.args[1]
  var₁  = dp.args[2]
  var₂  = dp.args[3]

  x = sdata[var₁]
  y = sdata[var₂]

  # fit average shifted histogram
  h = ash(x, y, mx=bandwidthx, my=bandwidthy)
  hx, hy, hz = h.rngx, h.rngy, h.z
  ls = quantile(vec(hz), quantiles)

  # 2D mean vector
  μ = mean([x y], dims=1)

  # plot scatter of points
  @series begin
    seriestype --> :scatter
    primary --> false
    seriescolor --> :black
    seriesalpha --> 0.5
    markersize --> 1
    x, y
  end

  # plot 2D mean
  @series begin
    primary --> false
    seriestype --> :vline
    seriescolor --> :green
    [μ[1]]
  end
  @series begin
    primary --> false
    seriestype --> :hline
    seriescolor --> :green
    [μ[2]]
  end
  @series begin
    primary --> false
    seriestype --> :scatter
    seriescolor --> :green
    marker --> :square
    markersize --> 4
    [μ[1]], [μ[2]]
  end

  seriestype --> :contour
  seriescolor --> :black
  aspect_ratio --> :equal
  framestyle --> :box
  legend --> false
  grid --> false
  levels --> ls
  xguide --> var₁
  yguide --> var₂

  hx, hy, hz
end


text:
# # [Implicit Nonlinear Dynamics : Cartpole](@id cartpole)
#
# The following is another example on how to use the [`ImplicitOptimizer`](@ref) that is taken from the [original paper](https://royalsocietypublishing.org/doi/10.1098/rspa.2020.0279).
# As always, we start by creating a corresponding dataset:

using DataDrivenDiffEq
using ModelingToolkit
using OrdinaryDiffEq
using LinearAlgebra
#md using Plots
#md gr()

function cart_pole(u, p, t)
    du = similar(u)
    F = -0.2 + 0.5*sin(6*t) # the input
    du[1] = u[3]
    du[2] = u[4]
    du[3] = -(19.62*sin(u[1])+sin(u[1])*cos(u[1])*u[3]^2+F*cos(u[1]))/(2-cos(u[1])^2)
    du[4] = -(sin(u[1])*u[3]^2 + 9.81*sin(u[1])*cos(u[1])+F)/(2-cos(u[1])^2)
    return du
end

u0 = [0.3; 0; 1.0; 0]
tspan = (0.0, 5.0)
dt = 0.05
cart_pole_prob = ODEProblem(cart_pole, u0, tspan)
solution = solve(cart_pole_prob, Tsit5(), saveat = dt)

X = solution[:,:]
DX = similar(X)
for (i, xi) in enumerate(eachcol(X))
    DX[:, i] = cart_pole(xi, [], solution.t[i])
end
t = solution.t

ddprob = ContinuousDataDrivenProblem(
    X , t, DX = DX[3:4, :], U = (u,p,t) -> [-0.2 + 0.5*sin(6*t)]
)

#md plot(ddprob)

# Note that we just included the third and forth time derivative, assuming that we already know that the velocity `x[3:4]` is equal to the time
# derivative of the position `x[1:2]`.
# Next, we define a sufficient [`Basis`](@ref). Again, we need to include `implicits` in the definition of
# our candidate functions and inform the [`Basis`](@ref) of it.

@parameters t
@variables u[1:4] du[1:2] x[1:1]
u, du, x = map(collect, [u, du, x])

polys = polynomial_basis(u, 2)
push!(polys, sin.(u[1]))
push!(polys, cos.(u[1]))
push!(polys, sin.(u[1])^2)
push!(polys, cos.(u[1])^2)
push!(polys, sin.(u[1]).*u[3:4]...)
push!(polys, sin.(u[1]).*u[3:4].^2...)
push!(polys, sin.(u[1]).*cos.(u[1])...)
push!(polys, sin.(u[1]).*cos.(u[1]).*u[3:4]...)
push!(polys, sin.(u[1]).*cos.(u[1]).*u[3:4].^2...)

implicits = [du;  du[1] .* u; du[2] .* u; du .* cos(u[1]);   du .* cos(u[1])^2; polys]
push!(implicits, x...)
push!(implicits, x[1]*cos(u[1]))
push!(implicits, x[1]*sin(u[1]))

basis= Basis(implicits, u, implicits = du, controls = x,  iv = t);
println(basis) # hide

# We solve the problem by varying over a sufficient set of thresholds for the associated optimizer.
# Additionally we activate the `scale_coefficients` option for the [`ImplicitOptimizer`](@ref), which helps to find sparse equations by normalizing the resulting coefficient matrix after each suboptimization.
# To evaluate the pareto optimal solution, we use the functions `f` and `g` which can be passed as keyword arguments into the `solve` function. `f` is a function with different signatures for different optimizers, but returns the ``L_0`` norm of the coefficients and the ``L_2`` error of the current model. `g` takes this vector and projects it down onto a scalar, using the ``L_2`` norm per default. However, here we want to use the `AIC`  of the output of `f`. A noteworthy exception is of course, that we want only results with two or more active coefficents. Hence, we modify `g` accordingly.

λ = [1e-4;5e-4;1e-3;2e-3;3e-3;4e-3;5e-3;6e-3;7e-3;8e-3;9e-3;1e-2;2e-2;3e-2;4e-2;5e-2]

opt = ImplicitOptimizer(λ)

g(x) = x[1] <= 1 ? Inf : 2*x[1]-2*log(x[2])

res = solve(ddprob, basis, opt, du, maxiter = 1000, g = g, show_progress = true)
system = result(res)
println(system) #hide

# We have recovered the correct equations of motion!
# Another visual check using the problem and the result yields

#md plot(
#md     plot(ddprob), plot(res), layout = (1,2)
#md )

#md # ## [Copy-Pasteable Code](@id cartpole_copy_paste)
#md #
#md # ```julia
#md # @__CODE__
#md # ```

## Test #src
for r_ in [res] #src
    @test all(l2error(r_) .< 0.5) #src
    @test all(aic(r_) .> 1e3) #src
    @test all(determination(r_) .>= 0.9) #src
end #src


text:
#!/usr/bin/env julia
using BinDeps2
using SHA

basedir = pwd()

# Our build products will go into ./products
out_path = joinpath(pwd(), "products")
rm(out_path; force=true, recursive=true)
mkpath(out_path)

products = Dict()

temp_prefix() do prefix
    src_path = joinpath(prefix, "src")
    build_path = joinpath(prefix, "build")
    try mkpath(src_path) end
    try mkpath(build_path) end

    # First, download the sources, store them into /src
    src_url = "https://computation.llnl.gov/projects/sundials/download/sundials-2.7.0.tar.gz"
    src_hash = "d39fcac7175d701398e4eb209f7e92a5b30a78358d4a0c0fcc23db23c11ba104"

    download_verify_unpack(src_url, src_hash, src_path; verbose=true)

    # Build for many platforms
    for platform in supported_platforms()
        cd(build_path) do
            mkdir(string(platform))
	    cd(string(platform)) do
                target = platform_triplet(platform)
                libsundials = LibraryResult(joinpath(libdir(prefix), target, "libsundials"))
                extras = platform in (:win64,:win32) ? `-DCMAKE_SYSTEM_NAME=Windows` : ``
                cmake_args = `-DCMAKE_INSTALL_PREFIX=/inst/$platform -DEXAMPLES_ENABLE=OFF $extras`
                steps = [`cmake $cmake_args ../../src/sundials-2.7.0/`,
                         `make -j4`,
                         `make install`]
                dep = Dependency("sundials", [libsundials], steps, platform, prefix)
                build(dep; verbose=true, force=true)
                rm("$basedir/libsundials_$(target).tar.gz"; force=true)
                tarball_path = package(BinDeps2.Prefix("$(prefix.path)/inst/$platform"),
			               "$basedir/products/libsundials",
                                       platform=platform, verbose=true)
                # Once we're built up, go ahead and package this prefix out
                tarball_hash = open(tarball_path, "r") do f
                    return bytes2hex(sha256(f))
                end
                products[target] = (basename(tarball_path), tarball_hash)
                info("Built and saved at $(tarball_path)")
            end
        end
    end
end

# In the end, dump an informative message telling the user how to download/install these
info("Hash/filename pairings:")
for target in keys(products)
    filename, hash = products[target]
    println("    \"$(target)\" => (\"\$prefix/$(filename)\", \"$(hash)\"),")
end


text:
using Plots
using Statistics
using LaTeXStrings
function std_calculater(arr,stdList)
    hight=[]
    for col in eachcol(arr)
        count=0
        for element in col
            if element!=0
                count+=1
            end
        end
        push!(hight, count)
    end
    std_num=std(hight)
    push!(stdList, std_num)
    return stdList
end
function BoundaryCondition(i)
    if i==L+1
        return 1
    elseif i==0
        return L
    else
        return i
    end
end
function height_cal(i,arr) #calculating the height of the layer in a column
    i=BoundaryCondition(i)
    row=L #starting from the top layer
    height=L
    for j in 1:L
        if arr[row,i]==0.0 #when the entry is empty
            height-=1 #decreasing height by one
            row-=1
        else
            break #exit the loop for the next column
        end
    end
    return height
end
function neighbor_checking(arr, col)
    h=[]
    for i in (col-1):1:(col+1)
        height= height_cal(i,arr)
        if i==col
            push!(h, (height+1))
        else
            push!(h, height)
        end
    end
    return maximum(h)
end
L=200
arr=zeros((L,L))
color=1
stdList=[]
t_interval=ceil.(Int, exp.(1:0.5:14))
for i in 1:26
    falling_n=t_interval[i+1]-t_interval[i]
    for each in 1:falling_n
        finalP=t_interval[18]
        col=rand(1:L)
        height=neighbor_checking(arr, col)
        if height>=L
            finalP=each
            break
        end
        arr[height, col]=color
        if i%(10*200*color)==0
            color+=1
        end
    end
    stdList=std_calculater(arr,stdList)
end
BigStdList=[]
run_number=100
for num in 1:run_number
    L=200
    arr=zeros((L,L))
    color=1
    stdList=[]
    time=ceil.(Int, exp.(0:0.5:14))
    for i in 1:26
        falling_n=time[i+1]-time[i]
        for each in 1:falling_n
            finalP=t_interval[18]
            col=rand(1:L)
            height=neighbor_checking(arr, col)
            if height>=L
                finalP=each
                break
            end
            arr[height, col]=color
            if i%(10*200*color)==0
                color+=1
            end
        end
        stdList=std_calculater(arr,stdList)
    end
    push!(BigStdList, stdList)
end
BigStdList
BiggerMeanList=[]
BiggerStdList=[]
for i in 1:26
    tempvalue=[]
    for j in 1:100
        push!(tempvalue,BigStdList[j][i])
    end
    push!(BiggerStdList, std(tempvalue))
    push!(BiggerMeanList, mean(tempvalue))
end
function LineFit(Time, Mean)
    a=[hcat(log.(Time)) reshape(ones(19), 19, 1)]
    b=hcat(log.(Mean))
    line=(a\b)
    x = 1:10
    y = x .* line[1] .+ line[2]
    return x, y, line
end
X,Y,Line=LineFit(t_interval[1:19], BiggerMeanList[1:19])
plot(log.(t_interval)[1:26], log.(BiggerMeanList)[1:26],yerr=BiggerStdList, xlabel="time", ylabel="w(t)")
plot!(X,Y,c= :black,label = L"y = %$(round(Line[1],digits= 2))x + %$(round(Line[2],digits= 2))", legend=:bottomright)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter3\\Fig\\3.3_3.png")


text:
module Day02
export part1, part2

using AdventOfCode.Computer

function part1()
    p = read_program("input02.dat")
    p.mem[2] = 12
    p.mem[3] = 02
    run_program!(p)
    p.mem[1]
end

function part2(output)
    p_orig = read_program("input02.dat")
    p = Program(similar(p_orig.mem.vect))

    for i in 0:99
        for j in 0:99
            reset!(p, p_orig)

            p.mem[2] = i
            p.mem[3] = j
            step_program!(p)
            if p.mem[1] == output
                return 100*i+j
            end
        end
    end
end

end


text:
using PMDPs

using StaticArrays
using POMDPs
using MCTS, DiscreteValueIteration
using StatsBase
using Plots
using POMDPSimulators
using D3Trees
using POMDPPolicies
using POMDPLinter
using Random
using DataFrames
using POMDPSimulators

using Traceur
using XLSX
using BSON

mdp = PMDPs.create_PMDP(PMDPe)

# POMDPLinter.@requirements_info SparseValueIterationSolver() mdp
# @requirements_info ValueIterationSolver() mdp

function value_heuristic(s::State, m::PMDP)
    α_t = 1.0/m.T
    α_c = 15.
    α_t*(m.T-s.t)*α_c*sum(s.c)
end
init_util = [value_heuristic(mdp.states[i], mdp) for i in 1:length(mdp.states)]

solver = SparseValueIterationSolver(max_iterations=100, belres=1e-6, verbose=true)#, init_util=init_util) # creates the solver
# POMDPs.@show_requirements POMDPs.solve(solver, mdp)

println("Solving...")
policy = POMDPs.solve(solver, mdp)
println("Done.")

# bson("policy.bson", Dict(:policy => policy))
# BSON.load("policy.bson")

# Get action counts
# df = DataFrame(pricei = policy.policy)
# df[:, :price] = [policy.action_map[i] for i in df.pricei]
# combine(groupby(df, :price), nrow)

qmat = policy.qmat
non_zero_Q_states = (sum(qmat, dims = 2).>0)[:,1]
qmat[non_zero_Q_states, :]

qdf = DataFrame(hcat( [repr(mdp.states[i].c) for i in 1:length(mdp.states)],
                      [mdp.states[i].t for i in 1:length(mdp.states)],
                      [repr(mdp.states[i].p) for i in 1:length(mdp.states)],policy.qmat),
                map(Symbol, vcat(["c", "t", "p"], policy.action_map)))

rm("q_mat.xlsx", force=true)
XLSX.writetable("q_mat.xlsx", qdf)

hr = HistoryRecorder(max_steps=100, capture_exception=true, rng=MersenneTwister(1234))
h = simulate(hr, mdp, policy)
collect(eachstep(h, "s, a, r"))
sum(h[:r])


text:
# Run tests on the fiber problem in NLEVP (bessel function nonlinearity)

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

push!(LOAD_PATH, string(@__DIR__, "/../src/gallery_extra"))
using GalleryNLEVP

@testset "NLEVP fiber" begin
    nep_org=nep_gallery(NLEVP_NEP,"fiber");
    n=size(nep_org,1);

    # An exact eigenvalue according (reported in NLEVP collection)
    sol_val= 7.139494306065948e-07;

    @bench @testset "basics" begin
        # Checking Mder
        λ=3;
        δ=0.0001;
        Ap=compute_Mder(nep_org,λ+δ)
        Am=compute_Mder(nep_org,λ-δ)
        Mp_approx=(Ap-Am)/(2δ);
        Mp=compute_Mder(nep_org,λ,1)
        @test opnorm(Mp-Mp_approx,1)<δ^2*100

        # Checking Mlincomb
        λ=2;
        x=randn(n)
        z_plus=compute_Mlincomb(nep_org,λ+δ,x)
        z_minus=compute_Mlincomb(nep_org,λ-δ,x)
        zp=compute_Mlincomb(nep_org,λ,x,[1],1)
        zp-(z_plus-z_minus)/(2δ)
        @test norm(zp-(z_plus-z_minus)/(2δ))<(δ^2*100)

        @info "Running Newton"
        (λstar,v)=newton(Float64,nep_org,λ=7e-7,v=ones(n),logger=displaylevel,
                         tol=1e-14)

        @test abs(sol_val-λstar)/abs(λ) < 1e-10
        if imag(λstar) != 0
            @warn "Newton switches to complex although it should be real: $λstar"
        end

        @info "Running quasi-newton w armijo (with eigval error as termination)"
        function myerrmeasure(λt,v)
            return abs((λt-λstar)/λstar)
        end

        tol=1e-11
        (λ,v)=quasinewton(nep_org,λ=7.1e-7,v=ones(n),
                          logger=displaylevel,errmeasure=myerrmeasure,
                          tol=tol,armijo_factor=0.5,armijo_max=10)

        @test abs(sol_val-λ)/abs(λ) < tol
        if imag(λ) != 0
            @warn "Quasinewton switches to complex although it should be real: $λ"
        end

        @info "Running resinv"
        tol=1e-8
        (λ,v)=resinv(nep_org,λ=7e-7,v=ones(n),
                     logger=displaylevel,errmeasure=myerrmeasure,
                     tol=tol)
        @test abs(sol_val-λ)/abs(λ) < tol
        if imag(λ) != 0
            @warn "resinv switches to complex although it should be real: $λ"
        end

        @info "Running MSLP"

        (λ,v)=mslp(Float64,nep_org,λ=7e-7,
                   logger=displaylevel,errmeasure=myerrmeasure,
                   tol=tol)
        @test abs(sol_val-λ)/abs(λ) < tol
        if imag(λ) != 0
            @warn "mslp switches to complex although it should be real: $λ"
        end
    end

    @bench @testset "PEP interpolation" begin
        # Create a PEP Which interpolates in a couple of points

        intpoints=[1e-8,1e-7,7e-7,9e-7]
        @info "Interpolating: $intpoints"
        pep=interpolate(nep_org,intpoints);
        @info "Running IAR"
        λ,v=iar(pep,σ=7e-7,maxit=100,logger=displaylevel,neigs=2)
        minerr1=minimum(abs.(sol_val.-λ)) ./ abs(sol_val)
        @info "Error: $minerr1"
        @test minerr1<1e-4

        # Create a new PEP which interpolates in the solution computed
        # in the previous iteration
        intpoints=vcat(intpoints,real(λ[1]))
        @info "Interpolating: $intpoints"
        pep=interpolate(nep_org,intpoints);
        @info "Running IAR"
        λ,v=iar(pep,σ=7e-7,maxit=100,logger=displaylevel,neigs=2)
        minerr2=minimum(abs.(sol_val.-λ)) ./ abs(sol_val)
        @info "Error: $minerr2"
        @test minerr2<minerr1
    end

end


text:
function matches(signal, chunk)
    l1 = length(signal)
    l2 = length(chunk)
    l = l1 - l2
    l < 0 && return Tuple{Int,Float64}[]
    res = zeros(l)
    for i=1:l+1
        s = @view signal[i:i+l2-1]
        res[i] = sum(ref.*s)
    end
end


text:
# ------------------------------------------------------------
# IsParallel -> Algorithm is run in parallel
# ------------------------------------------------------------
@define_trait IsParallel

@define_traitfn IsParallel init_subproblems!(lshaped::AbstractLShapedSolver{T,A,M,S},subsolver::AbstractMathProgSolver) where {T <: Real, A <: AbstractVector, M <: LQSolver, S <: LQSolver} = begin
    function init_subproblems!(lshaped::AbstractLShapedSolver{T,A,M,S},subsolver::AbstractMathProgSolver,!IsParallel) where {T <: Real, A <: AbstractVector, M <: LQSolver, S <: LQSolver}
        # Prepare the subproblems
        m = lshaped.structuredmodel
        load_subproblems!(lshaped,scenarioproblems(m),subsolver)
        return lshaped
    end

    function init_subproblems!(lshaped::AbstractLShapedSolver{T,A,M,S},subsolver::AbstractMathProgSolver,IsParallel) where {T <: Real, A <: AbstractVector, M <: LQSolver, S <: LQSolver}
        @unpack κ = lshaped.parameters
        # Partitioning
        (jobsize,extra) = divrem(lshaped.nscenarios,nworkers())
        # One extra to guarantee coverage
        if extra > 0
            jobsize += 1
        end
        # Load initial decision
        put!(lshaped.decisions,1,lshaped.x)
        # Create subproblems on worker processes
        m = lshaped.structuredmodel
        start = 1
        stop = jobsize
        active_workers = Vector{Future}(nworkers())
        for w in workers()
            lshaped.work[w-1] = RemoteChannel(() -> Channel{Int}(round(Int,10/κ)), w)
            put!(lshaped.work[w-1],1)
            lshaped.subworkers[w-1] = RemoteChannel(() -> Channel{Vector{SubProblem{T,A,S}}}(1), w)
            active_workers[w-1] = load_worker!(scenarioproblems(m),w,lshaped.subworkers[w-1],lshaped.x,start,stop,subsolver)
            if start > lshaped.nscenarios
                continue
            end
            start += jobsize
            stop += jobsize
            stop = min(stop,lshaped.nscenarios)
        end
        map(wait,active_workers)
        lshaped
    end
end

@define_traitfn IsParallel iterate!(lshaped::AbstractLShapedSolver) = begin
    function iterate!(lshaped::AbstractLShapedSolver,!IsParallel)
        iterate_nominal!(lshaped)
    end

    function iterate!(lshaped::AbstractLShapedSolver,IsParallel)
        iterate_parallel!(lshaped)
    end
end

@define_traitfn IsParallel init_workers!(lshaped::AbstractLShapedSolver) = begin
    function init_workers!(lshaped::AbstractLShapedSolver,IsParallel)
        active_workers = Vector{Future}(nworkers())
        for w in workers()
            active_workers[w-1] = remotecall(work_on_subproblems!,
                                             w,
                                             lshaped.subworkers[w-1],
                                             lshaped.work[w-1],
                                             lshaped.cutqueue,
                                             lshaped.decisions)
        end
        return active_workers
    end
end

@define_traitfn IsParallel close_workers!(lshaped::AbstractLShapedSolver,workers::Vector{Future}) = begin
    function close_workers!(lshaped::AbstractLShapedSolver,workers::Vector{Future},IsParallel)
        @async begin
            close(lshaped.cutqueue)
            map(wait,workers)
        end
    end
end

@define_traitfn IsParallel calculate_objective_value(lshaped::AbstractLShapedSolver,x::AbstractVector) = begin
    function calculate_objective_value(lshaped::AbstractLShapedSolver,x::AbstractVector,NullTrait)
        return lshaped.c⋅x + sum([subproblem.π*subproblem(x) for subproblem in lshaped.subproblems])
    end

    function calculate_objective_value(lshaped::AbstractLShapedSolver,x::AbstractVector,IsParallel)
        return lshaped.c⋅x + sum(fetch.([@spawnat w+1 calculate_subobjective(worker,x) for (w,worker) in enumerate(lshaped.subworkers)]))
    end
end

@define_traitfn IsParallel fill_submodels!(lshaped::AbstractLShapedSolver,scenarioproblems::StochasticPrograms.ScenarioProblems) = begin
    function fill_submodels!(lshaped::AbstractLShapedSolver,scenarioproblems::StochasticPrograms.ScenarioProblems,NullTrait)
        for (i,submodel) in enumerate(scenarioproblems.problems)
            fill_submodel!(submodel,lshaped.subproblems[i])
        end
    end

    function fill_submodels!(lshaped::AbstractLShapedSolver,scenarioproblems::StochasticPrograms.ScenarioProblems,IsParallel)
        j = 0
        for w = 1:length(lshaped.subworkers)
            n = remotecall_fetch((sw)->length(fetch(sw)),w+1,lshaped.subworkers[w])
            for i = 1:n
                fill_submodel!(scenarioproblems.problems[i+j],remotecall_fetch((sw,i)->get_solution(fetch(sw)[i]),w+1,lshaped.subworkers[w],i)...)
            end
            j += n
        end
    end
end

@define_traitfn IsParallel fill_submodels!(lshaped::AbstractLShapedSolver,scenarioproblems::StochasticPrograms.DScenarioProblems) = begin
    function fill_submodels!(lshaped::AbstractLShapedSolver,scenarioproblems::StochasticPrograms.DScenarioProblems,NullTrait)
        active_workers = Vector{Future}(length(scenarioproblems))
        j = 1
        for w = 1:length(scenarioproblems)
            n = remotecall_fetch((sp)->length(fetch(sp).problems),w+1,scenarioproblems[w])
            active_workers[w] = remotecall((subproblems,sp) -> begin
                                               scenarioproblems = fetch(sp)
                                               for (i,submodel) in enumerate(scenarioproblems.problems)
                                                 fill_submodel!(submodel,subproblems[i])
                                               end
                                             end,
                                             w+1,
                                             lshaped.subproblems[j:n],
                                             scenarioproblems[w])
            j += n
        end
    end

    function fill_submodels!(lshaped::AbstractLShapedSolver,scenarioproblems::StochasticPrograms.DScenarioProblems,IsParallel)
        active_workers = Vector{Future}(length(scenarioproblems))
        for w = 1:length(scenarioproblems)
            active_workers[w] = remotecall(fill_submodels!,
                                             w+1,
                                             lshaped.subworkers[w],
                                             scenarioproblems[w])
        end
        map(wait,active_workers)
    end
end

# Parallel routines #
# ======================================================================== #
mutable struct DecisionChannel{A <: AbstractArray} <: AbstractChannel
    decisions::Dict{Int,A}
    cond_take::Condition
    DecisionChannel(decisions::Dict{Int,A}) where A <: AbstractArray = new{A}(decisions, Condition())
end

function put!(channel::DecisionChannel, t, x)
    channel.decisions[t] = x
    notify(channel.cond_take)
    return channel
end

function take!(channel::DecisionChannel, t)
    x = fetch(channel,t)
    delete!(channel.decisions, t)
    return x
end

isready(channel::DecisionChannel) = length(channel.decisions) > 1
isready(channel::DecisionChannel, t) = haskey(channel.decisions,t)

function fetch(channel::DecisionChannel, t)
    wait(channel,t)
    return channel.decisions[t]
end

function wait(channel::DecisionChannel, t)
    while !isready(channel, t)
        wait(channel.cond_take)
    end
end

SubWorker{T,A,S} = RemoteChannel{Channel{Vector{SubProblem{T,A,S}}}}
ScenarioProblems{D,SD,S} = RemoteChannel{Channel{StochasticPrograms.ScenarioProblems{D,SD,S}}}
Work = RemoteChannel{Channel{Int}}
Decisions{A} = RemoteChannel{DecisionChannel{A}}
QCut{T} = Tuple{Int,T,SparseHyperPlane{T}}
CutQueue{T} = RemoteChannel{Channel{QCut{T}}}

function load_subproblems!(lshaped::AbstractLShapedSolver{T,A},scenarioproblems::StochasticPrograms.ScenarioProblems,subsolver::AbstractMathProgSolver) where {T <: Real, A <: AbstractVector}
    for i = 1:lshaped.nscenarios
        m = subproblem(scenarioproblems,i)
        y₀ = convert(A,rand(m.numCols))
        push!(lshaped.subproblems,SubProblem(m,
                                             parentmodel(scenarioproblems),
                                             i,
                                             probability(scenario(scenarioproblems,i)),
                                             copy(lshaped.x),
                                             y₀,
                                             subsolver))
    end
    return lshaped
end

function load_subproblems!(lshaped::AbstractLShapedSolver{T,A},scenarioproblems::StochasticPrograms.DScenarioProblems,subsolver::AbstractMathProgSolver) where {T <: Real, A <: AbstractVector}
    for i = 1:lshaped.nscenarios
        m = subproblem(scenarioproblems,i)
        y₀ = convert(A,rand(m.numCols))
        push!(lshaped.subproblems,SubProblem(m,
                                             i,
                                             probability(scenario(scenarioproblems,i)),
                                             copy(lshaped.x),
                                             y₀,
                                             masterterms(scenarioproblems,i),
                                             subsolver))
    end
    return lshaped
end

function load_worker!(sp::StochasticPrograms.ScenarioProblems,
                      w::Integer,
                      worker::SubWorker,
                      x::AbstractVector,
                      start::Integer,
                      stop::Integer,
                      subsolver::AbstractMathProgSolver)
    problems = [sp.problems[i] for i = start:stop]
    πs = [probability(sp.scenariodata[i]) for i = start:stop]
    return remotecall(init_subworker!,
                      w,
                      worker,
                      sp.parent,
                      problems,
                      πs,
                      x,
                      subsolver,
                      collect(start:stop))
end

function load_worker!(sp::StochasticPrograms.DScenarioProblems,
                      w::Integer,
                      worker::SubWorker,
                      x::AbstractVector,
                      start::Integer,
                      stop::Integer,
                      subsolver::AbstractMathProgSolver)
    return remotecall(init_subworker!,
                      w,
                      worker,
                      sp[w-1],
                      x,
                      subsolver,
                      collect(start:stop))
end

function init_subworker!(subworker::SubWorker{T,A,S},
                         parent::JuMP.Model,
                         submodels::Vector{JuMP.Model},
                         πs::A,
                         x::A,
                         subsolver::AbstractMathProgSolver,
                         ids::Vector{Int}) where {T <: Real, A <: AbstractArray, S <: LQSolver}
    subproblems = Vector{SubProblem{T,A,S}}(length(ids))
    for (i,id) = enumerate(ids)
        y₀ = convert(A,rand(submodels[i].numCols))
        subproblems[i] = SubProblem(submodels[i],parent,id,πs[i],x,y₀,subsolver)
    end
    put!(subworker,subproblems)
end

function init_subworker!(subworker::SubWorker{T,A,S},
                         scenarioproblems::ScenarioProblems,
                         x::A,
                         subsolver::AbstractMathProgSolver,
                         ids::Vector{Int}) where {T <: Real, A <: AbstractArray, S <: LQSolver}
    sp = fetch(scenarioproblems)
    subproblems = Vector{SubProblem{T,A,S}}(length(ids))
    for (i,id) = enumerate(ids)
        y₀ = convert(A,rand(sp.problems[i].numCols))
        subproblems[i] = SubProblem(sp.problems[i],sp.parent,id,probability(sp.scenariodata[i]),x,y₀,subsolver)
    end
    put!(subworker,subproblems)
end

function work_on_subproblems!(subworker::SubWorker{T,A,S},
                              work::Work,
                              cuts::CutQueue{T},
                              decisions::Decisions{A}) where {T <: Real, A <: AbstractArray, S <: LQSolver}
    subproblems::Vector{SubProblem{T,A,S}} = fetch(subworker)
    while true
        t::Int = take!(work)
        if t == -1
            # Worker finished
            return
        end
        x::A = fetch(decisions,t)
        update_subproblems!(subproblems,x)
        for subproblem in subproblems
            cut = subproblem()
            Q::T = cut(x)
            try
                put!(cuts,(t,Q,cut))
            catch err
                if err isa InvalidStateException
                    # Master closed the cut channel. Worker finished
                    return
                end
            end
        end
    end
end

function calculate_subobjective(subworker::SubWorker{T,A,S},
                                x::A) where {T <: Real, A <: AbstractArray, S <: LQSolver}
    subproblems::Vector{SubProblem{T,A,S}} = fetch(subworker)
    if length(subproblems) > 0
        return sum([subproblem.π*subproblem(x) for subproblem in subproblems])
    else
        return zero(T)
    end
end

function fill_submodels!(subworker::SubWorker{T,A,S},
                         scenarioproblems::ScenarioProblems) where {T <: Real, A <: AbstractArray, S <: LQSolver}
    sp = fetch(scenarioproblems)
    subproblems::Vector{SubProblem{T,A,S}} = fetch(subworker)
    for (i,submodel) in enumerate(sp.problems)
        fill_submodel!(submodel,subproblems[i])
    end
end

function fill_submodel!(submodel::JuMP.Model,x::AbstractVector,μ::AbstractVector,λ::AbstractVector,C::Real)
    submodel.colVal = x
    submodel.redCosts = μ
    submodel.linconstrDuals = λ
    submodel.objVal = C
end

function fill_submodel!(submodel::JuMP.Model,subproblem::SubProblem)
    fill_submodel!(submodel,get_solution(subproblem)...)
end

function iterate_parallel!(lshaped::AbstractLShapedSolver{T,A,M,S}) where {T <: Real, A <: AbstractVector, M <: LQSolver, S <: LQSolver}
    wait(lshaped.cutqueue)
    while isready(lshaped.cutqueue)
        # Add new cuts from subworkers
        t::Int,Q::T,cut::SparseHyperPlane{T} = take!(lshaped.cutqueue)
        if !bounded(cut)
            warn("Subproblem ",cut.id," is unbounded, aborting procedure.")
            return :Unbounded
        end
        addcut!(lshaped,cut,Q)
        lshaped.subobjectives[t][cut.id] = Q
        lshaped.finished[t] += 1
        if lshaped.finished[t] == lshaped.nscenarios
            lshaped.Q_history[t] = current_objective_value(lshaped,lshaped.subobjectives[t])
            if lshaped.Q_history[t] <= lshaped.solverdata.Q
                lshaped.solverdata.Q = lshaped.Q_history[t]
            end
            lshaped.x[:] = fetch(lshaped.decisions,t)
            take_step!(lshaped)
        end
    end
    # Resolve master
    t = lshaped.solverdata.timestamp
    if lshaped.finished[t] >= lshaped.parameters.κ*lshaped.nscenarios && length(lshaped.cuts) >= lshaped.nscenarios
        try
            lshaped.mastersolver(lshaped.mastervector)
        catch
            # Master problem could not be solved for some reason. Try to continue
            return :Valid
        end
        if status(lshaped.mastersolver) == :Infeasible
            warn("Master is infeasible, aborting procedure.")
            map(w->put!(w,-1),lshaped.work)
            return :Infeasible
        end
        # Update master solution
        update_solution!(lshaped)
        lshaped.solverdata.θ = calculate_estimate(lshaped)
        # Log progress at current timestamp
        log!(lshaped,t)
        # Check if optimal
        if check_optimality(lshaped)
            # Optimal
            map(w->put!(w,-1),lshaped.work)
            lshaped.solverdata.Q = calculate_objective_value(lshaped,lshaped.x)
            lshaped.Q_history[t] = lshaped.solverdata.Q
            return :Optimal
        end
        # Send new decision vector to workers
        put!(lshaped.decisions,t+1,lshaped.x)
        for w in lshaped.work
            put!(w,t+1)
        end
        # Prepare memory for next timestamp
        lshaped.solverdata.timestamp += 1
        push!(lshaped.subobjectives,zeros(lshaped.nscenarios))
        push!(lshaped.finished,0)
        # Log progress
        log!(lshaped)
    end
    # Just return a valid status for this iteration
    return :Valid
end


text:
#   Unit "pca.jl" of the Diagonalization.jl Package for Julia language
#
#   MIT License
#   Copyright (c) 2019-2021,
#   Marco Congedo, CNRS, Grenoble, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This unit implements Principal Component Analysys (PCA) and whitening
#   (sphering), that is, filters based on the eigendecomposition.

"""
```julia
(1)
function pca(C :: SorH;
             eVar     :: TeVaro = nothing,
             eVarMeth :: Function = searchsortedfirst,
             simple   :: Bool = false)

(2)
function pca(X::Mat;
             covEst   :: StatsBase.CovarianceEstimator = SCM,
             dims     :: Into = ○,
             meanX    :: Tmean = 0,
             wX       :: Tw = ○,
          eVar     :: TeVaro = ○,
          eVarMeth :: Function = searchsortedfirst,
          simple   :: Bool = false)

(3)
function pca(𝐗::VecMat;
             covEst   :: StatsBase.CovarianceEstimator = SCM,
             dims     :: Into = ○,
             meanX    :: Into = 0,
          eVar     :: TeVaro = ○,
          eVarMeth :: Function = searchsortedfirst,
          simple   :: Bool = false,
       metric   :: Metric = Euclidean,
       w        :: Vector = [],
       ✓w       :: Bool = true,
       init     :: SorHo = nothing,
       tol      :: Real = 0.,
       verbose  :: Bool = false)
```

Return a [LinearFilter](@ref) object:

**(1) Principal component analysis**
with real or complex covariance matrix `C` as input.

`C` must be flagged as `Symmetric`, if real, or `Hermitian`,
if either real or complex, see [data input](@ref).

`eVar` and `evarMeth` are keyword optional arguments for defining the
[subspace dimension](@ref) ``p`` using the `.arev` vector given by Eq. [pca.6].
The default values are:
- `eVar=0.999`
- `evarMeth=searchsortedfirst`

If `simple` is set to `true`, ``p`` is set equal to ``n`` and only the
fields `.F` and `.iF` are written in the constructed object.
This option is provided for low-level work when you don't need to define
a subspace dimension or you want to define it by your own methods.

**(2) Principal component analysis**
with a real or complex data matrix `X` as input.

`CovEst`, `dims`, `meanX`, `wX` are optional keyword arguments to
regulate the estimation of the covariance matrix of `X`.
See [covariance matrix estimations](@ref).

Once the covariance matrix estimated, method (1) is invoked
with optional keyword arguments `eVar`, `eVarMeth` and `simple`.

**(3) Principal component analysis**
with a vector of real or complex data matrices `𝐗` as input.

`CovEst`, `dims` and `meanX` are optional keyword arguments to
regulate the estimation of the covariance matrices for
all data matrices in `𝐗`. See [covariance matrix estimations](@ref).

A mean of these covariance matrices is computed using
optional keywords arguments `metric`, `w`, `✓w`, `init`, `tol` and `verbose`.
See [mean covariance matrix estimations](@ref).
By default, the arithmetic mean is computed.

Once the mean covariance matrix estimated, method (1) is invoked
with optional keyword arguments `eVar`, `eVarMeth` and `simple`.

**See also:** [Whitening](@ref), [CSP](@ref), [MCA](@ref), [AJD](@ref).

**Examples:**
```julia
using Diagonalizations, LinearAlgebra, PosDefManifold, Test

# Method (1) real
n, t=10, 100
X=genDataMatrix(n, t)
C=(X*X')/t
pC=pca(Hermitian(C); simple=true)
# or, shortly
pC=pca(ℍ(C); simple=true)

# Method (1) complex
Xc=genDataMatrix(ComplexF64, n, t)
Cc=(Xc*Xc')/t
pCc=pca(Hermitian(Cc); simple=true)

# Method (2) real
pX=pca(X; simple=true)
@test C≈pC.F*pC.D*pC.iF
@test C≈pC.F*pC.D*pC.F'
@test pX≈pC

# Method (2) complex
pXc=pca(Xc; simple=true)
@test Cc≈pCc.F*pCc.D*pCc.iF
@test Cc≈pCc.F*pCc.D*pCc.F'
@test pXc≈pCc

# Method (3) real
k=10
Xset=[genDataMatrix(n, t) for i=1:k]

# pca on the average covariance matrix
p=pca(Xset)

# ... selecting subspace dimension allowing an explained variance = 0.5
p=pca(Xset; eVar=0.5)

# ... averaging the covariance matrices using the logEuclidean metric
p=pca(Xset; metric=logEuclidean, eVar=0.5)

# ... giving weights `w` to the covariance matrices
p=pca(Xset; metric=logEuclidean, w=abs2.(randn(k)), eVar=0.5)

# ... subtracting the mean
p=pca(Xset; meanX=nothing, metric=logEuclidean, w=abs2.(randn(k)), eVar=0.5)

# pca on the average of the covariance matrices computed along dims 1
p=pca(Xset; dims=1)

# explained variance
p.eVar

# name of the filter
p.name

using Plots
# plot regularized accumulated eigenvalues
plot(p.arev)

# plot the original covariance matrix and the rotated covariance matrix
 Cmax=maximum(abs.(C));
 h1 = heatmap(C, clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="C");
 D=pC.F'*C*pC.F;
 Dmax=maximum(abs.(D));
 h2 = heatmap(D, clim=(0, Dmax), yflip=true, c=:amp, title="F'*C*F");
 📈=plot(h1, h2, size=(700, 300))
# savefig(📈, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigPCA.png")
```

 ![Figure PCA](assets/FigPCA.png)

```julia
# Method (3) complex
k=10
Xcset=[genDataMatrix(ComplexF64, n, t) for i=1:k]

# pca on the average covariance matrix
pc=pca(Xcset)

# ... selecting subspace dimension allowing an explained variance = 0.5
pc=pca(Xcset; eVar=0.5)

# ... averaging the covariance matrices using the logEuclidean metric
pc=pca(Xcset; metric=logEuclidean, eVar=0.5)

# ... giving weights `w` to the covariance matrices
pc=pca(Xcset; metric=logEuclidean, w=abs2.(randn(k)), eVar=0.5)

# ... subtracting the mean
pc=pca(Xcset; meanX=nothing, metric=logEuclidean, w=abs2.(randn(k)), eVar=0.5)

# pca on the average of the covariance matrices computed along dims 1
pc=pca(Xcset; dims=1)

# explained variance
pc.eVar

# name of the filter
pc.name
```

"""
function pca(C :: SorH;
             eVar     :: TeVaro = ○,
             eVarMeth :: Function = searchsortedfirst,
             simple   :: Bool = false)

   args=("Principal Component Analysis", false)

   λ, U = eig(C) # get evd
   λ = _checkλ(λ) # make sure no imaginary noise is present (for complex data)

   simple ? LF(U, Matrix(U'), Diagonal(λ), ○, ○, ○, args...) :
   begin
      eVar, D, U, p, arev=_ssd!(eVar, λ, U, _minDim(C), eVarMeth) # find subspace
      LF(U, Matrix(U'), D, eVar, λ, arev, args...)
   end
end

# Principal Component Analysis with data as input: API
function pca(X :: Mat;
             covEst   :: StatsBase.CovarianceEstimator = SCM,
             dims     :: Into = ○,
             meanX    :: Tmean = 0,
             wX       :: Tw = ○,
          eVar     :: TeVaro = ○,
          eVarMeth :: Function = searchsortedfirst,
          simple   :: Bool = false)

   dims===○ && (dims=_set_dims(X))
   _check_data(X, dims, covEst, meanX, wX)===○ && return
   args=("Principal Component Analysis", false)

   LF(_getEVD(X, covEst, dims, meanX, wX, eVar, eVarMeth, simple)..., args...)
end

function pca(𝐗 :: VecMat;
             covEst   :: StatsBase.CovarianceEstimator = SCM,
             dims     :: Into = ○,
             meanX    :: Into = 0,
          eVar     :: TeVaro = ○,
          eVarMeth :: Function = searchsortedfirst,
          simple   :: Bool = false,
       metric   :: Metric = Euclidean,
       w        :: Vector = [],
       ✓w       :: Bool = true,
       init     :: SorHo = nothing,
       tol      :: Real = 0.,
       verbose  :: Bool = false)

    Metric==VonNeumann && throw(ArgumentError(📌*", pca function: A solution for the mean is not available for the Von Neumann metric. Use another metric as `metric` argument"))
    dims===○ && (dims=_set_dims(𝐗))
    _check_data(𝐗, dims, covEst, meanX, ○)===○ && return

    𝐂=_cov(𝐗; covEst=covEst, dims = dims, meanX = meanX)

    pca(mean(metric, 𝐂;
             w = w, ✓w = ✓w,
             init = init===○ ? ○ : Hermitian(init), #just init here when you upfate PodDefManifold
             tol = tol, verbose = verbose),
        eVar = eVar, eVarMeth = eVarMeth, simple = simple)
end

"""
```julia
(1)
function whitening(C :: SorH;
                   eVar     :: TeVaro=○,
                   eVarMeth :: Function=searchsortedfirst,
                   simple   :: Bool=false)

(2)
function whitening(X::Mat;
                   covEst   :: StatsBase.CovarianceEstimator = SCM,
                   dims     :: Into = ○,
                   meanX    :: Tmean = 0,
                   wX       :: Tw = ○,
                eVar     :: TeVaro = ○,
                eVarMeth :: Function = searchsortedfirst,
                simple   :: Bool = false)

(3)
function whitening(𝐗::VecMat;
                   covEst   :: StatsBase.CovarianceEstimator = SCM,
                   dims     :: Into = ○,
                   meanX    :: Into = 0,
                eVar     :: TeVaro = ○,
                eVarMeth :: Function = searchsortedfirst,
                simple   :: Bool = false,
             metric   :: Metric = Euclidean,
             w        :: Vector = [],
             ✓w       :: Bool = true,
             init     :: SorHo = nothing,
             tol      :: Real = 0.,
             verbose  :: Bool = false)

```

Return a [LinearFilter](@ref) object:

**(1) Whitening**
with real or complex covariance matrix `C` as input.

`C` must be flagged as `Symmetric`, if real or `Hermitian`,
if real or complex, see [data input](@ref).

`eVar` and `evarMeth` are keyword optional arguments for defining the
[subspace dimension](@ref) ``p`` using the `.arev` vector given by Eq. [pca.6],
see [PCA](@ref).
The default values are:
- `eVar=0.999`
- `evarMeth=searchsortedfirst`

If `simple` is set to `true`, ``p`` is set equal to ``n`` and only the
fields `.F` and `.iF` are written in the constructed object.
This option is provided for low-level work when you don't need to define
a subspace dimension or you want to define it by your own methods.

**(2) Whitening**
with a real or complex data matrix `X` as input.

`CovEst`, `dims`, `meanX`, `wX` are optional keyword arguments to
regulate the estimation of the covariance matrix of `X`.
See [covariance matrix estimations](@ref).

Once the covariance matrix estimated, method (1) is invoked
with optional keyword arguments `eVar`, `eVarMeth` and `simple`.

**(3) Whitening**
with a vector of real or complex data matrices `𝐗` as input.

`CovEst`, `dims` and `meanX` are optional keyword arguments to
regulate the estimation of the covariance matrices for
all data matrices in `𝐗`. See [covariance matrix estimations](@ref).

A mean of these covariance matrices is computed using
optional keywords arguments `metric`, `w`, `✓w`, `init`, `tol` and `verbose`.
See [mean covariance matrix estimations](@ref).
By default, the arithmetic mean is computed.

Once the mean covariance matrix estimated, method (1) is invoked
with optional keyword arguments `eVar`, `eVarMeth` and `simple`.

**See also:** [PCA](@ref), [CCA](@ref).

**Examples:**
```julia
using Diagonalizations, LinearAlgebra, PosDefManifold, Test

# Method (1) real
n, t=10, 100
X=genDataMatrix(n, t)
C=(X*X')/t
wC=whitening(Hermitian(C); simple=true)
# or, shortly
wC=whitening(ℍ(C); simple=true)

# Method (1) complex
Xc=genDataMatrix(ComplexF64, n, t)
Cc=(Xc*Xc')/t
wCc=whitening(Hermitian(Cc); simple=true)

# Method (2) real
wX=whitening(X; simple=true)
@test wC.F'*C*wC.F≈I
@test wX.F'*C*wX.F≈I
@test wX≈wC

# Method (2) complex
wXc=whitening(Xc; simple=true)
@test wCc.F'*Cc*wCc.F≈I
@test wXc.F'*Cc*wXc.F≈I
@test wXc≈wCc

# Method (3) real
k=10
Xset=[genDataMatrix(n, t) for i=1:k]

# whitening on the average covariance matrix
w=whitening(Xset)

# ... selecting subspace dimension allowing an explained variance = 0.5
w=whitening(Xset; eVar=0.5)

# ... averaging the covariance matrices using the logEuclidean metric
w=whitening(Xset; metric=logEuclidean, eVar=0.5)

# ... giving weights `w` to the covariance matrices
w=whitening(Xset; metric=logEuclidean, w=abs2.(randn(k)), eVar=0.5)

# ... subtracting the mean
w=whitening(Xset; meanX=nothing, metric=logEuclidean, w=abs2.(randn(k)), eVar=0.5)

# whitening on the average of the covariance matrices computed along dims 1
w=whitening(Xset; dims=1)

# explained variance
w.eVar

# name of the filter
w.name

using Plots
# plot regularized accumulated eigenvalues
plot(w.arev)

# plot the original covariance matrix and the whitened covariance matrix
 Cmax=maximum(abs.(C));
 h1 = heatmap(C, clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="C");
 D=wC.F'*C*wC.F;
 h2 = heatmap(D, clim=(0, 1), yflip=true, c=:amp, title="F'*C*F");
 📈=plot(h1, h2, size=(700, 300))
# savefig(📈, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigWhitening.png")
```

 ![Figure Whitening](assets/FigWhitening.png)

```julia

# Method (3) complex
k=10
Xcset=[genDataMatrix(ComplexF64, n, t) for i=1:k]

# whitening on the average covariance matrix
wc=whitening(Xcset)

# ... selecting subspace dimension allowing an explained variance = 0.5
wc=whitening(Xcset; eVar=0.5)

# ... averaging the covariance matrices using the logEuclidean metric
wc=whitening(Xcset; metric=logEuclidean, eVar=0.5)

# ... giving weights `w` to the covariance matrices
wc=whitening(Xset; metric=logEuclidean, w=abs2.(randn(k)), eVar=0.5)

# ... subtracting the mean
wc=whitening(Xcset; meanX=nothing, metric=logEuclidean, w=abs2.(randn(k)), eVar=0.5)

# whitening on the average of the covariance matrices computed along dims 1
wc=whitening(Xcset; dims=1)

# explained variance
wc.eVar

# name of the filter
wc.name
```

"""
function whitening(C :: SorH;
                   eVar     :: TeVaro=○,
                   eVarMeth :: Function=searchsortedfirst,
                   simple   :: Bool=false)

   args=("Whitening", false)

   λ, U = eig(C) # get evd
   λ = _checkλ(λ) # make sure no imaginary noise is present (for complex data)
   if simple
     if eltype(C)<:Real
		  ispos(λ; tol=eps(eltype(C)), rev=true, 🔔=true,
		  msg="negative or almost zero eigenvalue") || throw(ArgumentError("A `simple` linear filter cannot be created. See the warning that has been printed in Julia's REPL"))
	 end
     D=Diagonal(λ)
     LF(U*D^-0.5, D^0.5*Matrix(U'), D, ○, ○, ○, args...)
   else
     eVar, D, U, p, arev=_ssd!(eVar, λ, U, _minDim(C), eVarMeth) # subspace dimension
     LF(U*D^-0.5, D^0.5*Matrix(U'), D, eVar, λ, arev, args...)
   end
end

function whitening(X::Mat;
                   covEst   :: StatsBase.CovarianceEstimator = SCM,
                   dims     :: Into = ○,
                   meanX    :: Tmean = 0,
                   wX       :: Tw = ○,
                eVar     :: TeVaro = ○,
                eVarMeth :: Function = searchsortedfirst,
                simple   :: Bool = false)

   dims===○ && (dims=_set_dims(X))
   _check_data(X, dims, covEst, meanX, wX)===○ && return
   args=("Whitening", false)

   LF(_getWhi(X, covEst, dims, meanX, wX, eVar, eVarMeth, simple)..., args...)
end

function whitening(𝐗::VecMat;
                   covEst   :: StatsBase.CovarianceEstimator = SCM,
                   dims     :: Into = ○,
                   meanX    :: Into = 0,
                eVar     :: TeVaro = ○,
                eVarMeth :: Function = searchsortedfirst,
                simple   :: Bool = false,
             metric   :: Metric = Euclidean,
             w        :: Vector = [],
             ✓w       :: Bool = true,
             init     :: SorHo = nothing,
             tol      :: Real = 0.,
             verbose  :: Bool = false)

   Metric==VonNeumann && throw(ArgumentError(📌*", whitening function: A solution for the mean is not available for the Von Neumann metric. Use another metric as `metric` argument"))
   dims===○ && (dims=_set_dims(𝐗))
   _check_data(𝐗, dims, covEst, meanX, ○)===○ && return

   whitening(mean(metric, _cov(𝐗; covEst=covEst, dims=dims, meanX=meanX);
                  w = w, ✓w = ✓w,
                  init = init===○ ? ○ : Hermitian(init), #jet init here when you upfate PodDefManifold
                  tol = tol, verbose = verbose),
             eVar = eVar, eVarMeth = eVarMeth, simple = simple)
end


text:
include("cairoWrapper.jl")

# get a cr and surface
y = initCairo("out.svg",400.0,400.0)

setLineWidthCairo(y,2.0)

drawRectCairo(y,140.0,140.0,50.0,50.0)

drawTextCairo(y,3.0, 40.0, 20.0,"Sans", "Julia wrapper for Cairo")

drawDiamondCairo(y,100.0,330.0,50.0,50.0)

drawRelLineCairo(y,40.0,200.0,200.0,100.0)

drawArrowCairo(y,60.0,130.0,30.0,50.0,270.0)

drawArrowLineCairo(y,100.0,300.0,-30.0,-30.0,10.0,10.0)

destroyCairo(y)

max_stars_repo_name: aklapatch/CodeVisualized

max_issues_repo_name: aklapatch/CodeVisualized

max_forks_repo_name: aklapatch/CodeVisualized


text:
# limits number of bytes written to an io stream (originally from https://github.com/JuliaDebug/Debugger.jl/blob/master/src/limitio.jl)
# useful to detect messages that would go over limit when converted to bytes.
mutable struct LimitIO{IO_t <: IO} <: IO
    io::IO_t
    maxbytes::Int
    n::Int # max bytes to write
end
LimitIO(io::IO, maxbytes) = LimitIO(io, maxbytes, 0)

function Base.write(io::LimitIO, v::UInt8)
    io.n > io.maxbytes && throw(gRPCMessageTooLargeException(io.maxbytes, io.n))
    nincr = write(io.io, v)
    io.n += nincr
    nincr
end

"""
Default maximum gRPC message size
"""
const DEFAULT_MAX_MESSAGE_LENGTH = 1024*1024*4

"""
    struct gRPCMessageTooLargeException
        limit::Int
        encountered::Int
    end

A `gRPMessageTooLargeException` exception is thrown when a message is
encountered that has a size greater than the limit configured.
Specifically, `max_recv_message_length` while receiving  and
`max_send_message_length` while sending.

A `gRPMessageTooLargeException` has the following members:

- `limit`: the limit value that was exceeded
- `encountered`: the amount of data that was actually received
    or sent before this error was triggered. Note that this may
    not correspond to the full size of the data, as error may be
    thrown before actually materializing the complete data.
"""
struct gRPCMessageTooLargeException <: gRPCException
    limit::Int
    encountered::Int
end

Base.show(io::IO, m::gRPCMessageTooLargeException) = print(io, "gRPMessageTooLargeException($(m.limit), $(m.encountered)) - Encountered message size $(m.encountered) > max configured $(m.limit)")


text:
export fund_args, norm_ang, ee_ct

"""
    fund_args(t)

Computes the fundamental arguments (mean elements) of the Sun and Moon.

# Arguments
- `t::Real`: TDB time in Julian centuries since J2000.0

# Returns
`[l,l′,F,D,Ω]` where
- `l`: Mean anomaly of the Moon
- `l′`: Mean anomaly of the Sun
- `F`: mean argument of the latitude of the Moon
- `D`: Mean elongation of the Moon from the Sun
- `Ω` Mean longitude of the Moon's ascending node
"""
function fund_args(t::T) where {T<:Real}
    l = 485868.249036 +
        t * (1717915923.2178 + t * (31.8792 + t * (0.051635 + t * (-0.00024470))))
    l′ = 1287104.79305 +
         t * (129596581.0481 + t * (-0.5532 + t * (0.000136 + t * (-0.00001149))))
    F = 335779.526232 +
        t * (1739527262.8478 + t * (-12.7512 + t * (-0.001037 + t * (0.00000417))))
    D = 1072260.70369 +
        t * (1602961601.2090 + t * (-6.3706 + t * (0.006593 + t * (-0.00003169))))
    Ω = 450160.398036 +
        t * (-6962890.5431 + t * (7.4722 + t * (0.007702 + t * (-0.00005939))))
    out = @SVector T[l, l′, F, D, Ω]
    return @. rem(out, ASEC360) * ASEC2RAD
end

function read_nu2000k()
    napl = readdlm(nu2000k("napl.csv"), ',')
    nals = readdlm(nu2000k("nals.csv"), ',')
    cpl = readdlm(nu2000k("cpl.csv"), ',')
    cls = readdlm(nu2000k("cls.csv"), ',')
    return nals, cls, napl, cpl
end

function read_iau2000a()
    lunisolar = readdlm(iers2003("tab5.3a.txt"); comment_char='*', comments=true)[1:678, :]
    planetary = readdlm(iers2003("tab5.3b.txt"); skipstart=5)
    nals = lunisolar[1:678, 1:5]                      # [:l,:l′,:F,:D,:Ω]
    cls = lunisolar[1:678, [7, 8, 11, 9, 10, 13]] .* 1e4  # [:A, :A′, :A″, :B, :B′, :B″]
    napl = planetary[1:687, 2:15]                    # [:l, :l′, :F, :D, :Ω, :Me, :Ve, :E, :Ma, :J, :Sa, :U, :Ne, :pA]
    cpl = planetary[1:687, 17:20] .* 1e4             # [:A, :A″, :B, :B″]
    return nals, cls, napl, cpl
end

function read_cterms()
    cterms = readdlm(iers2003("tab5.4.txt"); comment_char='j', comments=true, skipstart=50)
    ke0 = cterms[1:33, 4:17]
    ke1 = cterms[34, 4:17]
    se0 = cterms[1:33, 2:3] .* 1e-6
    se1 = cterms[34, 2:3] .* 1e-6
    return ke0, ke1, se0, se1
end

"""
    norm_ang(θ)

Normalizes an angle into range 0 <= θ <= 2π
"""
function norm_ang(angle::Real)
    return rem2pi(angle, RoundDown)
end

"""
    ee_ct(jd_high, jd_low)

Computes the "complementary terms" of the equation of the equinoxes.

# Arguments
- `jd_high::Real`: High-order part of UT1 Julian date
- `jd_low::Real`: Low-order part of UT1 Julian date

# Optional Arguments
- `accuracy::Symbol=:full`: Sets the accuracy level of `:full` or `:reduced`
"""
function ee_ct(jd_high::T, jd_low::T; accuracy::Symbol=:full) where {T<:Real}
    t = ((jd_high - T0) + jd_low) / 36525.0

    fa = zeros(T, 14)

    if accuracy == :full
        # Fundamental arguments
        # Mean anomaly of the moon
        fa[1] = norm_ang((485868.249036 +
                          (715923.2178 + (31.8792 + (0.051635 + (-0.00024470) * t) * t) * t) *
                          t) * ASEC2RAD + rem(1325.0 * t, 1.0) * 2π)
        # Mean anomaly of the Sun
        fa[2] = norm_ang((1287104.793048 +
                          (1292581.0481 +
                           (-0.5532 + (+0.000136 + (-0.00001149) * t) * t) * t) * t) *
                         ASEC2RAD + rem(99.0 * t, 1.0) * 2π)
        # Mean Longitude of the Moon minus Mean Longitude of the Ascending Node of the Moon.
        fa[3] = norm_ang((335779.526232 +
                          (295262.8478 +
                           (-12.7512 + (-0.001037 + (0.00000417) * t) * t) * t) * t) *
                         ASEC2RAD + rem(1342.0 * t, 1.0) * 2π)
        # Mean Elongation of the Moon from the Sun.
        fa[4] = norm_ang((1072260.703692 +
                          (1105601.2090 +
                           (-6.3706 + (0.006593 + (-0.00003169) * t) * t) * t) * t) *
                         ASEC2RAD + rem(1236.0 * t, 1.0) * 2π)
        # Mean Longitude of the Ascending Node of the Moon.
        fa[5] = norm_ang((450160.398036 +
                          (-482890.5431 + (7.4722 + (0.007702 + (-0.00005939) * t) * t) * t) *
                          t) * ASEC2RAD + rem(-5.0 * t, 1.0) * 2π)
        fa[6] = norm_ang(4.402608842 + 2608.7903141574 * t)
        fa[7] = norm_ang(3.176146697 + 1021.3285546211 * t)
        fa[8] = norm_ang(1.753470314 + 628.3075849991 * t)
        fa[9] = norm_ang(6.203480913 + 334.0612426700 * t)
        fa[10] = norm_ang(0.599546497 + 52.9690962641 * t)
        fa[11] = norm_ang(0.874016757 + 21.3299104960 * t)
        fa[12] = norm_ang(5.481293872 + 7.4781598567 * t)
        fa[13] = norm_ang(5.311886287 + 3.8133035638 * t)
        fa[14] = (0.024381750 + 0.00000538691 * t) * t
        # Evaluate the complementary terms.
        s0 = 0.0
        s1 = 0.0
        @inbounds @simd for i in 1:33
            a = 0.0
            for j in 1:14
                a += ke0[i, j] * fa[j]
            end
            s0 += se0[i, 1] * sin(a) + se0[i, 2] * cos(a)
        end
        a = 0.0
        @inbounds @simd for j in 1:14
            a += ke1[j] * fa[j]
        end
        s1 += se1[1] * sin(a) + se1[2] * cos(a)
        c_terms = s0 + s1 * t

    elseif accuracy == :reduced
        fa2 = fund_args(t)
        c_terms = 2640.96e-6 * sin(fa2[5]) +
                  63.52e-6 * sin(2.0 * fa2[5]) +
                  11.75e-6 * sin(2.0 * fa2[3] - 2.0 * fa2[4] + 3.0 * fa2[5]) +
                  11.21e-6 * sin(2.0 * fa2[3] - 2.0 * fa2[4] + fa2[5]) -
                  4.55e-6 * sin(2.0 * fa2[3] - 2.0 * fa2[4] + 2.0 * fa2[5]) +
                  2.02e-6 * sin(2.0 * fa2[3] + 3.0 * fa2[5]) +
                  1.98e-6 * sin(2.0 * fa2[3] + fa2[5]) - 1.72e-6 * sin(3.0 * fa2[5]) -
                  0.87e-6 * t * sin(fa2[5])
    end
    return c_terms * ASEC2RAD
end

# Preallocate constants
const ke0, ke1, se0, se1 = read_cterms()
const nals_k, cls_k, napl_k, cpl_k = read_nu2000k()
const nals_a, cls_a, napl_a, cpl_a = read_iau2000a()