text:
module Callbacks

using Photon
using Photon: getmetricvalue, history
using Photon.Metrics
using Printf

export Callback, EarlyStop, EpochSave, AutoSave, SilentMeter, ConsoleMeter, PlotMeter, FileMeter, TensorBoardMeter

abstract type Callback end

"""
A callback should implement the following function signature.
"""
function (::Callback)(workout::Workout, phase::Symbol) end

"""
Save the workout at the end of every epoch. Optionally provide a filename.

# Usage
```julia
train!(workout, data, cb=EpochSave())
```
"""
struct EpochSave <: Callback
    filename::Union{String, Nothing}
    EpochSave(filename=nothing) = new(filename)
end

function (c::EpochSave)(workout::Workout, phase::Symbol)
    if phase == :valid
        if c.filename !== nothing
            saveworkout(workout, c.filename)
        else
            saveworkout(workout)
        end
    end
end

"""
Save the Workout if a certain metric has improved since the last epoch

# Usage
```julia

# save as long as the validation loss is declining
train!(workout, data, cb=AutoSave(:val_loss))
```

"""
mutable struct AutoSave <: Callback
    value::Float64
    metric::Symbol
    filename::Union{String, Nothing}
    mode::Symbol
    function AutoSave(metric::Symbol, filename=nothing, mode=:min)
        init_value = mode == :min ? Inf : -Inf
        new(init_value, metric, filename, mode)
    end
end

function (c::AutoSave)(workout::Workout, phase::Symbol)
    if phase == :valid
        getmetricvalue(workout, c.metric) do x
            condition = c.mode == :min ? x < c.value : x > c.value
            if condition
                if c.filename !== nothing
                    saveworkout(workout, c.filename)
                else
                    saveworkout(workout)
                end
                c.value = x
            end
        end
    end
end

"""
Stop the training if a certain metric didn't improve
"""
mutable struct EarlyStop <: Callback
    value::Float64
    metric::Symbol
    mode::Symbol
    function EarlyStop(metric::Symbol, mode=:min)
        init_value = mode == :min ? Inf : -Inf
        new(init_value, metric, mode)
    end
end

function (c::EarlyStop)(workout::Workout, phase::Symbol)
    if phase == :valid
        getmetricvalue(workout, c.metric) do x
            condition = c.mode == :min ? x < c.value : x > c.value
            if condition
                stop(workout,"$(c.metric) didn't improve anymore")
            else
                c.value = x
            end
        end
    end
end

"""
A meter is responsible for presenting metric values. It can be used just as
a regular callback argument to the *train!* function.
A meter is not limited to printing results to the console output, it can also be
showing it on a TensorBoard or storing results in a database for example.
"""
abstract type Meter <: Callback end

"""
Use a SilentMeter in case no output is required. The default of
*train!* function is to use a ConsoleMeter and using the SilentMeter this
behavior can be overriden.

# Usage

```julia
train!(workout, data, val_data; epochs=5, cb=SilentMeter())
```
"""
struct SilentMeter <: Meter
end

(meter::SilentMeter)(workout::Workout, phase::Symbol) = ()

"""
Logs metrics to the console output. By default it will only log at the end of an epoch
and log the epoch, step, loss and validation loss. This is also the default configuration
when you run train! without specifying oter callbacks.

# Usage

```julia
meter = ConsoleMeter([:loss, :val_accuracy, :accuracy]; epochOnly=false)
train!(workout, data, val_data; epochs=5, cb=meter)
```
"""
mutable struct ConsoleMeter <: Meter
    throttle::Float64
    next::Float64
    metricnames::Vector{Symbol}
    epochOnly

    function ConsoleMeter(metrics=[:loss, :val_loss]; throttle=1.0, epochOnly=true)
        new(throttle, 0.0, metrics, epochOnly)
    end
end

function (meter::ConsoleMeter)(workout::Workout, phase::Symbol)
    meter.epochOnly && phase == :train && return
    now = time()
    if now > meter.next || phase == :valid
        result = ""
        for metricname in meter.metricnames
            getmetricvalue(workout, metricname) do value
                    result *= @sprintf  " - %s: %1.4f" metricname value
            end
        end
        if result != ""
            s = @sprintf "[%4d:%7d]" workout.epochs workout.steps
            print("\r", s, result)
            meter.next = now + meter.throttle
        end
        phase == :valid && println()
    end
end

"""
Logs metrics to a TensorBoard file so it can be viewed with TensorBoard. By default
it will log the metrics loss and val_loss at the end of each training step
and end of the validation phase.

This meter depends on the TensorBoardLogger to be installed. So if you didn't do so
already, please run:

    import Pkg; Pkg.add("TensorBoardLogger")
"""
mutable struct TensorBoardMeter <: Meter
    logger
    path::String
    metrics
    last_processed::IdDict{Symbol,Int}

    function TensorBoardMeter(path="./tensorboard_logs/runs", metrics=[:loss, :val_loss])
        try
            @eval import TensorBoardLogger
        catch
            @warn "Package TensorBoardLogger not installed"
        end
        new(nothing, path, metrics, IdDict())
    end
end

function (meter::TensorBoardMeter)(workout::Workout, phase::Symbol)
    if meter.logger === nothing
        meter.logger = TensorBoardLogger.TBLogger(meter.path)
    end

    for metric in meter.metrics
        getmetricvalue(workout, metric) do value
            last = get(meter.last_processed, metric, -1)
            workout.steps > last && TensorBoardLogger.log_value(meter.logger,
                string(metric), value, step=workout.steps)
            meter.last_processed[metric] = workout.steps
        end
    end
end

"""
Logs metrics to a delimeted text file. By default it will log the metrics loss and
val_loss at the end of each training step and end of the validation phase.

"""
mutable struct FileMeter <: Meter
    fileio
    filename::String
    metrics
    last_processed::IdDict{Symbol,Int}
    epochOnly::Bool

    function FileMeter(filename="traininglog.txt", metrics=[:loss, :val_loss], epochOnly=true)
        try
            @eval import DelimitedFiles
        catch
            @warn "Package DelimitedFiles not installed"
        end
        io = open(filename, "w")
        new(io, filename, metrics, IdDict(), epochOnly)
    end
end

function (meter::FileMeter)(workout::Workout, phase::Symbol)
    meter.epochOnly && phase == :train && return
    for metric in meter.metrics
        getmetricvalue(workout, metric) do value
            last = get(meter.last_processed, metric, -1)
            workout.steps > last && DelimitedFiles.writedlm(meter.fileio,
                [(metric, value, workout.steps)])
            meter.last_processed[metric] = workout.steps
        end
    end
    flush(meter.fileio)
end

"""
Plot metrics using the Plots module. At the end of each epoch the plot will
be updated with the values of the metrics. This works especially nice if you
are prototyping some model in an IDE like Juno.
"""
struct PlotMeter <: Meter
    plt
    metrics::Vector{Symbol}
    Plots::Module

    function PlotMeter(Plots::Module, metrics=[:loss, :val_loss])
        plt = Plots.plot(length(metrics),xlabel = "steps", ylabel="values", label=metrics)
        Plots.display(plt)
        new(plt, metrics, Plots)
    end
end

function (m::PlotMeter)(workout::Workout, phase::Symbol)
    phase == :train && return
    for (i, metric) in enumerate(m.metrics)
        if hasmetric(workout, metric)
            m.plt[i] = history(workout, metric)
        end
    end
    m.Plots.display(m.plt)
end

end


text:
function convert_common_kwargs(opt::NonconvexPercival.AugLag, opt_kwargs;
    cb=nothing,
    maxiters=nothing,
    maxtime=nothing,
    abstol=nothing,
    reltol=nothing)

    conv_opt_kwargs = (; opt_kwargs...)

    if !isnothing(cb)
        @warn "common callback argument is currently not used by $(opt)"
    end

    if !isnothing(maxiters)
        conv_opt_kwargs = (; conv_opt_kwargs..., max_iter=maxiters)
    end

    if !isnothing(maxtime)
        conv_opt_kwargs = (; conv_opt_kwargs..., max_time = maxtime)
    end

    if !isnothing(abstol)
        conv_opt_kwargs = (; conv_opt_kwargs..., atol = abstol)
    end

    if !isnothing(reltol)
        conv_opt_kwargs = (; conv_opt_kwargs..., rtol =reltol)
    end

    return conv_opt_kwargs
end

function __create_options(opt::NonconvexPercival.AugLag;
    opt_kwargs=nothing)

    options = !isnothing(opt_kwargs) ? NonconvexPercival.AugLagOptions(;opt_kwargs...) : NonconvexPercival.AugLagOptions()

    return options
end

function _create_options(opt::NonconvexPercival.AugLag;
    opt_kwargs=nothing,
    sub_options=nothing,
    convergence_criteria=nothing)

    options = (; options = __create_options(opt, opt_kwargs=opt_kwargs))

    return options
end

check_optimizer_backend(opt::NonconvexPercival.AugLag) = false

include("nonconvex.jl")

ext: jl


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Icestorm_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Icestorm")
JLLWrappers.@generate_main_file("Icestorm", UUID("d066d2d9-0d34-5c4c-ae28-f0a7a1fd9f53"))
end  # module Icestorm_jll


text:
module Hedge

import ..Clarus

function equivalents(;params...)
  return Clarus.api_request("Hedge","Equivalents";params...)
end

function equivalents01(;params...)
  return Clarus.api_request("Hedge","Equivalents01";params...)
end

end #MODULE-END


text:
export generate_ham_1d_hubbard

function generate_ham_1d_hubbard(t::Float64, U::Float64, nsite::Integer, μ::Float64)
    ham = FermionOperator()
    #斥力項
    for i in 1:nsite
        #up_index,down_indexの定義は、QC_materialを参照。
        up = up_index(i)
        down = down_index(i)
        ham += FermionOperator("$(up)^ $(down)^ $(up) $(down)", -U) #左側に生成演算子。右際に消滅演算子をもっていく過程で半交換関係が1回でマイナスをつける。
    end
    #ホッピング項
    for i in 1:nsite-1
        ham += FermionOperator("$(up_index(i+1))^ $(up_index(i))", t)
        ham += FermionOperator("$(up_index(i))^ $(up_index(i+1))", t)
        ham += FermionOperator("$(down_index(i+1))^ $(down_index(i))", t)
        ham += FermionOperator("$(down_index(i))^ $(down_index(i+1))", t)
    end

    #ケミカルポテンシャルの項
    for i in 1:nsite
        up = up_index(i)
        down = down_index(i)
        ham += FermionOperator("$(up)^  $(up) ", -μ)
        ham += FermionOperator("$(down)^ $(down)", -μ)
    end
    ham
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* File containg models of columns: distillation, stripping, absorbers
#* rectifier, ....
#*
#* The default nomenclature is:
#*		Type_Column_reboilertype_condensertyper
#*
#* where:
#*	Type = refluxed or reboiled or section
#*	Column = Stripping, Absorption, Rectifier, Distillation
#*	Reboiler type (if exists) = kettle or thermosyphon
#*	Condenser type (if exists) = with subccoling or without subcooling
#*
#*-----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id: column.mso 511 2008-05-12 17:25:33Z paula $
#*---------------------------------------------------------------------
type SectionColumnBasic
	SectionColumnBasic()=begin
		PP=outers.PP
		NComp=outers.NComp
		[

		],
		[

		],
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of components"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Feed Trays",
				:Default=>3,
				:Protected=>true
			)),
			fill(DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Feed tray Location",
				:Default=>2
			)),(NumberOfFeeds)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of trays",
				:Default=>8
			)),
			fill(DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Liquid Side Tray Index",
				:Default=>0,
				:Hidden=>true
			)),(NumberOfTrays)),
			fill(DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Vapour Side Tray Index",
				:Default=>0,
				:Hidden=>true
			)),(NumberOfTrays)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Liquid Side Stream Location",
				:Default=>2
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Vapour Side Stream Location",
				:Default=>2
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Valid=>["Reepmeyer", "Feehery_Fv", "Roffel_Fv", "Klingberg", "Wang_Fv", "Elgue"],
				:Default=>"Reepmeyer"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Valid=>["Francis", "Wang_Fl", "Olsen", "Feehery_Fl", "Roffel_Fl"],
				:Default=>"Francis"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Tray Diameter",
				:Default=>1.600
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Tray Spacing",
				:Default=>0.600
			)),
			fraction(Dict{Symbol,Any}(
				:Brief=>"Fraction of the active area that is occupied by the holes with respect to the total tray area",
				:Default=>0.10
			)),
			fraction(Dict{Symbol,Any}(
				:Brief=>"Fraction of the downcomer area with respect to the total tray area",
				:Default=>0.20
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Weir length",
				:Default=>1
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Weir height",
				:Default=>0.05
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Number of liquid passes in the tray",
				:Lower=>1,
				:Default=>1
			)),
			heat_rate(Dict{Symbol,Any}(
				:Brief=>"Rate of heat supply",
				:Default=>0
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Aeration fraction",
				:Default=>1
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Dry pressure drop coefficient",
				:Default=>0.60
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Murphree efficiency for All Trays",
				:Default=>1
			)),
			[
				:(TCI*"K" = TRAYS(CONTROL.Tindicator_TrayNumber).OutletVapour.T),
				:(PCI*"atm" = TRAYS(CONTROL.Pindicator_TrayNumber).OutletVapour.P),
				:(TRAYS.LFlowModel = 1),
				:(TRAYS.LFlowModel = 2),
				:(TRAYS.LFlowModel = 3),
				:(TRAYS.LFlowModel = 4),
				:(TRAYS.LFlowModel = 5),
				:(TRAYS.VFlowModel = 1),
				:(TRAYS.VFlowModel = 2),
				:(TRAYS.VFlowModel = 3),
				:(TRAYS.VFlowModel = 4),
				:(TRAYS.VFlowModel = 5),
				:(TRAYS.VFlowModel = 6),
				:(LiquidConnector.F = LiquidInlet.F),
				:(LiquidConnector.T = LiquidInlet.T),
				:(LiquidConnector.P = LiquidInlet.P),
				:(LiquidConnector.z = LiquidInlet.z),
				:(LiquidConnector.v = LiquidInlet.v),
				:(LiquidConnector.h = LiquidInlet.h),
				:(VapourConnector.F = VapourInlet.F),
				:(VapourConnector.T = VapourInlet.T),
				:(VapourConnector.P = VapourInlet.P),
				:(VapourConnector.z = VapourInlet.z),
				:(VapourConnector.v = VapourInlet.v),
				:(VapourConnector.h = VapourInlet.h),
				:(LiquidOutlet.F = TRAYS(NumberOfTrays).OutletLiquid.F),
				:(LiquidOutlet.T = TRAYS(NumberOfTrays).OutletLiquid.T),
				:(LiquidOutlet.P = TRAYS(NumberOfTrays).OutletLiquid.P),
				:(LiquidOutlet.z = TRAYS(NumberOfTrays).OutletLiquid.z),
				:(VapourOutlet.F = TRAYS(1).OutletVapour.F),
				:(VapourOutlet.T = TRAYS(1).OutletVapour.T),
				:(VapourOutlet.P = TRAYS(1).OutletVapour.P),
				:(VapourOutlet.z = TRAYS(1).OutletVapour.z),
				:(VapourDrawOff.F*VapSideTrayIndex= TRAYS.VapourSideStream.F),
				:(VapourDrawOff.T = TRAYS(VapourSideStreamLocation).VapourSideStream.T),
				:(VapourDrawOff.P = TRAYS(VapourSideStreamLocation).VapourSideStream.P),
				:(VapourDrawOff.z = TRAYS(VapourSideStreamLocation).VapourSideStream.z),
				:(VapourDrawOffFlow = VapourDrawOff.F),
				:(LiquidDrawOff.F*LiqSideTrayIndex= TRAYS.LiquidSideStream.F),
				:(LiquidDrawOff.T = TRAYS(LiquidSideStreamLocation).LiquidSideStream.T),
				:(LiquidDrawOff.P = TRAYS(LiquidSideStreamLocation).LiquidSideStream.P),
				:(LiquidDrawOff.z = TRAYS(LiquidSideStreamLocation).LiquidSideStream.z),
				:(LiquidDrawOffFlow = LiquidDrawOff.F),
			],
			[
				"Tray Temperature Indicator","Tray Pressure Indicator","","","","","","","","","","","","Liquid Inlet Flow","Liquid Inlet Temperature","Liquid Inlet Pressure","Liquid Inlet Composition","Liquid Inlet Vapour Fraction","Liquid Inlet Enthalpy","Vapour Inlet Flow","Vapour Inlet Temperature","Vapour Inlet Pressure","Vapour Inlet Composition","Vapour Inlet Vapour Fraction","Vapour Inlet Enthalpy","Liquid Outlet Flow","Liquid Outlet Temperature","Liquid Outlet Pressure","Liquid Outlet Composition","Vapour Outlet Flow","Vapour Outlet Temperature","Vapour Outlet Pressure","Vapour Outlet Composition","Vapour Draw Off Stream - Flow","Vapour Draw Off Stream - Temperature","Vapour Draw Off Stream - Pressure","Vapour Draw Off Stream - Composition","Vapour Draw Off Stream","Liquid Draw Off Stream - Flow","Liquid Draw Off Stream - Temperature","Liquid Draw Off Stream - Pressure","Liquid Draw Off Stream - Composition","Liquid Draw Off Stream",
			],
			[:PP,:NComp,:NumberOfFeeds,:FeedTrayLocation,:NumberOfTrays,:LiqSideTrayIndex,:VapSideTrayIndex,:LiquidSideStreamLocation,:VapourSideStreamLocation,:VapourFlowModel,:LiquidFlowModel,:TrayDiameter,:TraySpacing,:Fraction_HoleArea,:Fraction_DowncomerArea,:WeirLength,:WeirHeight,:TrayLiquidPasses,:HeatSupply,:AerationFraction,:DryPdropCoeff,:MurphreeEff,],
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	NumberOfFeeds::DanaInteger
	FeedTrayLocation::Array{DanaInteger}
	NumberOfTrays::DanaInteger
	LiqSideTrayIndex::Array{DanaInteger}
	VapSideTrayIndex::Array{DanaInteger}
	LiquidSideStreamLocation::DanaInteger
	VapourSideStreamLocation::DanaInteger
	VapourFlowModel::DanaSwitcher
	LiquidFlowModel::DanaSwitcher
	TrayDiameter::length
	TraySpacing::length
	Fraction_HoleArea::fraction
	Fraction_DowncomerArea::fraction
	WeirLength::length
	WeirHeight::length
	TrayLiquidPasses::positive
	HeatSupply::heat_rate
	AerationFraction::DanaReal
	DryPdropCoeff::DanaReal
	MurphreeEff::DanaReal
	initials::Array{Expr,1}
	initialNames::Array{String,1}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export SectionColumnBasic
function set(in::SectionColumnBasic)
	VapSideTrayIndex(VapourSideStreamLocation) =1
	 LiqSideTrayIndex(LiquidSideStreamLocation) =1
	 TRAYS.TrayDiameter_ = TrayDiameter
	 TRAYS.TraySpacing_ = TraySpacing
	 TRAYS.Fraction_HoleArea_ = Fraction_HoleArea
	 TRAYS.Fraction_DowncomerArea_ = Fraction_DowncomerArea
	 TRAYS.WeirLength_ = WeirLength
	 TRAYS.WeirHeight_ = WeirHeight
	 TRAYS.TrayLiquidPasses_ = TrayLiquidPasses
	 TRAYS.HeatSupply_ = HeatSupply
	 TRAYS.AerationFraction_ = AerationFraction
	 TRAYS.DryPdropCoeff_ = DryPdropCoeff
	 TRAYS.MurphreeEff_ = MurphreeEff

end
function initial(in::SectionColumnBasic)

end
function setEquationFlow(in::SectionColumnBasic)
	#TCI*'K' = TRAYS(min([NumberOfTrays, CONTROL.Tindicator_TrayNumber])).OutletVapour.T;
	addEquation(1)
	#PCI*'atm' = TRAYS(min([NumberOfTrays, CONTROL.Pindicator_TrayNumber])).OutletVapour.P;
	addEquation(2)
	let switch=LiquidFlowModel
		if switch=="Francis"
			addEquation(3)
		elseif switch=="Wang_Fl"
			addEquation(4)
		elseif switch=="Olsen"
			addEquation(5)
		elseif switch=="Feehery_Fl"
			addEquation(6)
		elseif switch=="Roffel_Fl"
			addEquation(7)
		end
	end
	let switch=VapourFlowModel
		if switch=="Reepmeyer"
			addEquation(8)
		elseif switch=="Feehery_Fv"
			addEquation(9)
		elseif switch=="Roffel_Fv"
			addEquation(10)
		elseif switch=="Klingberg"
			addEquation(11)
		elseif switch=="Wang_Fv"
			addEquation(12)
		elseif switch=="Elgue"
			addEquation(13)
		end
	end
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
	addEquation(19)
	addEquation(20)
	addEquation(21)
	addEquation(22)
	addEquation(23)
	addEquation(24)
	addEquation(25)
	addEquation(26)
	addEquation(27)
	addEquation(28)
	addEquation(29)
	addEquation(30)
	addEquation(31)
	addEquation(32)
	addEquation(33)
	addEquation(34)
	addEquation(35)
	addEquation(36)
	addEquation(37)
	addEquation(38)
	addEquation(39)
	addEquation(40)
	addEquation(41)
	addEquation(42)
	addEquation(43)
end
function atributes(in::SectionColumnBasic,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Model of a basic column section."
	fields[:Info]="Model of a basic column section containing a vetor of TRAYS numbered from the top-down."
	drive!(fields,_)
	return fields
end
SectionColumnBasic(_::Dict{Symbol,Any})=begin
	newModel=SectionColumnBasic()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# function circleshape(h, k, r)
#     theta = LinRange(0, 2*pi, 500)
#     return h .+ r*sin.(theta), k.+ r*cos.(theta)
# end

# function plotlayout!(p, turbinex, turbiney, rotordiameter; linecolor=:black, fillcolor=:blue, markeralpha=1, title="")
#     nturbines = length(turbinex)
#     for i in 1:nturbines
#         plot!(p, circleshape(turbinex[i],turbiney[i],rotordiameter[i]/2.0), seriestype=[:shape],
#             linecolor=linecolor, c=fillcolor, legend=false, aspect_ratio=1, alpha=markeralpha, title=title)
#     end
#     display(p)
# end
"""
    plotwindfarm!(ax, boundary_vertices, turbinex, turbiney, rotordiameter; nboundaries=1, aspect="equal", xlim=[], ylim=[], fill=false, color="k", markeralpha=1, title="")

    Convenience function for plotting wind farms

# Arguments
- `ax
- `boundary_vertices::Array{Float,1}(nvertices)`: an nx2 array of boundary vertices
- `turbinex::Array{Float,1}(nturbines)`: an array x coordinates of wind turbine locations
- `turbiney::Array{Float,1}(nturbines)`: an array y coordinates of wind turbine locations
- `rotordiameter::Array{Float,1}(nturbines)`: an array rotor diameters of wind turbines
- `nboundaries::Int`: number of discrete boundary regions
- `aspect::String`: set plot aspect ratio, default="equal"
- `xlim::Array`: limits in x coordinate. "[]" results in limits being automatically defined
- `ylim::Array`: limits in y coordinate. "[]" results in limits being automatically defined
- `fill::Bool`: determines whether turbine circle markers are filled or not
- `color::=String`: sets color for turbine markers
- `markeralpha::Int`: determines tranparancy of turbine markers
- `title::String`: optional title to include on the plot
"""
function plotwindfarm!(ax, turbinex, turbiney, rotordiameter; nboundaries=1,
    boundary_vertices=[], aspect="equal", xlim=[], ylim=[], fill=false, turbinecolor="k",
    boundarycolor="k", boundarylinestyle="-", turbinelinestyle="-", markeralpha=1, title="")

    # determine how many turbines are in the farm
    nturbines = length(turbinex)

    # add the wind turbines
    plotlayout!(ax, turbinex, turbiney, rotordiameter; aspect=aspect, fill=fill, color=turbinecolor, markeralpha=markeralpha, title=title, linestyle=turbinelinestyle)

    # add the bounary/ies
    if boundary_vertices != []
        plotboundary!(ax, boundary_vertices, nboundaries=nboundaries, color=boundarycolor, linestyle=boundarylinestyle)
    end

    # adjust plot x limits if not set
    if xlim !== nothing
        if xlim == [] && boundary_vertices != []
            if nboundaries == 1
                xlim = [minimum(boundary_vertices[:,1]) - rotordiameter[1], maximum(boundary_vertices[:,1]) + rotordiameter[1]]
            else
                for i = 1:nboundaries
                    xlim_tmp = [minimum(boundary_vertices[i][:,1]) - rotordiameter[1], maximum(boundary_vertices[i][:,1]) + rotordiameter[1]]
                    if i == 1
                        xlim = deepcopy(xlim_tmp)
                    end
                    xlim[1] = minimum([xlim[1], xlim_tmp[1]])
                    xlim[2] = maximum([xlim[2], xlim_tmp[2]])
                end
            end
        elseif boundary_vertices == []
                xlim = [minimum(turbinex)-sum(rotordiameter)/nturbines, maximum(turbinex)+sum(rotordiameter)/nturbines]
        end
        ax.set(xlim=xlim)
    end

    # adjust plot y limits if not set
    if ylim !== nothing
        if ylim == [] && boundary_vertices != []
            if nboundaries == 1
                ylim = [minimum(boundary_vertices[:,2]) - rotordiameter[1], maximum(boundary_vertices[:,2]) + rotordiameter[1]]
            else
                for i = 1:nboundaries
                    ylim_tmp = [minimum(boundary_vertices[i][:,2]) - rotordiameter[1], maximum(boundary_vertices[i][:,2]) + rotordiameter[1]]
                    if i == 1
                        ylim = deepcopy(ylim_tmp)
                    end
                    ylim[1] = minimum([ylim[1], ylim_tmp[1]])
                    ylim[2] = maximum([ylim[2], ylim_tmp[2]])
                end
            end
        elseif boundary_vertices == []
            ylim = [minimum(turbiney)-sum(rotordiameter)/nturbines, maximum(turbiney)+sum(rotordiameter)/nturbines]
        end
        ax.set(ylim=ylim)
    end

    # adjust aspect ratio
    ax.set(aspect=aspect)
end

"""
    plotlayout!(ax, turbinex, turbiney, rotordiameter; aspect="equal", xlim=[], ylim=[], fill=false, color="k", markeralpha=1, title="")

    Convenience function for plotting wind farm layouts

# Arguments
- `ax
- `turbinex::Array{Float,1}(nturbines)`: an array x coordinates of wind turbine locations
- `turbiney::Array{Float,1}(nturbines)`: an array y coordinates of wind turbine locations
- `rotordiameter::Array{Float,1}(nturbines)`: an array rotor diameters of wind turbines
- `aspect::String`: set plot aspect ratio, default="equal"
- `xlim::Array`: limits in x coordinate. "[]" results in limits being automatically defined
- `ylim::Array`: limits in y coordinate. "[]" results in limits being automatically defined
- `fill::Bool`: determines whether turbine circle markers are filled or not
- `color::=String`: sets color for turbine markers
- `markeralpha::Int`: determines tranparancy of turbine markers
- `itle::String`: optional title to include on the plot
"""
function plotlayout!(ax, turbinex, turbiney, rotordiameter; aspect="equal", fill=false, color="k", markeralpha=1, title="", linestyle="-")
    nturbines = length(turbinex)

    # add turbines
    for i in 1:nturbines
        circle = matplotlib.patches.Circle((turbinex[i], turbiney[i]), rotordiameter[i]/2.0, fill=fill, color=color, linestyle=linestyle)
        ax.add_patch(circle)
    end

end

"""
    plotsingleboundary!(ax, boundary_vertices; color="k", linestyle="-')

    Convenience function for plotting wind farm boundaries

# Arguments
- `ax
- `boundary_vertices::Array{Float,1}(nvertices)`: an nx2 array of boundary vertices for polygon or [[center_x, center_y], radius] for circle boundary
- `color::=String`: sets color for turbine markers
- `linestyle::String`: sets the line style for the boundary
"""
function plotsingleboundary!(ax, boundary_vertices; color="k", linestyle="-")

    # get the number of vertices defining the boundary
    nvertices = length(boundary_vertices[:,1])

    # add boundary
    if nvertices < 3
        # circle boundary
        farm_center = boundary_vertices[1]
        farm_radius = boundary_vertices[2]
        circle = matplotlib.patches.Circle((farm_center[1], farm_center[2]), farm_radius, color=color, linestyle=linestyle, fill=false)
        ax.add_patch(circle)
    else
        # polygon boundary
        x = push!(boundary_vertices[:,1], boundary_vertices[1,1])
        y = push!(boundary_vertices[:,2], boundary_vertices[1,2])
        ax.plot(x, y, color=color, linestyle=linestyle)
    end

end

"""
    plotboundary!(ax, boundary_vertices; color="k", linestyle="-", nboundaries=1)

    Convenience function for plotting wind farm boundaries

# Arguments
- `ax
- `boundary_vertices::Array{Float,1}(nvertices)`: an nx2 array of boundary vertices for polygon or [[center_x, center_y], radius] for circle boundary
- `color::=String`: sets color for turbine markers
- `linestyle::String`: sets the line style for the boundary
- `nboundaries::Int`: number of discrete boundary regions
"""
function plotboundary!(ax, boundary_vertices; color="k", linestyle="-", nboundaries=1)

    if nboundaries == 1
        plotsingleboundary!(ax, boundary_vertices, color=color, linestyle=linestyle)
    else
        for i = 1:nboundaries
            plotsingleboundary!(ax, boundary_vertices[i], color=color, linestyle=linestyle)
        end
    end

end

"""
    plotrotorsamplepoints!(ax, y, z; rotordiameter=2.0, aspect="equal", ylim=[], zlim=[], fill=false, color="k", markeralpha=1, title="")

    Convenience function for plotting where points are being sampled on the wind turbine rotor

# Arguments
- `ax
- `y::Array{Float,1}(nturbines)`: an array x coordinates of wind turbine locations
- `z::Array{Float,1}(nturbines)`: an array y coordinates of wind turbine locations
- `rotordiameter::Number`: rotor diameter of wind turbine
- `aspect::String`: set plot aspect ratio, default="equal"
- `ylim::Array`: limits in y coordinate. "[]" results in limits being automatically defined
- `zlim::Array`: limits in z coordinate. "[]" results in limits being automatically defined
- `fill::Bool`: determines whether turbine circle markers are filled or not
- `color::=String`: sets color for turbine markers
- `markeralpha::Int`: determines tranparancy of turbine markers between 0 (transparent) and 1 (opaque)
- `itle::String`: optional title to include on the plot
"""
function plotrotorsamplepoints!(ax, y, z; rotordiameter=2.0, aspect="equal", ylim=[], zlim=[], fill=false, color="k", markeralpha=1, title="")

    # set plot limits
    if ylim == []
        ylim = [-1.05, 1.05].*rotordiameter/2.0
    end
    if zlim == []
        zlim = [-1.05, 1.05].*rotordiameter/2.0
    end

    # add points
    ax.scatter(y, z, color=color, markeralpha)

    # add rotor swept area
    circle = matplotlib.patches.Circle((0.0, 0.0), rotordiameter/2.0, fill=fill, color=color)
    ax.add_patch(circle)
    println(ylim, zlim)
    ax.set(xlim=ylim, ylim=zlim, aspect=aspect, title=title)

end

"""
    plotwindresource!(ax::Array, windresource::ff.DiscretizedWindResource; roundingdigits=[1,3], fill=false, alpha=0.5, colors=["b", "b"], fontsize=8, edgecolor=nothing, rlabel_position=-45)

    Convenience function for visualizing the wind speed and wind frequency roses

# Arguments
- `ax::Array`: pre-initialized single dimension array of axes from pyplot with length at least 2
- `windresource::ff.DiscretizedWindResource`: wind rose information
- `roundingdigits::Array{Int, 1}`: how many significant digits to round to on each axis in ax
- `fill::Bool`: determines whether bars are filled or not
- `alpha::Number`: tranparancy of bars between 0 (transparent) and 1 (opaque)
- `colors::=Array{String, 1}`: sets color for turbine markers
- `fontsize::Int`: font size of text on figures
- `edgecolor`: color of edges of each bar in polar chart, nothing means no color
- `rlabel_position:Number`: Angle at which to draw the radial axes
"""
function plotwindresource!(windresource::ff.DiscretizedWindResource; ax::Array=[], roundingdigits=[1,3], fill=false, alpha=0.5, colors=["b", "b"], fontsize=8, edgecolor=nothing, rlabel_position=-45, titles=["Wind Speed", "Wind Probability"])

    axes_generated = false
    if length(ax) == 0
        fig, ax = plt.subplots(1, 2, subplot_kw=Dict("projection"=>"polar"))
        axes_generated = true
    end

    # extract wind resource elements for windrose plots
    d = windresource.wind_directions
    s = windresource.wind_speeds
    f = windresource.wind_probabilities

    # plot wind speed rose
    kwargs=(:edgecolor=>edgecolor, :alpha=>alpha, :color=>colors[1])
    plotwindrose!(ax[1], d, s, roundingdigit=roundingdigits[1], fontsize=fontsize, units="m/s", title=titles[1], kwargs=kwargs)

    # plot wind frequency rose
    kwargs=(:edgecolor=>edgecolor, :alpha=>alpha, :color=>colors[2])
    plotwindrose!(ax[2], d, f, roundingdigit=roundingdigits[2], fontsize=fontsize, units="%", title=titles[2], kwargs=kwargs)

    # return axis if newly generated
    axes_generated && return ax
end

"""
    plotwindrose!(ax, d, f; roundingdigit=1, color="C0",alpha=0.5,fontsize=8,
    dticks=(0,pi/4,pi/2,3*pi/4,pi,5*pi/4,3*pi/2,7*pi/4),
    dlabels=("E","NE","N","NW","W","SW","S","SW"),
    fticks=nothing, flabels=nothing, normalize=false, edgecolor=nothing, units="",
    rlabel_position=-45)

    Convenience function for creating a windrose from any polar data

# Arguments
- `ax::PyCall.PyObject`: pre-initialized axis from pyplot
- `d::Vector`: wind rose directions
- `f::Vector`: wind rose radial variable
- `roundingdigit::Int`: how many significant digits to round to
- `fontsize::Int`: font size of text on figures
- `dticks::Tuple`: contains angular tick locations
- `dlabels::Tuple`: contains angular tick labels
- `fticks::Tuple`: contains radial tick locations
- `flabels::Tuple`: contains radial tick labels
- `normalize::Bool`: choose whether or not to normalize by the sum of f
- `units::String`: Units to append to flabels
- `rlabel_position:Number`: Angle at which to draw the radial axis
- `plotcommand::String`: Type of plot. Can be ["bar", "plot"]
- `kwargs::Tuple`: tuple containing key word arguments to plotcommand in the form (key => "value")
"""
function plotwindrose!(ax, d, f; roundingdigit=1,fontsize=8,
    dticks=(0,pi/4,pi/2,3*pi/4,pi,5*pi/4,3*pi/2,7*pi/4),
    dlabels=("E","NE","N","NW","W","SW","S","SE"),
    fticks=nothing, flabels=nothing, normalize=false, units="",
    rlabel_position=-45, title="", plotcommand="bar", kwargs=(:edgecolor=>nothing, :alpha=>0.5, :color=>"C0"))

    # set up function ticks if not provided, scale and round appropriately
    if fticks === nothing
        fmin = floor(minimum(f), digits=roundingdigit)
        fmax = ceil(maximum(f), digits=roundingdigit)
        frange = fmax - fmin
        if frange == 0.0
            fticks = round.(collect((0:0.25:1.25).*fmax)[2:end-1], digits=roundingdigit)
        else
            fticks = round.(collect(fmin:frange/4.0:fmax)[2:end-1], digits=roundingdigit)
        end
    end
    # set up function labels if not provided
    if flabels === nothing
        flabels = string.(fticks).*" ".*units
    end

    # normalize if desired
    if normalize
        f = f./sum(f)
    end

    # adjust to radians if directions provided in degrees
    if maximum(d) > 10
        d = deg2rad.(d)
    end

    # get the number of wind directions
    ndirs = length(d)

    # plot wind rose
    if plotcommand == "bar"
        # specify bar width
        width = (2*pi/ndirs)*0.95
        ax.bar(pi/2 .-d, f; width, kwargs...)
    elseif plotcommand == "plot"
        ax.plot(pi/2 .-d, f; kwargs...)
    end

    # format polar plot
    ax.set_xticks(dticks)
    ax.set_xticklabels(dlabels,fontsize=fontsize)
    ax.set_rgrids(fticks,flabels,angle=rlabel_position,fontsize=fontsize, zorder=0)
    for tick in ax.yaxis.get_majorticklabels()
        tick.set_horizontalalignment("center")
    end
    ax.set_title(title, y=-0.25,fontsize=fontsize)

end

"""
    add_turbine!(ax; view="side", hubdiameter=0.1, hubheight=0.9, radius=0.5, chord=0.1,
        nacellewidth=0.3, nacelleheight=0.1, towerbottomdiam=0.1, towertopdiam=0.05,
        overhang=0.05, s=5)

    Convenience function for adding wind turbines to plots.

# Arguments
- `ax::PyCall.PyObject`: pre-initialized axis from pyplot
- `view::Number`: determines which turbine view to use "top" or "side" (default)
- `hubdiameter::Number`: hub diameter in axis coordinate frame
- `hubheight::Number`: hub height in axis coordinate frame
- `radius::Number`: full rotor radius in axis coordinate frame
- `chord::Number`: maximum chord in axis coordinate frame
- `nacellewidth::Number`: nacelle width in axis coordinate frame
- `nacelleheight::Number`: nacelle height in axis coordinate frame
- `towerbottomdiam::Number`: tower bottom diameter in axis coordinate frame
- `towertopdiam::Number`: tower top diameter in axis coordinate frame
- `overhang::Number`: overhang (distance from blade attachment to tower bottom in x axis) in axis coordinate frame
- `s::Number`: scales overhang and tower location in x direction to work with condensed x axis as in long contour plots
"""
function add_turbine!(ax; view="side", hubdiameter=0.1, hubheight=0.9, radius=0.5, chord=0.1, nacellewidth=0.3, nacelleheight=0.1, towerbottomdiam=0.1, towertopdiam=0.05, overhang=0.05, s=5, color="k")

    if view == "side"

        # create blade patches
        blade1 = plt.matplotlib.patches.Ellipse((0,hubheight+radius/2),chord, 0.5, color=color)
        blade2 = plt.matplotlib.patches.Ellipse((0,hubheight-radius/2),chord, 0.5, color=color)

        # create hub and nacelle patches
        hub = plt.matplotlib.patches.Ellipse((0,hubheight),3*hubdiameter, hubdiameter, color=color)
        nacelle = plt.matplotlib.patches.Rectangle((0,hubheight-hubdiameter/2),nacellewidth, nacelleheight, color=color)

        # scale if desired
        towerbottomdiam *= s
        towertopdiam *= s
        overhang *= s

        # get difference between top and bottom tower diameters
        ddiff = abs(towertopdiam-towerbottomdiam)

        # calculate polygon x points for tower
        p1x = overhang
        p2x = overhang+ddiff/2
        p3x = p2x + towertopdiam
        p4x = p1x + towerbottomdiam

        # create tower patch
        tower = plt.matplotlib.patches.Polygon([[p1x, 0.0],[p2x, hubheight],[p3x, hubheight],[p4x, 0.0]], closed=true, color=color)

        # add patches to axis
        ax.add_patch(blade1)
        ax.add_patch(blade2)
        ax.add_patch(hub)
        ax.add_patch(nacelle)
        ax.add_patch(tower)

    elseif view == "top"

        # create blade patches
        blade1 = plt.matplotlib.patches.Ellipse((0,radius/2),chord, 0.5, color=color)
        blade2 = plt.matplotlib.patches.Ellipse((0,radius/2),chord, 0.5, color=color)

        # create hub and nacelle patches
        hub = plt.matplotlib.patches.Ellipse((0.0,0.0),3*hubdiameter, hubdiameter, color=color)
        nacelle = plt.matplotlib.patches.Rectangle((0,-hubdiameter/2),nacellewidth, nacelleheight, color=color)

        # add patches to axis
        ax.add_patch(blade1)
        ax.add_patch(blade2)
        ax.add_patch(hub)
        ax.add_patch(nacelle)

    end
end


text:
mutable struct VecJacOperator{T,F,T1,T2,uType,P,tType,O} <:
               DiffEqBase.AbstractDiffEqLinearOperator{T}
    f::F
    cache1::T1
    cache2::T2
    u::uType
    p::P
    t::tType
    autodiff::Bool
    ishermitian::Bool
    opnorm::O

    function VecJacOperator{T}(
        f,
        p = nothing,
        t::Union{Nothing,Number} = nothing;
        autodiff = true,
        ishermitian = false,
        opnorm = true,
    ) where {T}
        p === nothing ? P = Any : P = typeof(p)
        t === nothing ? tType = Any : tType = typeof(t)
        new{T,typeof(f),Nothing,Nothing,Any,P,tType,typeof(opnorm)}(
            f,
            nothing,
            nothing,
            nothing,
            nothing,
            nothing,
            autodiff,
            ishermitian,
        )
    end
    function VecJacOperator{T}(
        f,
        u::AbstractArray,
        p = nothing,
        t::Union{Nothing,Number} = nothing;
        autodiff = true,
        ishermitian = false,
        opnorm = true,
    ) where {T}
        cache1 = similar(u)
        cache2 = similar(u)
        p === nothing ? P = Any : P = typeof(p)
        t === nothing ? tType = Any : tType = typeof(t)
        new{
            T,
            typeof(f),
            typeof(cache1),
            typeof(cache2),
            typeof(u),
            P,
            tType,
            typeof(opnorm),
        }(
            f,
            cache1,
            cache2,
            u,
            p,
            t,
            autodiff,
            ishermitian,
            opnorm,
        )
    end
    function VecJacOperator(f, u, args...; kwargs...)
        VecJacOperator{eltype(u)}(f, u, args...; kwargs...)
    end
end

LinearAlgebra.opnorm(L::VecJacOperator, p::Real = 2) = L.opnorm
LinearAlgebra.ishermitian(L::VecJacOperator) = L.ishermitian

Base.size(L::VecJacOperator) = (length(L.cache1), length(L.cache1))
Base.size(L::VecJacOperator, i::Int) = length(L.cache1)
function update_coefficients!(L::VecJacOperator, u, p, t)
    L.u = u
    L.p = p
    L.t = t
    !L.autodiff && L.cache1 !== nothing && L.f(L.cache1, L.u, L.p, L.t)
end

# Interpret the call as df/du' * u
function (L::VecJacOperator)(u, p, t::Number)
    update_coefficients!(L, u, p, t)
    L * u
end

function (L::VecJacOperator)(du, u, p, t::Number)
    update_coefficients!(L, u, p, t)
    mul!(du, L, u)
end

function Base.:*(L::VecJacOperator, x::AbstractVector)
    if hasmethod(L.f, typeof.((L.u, L.p, L.t)))
        return vec(
            L.autodiff ? auto_vecjac(_u -> L.f(_u, L.p, L.t), L.u, x) :
            num_vecjac(_u -> L.f(_u, L.p, L.t), L.u, x),
        )
    end
    return mul!(similar(vec(L.u)), L, x)
end

function LinearAlgebra.mul!(
    du::AbstractVector,
    L::VecJacOperator,
    x::AbstractVector,
)
    du = reshape(du, size(L.u))
    let p = L.p, t = L.t
        if L.cache1 === nothing
            if L.autodiff
                # For autodiff prefer non-inplace function
                if hasmethod(L.f, typeof.((L.u, L.p, L.t)))
                    auto_vecjac!(du, _u -> L.f(_u, p, t), L.u, x)
                else
                    auto_vecjac!(du, (_du, _u) -> L.f(_du, _u, p, t), L.u, x)
                end
            else
                if hasmethod(L.f, typeof.((du, L.u, L.p, L.t)))
                    num_vecjac!(
                        du,
                        (_du, _u) -> L.f(_du, _u, p, t),
                        L.u,
                        x;
                        compute_f0 = true,
                    )
                else
                    num_vecjac!(
                        du,
                        _u -> L.f(_u, p, t),
                        L.u,
                        x;
                        compute_f0 = true,
                    )
                end
            end
        else
            if L.autodiff
                # For autodiff prefer non-inplace function
                if hasmethod(L.f, typeof.((L.u, L.p, L.t)))
                    auto_vecjac!(
                        du,
                        _u -> L.f(_u, p, t),
                        L.u,
                        x,
                        L.cache1,
                        L.cache2,
                    )
                else
                    auto_vecjac!(
                        du,
                        (_du, _u) -> L.f(_du, _u, p, t),
                        L.u,
                        x,
                        L.cache1,
                        L.cache2,
                    )
                end
            else
                if hasmethod(L.f, typeof.((du, L.u, L.p, L.t)))
                    num_vecjac!(
                        du,
                        (_du, _u) -> L.f(_du, _u, p, t),
                        L.u,
                        x,
                        L.cache1,
                        L.cache2;
                        compute_f0 = true,
                    )
                else
                    num_vecjac!(
                        du,
                        _u -> L.f(_u, p, t),
                        L.u,
                        x,
                        L.cache1,
                        L.cache2;
                        compute_f0 = true,
                    )
                end
            end
        end
    end
    return vec(du)
end

function Base.resize!(J::VecJacOperator, i)
    resize!(J.cache1, i)
    resize!(J.cache2, i)
end


text:
module SlackThreads

using JSON3
using StructTypes
using FileIO
using Mocking

export AbstractSlackThread, SlackThread, slack_log_exception
export DummyThread, SlackCallRecord

abstract type AbstractSlackThread end

mutable struct SlackThread <: AbstractSlackThread
    channel::Union{String,Nothing}
    ts::Union{String,Nothing} # In Slack terminology, `ts_thread`: the ID of another un-threaded message to reply to (https://api.slack.com/reference/messaging/payload)
end

function Base.copyto!(thread::SlackThread, obj)
    thread.channel = obj.channel
    thread.ts = obj.ts
    return thread
end

StructTypes.StructType(::Type{SlackThread}) = StructTypes.Struct()

const CATCH_EXCEPTIONS = Ref(true)

# We turn off exception handling for our tests, to ensure we aren't throwing exceptions
# that we're missing. But we have it on by default, since in ordinary usage we want to
# be sure we are catching all exceptions.
# We wrap all public API methods in this, which should make it very difficult to throw
# an exception.
macro maybecatch(expr, log_str)
    quote
        if CATCH_EXCEPTIONS[]
            try
                $(esc(expr))
            catch e
                @error $(log_str) exception = (e, catch_backtrace())
                nothing
            end
        else
            let # introduce a local scope like `try` does
                $(esc(expr))
            end
        end
    end
end

"""
    SlackThread(channel=get(ENV, "SLACK_CHANNEL", nothing))

Constructs a `SlackThread`. A channel should be specified by it's ID (a number
like `C1H9RESGL` at the bottom of the "About" section of the channel).
"""
function SlackThread(channel=get(ENV, "SLACK_CHANNEL", nothing))
    thread = @maybecatch begin
        if channel === nothing
            #TODO: Use Preferences.jl for default channel?
            @warn "`channel` not passed, nor `SLACK_CHANNEL` environmental variable set, so will only emit logging statements."
        end
        SlackThread(channel, nothing)
    end "Error when constructing `SlackThread`"
    thread === nothing && return SlackThread(nothing, nothing)
    return thread
end

function format_slack_link(uri, msg=nothing)
    if msg === nothing
        return "<$(uri)>"
    else
        return "<$(uri)|$(msg)>"
    end
end

"""
    (thread::SlackThread)(text::AbstractString, uploads...;
                          combine_texts=combine_texts, options...)

Sends a message to the Slack thread with the contents `text`. If this is the
first message sent by `thread`, this starts a new thread (in `thread.channel`),
otherwise it updates the existing thread.

You can also include plots, images, and files by passing file paths or lists of
pairs as `uploads`.

Returns:

* If the request is successful, returns Slack's response for the message as a
  `JSON3.Object`. If more than one upload is present, or `options` are passed,
  the response will be for a text-only request, since the file uploads will
  be processed separately (using the strategy from
  https://stackoverflow.com/a/63391026/12486544).
* If the request is not successful, returns `nothing`.

## Uploads

Each item in `uploads` may be:

* a `Pair` of the form `name_with_extension::AbstractString => object`,
* or a path to a file (i.e. anything that supports `read` and `basename`)

Valid `object`s are:

* a vector of bytes (`Vector{UInt8}`)
* a string
* an object supporting `FileIO.save`

!!! note
    When using the pair syntax, including a file extension in the name
    helps Slack choose how to display the object, and helps FileIO choose how to
    save the object. E.g. `"my_plot.png"` instead of `"my_plot"`.

## Logging

* Emits `@debug` logs when sending requests and recieving responses from the Slack API.
* Emits `@warn` logs with the contents of requests to the Slack API when the channel or token is not configured correctly (in lieu of sending a request)
* Emits `@error` logs when an exception is encountered or Slack returns an error response.

## Message splitting

By default, [`SlackThreads.combine_texts`](@ref) is used to combine the message text and attachment link text
into a set of messages, which returns multiple messages in the case of many attachments.
One may pass any function to the `combine_texts` keyword argument which accepts and returns a vector of strings,
for example, `texts -> SlackThreads.combine_texts(texts; max_length=100, message_count_suffix=(i, n) -> "")`
to split messages after 100 characters, and to not add a `[\$i / \$n]` suffix to the messages.

## Options

One may pass any optional arguments supported by [`chat.postMessage`](https://api.slack.com/methods/chat.postMessage#args),
e.g. `link_names = true`, as keyword arguments.
"""
function (thread::SlackThread)(text::AbstractString, uploads...;
                               combine_texts=combine_texts, options...)
    return @maybecatch begin
        if thread.channel === nothing
            @warn "No Slack channel configured; message not sent." text uploads
            return nothing
        end
        if isempty(uploads)
            # send directly
            return send_message(thread, text; options...)
        end
        mktempdir() do dir
            if length(uploads) == 1 && thread.ts !== nothing && isempty(options)
                # special case: upload directly to thread
                # cannot do this if we don't have a `ts` already because
                # the response doesn't give us a `ts` to use.
                # i.e. we can post the message but don't have the `ts`
                # to thread from it. So in that case, we fallback
                # to the general case. We also can't do this if the user has passed
                # any options, since those are likely only valid for `chat.postMessage`.
                extra_args = ["-F", "initial_comment=$(text)", "-F",
                              "channels=$(thread.channel)", "-F", "thread_ts=$(thread.ts)"]
                return upload_file(local_file(only(uploads); dir); extra_args)
            end

            texts = String[text]

            for item in uploads
                r = upload_file(local_file(item; dir))
                r === nothing && continue
                push!(texts, format_slack_link(r.file.permalink, " "))
            end

            messages = combine_texts(texts)
            local r
            for msg in messages
                r = send_message(thread, msg; options...)
            end
            return r # return last response
        end
    end "Error when attempting to send message to Slack thread"
end

include("slack_api.jl")
include("slack_log_exception.jl")
include("utilities.jl")
include("dummy_slack_thread.jl")

end # module


text:
module Ex7_4c
using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

# Example 7.4c - Nullspace-based synthesis

println("Example 7.4c")

# define the state-space realizations of Gu, Gd and Gf
n = 5; p = 3; mu = 3; md = 1; mf = 5;       # enter dimensions
# define matrices of the state space realization
A = [ 0 0 1.132 0 -1;
0 -0.0538 -0.1712 0 0.0705;
0 0 0 1 0;
0 0.0485 0 -0.8556 -1.013;
0 -0.2909 0 1.0532 -0.6859];
Bu = [0 0 0;-0.12 1 0;0 0 0;4.419 0 -1.665;1.575 0 -0.0732];
Bd = Bu[:,mu]; Bf = [zeros(n,p) Bu[:,1:mu-1]];
C=eye(p,n); Du=zeros(p,mu); Dd=zeros(p,md); Df=eye(p,mf);
sys = dss(A,[Bu Bd Bf],C,[Du Dd Df]);        # define system

# setup the synthesis model
sysf = fdimodset(dss(A,[Bu Bd Bf],C,[Du Dd Df]),c =1:mu, d = mu .+ (1:md), f = (mu+md) .+ (1:mf));

# call of EFDSYN with the options for stability degree -1 and the synthesis
# of a scalar output filter using two design matrices
Q1, Rf1, info1 = efdsyn(sysf, HDesign = [0 1]; sdeg = -1, smarg = -1, rdim = 1);
Q2, Rf2, info2 = efdsyn(sysf, HDesign = [-1 1]; poles = [-1.0000 + 0.3992im, -1.0000 - 0.3992im], sdeg = -1, smarg = -1, rdim = 1);

# check admissibility, i.e., finite reciprocal sensitivity conditions
rscond1 = 1/fdiscond(Rf1,1)[1]
rscond2 = 1/fdiscond(Rf2,1)[1]
@test !isinf(rscond1) && !isinf(rscond1)
println("rscond1 = $rscond1 ")
println("rscond2 = $rscond2")

end # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
include("heat3d_model.jl")

# computes out <- exp(A * NSTEPS * dt) * b
function _expv(A, b, NSTEPS, dt; hermitian=false, m=min(30, size(A, 1)), tol=1e-7)

    # initialization of the krylov subspace
    TA, Tb = eltype(A), eltype(b)
    T = promote_type(TA, Tb)
    Ks = KrylovSubspace{T, real(T)}(length(b), m)
    arnoldi!(Ks, A, b; m=m, ishermitian=hermitian, tol=tol)

    out = similar(b)
    expv!(out, NSTEPS*dt, Ks)

    return out
end

function _phiv(A, b, NSTEPS, dt; hermitian=false, m=min(30, size(A, 1)), tol=1e-7)

    # initialization of the krylov subspace
    TA, Tb = eltype(A), eltype(b)
    T = promote_type(TA, Tb)
    Ks = KrylovSubspace{T, real(T)}(length(b), m)
    arnoldi!(Ks, A, b; m=m, ishermitian=hermitian, tol=tol)

    out = Matrix{Float64}(undef, size(A, 1), 3)
    phiv!(out, NSTEPS*dt, Ks, 2)

    return view(out, :, 3) .* dt^2
end

function _support_function_Phi_dir(ℓ, A, X0, δ; m=30, tol=1e-7, hermitian=false)
    At = transpose(A)
    vec = _expv(At, ℓ, 1, δ, m=m, tol=tol, hermitian=hermitian)
    ρ(vec, X0)
end

function _support_function_Eplus(ℓ::AbstractVector{N}, A::AbstractMatrix{N}, X0::Hyperrectangle, δ::Float64; cutoff=eps(N), m=min(30, size(A, 1)), tol=1e-7, hermitian=false) where {N}
    c = center(X0)
    r = radius_hyperrectangle(X0)
    A2 = A^2
    rin = abs.(A2 * c) + abs.(A2) * r
    rin = Vector(rin)
    Aabs = abs.(A)
    Pv = _phiv(Aabs, rin, 1, δ; m=m, tol=tol, hermitian=hermitian)
    return dot(abs.(ℓ), Pv)
end

# Non-Krylov version
function forward_non_krylov(instance=heat01; δ=0.02)
    A, Aᵀδ, X0, ℓ = instance(δ=δ)
    prob = @ivp(x' = A*x, x(0) ∈ X0)
    D = discretize(normalize(prob), δ, Forward())
    Ω0 = initial_state(D)
    res = ρ(ℓ, Ω0)
end

# Krylov version of Forward:
# ρ(ℓ, CH(X0, Φ*X0 ⊕ E⨥)) = ρ(ℓ, CH(X0, Φ*X0 ⊕ E⨥)) = max(ρ(ℓ, X0), ρ(Φ^T * ℓ, X0) + ρ(ℓ, E⨥))
function forward_krylov(instance=heat01; δ=0.02, m=30, tol=1e-7, hermitian=false)
    A, Aᵀδ, X0, ℓ = instance(δ=δ)
    aux1 = ρ(ℓ, X0)
    aux2 = _support_function_Phi_dir(ℓ, A, X0, δ, m=m, tol=tol, hermitian=hermitian)
    aux3 = _support_function_Eplus(ℓ, A, X0, δ, m=m, tol=tol, hermitian=hermitian)
    res = max(aux1, aux2 + aux3)
end

# Let Ω0 = CH(X0, Φ*X0 ⊕ E⨥), then we are interested in ρ(ℓ, Φ^k Ω0) for k = 0, 1, ..., NSTEPS-1
function forward_krylov_propagation(instance=heat01; NSTEPS=1, δ=0.02, m=30, tol=1e-7, hermitian=false)
    A, Aᵀδ, X0, ℓ = instance(δ=δ)
    At = copy(transpose(A))

    # initialization of the krylov subspace
    TA, Tb = eltype(At), eltype(ℓ)
    T = promote_type(TA, Tb)
    Ks = KrylovSubspace{T, real(T)}(length(ℓ), m)
    arnoldi!(Ks, At, ℓ; m=m, ishermitian=hermitian, tol=tol)

    ℓ′ = copy(ℓ)
    out = Vector{Float64}()
    for k in 0:NSTEPS-1
        if k > 0
            # ℓ′ = _expv(At, ℓ, k, δ; hermitian=false, m=min(30, size(A, 1)), tol=1e-7)

            # more efficient: use always the same precomputed Krylov subspace
            expv!(ℓ′, k*δ, Ks)
        end
        aux1 = ρ(ℓ′, X0)
        aux2 = _support_function_Phi_dir(ℓ′, A, X0, δ, m=m, tol=tol, hermitian=hermitian)
        aux3 = _support_function_Eplus(ℓ′, A, X0, δ, m=m, tol=tol, hermitian=hermitian)
        res = max(aux1, aux2 + aux3)
        push!(out, res)
    end
    out
end

function run_benchmarks()
    filename = "experiment_krylov_table_1.dat"
    open(filename, "w") do f
        forward_non_krylov(heat01) # warm-up runs
        forward_krylov(heat01)

        # ------------------------------------------------------------------------------
        write(f, "Case Heat01: 5x5x5 grid\n")
        t = @elapsed begin
                res = forward_non_krylov(heat01)
            end
        t = round(t, digits=3)
        write(f, "    Forward non-krylov: runtime: $t, result: $res\n")

        t = @elapsed begin
                res = forward_krylov(heat01, m=30, tol=1e-10, hermitian=true)
            end
        t = round(t, digits=3)
        write(f, "    Forward krylov: runtime: $t, result: $res\n")

        # ------------------------------------------------------------------------------
        write(f, "Case Heat02: 10x10x10 grid\n")
        t = @elapsed begin
                res = forward_non_krylov(heat02)
            end
        t = round(t, digits=3)
        write(f, "    Forward non-krylov: runtime: $t, result: $res\n")

        t = @elapsed begin
                res = forward_krylov(heat02, m=100, tol=1e-10, hermitian=true)
            end
        t = round(t, digits=3)
        write(f,"    Forward krylov: runtime: $t, result: $res\n")

        # ------------------------------------------------------------------------------
        write(f, "Case Heat03: 20x20x20 grid\n")
        write(f, "    Forward non-krylov: runtime: N/A, result: N/A\n")

        t = @elapsed begin
                res = forward_krylov(heat03, m=100, tol=1e-10, hermitian=true)
            end
        t = round(t, digits=3)
        write(f, "    Forward krylov: runtime: $t, result: $res\n")

        # ------------------------------------------------------------------------------
        write(f, "Case Heat04: 50x50x50 grid\n")
        write(f, "    Forward non-krylov: runtime: N/A, result: N/A\n")

        t = @elapsed begin
                res = forward_krylov(heat04, m=100, tol=1e-10, hermitian=true)
            end
        t = round(t, digits=3)
        write(f, "    Forward krylov: runtime: $t, result: $res\n")
    end
end

println("Running Heat3D model ...")
run_benchmarks()

nothing


text:
@testset "exports" begin
    # Ensure that every exported name is actually defined
    for name in names(Tri_Tet_Intersections)
        @test isdefined(Tri_Tet_Intersections, name)
        !isdefined(Tri_Tet_Intersections, name) && println("MISSING: ", name)
    end
end


text:
using Printf

const DEFAULT_ARMIJO_α₀ = 1.0
const DEFAULT_ARMIJO_σ₀ = 0.1
const DEFAULT_ARMIJO_σ₁ = 0.5
const DEFAULT_ARMIJO_p  = 0.5

struct BacktrackingState{OBJ,OPT,T} <: LinesearchState where {OBJ <: AbstractObjective, OPT <: Options, T <: Number}
    objective::OBJ
    config::OPT

    α₀::T
    ϵ::T
    p::T

    function BacktrackingState(objective; config = Options(),
                    α₀::T = DEFAULT_ARMIJO_α₀,
                    ϵ::T = DEFAULT_WOLFE_ϵ,
                    p::T = DEFAULT_ARMIJO_p) where {T}
        new{typeof(objective), typeof(config), T}(objective, config, α₀, ϵ, p)
    end
end

function BacktrackingState(F::Callable, x::Number = DEFAULT_ARMIJO_α₀; D = nothing, kwargs...)
    objective = UnivariateObjective(F, D, x)
    BacktrackingState(objective; kwargs...)
end

function BacktrackingState(F::Callable, x::AbstractVector; D = nothing, kwargs...)
    objective = MultivariateObjective(F, D, x)
    BacktrackingState(objective; kwargs...)
end

Base.show(io::IO, ls::BacktrackingState) = print(io, "Backtracking")

LinesearchState(algorithm::Backtracking, objective; kwargs...) = BacktrackingState(objective; kwargs...)

function (ls::BacktrackingState{<:UnivariateObjective})()
    local α = ls.α₀
    local y₀ = value!(ls.objective, zero(α))
    local d₀ = derivative!(ls.objective, zero(α))

    for _ in 1:ls.config.max_iterations
        if value!(ls.objective, α) ≥ y₀ + ls.ϵ * α * d₀
            α *= ls.p
        else
            break
        end
    end

    return α
end

backtracking(f, x; kwargs...) = BacktrackingState(f, x; kwargs...)()


text:
# Learning Rules.
# We are creating new data types here to hold each learning rule as a concept.
struct hebb end
struct hebbdecay end
struct oja end
struct willshaw end
struct antihebbdecay end

# Implementing different learning methods for different learning rule
# Here we are making use of the multiple dispatch feature in Julia
# Multiple Dispatch provides an unified framework for implementing different methods
learning(X,W,γ,r, rule::hebb;δ) = W + (γ .* X .* r)
learning(X,W,γ,r,rule::hebbdecay;δ) = W + (γ .* X .* r) - (γ .* δ .* W)
learning(X,W,γ,r,rule::oja;δ) = W + (γ .* X .* r) - (γ .* r^2 .* W)
learning(X,W,γ,r,rule::willshaw;δ) = W + (γ .* X .* r) .* (X .- W)
learning(X,W,γ,r,rule::antihebbdecay;δ) = W - (γ .* X .* r) - (γ * δ .* (W .- 1))

"""
USAGE:

`hebbianlearn(X,W,b, γ, r, hebb())` for hebbian learning \n
`hebbianlearn(X,W,b, γ, r, hebbdecay(), δ = delata_val)` for hebbian decay learning rule \n
`hebbianlearn(X,W,b, γ, r, oja(), δ = delata_val)` for oja learning rule \n
`hebbianlearn(X,W,b, γ, r, willshaw(), δ = delata_val)` for willshaw learning rule \n
`hebbianlearn(X,W,b, γ, r, antihebbdecay(), δ = delata_val)` for antihebbdecay learning rule

"""
function hebbianlearn(X, W, b, γ, r, rule;δ=0)
    Wb = [W b]
    Xnew = [X ones(size(X,1))]
    Wb = learning(Xnew,Wb, γ,r,rule, δ=δ)
    return Wb[:,1:size(W,2)], Wb[:,size(W,2)+1] # Unpacking W and b from Wb
end

hebbiandemo

max_stars_repo_name: coinslab/ComputationalCognitiveModeling

max_issues_repo_name: coinslab/ComputationalCognitiveModeling

max_forks_repo_name: coinslab/ComputationalCognitiveModeling


text:
module Semantics

export SemanticText

"""
An object representing a string of text but with additional semantic
information. These objects convert to `String`s through the `string` function,
but also typically support other operations.
"""
abstract type SemanticText end

end

using .Semantics: SemanticText
export SemanticText


text:
export link!

function link!(dst::Module, src::Module)
    status = convert(Core.Bool, API.LLVMLinkModules2(ref(dst), ref(src)))
    @assert !status # caught by diagnostics handler

    return nothing
end


text:
using Displaz
using Test

function displaz_args(f::Function)
    old_binding = Displaz.run_displaz
    try
        arglist = []
        capture_args = (args)->push!(arglist, args)
        Displaz.eval(:(run_displaz = $capture_args))
        f()
        @test length(arglist) == 1
        return arglist[1]
    finally
        Displaz.eval(:(run_displaz = $old_binding))
    end
end

macro displaz_args(ex)
    quote
        displaz_args(()->$(esc(ex)))
    end
end

@testset "clearplot" begin
    @test @displaz_args(clearplot()) == `-script -server default -clear`
end

@testset "viewplot" begin
    @test @displaz_args(viewplot(center=[1,2,3])) ==
          `-script -server default -viewposition 1 2 3`
    @test @displaz_args(viewplot(center=[1.1,2.2,3.3])) ==
          `-script -server default -viewposition 1.1 2.2 3.3`

    @test @displaz_args(viewplot(radius=10)) ==
          `-script -server default -viewradius 10`
    @test @displaz_args(viewplot(radius=10.1)) ==
          `-script -server default -viewradius 10.1`

    @test @displaz_args(viewplot(rotation=(1,2,3))) ==
          `-script -server default -viewangles 1 2 3`
    @test @displaz_args(viewplot(rotation=[0 1 0; 1 0 0; 0 0 -1])) ==
          `-script -server default -viewrotation 0 1 0 1 0 0 0 0 -1`

    # All args together
    @test @displaz_args(viewplot(center=[1,2,3], radius=4, rotation=(5,6,7))) ==
          `-script -server default -viewposition 1 2 3 -viewradius 4 -viewangles 5 6 7`
end

@testset "annotation" begin
    @test @displaz_args(annotation([0, 1, 2], "Hello world")) == `-script -server default -annotation "Hello world" 0 1 2 -label "Hello world"`
    @test @displaz_args(annotation([0, 1, 2], "Hello world", "A label")) == `-script -server default -annotation "Hello world" 0 1 2 -label "A label"`
end

@testset "hook" begin

    # first make sure the command is OK:
    command = Displaz.hook_command("someserver",
                                KeyEvent("c")=>Nothing,
                                KeyEvent("p")=>CursorPosition)

    @test command == `$(Displaz._displaz_cmd) -server someserver -hook key:c null -hook key:p cursor`

    # prepare some fake Displaz outputs:
    lines = ["key:c null",
             "key:p cursor 0 0 0"]
    event_stream = IOBuffer()
    foreach(l -> write(event_stream, l * "\n"), lines)
    seek(event_stream, 0)

    # prepare our received events buffer and a callback function:
    received = Vector{Any}()
    callback(x...) = push!(received, x)

    # this should plow through our fake output lines and return:
    Displaz.handle_events(callback, event_stream)

    # see if we got what we expected:
    @test length(received) == 2
    @test received[1] == (KeyEvent("c"), nothing)
    @test received[2] == (KeyEvent("p"), CursorPosition([0.0, 0.0, 0.0]))

end


text:
module exdiff
using SDE
using Diffusion
	println("euler")
	println(quvar(euler(0, 1, (t,x)-> -5x, (t,x) -> 1, diff(linspace(0., 1.,1000)))))

	dt = diff(linspace(0., 0.5, 100))
	dw = dW1(dt)

	println("Brownian motion with drift")
	println(" sint(4 .. dt) + ito (8 .. dw)")
	println((ito(4 .. dt) + ito(8 .. dw))[1:10])

	println(bracket(brown1(0, 2, 100))[1:10])

	#%discretization error of quadratic variation ?
	tic()
	println("Quadratic variation of B(0<t<2) = 2 \u2248 ", mean([quvar(brown1(0, 2, 5)) for i in 1:100000]))
	toc()

	tic()
	dt2 = diff(linspace(0., 2,5))
	println("Quadratic variation of B(0<t<2) = 2 \u2248 ", mean([quvar(ito(dW1(dt2))) for i in 1:100000]))
	toc()

	println("W_aug")
	adw = aug(dw,dt,10)
	println(adw[1:10])
	println("W(t) - W_aug(t):", sum(dw) - sum(adw))
	println("quvar(W, W_aug):", quvar(ito(1 .. dw)), " ", quvar(ito(1 .. adw)))

	dx = dW1(5,3)
	println(ito(dx .. dx)[end])
end


text:
using HTTP, HTTP.IOExtras, JSON, MusicProcessing
HTTP.open("POST", "http://music.com/play") do io
    write(io, JSON.json([
        "auth" => "12345XXXX",
        "song_id" => 7,
    ]))
    r = startread(io)
    @show r.status
    while !eof(io)
        bytes = readavailable(io)
        play(bytes)
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using StatsBase, Turing, ContinuousBenchmarks
using Gadfly
using Mamba: describe

# include("ASCIIPlot.jl")
import Gadfly.ElementOrFunction

# First add a method to the basic Gadfly.plot function for QQPair types (generated by Distributions.qqbuild())
Gadfly.plot(qq::QQPair, elements::ElementOrFunction...) = Gadfly.plot(x=qq.qx, y=qq.qy, Geom.point, Theme(highlight_width=0px), elements...)

# Now some shorthand functions
qqplot(x, y, elements::ElementOrFunction...) = Gadfly.plot(qqbuild(x, y), elements...)
qqnorm(x, elements::ElementOrFunction...) = qqplot(Normal(), x, Guide.xlabel("Theoretical Normal quantiles"), Guide.ylabel("Observed quantiles"), elements...)

NSamples = 30000

@model gdemo2(x, bkstep) = begin
    if bkstep == false
        # Forward Sample
        s ~ InverseGamma(2,3)
        m ~ Normal(0,sqrt(s))
        y ~ MvNormal([m; m], [sqrt(s) 0; 0 sqrt(s)])
    elseif bkstep == true
        # Backward Step 1: theta ~ theta | x
        s ~ InverseGamma(2,3)
        m ~ Normal(0,sqrt(s))
        x ~ MvNormal([m; m], [sqrt(s) 0; 0 sqrt(s)])
        # Backward Step 2: x ~ x | theta
        y ~ MvNormal([m; m], [sqrt(s) 0; 0 sqrt(s)])
    end
    return s, m, y
end

fw = PG(20, NSamples)
# bk = Gibbs(10, PG(10,10, :s, :y), HMC(1, 0.25, 5, :m));
bk = PG(20,50);

s = sample(gdemo2([1.5, 2], false), fw);
describe(s)

N = div(NSamples, 50)

x = [s[:y][1]...]
s_bk = Array{Turing.Chain}(undef, N)

#set_verbosity(0)
for i = 1:N
    global x, bk
    s_bk[i] = sample(gdemo2(x, true), bk);
    x = [s_bk[i][:y][end]...];
end
#set_verbosity(1)

s2 = vcat(s_bk...);
describe(s2)

# qqplot(s[:m], s2[:m])
# qqplot(s[:s], s2[:s])

using UnicodePlots

qqs = qqbuild(s[:s], s2[:s])
show(scatterplot(qqs.qx, qqs.qy, title = "QQ plot for s", canvas = DotCanvas))
show(scatterplot(qqs.qx[51:end-50], qqs.qy[51:end-50], title = "QQ plot for s (removing last 50 quantiles)", canvas = DotCanvas))

qqm = qqbuild(s[:m], s2[:m])
show(scatterplot(qqm.qx, qqm.qy, title = "QQ plot for m", canvas = DotCanvas))
show(scatterplot(qqm.qx[51:end-50], qqm.qy[51:end-50], title = "QQ plot for m (removing first and last 50 quantiles):", canvas = DotCanvas))


text:
"""
Build a q network given an architecture of the form
[CONV]-[FC]
with ReLU activation apart from the output layer
"""
function build_q(input::Tensor,
                arch::QNetworkArchitecture,
                env::AbstractEnvironment;
                scope::String="",
                reuse::Bool=false,
                dueling::Bool=false)
    return cnn_to_mlp(input, arch.conv, arch.fc, n_actions(env), scope=scope, reuse=reuse, dueling=dueling)
end

"""
Build a recurrent q network given an architecture of the form
[CONV]-[FC]-[LSTM]-[FC]
with ReLU activation apart from the output layer
"""
function build_recurrent_q_network(inputs::Tensor,
                                   convs::Vector{Tuple{Int64, Vector{Int64}, Int64}},
                                   hiddens_in::Vector{Int64},
                                   hiddens_out::Vector{Int64},
                                   lstm_size::Int64,
                                   num_output::Int64;
                                   final_activation = identity,
                                   scope::String= "drqn",
                                   reuse::Bool = false,
                                   dueling::Bool = false)

    # retrieve static dims
    trace_length = get(get_shape(inputs).dims[2])
    obs_dim = obs_dim = [get(d) for d in get_shape(inputs).dims[3:end]]
    # flatten time dim
    out = reshape(inputs, (-1, obs_dim...)) # should not do anything but get the tensor shape not unknown
    # feed into conv_to_mlp
    out = cnn_to_mlp(out, convs, hiddens_in, 0, scope=scope, reuse=reuse, dueling=false, final_activation=nn.relu)
    #retrieve time dimension
    flat_dim = get(get_shape(out).dims[end])
    out = reshape(out, (-1, trace_length, flat_dim))

    # build RNN
    rnn_cell = nn.rnn_cell.LSTMCell(lstm_size)
    c = placeholder(Float32, shape=[-1, lstm_size])
    h =  placeholder(Float32, shape=[-1, lstm_size])
    state_in = LSTMStateTuple(c, h)
    out, state_out = tf.nn.rnn(rnn_cell,
                               out,
                               initial_state=state_in,
                               reuse=reuse,
                               scope=scope)
    out = stack(out, axis=2)
    # output with dueling
    out = reshape(out, (-1, lstm_size))
    out = cnn_to_mlp(out, [], hiddens_out, num_output, scope=scope, reuse=reuse, dueling=dueling, scope=scope)
    #END OF Q NETWORK GRAPH
    return out, state_in, state_out
end

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
# # Adding additional fields
#
# All components currently implemented in `PowerSystems.jl` have an `ext` field that contains
# an empty `Dictionary`. This empty dictionary can be accessed with the method `get_ext`.
# shows an example of adding additional data to the [`ThermalStandard`](@ref) structs for other
# modeling requirements.

# This method is useful if there is no requirement to create new behaviours for components of
# a particular type but additional fields are required. A simple example is the addition of geographic information if needed.

using PowerSystems
const PSY = PowerSystems

DATA_DIR = "../../../data" #hide
system = System(joinpath(DATA_DIR, "matpower/case5_re.m"))

for g in get_components(ThermalStandard, system)
    external_field = get_ext(g)
    external_field["my_data"] = 1.0
end

gen_alta = get_component(ThermalStandard, system, "Alta")

my_data = get_ext(gen_alta)["my_data"]


text:
using jInv.Mesh
using Test

# setup 3D problem
nc = [5, 7, 8]
x0 = [0.1,0.2, 0.4]
domain = [x0[1], 4, x0[2], 2, x0[3], 6]
h   = (domain[2:2:end]-domain[1:2:end])./nc
h1  = h[1]*ones(nc[1])
h2  = h[2]*ones(nc[2])
h3  = h[3]*ones(nc[3])

Mt = getTensorMesh3D(h1,h2,h3,x0)
Mr = getRegularMesh(domain,nc)

print("\ttest cell-centered axes...")
xt = getCellCenteredAxes(Mt)
xr = getCellCenteredAxes(Mr)
for d=1:3
	@test norm(xt[d]-xr[d],Inf)<1e-14
end
print("passed\n")

print("\ttest nodal axes...")
xt = getNodalAxes(Mt)
xr = getNodalAxes(Mr)
for d=1:3
	@test norm(xt[d]-xr[d],Inf)<1e-14
end
print("passed\n")

print("\ttest cell-centered grid...")
xt = getCellCenteredGrid(Mt)
xr = getCellCenteredGrid(Mr)
@test norm(xt-xr,Inf)<1e-14
print("passed\n")

print("\ttest nodal grid...")
xt = getNodalGrid(Mt)
xr = getNodalGrid(Mr)
@test norm(xt-xr,Inf)<1e-14
print("passed\n")

print("\ttest face grid...")
xt = getFaceGrids(Mt)
xr = getFaceGrids(Mr)
for d=1:3
	@test norm(xt[d]-xr[d],Inf)<1e-14
end
print("passed\n")

print("\ttest edge grid...")
xt = getEdgeGrids(Mt)
xr = getEdgeGrids(Mr)
for d=1:3
	@test norm(xt[d]-xr[d],Inf)<1e-14
end
print("passed\n")

print("\ttest volume...")
xt = getVolume(Mt)
xr = getVolume(Mr)
@test norm(xt-xr,Inf)<1e-14
xt = getVolumeInv(Mt)
xr = getVolumeInv(Mr)
@test norm(xt-xr,Inf)<1e-14
print("passed\n")

print("\ttest area...")
xt = getFaceArea(Mt)
xr = getFaceArea(Mr)
@test norm(xt-xr,Inf)<1e-14
xt = getFaceAreaInv(Mt)
xr = getFaceAreaInv(Mr)
@test norm(xt-xr,Inf)<1e-14
print("passed\n")

print("\ttest length...")
# test length
xt = getLength(Mt)
xr = getLength(Mr)
@test norm(xt-xr,Inf)<1e-14
# test length inverse
xt = getLengthInv(Mt)
xr = getLengthInv(Mr)
@test norm(xt-xr,Inf)<1e-14
print("passed\n")

print("\ttest counting...")
@test Mt.nc==Mr.nc
@test Mt.ne==Mr.ne
@test Mt.nf==Mr.nf
@test Mt.n ==Mr.n
print("passed\n")

print("\ttest nodal gradient matrix...")
Gr  = getNodalGradientMatrix(Mr)
Gt  = getNodalGradientMatrix(Mt)
@test norm(Gr-Gt,1)/norm(Gr,1) < 1e-12
print("passed\n")

print("\ttest divergence matrix...")
Dr  = getDivergenceMatrix(Mr)
Dt  = getDivergenceMatrix(Mt)
@test norm(Dr-Dt,1)/norm(Dr,1) < 1e-12
print("passed\n")

print("\ttest curl matrix...")
Cr  = getCurlMatrix(Mr)
Ct  = getCurlMatrix(Mt)
@test norm(Cr-Ct,1)/norm(Cr,1) < 1e-12
print("passed\n")

print("\ttest face average matrix...")
Afr = getFaceAverageMatrix(Mr)
Aft = getFaceAverageMatrix(Mt)
@test norm(Afr-Aft,1)/norm(Afr,1) < 1e-12
print("passed\n")

print("\ttest edge average matrix...")
Aer = getEdgeAverageMatrix(Mr)
Aet = getEdgeAverageMatrix(Mt)
@test norm(Aer-Aet,1)/norm(Aer,1) < 1e-12
print("passed\n")


text:
### Generic rand methods

"""
    rand([rng::AbstractRNG,] s::Sampleable)

Generate one sample for `s`.

    rand([rng::AbstractRNG,] s::Sampleable, n::Int)

Generate `n` samples from `s`. The form of the returned object depends on the variate form of `s`:

- When `s` is univariate, it returns a vector of length `n`.
- When `s` is multivariate, it returns a matrix with `n` columns.
- When `s` is matrix-variate, it returns an array, where each element is a sample matrix.

    rand([rng::AbstractRNG,] s::Sampleable, dim1::Int, dim2::Int...)
    rand([rng::AbstractRNG,] s::Sampleable, dims::Dims)

Generate an array of samples from `s` whose shape is determined by the given
dimensions.
"""
rand(s::Sampleable) = rand(GLOBAL_RNG, s)
rand(s::Sampleable, dims::Dims) = rand(GLOBAL_RNG, s, dims)
rand(s::Sampleable, dim1::Int, moredims::Int...) =
    rand(GLOBAL_RNG, s, (dim1, moredims...))
rand(rng::AbstractRNG, s::Sampleable, dim1::Int, moredims::Int...) =
    rand(rng, s, (dim1, moredims...))

"""
    rand!([rng::AbstractRNG,] s::Sampleable, A::AbstractArray)

Generate one or multiple samples from `s` to a pre-allocated array `A`. `A` should be in the
form as specified above. The rules are summarized as below:

- When `s` is univariate, `A` can be an array of arbitrary shape. Each element of `A` will
  be overriden by one sample.
- When `s` is multivariate, `A` can be a vector to store one sample, or a matrix with each
  column for a sample.
- When `s` is matrix-variate, `A` can be a matrix to store one sample, or an array of
  matrices with each element for a sample matrix.
"""
function rand! end
rand!(s::Sampleable, X::AbstractArray{<:AbstractArray}, allocate::Bool) =
    rand!(GLOBAL_RNG, s, X, allocate)
rand!(s::Sampleable, X::AbstractArray) = rand!(GLOBAL_RNG, s, X)
rand!(rng::AbstractRNG, s::Sampleable, X::AbstractArray) = _rand!(rng, s, X)

"""
    sampler(d::Distribution) -> Sampleable
    sampler(s::Sampleable) -> s

Samplers can often rely on pre-computed quantities (that are not parameters
themselves) to improve efficiency. If such a sampler exists, it can be provided
with this `sampler` method, which would be used for batch sampling.
The general fallback is `sampler(d::Distribution) = d`.
"""
sampler(s::Sampleable) = s

# Random API
Random.Sampler(::Type{<:AbstractRNG}, s::Sampleable, ::Val{1}) = s
Random.Sampler(::Type{<:AbstractRNG}, s::Sampleable, ::Val{Inf}) = sampler(s)


text:
# FIXME: Add some condition to guarantee that: TR <: Union(T, complex(T)) !

# T the numeric type of the transfer function
# TR the type of the roots
struct SisoZpk{T,TR<:Number} <: SisoTf{T}
    z::Vector{TR}
    p::Vector{TR}
    k::T
    function SisoZpk{T,TR}(z::Vector{TR}, p::Vector{TR}, k::T) where  {T<:Number, TR<:Number}
        if k == zero(T)
            p = TR[]
            z = TR[]
        end
        if TR <: Complex && T <: Real
            @assert check_real(z) "zpk model should be real-valued, but zeros do not come in conjugate pairs."
            @assert check_real(p) "zpk model should be real-valued, but poles do not come in conjugate pairs."
        end
        new{T,TR}(z, p, k)
    end
end
function SisoZpk{T,TR}(z::Vector, p::Vector, k::Number) where  {T<:Number, TR<:Number}
    SisoZpk{T,TR}(Vector{TR}(z), Vector{TR}(p), T(k))
end
function SisoZpk{T}(z::Vector, p::Vector, k::Number) where T
    TR = complex(T)
    SisoZpk{T,TR}(Vector{TR}(z), Vector{TR}(p), T(k))
end
function SisoZpk(z::AbstractVector{TZ}, p::AbstractVector{TP}, k::T) where {T<:Number, TZ<:Number, TP<:Number} # NOTE: is this constructor really needed?
    TR = promote_type(TZ,TP)
    # Could check if complex roots come with their conjugates,
    # i.e., if the SisoZpk corresponds to a real-valued system

    if TR <: Complex && T <: Real
        @assert check_real(z) "zpk model should be real-valued, but zeros do not come in conjugate pairs."
        @assert check_real(p) "zpk model should be real-valued, but poles do not come in conjugate pairs."
    end
    SisoZpk{T,TR}(Vector{TR}(z), Vector{TR}(p), k)
end

Base.zero(::Type{SisoZpk{T}}) where T = SisoZpk{T}(T[], T[], zero(T))
Base.one(::Type{SisoZpk{T}}) where T = SisoZpk{T}(T[], T[], one(T))

Base.zero(::Type{SisoZpk{T,TR}}) where {T, TR} = SisoZpk{T,TR}(TR[], TR[], zero(T))
Base.one(::Type{SisoZpk{T,TR}}) where {T, TR} = SisoZpk{T,TR}(TR[], TR[], one(T))

Base.one(f::SisoZpk) = one(typeof(f))
Base.zero(f::SisoZpk) = zero(typeof(f))

# tzero is not meaningful for transfer function element? But both zero and zeros are taken...
tzero(f::SisoZpk) = f.z # Do minreal first?,
pole(f::SisoZpk) = f.p # Do minreal first?

numpoly(f::SisoZpk{<:Real}) = f.k*prod(roots2real_poly_factors(f.z))
denpoly(f::SisoZpk{<:Real}) = prod(roots2real_poly_factors(f.p))

numpoly(f::SisoZpk) = f.k*prod(roots2poly_factors(f.z))
denpoly(f::SisoZpk) = prod(roots2poly_factors(f.p))

numvec(f::SisoZpk) = reverse(coeffs(numpoly(f))) # FIXME: reverse?!
denvec(f::SisoZpk) = reverse(coeffs(denpoly(f))) # FIXME: reverse?!

isproper(f::SisoZpk) = (length(f.z) <= length(f.p))

# Will not be able to create complex polynomials without type instability

function minreal(sys::SisoZpk{T,TR}, eps::Real) where {T, TR}
    if length(sys.p) == 0
        return sys
    end

    newZ = copy(sys.z)
    newP = Vector{TR}()

    pidx = 1
    p = sys.p[pidx]
    while true
        if isempty(newZ)
            push!(newP, p)
        else
            distance, zidx = findmin(abs.(p .- newZ))

            if distance < eps
                if imag(p) == 0 && imag(newZ[zidx]) != 0
                    newZ[zidx+1] = real(newZ[zidx+1])
                end
                if imag(newZ[zidx]) == 0 && imag(p) != 0
                    pidx += 1
                    p = real(sys.p[pidx])
                    deleteat!(newZ, zidx)
                    continue
                end
                deleteat!(newZ, zidx)
            else
                push!(newP, p)
            end
        end

        pidx += 1
        if pidx > length(sys.p)
            break
        end
        p = sys.p[pidx]
    end
    SisoZpk{T, TR}(newZ, newP, sys.k)
end

# FIXME: Perhaps move to some other file with auxilliary functions,
# the name could also be imporoved. Perhaps this functionality can be found in some other package.
""" If TR is Complex and T is Real, check that every pole is matched to its conjugate
    this assumes that the compelx poles are ordered as they are output by the LAPACK
    routines that return complex-conjugated values, i.e., (x+iy) is followed by (x-iy)"""
function check_real(r_vec::AbstractVector{<:Complex})
    k = 1
    while k <= length(r_vec)
        if isreal(r_vec[k])
            # Do nothing
        else
            if k == length(r_vec) # Element is complex and there is no more elements to match up with
                return false
            elseif conj(r_vec[k]) == r_vec[k+1] # Element is complex and succeeding element is its conjugate
                k=k+1
            else
                return false
            end
        end
        k += 1
    end
    return true
end

function evalfr(f::SisoZpk{T1,TR}, s::T2) where {T1<:Number, TR<:Number, T2<:Number}
    T0 = promote_type(T2, TR)
    T = promote_type(T1, Base.promote_op(/, T0, T0))
    den = mapreduce(a -> s-a, *, f.p, init = one(T0)) # Default to one

    if den == zero(T0)
        return convert(T, Inf)
    else
        num = mapreduce(b -> s-b, *, f.z, init=one(T0))
        return f.k*num/den
    end
end

function poly_factors2string(poly_factors::AbstractArray{<:Poly{T}}, var) where T
    if length(poly_factors) == 0
        str = sprint(printpolyfun(var), Poly(one(T)))
    elseif length(poly_factors) == 1
        str = sprint(printpolyfun(var), poly_factors[1])
    else
        str = prod(z -> "("*sprint(printpolyfun(var), z)*")", poly_factors)
    end
end

""" Heurisitc function that tries to add parentheses when printing the coeffcient
    for systems on zpk form. Should at least handle the following types
    Measurment, Dual, Sym. """
function _printcoefficient(nbr::Number)
    # Print type information as in 1.0f0 for Float32
    # showcompact might be better, but is not consistent with polynomials
    nbr_string = sprint(show,nbr)
    if occursin(" + ", nbr_string) || occursin(" - ", nbr_string) || occursin(" ± ", nbr_string)
        return "(" * nbr_string * ")" # Add parens
    else
        return nbr_string
    end
end

function print_siso(io::IO, t::SisoZpk, var=:s)
    if typeof(t.k) <: Real
        numstr = poly_factors2string(roots2real_poly_factors(t.z), var)
        denstr = poly_factors2string(roots2real_poly_factors(t.p), var)
    else
        numstr = poly_factors2string(roots2poly_factors(t.z), var)
        denstr = poly_factors2string(roots2poly_factors(t.p), var)
    end
    # Figure out the length of the separating line
    len_num = length(numstr)
    len_den = length(denstr)
    dashcount = max(len_num, len_den)

    # Center the numerator or denominator
    if len_num < dashcount
        numstr = "$(repeat(" ", div(dashcount - len_num, 2)))$numstr"
    else
        denstr = "$(repeat(" ", div(dashcount - len_den, 2)))$denstr"
    end

    gainstr = _printcoefficient(t.k)
    #Add spaces to account for gain string
    numstr = " "^(length(gainstr))*numstr
    denstr = " "^(length(gainstr))*denstr
    println(io, numstr)
    println(io, gainstr*repeat("-", dashcount))
    println(io, denstr)
end

==(f1::SisoZpk, f2::SisoZpk) = (f1-f2).k == 0
function isapprox(f1::SisoZpk, f2::SisoZpk; rtol::Real=sqrt(eps()), atol::Real=sqrt(eps()))
    fdiff = f1 - f2
    isapprox(fdiff.k, 0, atol=atol, rtol=rtol)
end

function +(f1::SisoZpk{T1,TR1}, f2::SisoZpk{T2,TR2}) where {T1<:Number,T2<:Number,TR1<:Number,TR2<:Number}
    numPoly = numpoly(f1)*denpoly(f2) + numpoly(f2)*denpoly(f1)

    TRtmp = promote_type(TR1, TR2)
    # Calculating roots can make integers floats
    TRnew = Base.promote_op(/, TRtmp, TRtmp)
    z = convert(Vector{TRnew}, roots(numPoly))
    #TODO gains could remain integers, but numerical precision inhibits this
    Ttmp = promote_type(T1,T2)
    Tnew = Base.promote_op(/, Ttmp, Ttmp)
    if length(numPoly) > 0
        k = convert(Tnew, numPoly[end])
        p = convert(Vector{TRnew}, [f1.p;f2.p])
    else
        k = zero(Tnew)
        p = TRnew[]
    end

    # FIXME:
    # Threshold for pole-zero cancellation should depend on the roots of the system
    # Note the difference between continuous and discrete-time systems...
    minreal(SisoZpk(z,p,k), sqrt(eps()))
end

+(f::SisoZpk, n::T) where {T<:Number} = f + SisoZpk{T,T}(T[],T[],n)
+(n::Number, f::SisoZpk) = f + n
#.+(t::SisoZpk, n::Number) = t + n
#.+(n::Number, t::SisoZpk) = t + n

-(t1::SisoZpk, t2::SisoZpk) = +(t1,-t2)
-(n::T, f::SisoZpk) where {T<:Number} = SisoZpk{T,T}(T[],T[],n) - f
-(t::SisoZpk, n::Number) = +(t, -n)
#.-(t::SisoZpk, n::Number) = -(t, n)
#.-(n::Number, t::SisoZpk) = -(n, t)

-(f::SisoZpk) = SisoZpk(f.z, f.p, -f.k)

*(f1::SisoZpk, f2::SisoZpk) = SisoZpk([f1.z;f2.z], [f1.p;f2.p], f1.k*f2.k)

*(f::SisoZpk, n::Number) = SisoZpk(f.z, f.p, f.k*n)
*(n::Number, f::SisoZpk) = *(f, n)
#.*(f1::SisoZpk, f2::SisoZpk) = *(f1, f2)
#.*(t::SisoZpk, n::Number) = *(t, n)
#.*(n::Number, t::SisoZpk) = *(t, n)

/(n::Number, f::SisoZpk) = SisoZpk(f.p, f.z, n/f.k)
/(f::SisoZpk, n::Number) = SisoZpk(f.z, f.p, f.k/n)
/(f1::SisoZpk, f2::SisoZpk) = f1*(1/f2)


text:
using Turing, Random, Test

dir = splitdir(splitdir(pathof(Turing))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")

@stage_testset "stan" "stan-interface.jl" begin
    using CmdStan

    chn = sample(gdemo_default, 2000, 1000, false, 1,
        CmdStan.Adapt(), CmdStan.Hmc(stepsize_jitter=0))
    check_gdemo(chn)
end


text:
export SingleMissingPredictive, PosteriorPredictive

struct SingleMissingPredictive
    S::InteractionSequence{Int}
    ind::Tuple{Int,Int}
    p::Vector{Float64}
end

function Base.show(io::IO, pred::SingleMissingPredictive)
    title = "Missing Entry Predictive Distribution"
    println(io, title)
    println(io, "-"^length(title))
    println(io, "Observation: $(pred.S)")
    println(io, "Missing entry: $(pred.ind)")
end

struct PosteriorPredictive
    posterior::Union{SisPosteriorMcmcOutput, SimPosteriorMcmcOutput}
    model_type::DataType
    function PosteriorPredictive(p::Union{SisPosteriorMcmcOutput, SimPosteriorMcmcOutput})
        T = typeof(p.posterior.S_prior)
        new(p, T)
    end
end


text:
# # Even reduction

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Symmetry/even_reduction.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Symmetry/even_reduction.ipynb)

using Test #src
using DynamicPolynomials
@polyvar x

# We would like to find the minimum value of the following polynomial:

poly = x^4 - 2x^2

using SumOfSquares

# We define the custom action as follows:

struct OnSign <: Symmetry.OnMonomials end
using PermutationGroups
import SymbolicWedderburn
SymbolicWedderburn.coeff_type(::OnSign) = Float64
function SymbolicWedderburn.action(::OnSign, p::Permutation, mono::AbstractMonomial)
    if isone(p) || iseven(DynamicPolynomials.degree(mono))
        return 1 * mono
    else
        @assert p.perm == perm"(1,2)"
        return -1 * mono
    end
end
G = PermGroup([perm"(1,2)"])

# We can exploit the symmetry as follows:

import CSDP
solver = CSDP.Optimizer
model = Model(solver)
@variable(model, t)
@objective(model, Max, t)
pattern = Symmetry.Pattern(G, OnSign())
con_ref = @constraint(model, poly - t in SOSCone(), symmetry = pattern)
optimize!(model)
@test value(t) ≈ -1 #src
value(t)

# We indeed find `-1`, let's verify that symmetry was exploited:

@test length(gram_matrix(con_ref).sub_gram_matrices) == 2 #src
@test gram_matrix(con_ref).sub_gram_matrices[1].basis.polynomials == [x^2, 1] #src
@test gram_matrix(con_ref).sub_gram_matrices[2].basis.polynomials == [x] #src
gram_matrix(con_ref)

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
### A Pluto.jl notebook ###
# v0.18.4

using Markdown
using InteractiveUtils

# ╔═╡ 4fc964c9-4690-412f-86ca-0ad975e5dc44
struct Squares
    count::Int
end

# ╔═╡ 2a1988c7-e51f-4a8f-b8c0-f71a744b3c68
begin
	Base.iterate(S::Squares, state=1) =
    	state > S.count ? nothing : (state*state, state+1)
			Base.eltype(::Type{Squares}) = Int # Note that this is defined for the type
	Base.length(S::Squares) = S.count
	Base.sum(S::Squares) = (n = S.count; return n*(n+1)*(2n+1)÷6)
	Base.iterate(rS::Iterators.Reverse{Squares}, state=rS.itr.count) = state < 1 ? 		nothing : (state*state, state-1)

	function Base.getindex(S::Squares, i::Int)
		1 <= i <= S.count || throw(BoundsError(S, i))
		return i*i
	end

	Base.firstindex(S::Squares) = 1
	Base.lastindex(S::Squares) = length(S)

	Base.getindex(S::Squares, i::Number) = S[convert(Int, i)]
	Base.getindex(S::Squares, I) = [S[i] for i in I]
end

# ╔═╡ 99e4bdcd-696c-4a87-bd2b-9445b3f32bff
begin
	items = Int[]
	for item in Squares(9) # or "for item = Squares(9)"
	    # body
	    append!(items, item)
	end
	items
end

# ╔═╡ ba0a922b-37fd-4608-a855-0b2fb52e1e3b
25 in Squares(10)

# ╔═╡ fe592470-cca8-4f18-a7c1-ae72ea478067
collect(Squares(8))

# ╔═╡ 13631280-0ee6-466d-b46e-af49d8478702
sum(Squares(3))

# ╔═╡ 542e5321-e4f4-4991-b888-c2b373e07be1
collect(Iterators.reverse(Squares(4)))

# ╔═╡ e0d63466-3c0d-43f1-bc72-09c15e361edd
Squares(100)[23]

# ╔═╡ 41f90116-82f9-4a64-a300-8e9f119cb6f8
Squares(23)[end]

# ╔═╡ 0a09c030-25f0-47ca-ad3b-ad745c1e3052
Squares(10)[[3, 4., 5]]

# ╔═╡ ba97acaf-8b4c-4743-9277-74de9db40ac4
begin
	struct SquaresVector <: AbstractArray{Int, 1}
		count::Int
	end
	Base.size(S::SquaresVector) = (S.count,)
	Base.IndexStyle(::Type{<:SquaresVector}) = IndexLinear()
	Base.getindex(S::SquaresVector, i::Int) = i*i
end

# ╔═╡ fd5337b3-9ef4-4ace-b69f-cfbf036c58d4
s = SquaresVector(4)

# ╔═╡ 320f9fb8-be12-4287-aed7-1847329cc679
s[s .> 8]

# ╔═╡ eb629fdf-7810-4ce4-9133-6825f669ac2b
s + s

# ╔═╡ e50a56db-3eaa-4323-8ae6-b36b716c96e6
sin.(s)

# ╔═╡ a702d978-462c-4879-859d-5a92ba4b6c83
IndexStyle(SquaresVector)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.9.0-DEV"
manifest_format = "2.0"
project_hash = "da39a3ee5e6b4b0d3255bfef95601890afd80709"

[deps]
"""

# ╔═╡ Cell order:
# ╠═4fc964c9-4690-412f-86ca-0ad975e5dc44
# ╠═2a1988c7-e51f-4a8f-b8c0-f71a744b3c68
# ╠═99e4bdcd-696c-4a87-bd2b-9445b3f32bff
# ╠═ba0a922b-37fd-4608-a855-0b2fb52e1e3b
# ╠═fe592470-cca8-4f18-a7c1-ae72ea478067
# ╠═13631280-0ee6-466d-b46e-af49d8478702
# ╠═542e5321-e4f4-4991-b888-c2b373e07be1
# ╠═e0d63466-3c0d-43f1-bc72-09c15e361edd
# ╠═41f90116-82f9-4a64-a300-8e9f119cb6f8
# ╠═0a09c030-25f0-47ca-ad3b-ad745c1e3052
# ╠═ba97acaf-8b4c-4743-9277-74de9db40ac4
# ╠═fd5337b3-9ef4-4ace-b69f-cfbf036c58d4
# ╠═320f9fb8-be12-4287-aed7-1847329cc679
# ╠═eb629fdf-7810-4ce4-9133-6825f669ac2b
# ╠═e50a56db-3eaa-4323-8ae6-b36b716c96e6
# ╠═a702d978-462c-4879-859d-5a92ba4b6c83
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using MakieLayout
using Makie

begin
    scene = Scene(resolution = (1000, 1000));
    screen = display(scene)
    campixel!(scene);

    maingl = GridLayout(scene, alignmode = Outside(30))

    ni = 2
    nj = 2

    las = maingl[1:ni, 1:nj] = [LAxis(scene, titlesize=20) for i in 1:ni, j in 1:nj]

    lines!(las[1, 1], 0:0.1:100, sin.((0:0.1:100) ./ 5) .+ 1)
    lines!(las[1, 2], (0:0.1:100) .+ 50, -sin.((0:0.1:100) ./ 5), color=:blue)

    scatter!(las[2, 1], rand(100, 2) .* 200, markersize=4)

    las[2, 2].attributes.xautolimitmargin = (0, 0)
    las[2, 2].attributes.yautolimitmargin = (0, 0)
    las[2, 2].attributes.aspect = DataAspect()
    image!(las[2, 2], rand(500, 400))

    linkxaxes!(las[1, :]...)
    linkyaxes!(las[1, :]...)

    autolimits!(las[1, 1])
end


text:
# Add fallbacks for predict_* which require mean, mode, median.

const BadModeTypes   = Union{AbstractArray{Continuous},Table(Continuous)}
const BadMeanTypes   = Union{AbstractArray{<:Finite},Table(Finite)}
const BadMedianTypes = Union{AbstractArray{<:Finite},Table(Finite)}

# mode:
predict_mode(m::Probabilistic, fitres, Xnew) =
    predict_mode(m, fitres, Xnew, Val(target_scitype(m)))
predict_mode(m, fitresult, Xnew, ::Any) =
    mode.(predict(m, fitresult, Xnew))
predict_mode(m, fitresult, Xnew, ::Val{<:BadModeTypes}) =
    throw(ArgumentError("Attempting to compute mode of predictions made "*
                        "by a model expecting `Continuous` targets. "))

# mean:
predict_mean(m::Probabilistic, fitresult, Xnew) =
    predict_mean(m, fitresult, Xnew, Val(target_scitype(m)))
predict_mean(m, fitresult, Xnew, ::Any) =
    mean.(predict(m, fitresult, Xnew))
predict_mean(m, fitresult, Xnew, ::Val{<:BadMeanTypes}) =
    throw(ArgumentError("Attempting to compute mean of predictions made "*
                        "by a model expecting `Finite` targets. "))

# median:
predict_median(m::Probabilistic, fitresult, Xnew) =
    predict_median(m, fitresult, Xnew, Val(target_scitype(m)))
predict_median(m, fitresult, Xnew, ::Any) =
    median.(predict(m, fitresult, Xnew))
predict_median(m, fitresult, Xnew, ::Val{<:BadMedianTypes}) =
    throw(ArgumentError("Attempting to compute median of predictions made "*
                        "by a model expecting `Finite` targets. "))

# # operations implemented by some meta-models:
# function evaluate end

# not in MLJModelInterface as methodswith requires InteractiveUtils
MMI.implemented_methods(::FI, M::Type{<:MLJType}) =
    getfield.(methodswith(M), :name)


text:
using Images
using TestImages
using StaticArrays
using OffsetArrays
using Random
using CoordinateTransformations

function evaluate_error(dims, flow::Array{SVector{2, Float64}, 1}, Δ, tol)
    error_count = 0
    maximum_error = 0.0
    for i in eachindex(flow)
        δ = flow[i]
        if abs(δ[1] - Δ[1]) > tol || abs(δ[2] - Δ[2]) > tol
            error_count += 1
        end
        error = sqrt( (δ[1] - Δ[1])^2 + (δ[2] - Δ[2])^2 )
        if error > maximum_error
            maximum_error = error
        end
    end
    return error_count, maximum_error
end

function get_keypoints(image, n_points)
    corners = imcorner(image, method=shi_tomasi)
    I = CartesianIndices(corners)[corners .!= 0]
    r, c = (getindex.(I, 1), getindex.(I, 2))
    points = map((ri, ci) -> SVector{2, Float64}(ri, ci), r, c)
    rand(points, n_points)
end

@testset "Lucas-Kanade" begin
    @testset "Horizontal Motion" begin
        @info "Running Horizontal Motion test."
        Random.seed!(9876)

        tol = 1e-1
        maximum_percentage_error = 7
        number_test_pts = 500
        img1 = Gray{Float64}.(testimage("mandrill"))
        algorithm = LucasKanade()

        Δ = (0.0, 3.0)
        trans = Translation(-Δ[1], -Δ[2])
        img2 = warp(img1, trans, axes(img1))

        points = get_keypoints(img1, number_test_pts)

        flow, status = optical_flow(img1, img2, points, algorithm)
        error_count, maximum_error = evaluate_error(size(img1), flow[status], Δ, tol)
        percentage_error = (error_count / sum(status)) * 100
        @info "Case: Horizontal Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error < maximum_percentage_error

        # Same as above, except that we pass in an initial displacement field.
        displacement = fill(SVector{2, Float64}(0.0, 0.0), length(points))
        flow, status = optical_flow(img1, img2, points, displacement, algorithm)
        error_count, maximum_error = evaluate_error(size(img1), flow[status], Δ, tol)
        percentage_error = (error_count / sum(status)) * 100
        @info "Case: Horizontal Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error < maximum_percentage_error
    end

    @testset "Vertical Motion" begin
        @info "Running Vertical Motion test."
        Random.seed!(9876)

        tol = 1e-1
        maximum_percentage_error = 7
        number_test_pts = 500
        img1 = Gray{Float64}.(testimage("mandrill"))
        algorithm = LucasKanade()

        Δ = (3.0, 0.0)
        trans = Translation(-Δ[1], -Δ[2])
        img2 = warp(img1, trans, axes(img1))

        points = get_keypoints(img1, number_test_pts)

        flow, status_array = optical_flow(img1, img2, points, algorithm)
        error_count, maximum_error = evaluate_error(size(img1), flow[status_array], Δ, tol)
        percentage_error = (error_count / sum(status_array)) * 100
        @info "Case: Vertical Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error < maximum_percentage_error

        # Same as above, except that we pass in an initial displacement field.
        displacement = fill(SVector{2, Float64}(0.0, 0.0), length(points))
        flow, status_array = optical_flow(img1, img2, points, displacement, algorithm)
        error_count, maximum_error = evaluate_error(size(img1), flow[status_array], Δ, tol)
        percentage_error = (error_count / sum(status_array)) * 100
        @info "Case: Vertical Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error < maximum_percentage_error
    end

    @testset "Combined Motion" begin
        @info "Running Combined Motion test."
        Random.seed!(9876)

        tol = 1e-1
        maximum_percentage_error = 7
        number_test_pts = 500
        img1 = Gray{Float64}.(testimage("mandrill"))
        algorithm = LucasKanade()

        Δ = (3.0, -1.0)
        trans = Translation(-Δ[1], -Δ[2])
        img2 = warp(img1, trans, axes(img1))

        points = get_keypoints(img1, number_test_pts)

        flow, status_array = optical_flow(img1, img2, points, algorithm)
        error_count, maximum_error = evaluate_error(size(img1), flow[status_array], Δ, tol)
        percentage_error = (error_count / sum(status_array)) * 100
        @info "Case: Combined Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error < maximum_percentage_error

        # Same as above, except that we pass in an initial displacement field.
        displacement = fill(SVector{2, Float64}(0.0, 0.0), length(points))
        flow, status_array = optical_flow(img1, img2, points, displacement, algorithm)
        error_count, maximum_error = evaluate_error(size(img1), flow[status_array], Δ, tol)
        percentage_error = (error_count / sum(status_array)) * 100
        @info "Case: Combined Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error < maximum_percentage_error
    end
end