text:
"""
A Dictionary to store contract attributes.
"""
struct ContractAttributes
    attributes::Dict{String, Any}
end

#
# Observables
#

"""
Represents a time-varying quantity of type {T}.
So, its value is indexed by a date or datetime.
"""
abstract type Observable{T} end

#=
"""
The value of `CurrentDate()` at date `dt` is just `dt`.
"""
struct CurrentDate <: Observable{Date} end
=#

"""
    PricingDate <: Observable{Date}

Represents the pricing date
"""
struct PricingDate <: Observable{Date} end

"""
    AdjustedDate{D<:Observable{Date}} <: Observable{Date}

Adjusts to last business day if `date` is not a bday
following `calendar` holiday calendar.

# Fields

* `date`

* `calendar`
"""
struct AdjustedDate{D<:Observable{Date}} <: Observable{Date}
    reference_date::D
    calendar::BusinessDays.HolidayCalendar
end

"""
The value of `Konst(x)` is always `x` at any date.
"""
struct Konst{T} <: Observable{T}
    val::T
end

struct ObservableAt{D<:Observable{Date},T} <: Observable{T}
    at::D
    o::Observable{T}
end

"""
    at(date, observable)

Fixes observation date of `observable`
at `date`.

`date` might be a `Date` or `Observable{Date}`.
"""
at(date::Date, k::Konst) = k
at(date::Date, o::Observable) = ObservableAt(Konst(date), o)

function at(d::D, o::Observable) where {D<:Observable{Date}}
    return ObservableAt(d, o)
end

function at(date::Date, oat::ObservableAt)
    @assert date == oat.at
    return oat
end

"""
MarketData historical series value observed at
date infered by `at` observable,
asof `PricingDate`.
"""
struct HistoricalValue{D<:Observable{Date}} <: Observable{Float64}
    serie_name::String
    locf::Bool
    at::D
end

"""
First order risk factor.
"""
abstract type RiskFactor <: Observable{Float64} end

"""
Second order risk factor.

See [`DeltaGammaApproxExposuresMethod`](@ref).
"""
struct SecondOrderRiskFactor{T<:RiskFactor} <: RiskFactor
    rf::T
end

"""
Identifies a Risk Factor that has no currency value.
"""
abstract type NonCashRiskFactor <: RiskFactor end

"""
Identifies a Risk Factor that has a currency value.
"""
abstract type CashRiskFactor <: RiskFactor end

"""
Spot price for a currency.
"""
struct SpotCurrency{C<:Currencies.Currency} <: CashRiskFactor
    currency::C
end

"""
Discount factor from an interest rate curve.
"""
struct DiscountFactor <: NonCashRiskFactor
    sym::Symbol # curve symbol (:cpUSD, :PRE, etc.)
    maturity::Date
end

"""
Discount Factor from `start_date` maturing at `end_date`.

Given dates *a*, *b*, *c*, where ``a < b < c``,
we have the following relationship between effective factors:

`` f^{a \\to b} f_{FWD}^{ b \\to c } = f^{ a \\to c} ``

Rewriting as discount factors, we have:

`` DF_{FWD}^{ b \\to c } = \\frac{ DF^{ a \\to c} }{ DF^{ a \\to b } } ``
"""
struct DiscountFactorForward <: Observable{Float64}
    sym::Symbol
    start_date::Date
    end_date::Date
end

"Stock price."
struct Stock <: CashRiskFactor
    ticker::Symbol
end

"Risk Factor volatility."
struct Volatility{T<:RiskFactor} <: NonCashRiskFactor
    underlying::T
end

Currencies.currency_symbol(::SpotCurrency{Currencies.Currency{C}}) where {C} = C # see Currencies.currency_symbol(::Currencies.Currency{C})

"""
    risk_factor_symbol(rf::RiskFactor) :: Symbol

Returns a unique symbol used to identify the Risk Factor.

This is used to map the arguments of the pricing function
to the Risk Factors.

In the case of Discount Factors, the risk factor symbol
also carries the maturity date.
For instance, `:cpUSD_20200226` represents `DiscountFactor(:cpUSD, Date(2020, 2, 26))`.

See also [`market_data_symbol`](@ref).
"""
risk_factor_symbol(::SpotCurrency{Currencies.Currency{C}}) where {C} = C
risk_factor_symbol(rf::Stock) = rf.ticker

# :cpUSD_20200226
risk_factor_symbol(rf::DiscountFactor) = Symbol(rf.sym, :_, Dates.format(rf.maturity, dateformat"yyyymmdd"))

"""
    market_data_symbol(rf::RiskFactor) :: Symbol

Returns the name for the historical series for a Risk Factor.
This is used to search series in the `MarketData` module.

Discount Factors are identified by the curve name.
Therefore, this function differs from `risk_factor_symbol`
for Discount Factors, since it does not carry information
about the maturity of the Discount Factor.

See also [`risk_factor_symbol`](@ref).
"""
market_data_symbol(rf::RiskFactor) = risk_factor_symbol(rf)
market_data_symbol(rf::DiscountFactor) = rf.sym

"""
Models a pricing model input that is given as a fixed constant,
but that we would like to consider as a Risk Factor.
"""
struct FixedCashRiskFactor{R<:CashRiskFactor} <: CashRiskFactor
    rf::R
    val::Currencies.Cash
end

struct FixedNonCashRiskFactor{R<:NonCashRiskFactor} <: NonCashRiskFactor
    rf::R
    val::Float64
end

for T in (:FixedCashRiskFactor, :FixedNonCashRiskFactor)
    @eval begin
        risk_factor_symbol(rf::$T) = risk_factor_symbol(rf.rf)
        Base.isless(x, rf::$T) = isless(x, rf.val)
        Base.isless(rf::$T, x) = isless(rf.val, x)
    end
end

#=
"""
The value of `At(dt)` is `true` when `dt == value(CurrentDate())`.
"""
struct At <: Observable{Bool}
    date::Date
end
=#

@inline obstype(::Type{P}) where {T,P<:Observable{T}} = T
@inline obstype(obs::Observable) = obstype(typeof(obs))

struct LiftObs{F<:Function, R} <: Observable{R}
    f::F
    o::Observable
end

function LiftObs(f::Function, a::Observable)
    RR = Base.return_types(f, Tuple{obstype(a)})
    R = length(RR) > 1 ? Any : RR[1]
    return LiftObs{typeof(f), R}(f, a)
end

struct LiftObs2{F<:Function, R} <: Observable{R}
    f::F
    o1::Observable
    o2::Observable
end

function LiftObs2(f::Function, o1::Observable, o2::Observable)
    RR = Base.return_types(f, Tuple{obstype(o1), obstype(o2)})
    R = length(RR) > 1 ? Any : RR[1]
    return LiftObs2{typeof(f), R}(f, o1, o2)
end

struct ReduceObs{F<:Function, T} <: Observable{T}
    f::F
    observables::Vector{Observable{T}}
end

function ReduceObs(f::Function, ::Type{T}) where {T}
    return ReduceObs{typeof(f), T}(f, Vector{Observable{T}}())
end

function push_reduce_obs!(red::ReduceObs{F, T}, obs::Observable{T}) where {F, T}
    push!(red.observables, obs)
    nothing
end

# redefining push_obs! for Konst to avoid stackoverflow at convert Konst method
function push_reduce_obs!(red::ReduceObs{F, T}, obs::Konst{T}) where {F, T}
    push!(red.observables, obs)
    nothing
end

function push_reduce_obs!(red::ReduceObs{F, T}, obs::Konst{R}) where {F, T, R}
    push_reduce_obs!(red, Konst(convert(T, obs.val)))
end

"""
This struct is used to allow the reuse of
the result of `Compiler.lower!` for observables
that get repeated in the contract definition.

See also [`buffer`](@ref).
"""
mutable struct BufferedObservable{T} <: Observable{T}
    o::Observable{T}
end

"""
    buffer(observable)

Marks an `observable` as buffered.
For constants, this is a *no-op*.

See [`BufferedObservable`](@ref).
"""
@inline buffer(o::Observable) = BufferedObservable(o)
@inline buffer(k::Konst) = k # no-op for constants

"""
A Financial Contract.
"""
abstract type Contract end

#
# Contract primitives
#

"""
`Worthless` is a contract that has no rights and no obligations.
"""
struct Worthless <: Contract
end

"""
    Unit(o::SpotCurrency; carryless::Bool=false)

If you acquire `Unit(k)` you immediately receive one unit of the currency `k`.

`carryless::Bool`: if true, the unit does not use carry curves when it is being valued (discounted).
"""
struct Unit{RF<:CashRiskFactor} <: Contract
    o::RF
    carryless::Bool # controls wether we should use the carry curve when discounting

    function Unit(o::RF; carryless::Bool=false) where {RF<:CashRiskFactor}
        return new{RF}(o, carryless)
    end
end

"""
`Give(c)` inverts the position (right or obligation) of a contract `c`.
"""
struct Give <: Contract
    c::Contract
end

"""
Acquiring `Scale(o, c)` is equivalent to acquiring `c` with all
payments multiplied by `o` sampled at the time of acquisition.
"""
struct Scale <: Contract
    s::Observable
    c::Contract
end

# avoids Scale(s, scale) chain
Scale(s::Observable, c::Scale) = Scale(s * c.s, c.c)

"""
Acquiring `Both(c1, c2)` implies in acquiring immediately both contracts `c1` and `c2`.
"""
struct Both <: Contract
    c1::Contract
    c2::Contract
end

"""
`Either(c1, c2)` will choose one contract to aquire: either `c1` or `c2`.
"""
struct Either{C1<:Contract,C2<:Contract} <: Contract
    c1::C1
    c2::C2
end

#=
"""
`When(o, c)` acquires `c` as soon as `o` becomes `true`.
It's worthless after
"""
struct When{P<:Observable{Bool}, C<:Contract} <: Contract
    p::P
    c::C
end
=#

"""
`WhenAt(mat::Date, c::Contract; expires_at_maturity::Bool=false)` will acquire `c` at date `mat`.

Optional kwarg `expires_at_maturity` controls the behavior of `get_horizon(c)` for this contract:

* if `expires_at_maturity = false` (default), `get_horizon(c)` will return the maximum horizon found
in `c.maturity` and the result of `get_horizon(c.c)`.

* if `expires_at_maturity = true`, `get_horizon(c)` will always return `c.maturity`.

This is useful for Forward contracts, where the underlying is another fixed income contract,
but the settlement of the Forward contract occurs using cash at `maturity`.
"""
struct WhenAt{C<:Contract} <: Contract
    expires_at_maturity::Bool
    maturity::Date
    c::C

    function WhenAt(maturity::Date, c::C; expires_at_maturity::Bool=false) where {C<:Contract}
        new{C}(expires_at_maturity, maturity, c)
    end
end

"""
`Anytime(expiry::Date, c::Contract)` is an option to acquire `c` anytime until `expiry`.
"""
struct Anytime{C<:Contract} <: Contract
    expiry::Date
    c::C
end

#
# Fixed Income
#

"""
`AbstractFixedIncomeContract` identifies fixed income contracts.

These are contracts that are composed of a set of fixed income events,
where each event wraps a `WhenAt` primitive contract.

# Concrete Types:

* `FixedIncomeEvent`

* `FixedIncomeContract`

* `SwapContract`

"""
abstract type AbstractFixedIncomeContract <: Contract end

"""
Defines a financial event of a `FixedIncomeContract`.

# Constructors

```julia
evn = FixedIncomeEvent(:EVENT_NAME, c)
```

where `c` is a `WhenAt` contract.

The event name can be queried by `event_symbol` function.

```julia
evn = FixedIncomeEvent(:EVENT_NAME, Date(2019, 2, 1), 1.0, BRL)
```
"""
struct FixedIncomeEvent{E,C<:Union{Scale, Unit}} <: AbstractFixedIncomeContract
    c::WhenAt{C}

    function FixedIncomeEvent(evn::Symbol, contract::WhenAt{C}) where {C<:Unit}
        @assert Events.is_valid(evn) "$evn is not a valid event."
        return new{evn, C}(contract)
    end

    function FixedIncomeEvent(evn::Symbol, contract::WhenAt{C}) where {C<:Scale}
        @assert contract.maturity == get_horizon(contract) "contract maturity ($(contract.maturity)) should be equal to its horizon ($(get_horizon(contract)))."
        scale = contract.c
        @assert isa(scale.c, Unit) "FixedIncomeEvent WhenAt(scale) should have a Unit as underlying. Got type$(typeof(scale.c)): $(scale.c)."
        return new{evn, C}(contract)
    end
end

"""
A `FixedIncomeContract` is defined as a set of `FixedIncomeEvent`s.
"""
struct FixedIncomeContract <: AbstractFixedIncomeContract
    events::Vector{FixedIncomeEvent}
end

"""
A `SwapContract` has two legs: `asset::FixedIncomeContract`, `liability::FixedIncomeContract`.

The `liability` leg will be valued as `Give(liability)`.
"""
struct SwapContract <: AbstractFixedIncomeContract
    asset::FixedIncomeContract
    liability::FixedIncomeContract
end

include("types_models.jl")

#
# Scenario
#

"""
Provides values for the risk factors.

A Scenario is not required to be linked to a pricing date.

# Interface

* `Base.haskey(s, rf::RiskFactor) :: Bool`

* `Base.getindex(s, rf::CashRiskFactor) :: Cash`

* `Base.getindex(s, rf::NonCashRiskFactor) :: Float64`

Se *writable*:

* `Base.setindex!(s, val, rf::RiskFactor)`

# Provided Methods

* `Base.haskey(scenario::Scenario, c::Currencies.Currency) :: Bool`

* `Base.getindex(scenario::Scenario, c::Currencies.Currency) :: Currencies.Cash`

* `Base.setindex!(scenario::Scenario, val::Currencies.Cash, c::Currencies.Currency)`
"""
abstract type Scenario end

"""
    FixedScenario

A scenario implementation based on a dictionary of fixed values for the risk factors.
"""
struct FixedScenario <: Scenario
    cash_risk_factors::Dict{CashRiskFactor, Currencies.Cash}
    non_cash_risk_factors::Dict{NonCashRiskFactor, Float64}

    function FixedScenario()
        return new(Dict{CashRiskFactor, Currencies.Cash}(), Dict{NonCashRiskFactor, Float64}())
    end
end

"""
    ActualScenario

A scenario implementation that is linked to a `MarketData.AbstractMarketDataProvider`.
"""
struct ActualScenario{P<:MarketData.AbstractMarketDataProvider, F<:Function} <: Scenario
    provider::P
    date::Date
    locf_option::F
end

"""
    CompositeScenario

Builds a new scenario from a list of existing scenarios.

When a risk factor value is requested from this scenario,
it returns the first occurence of the risk factor
in the list of underlying scenarios.

This is useful to override some risk factors from the `ActualScenario`.
"""
struct CompositeScenario <: Scenario
    scenarios::Vector{Scenario}
end

"""
    DebugScenario

Wrapps an existing scenario, and tracks which risk factors
are asked from this scenario.
"""
struct DebugScenario <: Scenario
    scenario::Scenario
    record::Dict{RiskFactor, Union{obstype(RiskFactor), Currencies.Cash}}

    function DebugScenario(scenario::Scenario)
        return new(scenario, Dict{RiskFactor, Union{obstype(RiskFactor), Currencies.Cash}}())
    end
end

"""
    ScenarioMap{F<:Function, S<:Scenario} <: Scenario

Derives a new scenario by mapping a function
to an existing scenario.

For a `ScenarioMap`, `getindex` is implemented as :

```julia
Base.getindex(s::ScenarioMap, rf::RiskFactor) = s.f(rf, s.scenario[rf])
```

# Fields

* `f` mapping function: `f(risk_factor, risk_factor_value) -> new_risk_factor_value`.

* `scenario`: base scenario.

# Constructor

```julia
scenario_map = ScenarioMap(f, scenario)
```
"""
struct ScenarioMap{F<:Function, S<:Scenario} <: Scenario
    f::F
    scenario::S
end

#
# Pricer
#

"""
    AbstractPricer

Represents a pricing function along with
the Risk Factors that are arguments of the
pricing function.

# Interface

* `price(p::AbstractPricer, scenario::Scenario) :: Number`

* `exposures(method::AbstractExposuresMethod, p::AbstractPricer, scenario::Scenario) :: ExposureResult`

* `riskfactors(pricer::AbstractPricer) -> itr`

* `get_functional_currency(pricer::AbstractPricer) :: Currencies.Currency`

"""
abstract type AbstractPricer end

"""
    AbstractCashflowPricer

Used for Cashflow projection.
Has the same interface as `AbstractPricer`.

See also [`foreachcashflow`](@ref), [`collect_cashflows`](@ref).
"""
abstract type AbstractCashflowPricer <: AbstractPricer end

#
# Cashflow projection
#

struct Cashflow
    event::Symbol
    maturity::Date
    value::Float64
    currency::Currencies.Currency
end

#
# Methods for Risk Factors exposures
#

"""
Identifies the algorithm to calculate exposures to risk factors.
See [`exposures`](@ref).
"""
abstract type AbstractExposuresMethod end

"""
Linear exposures to risk factors.

``S \\frac{\\partial P}{\\partial S}``
"""
struct DeltaNormalExposuresMethod <: AbstractExposuresMethod
end

"""
Adds second order exposures approximation
to the standard linear exposures.

Known as delta-gamma approximation method.

As an example, for the exposure to Risk Factor `Stock(:PETR4)`,
the gamma exposure (2nd order) is approximated as:

``\\frac{1}{2} \\frac{ \\partial^2 P }{ \\partial S^2 } S^2``

Second order risk factors are wrapped by [`SecondOrderRiskFactor`](@ref).
"""
struct DeltaGammaApproxExposuresMethod <: AbstractExposuresMethod
end


text:
using SeisPlot, SeisAcoustic

# ========================================================
#        testing the code to get source-side wave field
# ========================================================
# 2-layers velocity model
vel = 3000 * ones(101, 301);  # m/s
vel[51:end,:] .= 3500;

# constant density model
rho = 2000 * ones(101, 301);  # kg/m^3

# number of PML layers
npml = 20;

# top boundary condition
free_surface = true;

# vertical and horizontal grid size
dz = 10; dx = 10;

# time step size and maximum modelling length
dt = 0.001; tmax = 2.0;

# organize these parameters into a structure
params = TdParams(rho, vel, free_surface, dz, dx, dt, tmax;
         data_format=Float64, fd_flag="taylor", order=2, npml=20, apml=900.);

# initialize a single source
# src = Source(27, 49, params; ot=0.0, fdom=20.0,
#       type_flag="ricker", amp=100000, location_flag="index");

# initialize multi-sources
isx = collect(5:60:295); ns=length(isx); isz = 2*ones(ns);
ot  = 0.5*rand(ns);
srcs = get_multi_sources(isz, isx, params; amp=100000, ot=ot, fdom=15);

# initialize recordings
irx = collect(1:2:params.nx);
irz = 2 * ones(length(irx));
rec = Recordings(irz, irx, params);

multi_step_forward!(rec, src, params)
SeisPlotTX(rec.p)

# generate recordings and boundary value
path_bnd = joinpath(homedir(), "Desktop/bnd.rsb");
path_wfd = joinpath(homedir(), "Desktop/wfd.rsb");

# multi_step_forward!(rec, src, params; path_bnd=path_bnd, path_wfd=path_wfd)
multi_step_forward!(rec, src, params; path_bnd=path_bnd, path_wfd=path_wfd)
SeisPlotTX(rec.p)

# compute source-side wavefield
pre1 = get_sourceside_wavefield(src, params; iflag=1);
pre2 = get_sourceside_wavefield(src, params; iflag=2);
pre3 = sourceside_reconstruct_forward(path_bnd, src, params);
pre4 = sourceside_reconstruct_backward(path_bnd, path_wfd, src, params);

# check the consistent of the value
norm(pre1 - pre2) / norm(pre1)
norm(pre2 - pre3) / norm(pre2)
norm(pre2 - pre4) / norm(pre2)

it = 399
SeisPlotTX(pre1[:,:,it], cmap="seismic", wbox=9, hbox=3)
SeisPlotTX(pre2[:,:,it], cmap="seismic", wbox=9, hbox=3)
SeisPlotTX(pre3[:,:,it], cmap="seismic", wbox=9, hbox=3)

# ==============================================================================
#            check the analytical gradient against numerical gradient
# ==============================================================================
using SeisPlot, SeisAcoustic, PyPlot, LinearAlgebra

# size of model
nz = 101; nx = 301;

# 2-layers velocity model
vel = 3000 * ones(nz, nx);  # m/s
vel[51:end,:] .= 3500;

# constant density model
rho = 2000 * ones(nz, nx);  # kg/m^3

# top boundary condition
free_surface = false;

# vertical and horizontal grid size
dz = 10; dx = 10;

# time step size and maximum modelling length
dt = 0.001; tmax = 2.0;

# precision
data_format=Float64;
order=5

# organize these parameters into a structure
params = TdParams(rho, vel, free_surface, dz, dx, dt, tmax;
                  data_format=data_format, order=order);

# initialize a source
isz = 2; isx = 150;
src = Source(isz, isx, params; amp=100000, fdom=5);

# isx = collect(5:10:295); ns=length(isx);
# isz = 2   * ones(ns);
# ot  = 0.5 * rand(ns);
# pol = rand(Bool, ns);
# amp = 100000 * ones(ns);
# for i = 1 : ns
#     if !pol[i]
#        amp[i] = - amp[i]
#     end
# end
# srcs = get_multi_sources(isz, isx, params; amp=amp, ot=ot, fdom=5);

# generate observed data
irx = collect(1:2:params.nx);
irz = 2 * ones(length(irx));
dobs= Recordings(irz, irx, params);

# generate recordings and boundary value
multi_step_forward!(dobs, src, params)
SeisPlotTX(dobs.p, cmap="seismic", pclip=98,
           xlabel="Traces", ylabel="Time (s)", dy=0.001)

# ==============================================================================
#                      the gradient based on adjoint method
# initial velocity model
vel0 = 3000 * ones(nz, nx);
params0 = TdParams(rho, vel0, free_surface, dz, dx, dt, tmax;
                   data_format=data_format, order=order);

# generate synthetic data and save the boundary value
path_bnd = joinpath(homedir(), "Desktop/bnd.rsb");
path_wfd = joinpath(homedir(), "Desktop/wfd.rsb");
dsyn = Recordings(irz, irx, params0);

multi_step_forward!(dsyn, src, params0, path_bnd=path_bnd, path_wfd=path_wfd);
SeisPlotTX(dsyn.p, cmap="seismic", pclip=98,
           xlabel="Traces", ylabel="Time (s)", dy=0.001)

# compute the residue
dres = get_residue(dsyn, dobs);
SeisPlotTX(dres.p, cmap="seismic", pclip=98,
           xlabel="Traces", ylabel="Time (s)", dy=0.001)

# get the gradient
@time g_adj = velocity_gradient(dres, path_bnd, path_wfd, src, params0);
g_adj = reshape(g_adj, nz, nx);
SeisPlotTX(g_adj, cmap="seismic", wbox=9, hbox=3)

# ==============================================================================
#                      compute gradient numerically
# initial velocity model as homogenous one
vel0    = 3000 * ones(nz, nx);
delta_m = 1e-7;

# one model parameter's gradient
iz = 51; ix = 151;

# positive velocity model partubation
vel0[iz,ix] = vel0[iz,ix] + delta_m

# model parameter
params0 = TdParams(rho, vel0, free_surface, dz, dx, dt, tmax;
                   data_format=data_format, order=order);

# get synthetic data
dsyn = Recordings(irz, irx, params0);
multi_step_forward!(dsyn, srcs, params0);

# compute the residue
dres = get_residue(dsyn, dobs);
J_plus = 1/2.0 * (norm(dres.p))^2

# negative velocity model partubation
vel0[iz,ix] = vel0[iz,ix] - 2*delta_m

# model parameter
params0 = TdParams(rho, vel0, free_surface, dz, dx, dt, tmax;
                   data_format=data_format, order=order);

# synthetic data
dsyn = Recordings(irz, irx, params0);
multi_step_forward!(dsyn, srcs, params0);

# compute the residue
dres = get_residue(dsyn, dobs);
J_minus = 1/2.0 * (norm(dres.p))^2;

# compute gradient numerically
g_num = (J_plus - J_minus) / (2 * delta_m)
g_adj[iz,ix]
g_adj[iz,ix] - g_num

# some codes used for testing
# multi_step_adjoint!(path_adj, dres, params0, save_flag="pressure")
# (hdr, adj) = read_RSdata(path_adj);
#
# # get the source side wave field
# pre = get_sourceside_wavefield(src, params0);
#
# # compute the gradient analytically
# g_ana = zeros(params.data_format, params.nz, params.nx);
# for it = 1 : params0.nt-1
#     for i2 = 1 : params.nx
#         for i1 = 1 : params.nz
#             g_ana[i1,i2] = g_ana[i1,i2] + pre[i1,i2,it] * adj[i1,i2,it+1] / 2.0
#         end
#     end
# end
#
# g_ana[iz,ix]
# g_num
# SeisPlotTX(g_ana, cmap="seismic", hbox=4, wbox=9)

ext: jl


text:
#=
Contains the pieces to build the stochastic Neoclassical growth model
described in the Arellano, Maliar, Maliar, Tsyrennikov 2016 paper.
=#

#
# The Model
#
"""
The stochastic Neoclassical growth model type contains parameters
which define the model

* α: Capital share in output
* β: Discount factor
* δ: Depreciation rate
* γ: Risk aversion
* ρ: Persistence of the log of the productivity level
* σ: Standard deviation of shocks to log productivity level
* A: Coefficient on C-D production function

* kgrid: Grid over capital
* zgrid: Grid over productivity
* grid: Grid of (k, z) pairs
* eps_nodes: Nodes used to integrate
* weights: Weights used to integrate
* z1: A grid of the possible z1s tomorrow given eps_nodes and zgrid
"""
@with_kw struct NeoclassicalGrowth
    # Parameters
    α::Float64 = 0.36
    β::Float64 = 0.99
    δ::Float64 = 0.02
    γ::Float64 = 2.0
    ρ::Float64 = 0.95
    σ::Float64 = 0.01
    A::Float64 = (1.0/β - (1 - δ)) / α

    # Grids
    kgrid::Vector{Float64} = collect(range(0.9, stop=1.1, length=10))
    zgrid::Vector{Float64} = collect(range(0.9, stop=1.1, length=10))
    grid::Matrix{Float64} = gridmake(kgrid, zgrid)
    eps_nodes::Vector{Float64} = qnwnorm(5, 0.0, σ^2)[1]
    weights::Vector{Float64} = qnwnorm(5, 0.0, σ^2)[2]
    z1::Matrix{Float64} = (zgrid.^(ρ))' .* exp.(eps_nodes)
end

# Helper functions
f(ncgm::NeoclassicalGrowth, k, z) = @. z * (ncgm.A * k^ncgm.α)
df(ncgm::NeoclassicalGrowth, k, z) = @. ncgm.α * z * (ncgm.A * k^(ncgm.α - 1.0))

u(ncgm::NeoclassicalGrowth, c) = c > 1e-10 ? @.(c^(1-ncgm.γ)-1)/(1-ncgm.γ) : -1e10
du(ncgm::NeoclassicalGrowth, c) = c > 1e-10 ? c.^(-ncgm.γ) : 1e10
duinv(ncgm::NeoclassicalGrowth, u) = u .^ (-1 / ncgm.γ)

expendables_t(ncgm::NeoclassicalGrowth, k, z) = (1-ncgm.δ)*k + f(ncgm, k, z)

# Types for solution methods
abstract type SolutionMethod end

struct IterateOnPolicy <: SolutionMethod end
struct VFI_ECM <: SolutionMethod end
struct VFI_EGM <: SolutionMethod end
struct VFI <: SolutionMethod end
struct PFI_ECM <: SolutionMethod end
struct PFI <: SolutionMethod end
struct dVFI_ECM <: SolutionMethod end
struct EulEq <: SolutionMethod end

#
# Type for Approximating Value and Policy
#
mutable struct ValueCoeffs{T <: SolutionMethod,D <: Degree}
    d::D
    v_coeffs::Vector{Float64}
    k_coeffs::Vector{Float64}
end

function ValueCoeffs(::Type{Val{d}}, method::T) where T <: SolutionMethod where d
    # Initialize two vectors of zeros
    deg = Degree{d}()
    n = n_complete(2, deg)
    v_coeffs = zeros(n)
    k_coeffs = zeros(n)

    return ValueCoeffs{T,Degree{d}}(deg, v_coeffs, k_coeffs)
end

function ValueCoeffs(
        ncgm::NeoclassicalGrowth, ::Type{Val{d}}, method::T
    ) where T <: SolutionMethod where d
    # Initialize with vector of zeros
    deg = Degree{d}()
    n = n_complete(2, deg)
    v_coeffs = zeros(n)

    # Policy guesses based on k and z
    k, z = ncgm.grid[:, 1], ncgm.grid[:, 2]
    css = ncgm.A - ncgm.δ
    yss = ncgm.A
    c_pol = f(ncgm, k, z) * (css/yss)

    # Figure out what kp is
    k_pol = expendables_t(ncgm, k, z) - c_pol
    k_coeffs = complete_polynomial(ncgm.grid, d) \ k_pol

    return ValueCoeffs{T,Degree{d}}(deg, v_coeffs, k_coeffs)
end

solutionmethod(::ValueCoeffs{T}) where T <:SolutionMethod = T

# A few copy methods to make life easier
Base.copy(vp::ValueCoeffs{T,D}) where T where D =
    ValueCoeffs{T,D}(vp.d, vp.v_coeffs, vp.k_coeffs)

function Base.copy(vp::ValueCoeffs{T1,D}, ::T2) where T1 where D where T2 <: SolutionMethod
    ValueCoeffs{T2,D}(vp.d, vp.v_coeffs, vp.k_coeffs)
end

function Base.copy(
        ncgm::NeoclassicalGrowth, vp::ValueCoeffs{T}, ::Type{Val{new_degree}}
    ) where T where new_degree
    # Build Value and policy matrix
    deg = Degree{new_degree}()
    V = build_V(ncgm, vp)
    k = build_k(ncgm, vp)

    # Build new Phi
    Phi = complete_polynomial(ncgm.grid, deg)
    v_coeffs = Phi \ V
    k_coeffs = Phi \ k

    return ValueCoeffs{T,Degree{new_degree}}(deg, v_coeffs, k_coeffs)
end

"""
Updates the coefficients for the value function inplace in `vp`
"""
function update_v!(vp::ValueCoeffs, new_coeffs::Vector{Float64}, dampen::Float64)
    vp.v_coeffs = (1-dampen)*vp.v_coeffs + dampen*new_coeffs
end

"""
Updates the coefficients for the policy function inplace in `vp`
"""
function update_k!(vp::ValueCoeffs, new_coeffs::Vector{Float64}, dampen::Float64)
    vp.k_coeffs = (1-dampen)*vp.k_coeffs + dampen*new_coeffs
end

"""
Builds either V or dV depending on the solution method that is given. If it
is a solution method that iterates on the derivative of the value function
then it will return derivative of the value function, otherwise the
value function itself
"""
build_V_or_dV(ncgm::NeoclassicalGrowth, vp::ValueCoeffs) =
    build_V_or_dV(ncgm, vp, solutionmethod(vp)())

build_V_or_dV(ncgm, vp::ValueCoeffs, ::SolutionMethod) = build_V(ncgm, vp)
build_V_or_dV(ncgm, vp::ValueCoeffs, T::dVFI_ECM) = build_dV(ncgm, vp)

function build_dV(ncgm::NeoclassicalGrowth, vp::ValueCoeffs)
    Φ = complete_polynomial(ncgm.grid, vp.d, Derivative{1}())
    Φ*vp.v_coeffs
end

function build_V(ncgm::NeoclassicalGrowth, vp::ValueCoeffs)
    Φ = complete_polynomial(ncgm.grid, vp.d)
    Φ*vp.v_coeffs
end

function build_k(ncgm::NeoclassicalGrowth, vp::ValueCoeffs)
    Φ = complete_polynomial(ncgm.grid, vp.d)
    Φ*vp.k_coeffs
end

#
# Helper functions
#
function compute_EV!(cp_kpzp::Vector{Float64}, ncgm::NeoclassicalGrowth,
                     vp::ValueCoeffs, kp, iz)
    # Pull out information from types
    z1, weightsz = ncgm.z1, ncgm.weights

    # Get number nodes
    nzp = length(weightsz)

    EV = 0.0
    for izp in 1:nzp
        zp = z1[izp, iz]
        complete_polynomial!(cp_kpzp, [kp, zp], vp.d)
        EV += weightsz[izp] * dot(vp.v_coeffs, cp_kpzp)
    end

    return EV
end

function compute_EV(ncgm::NeoclassicalGrowth, vp::ValueCoeffs, kp, iz)
    cp_kpzp = Array{Float64}(undef, n_complete(2, vp.d))

    return compute_EV!(cp_kpzp, ncgm, vp, kp, iz)
end

function compute_EV(ncgm::NeoclassicalGrowth, vp::ValueCoeffs)
    # Get length of k and z grids
    kgrid, zgrid = ncgm.kgrid, ncgm.zgrid
    nk, nz = length(kgrid), length(zgrid)
    temp = Array{Float64}(undef, n_complete(2, vp.d))

    # Allocate space to store EV
    EV = Array{Float64}(undef, nk*nz)

    _inds = LinearIndices((nk, nz))

    for ik in 1:nk, iz in 1:nz
        # Pull out states
        k = kgrid[ik]
        z = zgrid[iz]
        ikiz_index = _inds[ik, iz]

        # Pass to scalar EV
        complete_polynomial!(temp, [k, z], vp.d)
        kp = dot(vp.k_coeffs, temp)
        EV[ikiz_index] = compute_EV!(temp, ncgm, vp, kp, iz)
    end

    return EV
end

function compute_dEV!(cp_dkpzp::Vector,  ncgm::NeoclassicalGrowth,
                      vp::ValueCoeffs, kp, iz)
    # Pull out information from types
    z1, weightsz = ncgm.z1, ncgm.weights

    # Get number nodes
    nzp = length(weightsz)

    dEV = 0.0
    for izp in 1:nzp
        zp = z1[izp, iz]
        complete_polynomial!(cp_dkpzp, [kp, zp], vp.d, Derivative{1}())
        dEV += weightsz[izp] * dot(vp.v_coeffs, cp_dkpzp)
    end

    return dEV
end

function compute_dEV(ncgm::NeoclassicalGrowth, vp::ValueCoeffs, kp, iz)
    compute_dEV!(Array{Float64}(undef, n_complete(2, vp.d)), ncgm, vp, kp, iz)
end

function env_condition_kp!(cp_out::Vector{Float64}, ncgm::NeoclassicalGrowth,
                           vp::ValueCoeffs, k::Float64, z::Float64)
    # Compute derivative of VF
    dV = dot(vp.v_coeffs, complete_polynomial!(cp_out, [k, z], vp.d, Derivative{1}()))

    # Consumption is then computed as
    c = duinv(ncgm, dV / (1 - ncgm.δ .+ df(ncgm, k, z)))

    return expendables_t(ncgm, k, z) - c
end

function env_condition_kp(ncgm::NeoclassicalGrowth, vp::ValueCoeffs,
                          k::Float64, z::Float64)
    cp_out = Array{Float64}(undef, n_complete(2, vp.d))
    env_condition_kp!(cp_out, ncgm, vp, k, z)
end

function env_condition_kp(ncgm::NeoclassicalGrowth, vp::ValueCoeffs)
    # Pull out k and z from grid
    k = ncgm.grid[:, 1]
    z = ncgm.grid[:, 2]

    # Create basis matrix for entire grid
    dPhi = complete_polynomial(ncgm.grid, vp.d, Derivative{1}())

    # Compute consumption
    c = duinv(ncgm, (dPhi*vp.v_coeffs) ./ (1-ncgm.δ.+df(ncgm, k, z)))

    return expendables_t(ncgm, k, z) .- c
end

max_stars_repo_name: sglyon/CLMMJuliaPythonMatlab

max_issues_repo_name: sglyon/CLMMJuliaPythonMatlab

max_forks_repo_name: sglyon/CLMMJuliaPythonMatlab


text:
module MyInterpolations

export my_lin_interp

function my_lin_interp(grid, vals)
  function func(x)
    for i in 1:(length(grid)-1)
      if (grid[i]<=x)&&(x<=grid[i+1])
        return ((vals[i+1]-vals[i])/(grid[i+1]-grid[i]))*(x-grid[i])+vals[i]
      end
    end
  end

  return func
end

end


text:
using kRPC
using Base.Test

# It's an RPC library. Testing this without a Kerbal instance around is not really possible....

# write your own tests here
@test 1 == 1


text:
"""
Maps from equation number to a Dict: variable -> Set(time_periods)
"""
struct IncidenceTable
    by_eq::Dict{Int,Dict{Symbol,Set{Int}}}
    by_var::Dict{Symbol,Set{Int}}
    by_date::Dict{Int,Set{Symbol}}
end

IncidenceTable() = IncidenceTable(Dict(), Dict(), Dict())

function IncidenceTable(eqs::AbstractVector, skip::AbstractVector=Symbol[])
    # create incidence
    it = IncidenceTable()
    for (i, eq) in enumerate(eqs)
        visit!(it, eq, i, 0, skip)
        # make sure the Dict by_eq has an entry for this equation, even if it
        # is empty.
        get!(it.by_eq, i, Dict{Symbol,Set{Int}}())
    end
    it
end

function IncidenceTable(eq::Expr, skip::AbstractVector=Symbol[])
    it = IncidenceTable([eq], skip)
end

Base.getindex(it::IncidenceTable, i::Int) = it.by_date[i]
Base.getindex(it::IncidenceTable, s::Symbol) = it.by_var[s]

function visit!(it::IncidenceTable, s::Symbol, n::Int, shift::Int,
                skip::AbstractVector=Symbol[])

    if s in skip
        return nothing
    end

    for_eq = get!(it.by_eq, n, Dict{Symbol,Set{Int}}())
    for_sym = get!(for_eq, s, Set{Int}())
    push!(for_sym, shift)

    # update by_var
    push!(get!(it.by_var, s, Set{Int}()), shift)

    # upate by_date
    push!(get!(it.by_date, shift, Set{Symbol}()), s)

    nothing
end

# don't worry about anything else
function visit!(it::IncidenceTable, s::Any, n::Int, shift::Int,
                skip::AbstractVector=Symbol[])
    nothing
end

function visit!(it::IncidenceTable, ex::Expr, n::Int, shift::Int,
                skip::AbstractVector=Symbol[])

    if ex.head == :line
        # this is a line number node in a quoted expression. There is nothing
        # for us to do here.
        return nothing
    end
    if is_time_shift(ex)
        var = ex.args[1]
        i = ex.args[2]
        visit!(it, var, n, i+shift, skip)
        return nothing
    end

    # don't visit the function name
    if ex.head == :call
        for an_arg in ex.args[2:end]
            visit!(it, an_arg, n, shift, skip)
        end
        return nothing
    end

    # Otherwise just visit everything
    for an_arg in ex.args
        visit!(it, an_arg, n, shift, skip)
    end
    nothing
end

max_issues_repo_name: EconForge/Dolang


text:
module MaxHelpingHandMadelineSilhouetteTrigger

using ..Ahorn, Maple

@mapdef Trigger "MaxHelpingHand/MadelineSilhouetteTrigger" MadelineSilhouetteTrigger(x::Integer, y::Integer,
    width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight, enable::Bool=true)

const placements = Ahorn.PlacementDict(
    "Madeline Silhouette (max480's Helping Hand)" => Ahorn.EntityPlacement(
        MadelineSilhouetteTrigger,
        "rectangle"
    )
)

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
@inline function discretizestate!(body,x1,q1,v1,v2,ω1,ω2,Δt)
    state = body.state

    state.xc = x1
    state.qc = q1
    state.vc = v1
    state.ωc = ω1
    state.vsol[2] = v2
    state.ωsol[2] = ω2

    state.xk[1] = x1
    state.xk[2] = x1 + v1*Δt
    state.qk[1] = q1
    state.qk[2] = q1 * ωbar(ω1,Δt) * Δt / 2

    state.xsol[2] = state.xk[2]
    state.qsol[2] = state.qk[2]

    return
end


text:
#=
Only Whiskers are identified that are inside a region of interest (ROI)
That ROI is a circle centered on the whisker pad of radius 100
=#

function apply_roi(whiskers::Array{Whisker1,1},pad_pos::Tuple{Float32,Float32})

    remove_whiskers=Array{Int64,1}()

    #for i=1:length(whiskers)

        #if sqrt((pad_pos[1]-whiskers[i].x[end])^2+(pad_pos[2]-whiskers[i].y[end])^2)>100.0
            #push!(remove_whiskers,i)
        #end
    #end

    #deleteat!(whiskers,remove_whiskers)

    nothing
end

function total_frames(tt,fps)
    h=Base.Dates.hour(tt)
    m=Base.Dates.minute(tt)
    s=Base.Dates.second(tt)
    (h*3600+m*60+s)*fps + 1
end

function frames_between(tt1,tt2,fps)
    total_frames(tt2,fps)-total_frames(tt1,fps)
end

function get_follicle(han::Tracker_Handles)
    x=0.0
    y=0.0
    count=0

    for i=1:length(han.tracked)
        if han.tracked[i]
            x+=han.woi[i].x[end]
            y+=han.woi[i].y[end]
            count+=1
        end
    end
    x=x/count
    y=y/count
    (x,y)
end

#Need to optimize this significantly
#Really we want this to be such that the length of each line is equal.
#Probably can do this by interpolating along each segment or something
#I also need to reverse these

whisker_similarity(han) = whisker_similarity(han,1)

function whisker_similarity(han::Tracker_Handles,prev)

    cor_length=20

    w2_x=han.woi[han.frame-prev].x
    w2_y=han.woi[han.frame-prev].y
    mincor=10000.0.*ones(length(han.wt.whiskers))
    w_id = 0;
    for i=1:length(han.wt.whiskers)
        w1_x=han.wt.whiskers[i].x
        w1_y=han.wt.whiskers[i].y

        for j=1:(size(w1_x,1)-cor_length)
            for k=1:(size(w2_x,1)-cor_length)
                mycor=w_dist(w1_x,w2_x,w1_y,w2_y,j,k,cor_length)
                if mycor < mincor[i]
                    mincor[i]=mycor
                end
            end
        end
    end
    findmin(mincor)
end

function w_dist(x1,x2,y1,y2,i1,i2,cor_length)

    mysum=0.0
    for jj=1:cor_length
        mysum += sqrt((x2[i2+jj] - x1[i1+jj]) ^ 2 + (y2[i2+jj] - y1[i1+jj]) ^2)
    end
    mysum
end

function smooth(x, window_len=7)
    w = getfield(DSP.Windows, :lanczos)(window_len)
    DSP.filtfilt(w ./ sum(w), [1.0], x)
end

function assign_woi(han::Tracker_Handles)

    han.woi[han.displayed_frame] = deepcopy(han.wt.whiskers[han.woi_id])

    nothing
end

function get_frame_list(han::Tracker_Handles)
    sort(collect(keys(han.woi)))
end

function get_woi_array(han::Tracker_Handles)
    [han.woi[i] for i in get_frame_list(han)]
end

function get_frame_index(woi,frame_num)
    findfirst(sort(collect(keys(woi))).==frame_num)
end

function WT_length_constraint(whiskers::Array{Whisker1,1},min_length::Int)

    remove_whiskers=Array{Int64,1}()

    #length constraint
    for i=1:length(whiskers)
        if whiskers[i].len<min_length
            push!(remove_whiskers,i)
        end
    end

    deleteat!(whiskers,remove_whiskers)

    nothing
end

#Order anterior to posterior
function reorder_whiskers(wt::Tracker)

    for i=1:length(wt.all_whiskers)

        xpos=[wt.all_whiskers[i][j].x[end] for j=1:length(wt.all_whiskers[i])]
        wt.all_whiskers[i]=wt.all_whiskers[i][sortperm(xpos)]
    end

    nothing
end

function WT_reorder_whisker(whiskers::Array{Whisker1,1},pad_pos::Tuple{Float32,Float32})

    #order whiskers so that the last index is closest to the whisker pad
    for i=1:length(whiskers)
        WT_reorder_whisker(whiskers[i],pad_pos)
    end

    nothing
end

function WT_reorder_whisker(whisker::Whisker1,pad_pos::Tuple{Float32,Float32})

    front_dist = (whisker.x[1]-pad_pos[1])^2+(whisker.y[1]-pad_pos[2])^2
    end_dist = (whisker.x[end]-pad_pos[1])^2+(whisker.y[end]-pad_pos[2])^2

    if front_dist < end_dist #
        reverse!(whisker.x,1)
        reverse!(whisker.y,1)
        reverse!(whisker.scores,1)
        reverse!(whisker.thick,1)
    end
    nothing
end

function eliminate_redundant(whiskers::Array{Whisker1,1},keep_thres=20.0)

    i=1

    while i<length(whiskers)
        w2=[whiskers[i].x[(end-19):end] whiskers[i].y[(end-19):end]]

        mincor=10000.0
        w_id = 0;
        for j=1:length(whiskers)

            if j != i
                w1=[whiskers[j].x[(end-19):end] whiskers[j].y[(end-19):end]]

                mycor=euclidean(w1,w2)
                if mycor < mincor
                    mincor=mycor
                    w_id = i
                end
            end
            if mincor<keep_thres
                w1_score=mean(whiskers[j].scores)
                w2_score=mean(whiskers[i].scores)

                if w1_score > w2_score
                    deleteat!(whiskers,i)
                else
                    deleteat!(whiskers,j)
                end

                i=1
                break
            end
        end
        i+=1
    end

    nothing
end


text:
using Random

struct Gaussianθ
   A::Float64
   B::Float64
   C::Float64
   D::Float64
end

mutable struct Float64Particle
  x::Float64
  Float64Particle() = new()
end

function makeGaussianModel(θ::Gaussianθ, ys::Vector{Float64} = Vector{Float64}(undef, 0))
  n::Int64 = length(ys)
  sqrtB::Float64 = sqrt(θ.B)
  invDover2::Float64 = 0.5/θ.D
  logncG::Float64 = -0.5 * log(2 * π * θ.D)
  @inline function G(p::Int64, px::Float64Particle)
    @inbounds v::Float64 = θ.C * px.x - ys[p]
    lG::Float64 = logncG - v * invDover2 * v
    return exp(lG)
  end
  @inline function μ!(px::Float64Particle, rng::RNG) where RNG <: AbstractRNG
    px.x = θ.A + sqrtB*randn(rng)
  end
  @inline function simulate(px::Float64Particle, rng::RNG) where RNG <: AbstractRNG
    return θ.C*px.x + sqrt(θ.D)*randn(rng)
  end
  return Model(μ!, G, n, simulate, Float64Particle, Float64)
end

@inline function ldnorm(x, mu, var)
  lnc = -0.5*log(2*π*var)
  return lnc - 0.5/var*(x-mu)^2
end

function gaussianExactLL(θ::Gaussianθ, ys::Vector{Float64})
  n::Int64 = length(ys)
  v::Float64 = 0.0
  for p in 1:n
    v += ldnorm(ys[p], θ.A*θ.C, θ.C^2*θ.B+θ.D)
  end
  return v
end

# function defaultGaussianModel(n::Int64)
#   θ = Gaussianθ(0.0, 1.5, 2.5, 0.5)
#   model = makeGaussianModel(θ)
#   ys = simulateModel(model, n)
#   model = makeGaussianModel(θ, ys)
#   exactLL = gaussianExactLL(θ, ys)
#   return model, θ, ys, exactLL
# end


text:
using GeoPhyInv
using Calculus
using Test
#include("core.jl")

nx=5
nz=7
nznx=nz*nx

p=randn(nznx)
σ=randn(nz,nx)
Q=randn(nz,nx)

function func_test(x)
	x=reshape(x,nz,nx)
	paA=J.Poisson.Param(nx,nz,x)

	global p
	psi=(paA.A*p)

	# some dummy functional
	f=sum(abs2,psi)
	return f
end

function grad_test(x)
	global nznx, nz, nx,p
	x=reshape(x,nz,nx)
	paA=J.Poisson.Param(nx,nz,x)

	psi=paA.A*p

	# gradient, that involves dAdx for testing
	g1= vec(2.0*(paA.A)*p*transpose(p))'*hcat(vec.(paA.dAdx)...)
	return reshape(g1,nz,nx)
end

f = x->func_test(x);

gcalc = Calculus.gradient(f);
g1=reshape(gcalc(vec(Q)),nz,nx);

g2=grad_test(Q)

@test g1≈g2


text:
"""
    BatchView(data, batchsize; partial=true, collate=nothing)
    BatchView(data; batchsize=1, partial=true, collate=nothing)

Create a view of the given `data` that represents it as a vector
of batches. Each batch will contain an equal amount of
observations in them. The batch-size
can be specified using the  parameter `batchsize`.
In the case that the size of the dataset is not dividable by the
specified `batchsize`, the remaining observations will
be ignored if `partial=false`. If  `partial=true` instead
the last batch-size can be slightly smaller.

Note that any data access is delayed until `getindex` is called.

If used as an iterator, the object will iterate over the dataset
once, effectively denoting an epoch.

For `BatchView` to work on some data structure, the type of the
given variable `data` must implement the data container
interface. See [`ObsView`](@ref) for more info.

# Arguments

- **`data`** : The object describing the dataset. Can be of any
    type as long as it implements [`getobs`](@ref) and
    [`numobs`](@ref) (see Details for more information).

- **`batchsize`** : The batch-size of each batch.
    It is the number of observations that each batch must contain
    (except possibly for the last one).

- **`partial`** : If `partial=false` and the number of observations is
    not divisible by the batch-size, then the last mini-batch is dropped.

- **`collate`**: Batching behavior. If `nothing` (default), a batch
    is `getobs(data, indices)`. If `false`, each batch is
    `[getobs(data, i) for i in indices]`. When `true`, applies [`batch`](@ref)
    to the vector of observations in a batch, recursively collating
    arrays in the last dimensions. See [`batch`](@ref) for more information
    and examples.

# Examples

```julia
using MLUtils
X, Y = MLUtils.load_iris()

A = BatchView(X, batchsize=30)
@assert typeof(A) <: BatchView <: AbstractVector
@assert eltype(A) <: SubArray{Float64,2}
@assert length(A) == 5 # Iris has 150 observations
@assert size(A[1]) == (4,30) # Iris has 4 features

# 5 batches of size 30 observations
for x in BatchView(X, batchsize=30)
    @assert typeof(x) <: SubArray{Float64,2}
    @assert numobs(x) === 30
end

# 7 batches of size 20 observations
# Note that the iris dataset has 150 observations,
# which means that with a batchsize of 20, the last
# 10 observations will be ignored
for (x, y) in BatchView((X, Y), batchsize=20, partial=false)
    @assert typeof(x) <: SubArray{Float64,2}
    @assert typeof(y) <: SubArray{String,1}
    @assert numobs(x) == numobs(y) == 20
end

# collate tuple observations
for (x, y) in BatchView((rand(10, 3), ["a", "b", "c"]), batchsize=2, collate=true, partial=false)
    @assert size(x) == (10, 2)
    @assert size(y) == (2,)
end

# randomly assign observations to one and only one batch.
for (x, y) in BatchView(shuffleobs((X, Y)), batchsize=20)
    @assert typeof(x) <: SubArray{Float64,2}
    @assert typeof(y) <: SubArray{String,1}
end
```
"""
struct BatchView{TElem,TData,TCollate} <: AbstractDataContainer
    data::TData
    batchsize::Int
    count::Int
    partial::Bool
    imax::Int
end

function BatchView(data::T; batchsize::Int=1, partial::Bool=true, collate=Val(nothing)) where {T}
    n = numobs(data)
    if n < batchsize
        @warn "Number of observations less than batch-size, decreasing the batch-size to $n"
        batchsize = n
    end
    collate = collate isa Val ? collate : Val(collate)
    if !(collate ∈ (Val(nothing), Val(true), Val(false)))
        throw(ArgumentError("`collate` must be one of `nothing`, `true` or `false`."))
    end
    E = _batchviewelemtype(data, collate)
    imax = partial ? n : n - batchsize + 1
    count = partial ? ceil(Int, n / batchsize) : floor(Int, n / batchsize)
    BatchView{E,T,typeof(collate)}(data, batchsize, count, partial, imax)
end

_batchviewelemtype(::TData, ::Val{nothing}) where TData =
    Core.Compiler.return_type(getobs, Tuple{TData, UnitRange{Int}})
_batchviewelemtype(::TData, ::Val{false}) where TData =
    Vector{Core.Compiler.return_type(getobs, Tuple{TData, Int})}
_batchviewelemtype(data, ::Val{true}) =
    Core.Compiler.return_type(batch, Tuple{_batchviewelemtype(data, Val(false))})

"""
    batchsize(data) -> Int

Return the fixed size of each batch in `data`.
"""
batchsize(A::BatchView) = A.batchsize

Base.length(A::BatchView) = A.count

Base.@propagate_inbounds function getobs(A::BatchView)
    return _getbatch(A, 1:numobs(A.data))
end

Base.@propagate_inbounds function Base.getindex(A::BatchView, i::Int)
    obsindices = _batchrange(A, i)
    _getbatch(A, obsindices)
end

Base.@propagate_inbounds function Base.getindex(A::BatchView, is::AbstractVector)
    obsindices = union((_batchrange(A, i) for i in is)...)::Vector{Int}
    _getbatch(A, obsindices)
end

function _getbatch(A::BatchView{TElem, TData, Val{true}}, obsindices) where {TElem, TData}
    batch([getobs(A.data, i) for i in obsindices])
end
function _getbatch(A::BatchView{TElem, TData, Val{false}}, obsindices) where {TElem, TData}
    return [getobs(A.data, i) for i in obsindices]
end
function _getbatch(A::BatchView{TElem, TData, Val{nothing}}, obsindices) where {TElem, TData}
    getobs(A.data, obsindices)
end

Base.parent(A::BatchView) = A.data
Base.eltype(::BatchView{Tel}) where Tel = Tel

# override AbstractDataContainer default
Base.iterate(A::BatchView, state = 1) =
    (state > numobs(A)) ? nothing : (A[state], state + 1)

# Helper function to translate a batch-index into a range of observations.
@inline function _batchrange(a::BatchView, batchindex::Int)
    @boundscheck (batchindex > a.count || batchindex < 0) && throw(BoundsError())
    startidx = (batchindex - 1) * a.batchsize + 1
    endidx = min(a.imax, startidx + a.batchsize -1)
    return startidx:endidx
end

function Base.showarg(io::IO, A::BatchView, toplevel)
    print(io, "BatchView(")
    Base.showarg(io, parent(A), false)
    print(io, ", ")
    print(io, "batchsize=$(A.batchsize), ")
    print(io, "partial=$(A.partial)")
    print(io, ')')
    toplevel && print(io, " with eltype ", nameof(eltype(A))) # simplify
end

# --------------------------------------------------------------------

max_stars_repo_name: JuliaML/MLBase

max_issues_repo_name: JuliaML/MLBase

max_forks_repo_name: JuliaML/MLBase


text:
#all types shoud inheret from AbstractAmphitheatre and provide a visualize! function
abstract type AbstractAmphitheatre end

#pose landmark draw property #TODO Move to core?
mutable struct plDrawProp
	scale::Float64
	sphereScale::Float64
	color::RGBA
end

plDrawProp() = plDrawProp(0.3, 0.1, RGBA())

# ============================================================
# --------------------BasicFactorGraphPose--------------------
# ============================================================

#TODO All DFGs now has robotId and sessionId, change to use that
struct BasicFactorGraphPose <: AbstractAmphitheatre
	robotId::String
	sessionId::String
	fg::AbstractDFG

	nodes::Dict{Symbol, AbstractPointPose} #poseId, Softtype
	meanmax::Symbol
	zoffset::Float64
	drawPath::Bool
	#pose drawing propeties
	poseProp::plDrawProp
	#landmark drawing propeties
	landmarkProp::plDrawProp
end

"""
    $(SIGNATURES)
Basic visualizer object to draw poses and landmarks.
"""
function BasicFactorGraphPose(robotId::String, sessionId::String, fg::AbstractDFG;
							  meanmax::Symbol=:max,
						      zoffset::Float64=0.0,
							  drawPath::Bool=false,
							  poseProp::plDrawProp = plDrawProp(0.15, 0.05, RGBA(1,1,0,0.5)),
							  landmarkProp::plDrawProp = plDrawProp(0.2, 0.1, RGBA(0,1,0,0.5)))

    return BasicFactorGraphPose(robotId, sessionId, fg, Dict{Symbol,AbstractPointPose}(),
							    meanmax,
							    zoffset,
							    drawPath,
							    poseProp,
							    landmarkProp)
end

function BasicFactorGraphPose(fg::AbstractDFG;
							  meanmax::Symbol=:max,
						      zoffset::Float64=0.0,
							  drawPath::Bool=false,
							  poseProp::plDrawProp = plDrawProp(0.15, 0.05, RGBA(1,1,0,0.5)),
							  landmarkProp::plDrawProp = plDrawProp(0.2, 0.1, RGBA(0,1,0,0.5)))

    return BasicFactorGraphPose(fg.robotId, fg.sessionId, fg, Dict{Symbol,AbstractPointPose}(),
							    meanmax,
							    zoffset,
							    drawPath,
							    poseProp,
							    landmarkProp)
end

"""
    $(SIGNATURES)
Basic visualizer object visualize! function.
"""
function visualize!(vis::Visualizer, basicfg::BasicFactorGraphPose)::Nothing
	#TODO maybe improve this function to lower memmory allocations

	# TODO
	solveKey = :default

	fg = basicfg.fg

	robotId = basicfg.robotId
	sessionId = basicfg.sessionId

	# get all variables
    vars = getVariables(fg)

	basicfg.drawPath && (trackPoints = Point3f0[])
    # update the variable point-estimate cache
    for v in vars
		vsym = v.label

		ppeDict = DistributedFactorGraphs.getPPEDict(v)
		if haskey(ppeDict, solveKey)
			#TODO find a better way to use PPE
			# with this RoME pose and ampi pose has to be kept consistent.
			#TODO maybe use timestamp?
			# varEst = ppe[basicfg.meanmax]
			#TODO meanmax is ignored an suggested used for now
			sugPPE = getSuggestedPPE(ppeDict[solveKey])

			typestr = split(typeof(getSofttype(v)) |> string, ".")[end]
			typesym = Symbol(typestr)
			nodef = getfield(Amphitheatre, typesym)

			nodestruct = nodef(sugPPE...)
		else
			continue
		end
		# if ppe === nothing
	    #     # get vertex and estimate from the factor graph object
	    #     X = getKDE(v)

	    #     xmx = basicfg.meanmax == :max ? getKDEMax(X) : getKDEMean(X)
		# 	# get the variable type
	    #     typestr = split(solverData(v).softtype |> typeof |> string, ".")[end]
		# 	typesym = Symbol(typestr)

		# 	nodef = getfield(Amphitheatre, typesym)

		# 	#NOTE make sure storage order and softtypes are always the same
		# 	nodestruct = nodef(xmx...)

		# else
		# 	#TODO find a better way to use VariableEstimate
		# 	# with this RoME pose and ampi pose has to be kept consistent.
		# 	#TODO maybe use timestamp?
		# 	varEst = ppe[basicfg.meanmax]

		# 	typestr = split(softtype(v) |> string, ".")[end]
		# 	typesym = Symbol(typestr)
		# 	nodef = getfield(Amphitheatre, typesym)

		# 	nodestruct = nodef(varEst.estimate...)
		# end

		if in(:POSE, v.tags)
			groupsym = :poses
			drawProp = basicfg.poseProp
		elseif in(:LANDMARK, v.tags)
			groupsym = :landmarks
			drawProp = basicfg.landmarkProp

		# TODO remove? assume l.. is landmark and x.. is pose
		elseif string(vsym)[1] == 'l'
			groupsym = :landmarks
			drawProp = basicfg.landmarkProp
		elseif string(vsym)[1] == 'x'
			groupsym = :poses
			drawProp = basicfg.poseProp
		elseif in(:VARIABLE, v.tags)
			@warn "Unknown symbol encountered $vsym"
			groupsym = :variables
		end

		isnewnode = !haskey(basicfg.nodes, vsym)
		if isnewnode
			push!(basicfg.nodes, vsym=>nodestruct)
		else
			basicfg.nodes[vsym] = nodestruct
		end

		visNode!(vis[robotId][sessionId][groupsym][vsym], nodestruct, isnewnode,
				 zoffset = basicfg.zoffset,
				 scale = drawProp.scale,
				 sphereScale = drawProp.sphereScale,
				 color = drawProp.color)

		#FIXME this is bad, but just testing feature TODO figure out how to do it properly
		basicfg.drawPath && groupsym == :poses && push!(trackPoints, Point3f0(xmx[1],xmx[2],0.0))

	# sleep(0.1)

    end

	basicfg.drawPath && setobject!(vis[robotId][sessionId][:track], Object(PointCloud(trackPoints), LineBasicMaterial(), "Line"))

    return nothing
end

max_stars_repo_name: Affie/Amphitheatre

max_issues_repo_name: Affie/Amphitheatre


text:
#inspired by flatwhite syntax theme for Atom
const sand_palette = [colorant"#6494ED", # blue
                      colorant"#73C990", # green
                      colorant"#E2C08D", # brown
                      colorant"#FF6347", # red
                      colorant"#2E2C29", # dark,
                      colorant"#4B4844"] # medium

sand_bg = colorant"#F7F3EE"

sand = PlotTheme(sand_bg,
                    colorant"#E2DCD4",
                    colorant"#CBBFAF",
                    colorant"#725B61",
                    expand_palette(sand_bg, sand_palette),
                    sand_palette[[1,4]])

add_plots_theme(:sand, sand)


text:
"""
    optimal_cell_size(Δ_min, Δ_max, x_c; n_particles = 100, p_norm = Inf)

Uses particle swarm optimization to return a cell size in the range
[Δ_min,Δ_max] such that the nodes `x_c` are nearly integer multiples of the optimized cell size.

# Arguments
- `Δ_min::Float64`: minimum cell size
- `Δ_max::Float64`: maximum cell size
- `x_c::Vector{Float64}`: constraint node coordinates
- `n_particles::Int`: number of particles in the swarm (default: 100)
- `p_norm::Int`: norm used to calculate the misfit objective function (default: `Inf`)
- `range::Bool`: set `true` to return a `LinRange` instead  of cell spacing.

# Examples
```julia-repl
julia> optimal_cell_size(0.1,.3,[-1,-.5,1,1.25,2])
[ Info: finding optimal cell size...
[ Info: optimal cell size found is 0.125 with loss = 0.0
0.125
```

# References
The details of the algorithm are described in:
[Structured Mesh Generation: Open-source automatic nonuniform mesh generation for FDTD simulation](https://ieeexplore.ieee.org/document/7458133)
"""
function optimal_cell_size(Δ_min, Δ_max, x_c; n_particles = 100, p_norm = Inf,range = false)
    x0 = min(x_c...)
    xL = max(x_c...)
    Lx = xL - x0
    M_min = Lx / Δ_max
    M_max = Lx / Δ_min
    objective(M) = norm(M[1] / Lx * x_c - round.(M[1] / Lx * x_c), p_norm)
    @info "finding optimal cell size..."
    M_initial = 0.5*(M_max + M_min)
    op = optimize(objective, [M_initial], ParticleSwarm(; lower = [M_min], upper = [M_max], n_particles = n_particles))
    @info "optimal cell size found is $(Lx/round.(Int,op.minimizer[1])) with loss = $(op.minimum)"
    return range ? LinRange(x0,xL,round.(Int, op.minimizer[1] + 1))  :  Lx / round.(Int, op.minimizer[1])
end


text:
## Cross-layer IP-optical restoration-aware TE
function OPTIMAL(GRB_ENV, edges, capacity, flows, demand, scenarios, fiber_scenarios, T, Tf, OpticalTopo, rerouting_K, ILP_LP, time_limit, solve_or_not)
    printstyled("\n** solving cross-layer RWA-TE Optimal super ILP...\n", color=:yellow)
    nedges = size(edges,1)
    nflows = size(flows,1)
    ntunnels = size(T,1)  # all tunnels routing info on IP edge index
    nscenarios = size(scenarios,1)
    nFibers = length(fiber_scenarios[1])
    FibercapacityCode = OpticalTopo["capacityCode"]
    nwavelength = size(FibercapacityCode,2)

    tunnel_num = 0 # tunnel number for each flow, some flow may have less number of tunnels, take the largest
    for x in 1:length(Tf)
        if size(Tf[x],1) > tunnel_num
            tunnel_num = size(Tf[x],1)
        end
    end

    # fiber links are bidirectional
    uni_IPedges = []
    reverse_IPedges = []
    for edge_index in 1:nedges
        e = findfirst(x -> x == (edges[edge_index][2], edges[edge_index][1], edges[edge_index][3]), edges)
        if edge_index < e
            push!(uni_IPedges, edge_index)
            push!(reverse_IPedges, e)
        end
    end

    AllIPedge = []
    for e in 1:nedges
        push!(AllIPedge, edges[e])
        push!(AllIPedge, (edges[e][2], edges[e][1], edges[e][3]))
    end

    IPBranchRoutingFiber, _, IPbranchIndexAll, IPbrachIndexGroup = AllIPWaveRerouting(OpticalTopo, AllIPedge, rerouting_K)
    nIPedgeBranchAll = length(IPbranchIndexAll)

    # for each failure scenario, generate the IP link health state matrix
    IPhealth = ones(nscenarios, nedges)
    for q in 1:nscenarios
        for e in 1:nedges
            IPhealth[q,e] = scenarios[q][e]
        end
    end

    # for each failure scenario, generate the Fiber health state matrix
    Fiberhealth = ones(nscenarios, nFibers)
    for q in 1:nscenarios
        for b in 1:nFibers
            Fiberhealth[q,b] = fiber_scenarios[q][b]
        end
    end

    # how tunnels are routed on IP links
    L = zeros(ntunnels, nedges)
    for t in 1:ntunnels
        for e in 1:nedges
            if in(e, T[t])
                L[t,e] = 1
            end
        end
    end

    # how each IP branch is routed on fibers
    pai = zeros(nIPedgeBranchAll, nFibers)
    for e in 1:nIPedgeBranchAll  # nIPedgeBranchAll is global indexed
        for f in 1:nFibers
            if in(f, IPBranchRoutingFiber[e])
                pai[e,f] = 1  # this IP branch is routed, but take care it this fiber may break in some scenarios
            end
        end
    end

    # CREATE RESIDUAL TUNNELS BY SCENARIO BY FLOW (References Tf): fiber cut scenarios
    Tsf = []  # residual tunnels after failure
    dTaf = []  # affected / restorable tunnels after failure
    for s in 1:size(scenarios,1)
        sft = []
        aft = []
        for f in 1:size(Tf,1)
            ft = []
            at = []
            for t in 1:size(Tf[f],1)
                up = true
                if (length(T[Tf[f][t]]) == 0)
                    up = false
                end
                for e in T[Tf[f][t]]
                    if scenarios[s][e] == 0
                        up = false
                    end
                end
                if up
                    push!(ft, t)  # this tunnel is not affected by failures
                else
                    push!(at, t)  # this tunnel is affected by failures
                end
            end
            push!(sft, ft)
            push!(aft, at)
        end
        push!(Tsf, sft)
        push!(dTaf, aft)
    end

    println("nscenarios ", nscenarios)
    println("nedges ", nedges)
    println("rerouting_K ", rerouting_K)
    println("nIPedgeBranchAll ", nIPedgeBranchAll)
    println("nFibers ", nFibers)
    println("nwavelength ", nwavelength)
    println("nflows ", nflows)
    println("tunnel_num ", tunnel_num)

    model = Model(() -> Gurobi.Optimizer(GRB_ENV))
    set_optimizer_attribute(model, "OutputFlag", 1)
    set_optimizer_attribute(model, "Threads", 32)
    set_time_limit_sec(model, time_limit)

    solve_runtime = 0
    model_def_time = @elapsed begin
    ## The super ILP formulation

    ## optical variables for the optical layer
    if ILP_LP == 1
        @variable(model, IPBranch_bw[1:nscenarios, 1:nedges, 1:rerouting_K] >= 0, Int)  # bandwidth allocation for all IP multipath branches during restoration
        @variable(model, lambda[1:nscenarios, 1:nIPedgeBranchAll, 1:nFibers, 1:nwavelength] >=0, Bin)  # if IP link's branch use fiber and wavelength
    else
        printstyled("** [super ILP model relaxed to LP]\n", color=:yellow)
        @variable(model, IPBranch_bw[1:nscenarios, 1:nedges, 1:rerouting_K] >= 0)  # bandwidth allocation for all IP multipath branches during restoration
        @variable(model, 0<=lambda[1:nscenarios, 1:nIPedgeBranchAll, 1:nFibers, 1:nwavelength] <=1)  # if IP link's branch use fiber and wavelength
    end

    ## TE variables for the entire IP layer network
    @variable(model, restored_capacity[1:nscenarios, 1:nedges] >= 0)   # the bw after restoration (if no failure original bw) for each IP links
    @variable(model, b[1:nflows] >= 0)  # flow bw
    @variable(model, a[1:nflows,1:tunnel_num] >= 0)  # tunnels bw, tunnels are per flow indexed

    # Equation 18, [IP] sum tunnel bandwidth larger than flow bandwidth
    for f in 1:nflows
        @constraint(model, sum(a[f,t] for t in 1:size(Tf[f],1)) >= b[f])
    end

    # Equation 19, [IP] overlapping flows cannot add up to the capacity of that link
    for e in 1:nedges
        @constraint(model, sum(a[f,t] * L[Tf[f][t],e] for f in 1:nflows for t in 1:size(Tf[f],1)) <= capacity[e])
    end

    # Equation 20, [IP] flow bw no larger than demand
    for f in 1:nflows
        @constraint(model, b[f] <= demand[f])
    end

    # Equation 21, [IP] residual tunnels + dynamic restorable tunnels must be able to carry bandwidth for all scenarios
    for f in 1:nflows
        for q in 1:nscenarios
            @constraint(model, sum(a[f,t] for t in Tsf[q][f]) + sum(a[f,t] for t in dTaf[q][f]) >= b[f])
        end
    end

    # Equation 22, [IP] restorable tunnel bandwidth no larger than restored capacity of links
    for q in 1:nscenarios
        for e in 1:nedges
            @constraint(model, sum(a[f,t] * L[Tf[f][t],e] for f in 1:nflows for t in dTaf[q][f]) <= restored_capacity[q,e])
        end
    end

    # Equation 23, [optical] wavelength resource used only once if the resource is usable
    for q in 1:nscenarios
        for w in 1:nwavelength
            for b in 1:nFibers
                @constraint(model, sum(lambda[q,t,b,w] for t in 1:nIPedgeBranchAll) <= FibercapacityCode[b,w])
            end
        end
    end

    # Equation 24, [optical] translate wavelength usage to IPBranch_bw
    for q in 1:nscenarios
        for e in 1:nedges
            for t in 1:rerouting_K  # t is the index for branches of the failIP link, not global branch index
                if t <= length(IPbrachIndexGroup[e])
                    for f in 1:nFibers
                        @constraint(model, IPBranch_bw[q,e,t]*pai[IPbrachIndexGroup[e][t],f] == sum(lambda[q,IPbrachIndexGroup[e][t],f,w]*Fiberhealth[q,f] for w in 1:nwavelength))
                    end
                end
            end
        end
    end

    # Equation 25, [optical] wavelength continuity
    for q in 1:nscenarios
        for t in 1:nIPedgeBranchAll
            for b in IPBranchRoutingFiber[t]
                for bb in IPBranchRoutingFiber[t]
                    for w in 1:nwavelength
                        @constraint(model, lambda[q,t,b,w]*pai[t,b] == lambda[q,t,bb,w]*pai[t,bb])
                    end
                end
            end
        end
    end

    # Equation 26, [optical] restored capacity should not exceed the initial capacity
    for q in 1:nscenarios
        for e in 1:nedges
            @constraint(model, restored_capacity[q,e] <= capacity[e])
            @constraint(model, restored_capacity[q,e] >= capacity[e]*IPhealth[q,e])
        end
    end

    # Equation 27, [optical] restored bw is the sum of all branches with modulation considered
    for q in 1:nscenarios
        for e in 1:nedges
            @constraint(model, restored_capacity[q,e] == 100*sum(IPBranch_bw[q,e,t] for t in 1:rerouting_K)) # modulation is 100 G per wave
        end
    end

    # Auxiliary: IP links are bidirectional (bandwidth equal)
    for q in 1:nscenarios
        for e in 1:length(uni_IPedges)
            @constraint(model, restored_capacity[q,uni_IPedges[e]] == restored_capacity[q,reverse_IPedges[e]])
        end
    end

    # Auxiliary: if surrogate fiber path is less than k
    for q in 1:nscenarios
        for e in 1:nedges
            for t in 1:rerouting_K
                if t > length(IPbrachIndexGroup[e])
                    @constraint(model, IPBranch_bw[q,e,t] == 0)
                end
            end
        end
    end

    @objective(model, Max, sum(b[i] for i in 1:size(b,1)))
    end  # model_def_time end

    # get problem size in terms of number of variables and constraints
    var_num = num_variables(model)
    println("var_num ", var_num)

    sum_constraint = 0
    xt = list_of_constraint_types(model)
    for x in xt
        println(x)
        cons_num = num_constraints(model, x[1], x[2])
        sum_constraint += cons_num
        println("cons_num ", cons_num)
    end
    println("sum_constraint ", sum_constraint)

    # solve the model or not
    if solve_or_not
        optimize!(model)
        solve_runtime = solve_time(model)
        opt_runtime = solve_runtime + model_def_time
        return value.(a), value.(b), objective_value(model), solve_runtime, opt_runtime, var_num, sum_constraint
    else
        return var_num, sum_constraint
    end
end

max_stars_repo_name: hipersys-team/arrow

max_issues_repo_name: arrow-bow/arrow

max_forks_repo_name: arrow-bow/arrow


text:
using Distributions, SearchDiscoveryReplication, GalacticOptim, NLopt
using BSON: @save, @load
using PDMats # required for BSON to load data again

##############################################
## Simulation for smoothing Weitzman1 (with ranking effects)

# NOTE: Throughout uses default estimation options (see definitions_struct.jl)

# Set parameter values (same as baseline)
β0 = [1.0 , -1.0,3.5]
cs0 = 0.03
cd0 = 0.08
dChars = MvNormal([2., 3.5], [3.0, 1.0])
outDummy = true
dV = Normal(0, 0)
nCons = 2000
# By having many products in total, it is unlikely that the same product is shown to two
# different consumers
nProdTot = 1000000
nProdCons = 30
nA0 = 1
seed = 23

# Create DS model
W = Weitzman(
	β = β0,
	cs = cs0,
	cd = cd0,
	nd = 1,
	nA0 = nA0,
	dChars = dChars,
	outDummy = outDummy,
	seed = seed ,
	dV = dV,
	dE = Normal(),
	cfun = "linear",
	smo = 10,
	options = estimOptions(algo=Opt(:LN_BOBYQA,length(β0)+2),maxIter = 1e8,f_tol = 1e-8, x_tol = 1e-8)
	)

# Generate data
d = genData(W,nCons,nProdTot,nProdCons)

# Estimate and check results at end
res = estimateModel(W,d; startvals = vcat(β0,cs0,cd0),  pars_true = vcat(β0,cs0,cd0))

##############################################
## Simulation for smoothing Weitzman2 (no ranking)

# Create DS model
Wf = Weitzman(
	β = β0,
	cs = cs0,
	cd = 0,
	nd = 1,
	nA0 = nA0,
	dChars = dChars,
	outDummy = outDummy,
	seed = seed ,
	dV = dV,
	dE = Normal(),
	cfun = "linear",
	smo = 10,
	options = estimOptions(algo=Opt(:LN_BOBYQA,length(β0)+1),maxIter = 1e8,f_tol = 1e-8, x_tol = 1e-8)
	)

# Generate data
d = genData(Wf,nCons,nProdTot,nProdCons)

# Estimate and check results at end
res = estimateModel(Wf,d; startvals = vcat(β0,cs0),  pars_true = vcat(β0,cs0))

##############################################
## Simulation for smoothing RS
RS = RandomSearch(
	β = β0,
	cs = cs0+cd0,
	dChars = dChars,
	outDummy = outDummy,
	seed = seed,
	dV = dV,
	dE = Normal(),
	cfun = "linear",
	smo = 10,
	options = estimOptions(algo=Opt(:LN_BOBYQA,length(β0)+1))
	)
d = genData(RS,nCons,nProdTot,nProdCons)

# Estimate and check results at end (fct uses random starting values)
res = estimateModel(RS,d; startvals = vcat(β0,cs0+cd0),  pars_true = vcat(β0,cs0+cd0))

##############################################
## Simulation for smoothing FI
FI = FullInfo(
	β = β0,
	dChars = dChars,
	outDummy = outDummy,
	seed = seed + 1,
	dV = dV,
	dE = Normal(),
	smo = 10,
	options = estimOptions(algo=Opt(:LN_BOBYQA,length(β0)))
	)
d = genData(FI,nCons,nProdTot,nProdCons)
res = estimateModel(FI,d; startvals = β0,  pars_true = β0)


text:
immutable InvFourierLayer <: Mocha.Layer
    bottoms::Vector{Symbol}
    flags::UInt32
    name::AbstractString
    pathkeys::Vector{PathKey}
    timelimit::Float64
    tops::Vector{Symbol}
end

function InvFourierLayer( ;
        bottoms::Vector{Symbol} = Symbol[],
        flags::UInt32 = FFTW.ESTIMATE,
        name::AbstractString = "invfourier",
        pathkeys::Vector{PathKey} = PathKey[],
        timelimit::Float64 = Inf,
        tops::Vector{Symbol} = Symbol[])
    @assert length(bottoms) == 1
    @assert length(tops) == 1
    InvFourierLayer(bottoms, flags, name, pathkeys, timelimit, tops)
end

Mocha.can_do_bp(::InvFourierLayer) = false
Mocha.has_neuron(::InvFourierLayer) = false
Mocha.has_param(::InvFourierLayer) = false


text:
module JuMPExtension

# Simple example of JuMP extension used in the tests to check that JuMP works well with extensions
# The main difference between `JuMP.Model` and `JuMPExtension.MyModel` is the fact that in `add_variable` (resp. `add_constraint`),
# `JuMP.Model` applies the modification to its `moi_backend` field while
# `JuMPExtension.MyModel` stores the `AbstractVariable` (resp. `AbstractConstraint`) in a list.

using JuMP

struct ConstraintIndex
    value::Int # Index in `model.constraints`
end

mutable struct MyModel <: JuMP.AbstractModel
    nextvaridx::Int                                 # Next variable index is nextvaridx+1
    variables::Dict{Int,JuMP.ScalarVariable}       # Map varidx -> variable
    var_to_name::Dict{Int,String}                  # Map varidx -> name
    name_to_var::Union{Dict{String,Int},Nothing}  # Map varidx -> name
    nextconidx::Int                                 # Next constraint index is nextconidx+1
    constraints::Dict{ConstraintIndex,JuMP.AbstractConstraint}      # Map conidx -> variable
    con_to_name::Dict{ConstraintIndex,String}      # Map conidx -> name
    name_to_con::Union{Dict{String,ConstraintIndex},Nothing}                     # Map name -> conidx
    objectivesense::MOI.OptimizationSense
    objective_function::JuMP.AbstractJuMPScalar
    obj_dict::Dict{Symbol,Any}                     # Same that JuMP.Model's field `obj_dict`
    function MyModel()
        return new(
            0,
            Dict{Int,JuMP.AbstractVariable}(),
            Dict{Int,String}(),
            nothing,                        # Variables
            0,
            Dict{ConstraintIndex,JuMP.AbstractConstraint}(),
            Dict{ConstraintIndex,String}(),
            nothing,            # Constraints
            MOI.FEASIBILITY_SENSE,
            zero(JuMP.GenericAffExpr{Float64,MyVariableRef}),
            Dict{Symbol,Any}(),
        )
    end
end
Base.broadcastable(model::MyModel) = Ref(model)

JuMP.object_dictionary(model::MyModel) = model.obj_dict

# Variables
struct MyVariableRef <: JuMP.AbstractVariableRef
    model::MyModel # `model` owning the variable
    idx::Int       # Index in `model.variables`
end
Base.copy(v::MyVariableRef) = v

function Base.:(==)(v::MyVariableRef, w::MyVariableRef)
    return v.model === w.model && v.idx == w.idx
end
Base.broadcastable(v::MyVariableRef) = Ref(v)
JuMP.isequal_canonical(v::MyVariableRef, w::MyVariableRef) = v == w
JuMP.variable_type(::MyModel) = MyVariableRef
function JuMP.add_variable(
    m::MyModel,
    v::JuMP.AbstractVariable,
    name::String = "",
)
    m.nextvaridx += 1
    vref = MyVariableRef(m, m.nextvaridx)
    m.variables[vref.idx] = v
    JuMP.set_name(vref, name)
    return vref
end
function JuMP.add_variable(
    model::MyModel,
    variable::JuMP.VariableConstrainedOnCreation,
    name::String,
)
    var_ref = JuMP.add_variable(model, variable.scalar_variable, name)
    JuMP.add_constraint(model, JuMP.ScalarConstraint(var_ref, variable.set))
    return var_ref
end
function JuMP.add_variable(
    model::MyModel,
    variable::JuMP.VariablesConstrainedOnCreation,
    names,
)
    var_refs =
        JuMP.add_variable.(
            model,
            variable.scalar_variables,
            JuMP.vectorize(names, variable.shape),
        )
    JuMP.add_constraint(model, JuMP.VectorConstraint(var_refs, variable.set))
    return JuMP.reshape_vector(var_refs, variable.shape)
end
function JuMP.delete(model::MyModel, vref::MyVariableRef)
    @assert JuMP.is_valid(model, vref)
    delete!(model.variables, vref.idx)
    return delete!(model.var_to_name, vref.idx)
end
function JuMP.delete(model::MyModel, vrefs::Vector{MyVariableRef})
    return JuMP.delete.(model, vrefs)
end
function JuMP.is_valid(model::MyModel, vref::MyVariableRef)
    return (model === vref.model && vref.idx in keys(model.variables))
end
JuMP.num_variables(m::MyModel) = length(m.variables)

# Internal functions
variable_info(vref::MyVariableRef) = vref.model.variables[vref.idx].info
function update_variable_info(vref::MyVariableRef, info::JuMP.VariableInfo)
    return vref.model.variables[vref.idx] = JuMP.ScalarVariable(info)
end

JuMP.has_lower_bound(vref::MyVariableRef) = variable_info(vref).has_lb
function JuMP.lower_bound(vref::MyVariableRef)::Float64
    @assert !JuMP.is_fixed(vref)
    return variable_info(vref).lower_bound
end
function JuMP.set_lower_bound(vref::MyVariableRef, lower)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            true,
            lower,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            info.integer,
        ),
    )
end
function JuMP.delete_lower_bound(vref::MyVariableRef)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            false,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            info.integer,
        ),
    )
end
JuMP.has_upper_bound(vref::MyVariableRef) = variable_info(vref).has_ub
function JuMP.upper_bound(vref::MyVariableRef)::Float64
    @assert !JuMP.is_fixed(vref)
    return variable_info(vref).upper_bound
end
function JuMP.set_upper_bound(vref::MyVariableRef, upper)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            true,
            upper,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            info.integer,
        ),
    )
end
function JuMP.delete_upper_bound(vref::MyVariableRef)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            false,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            info.integer,
        ),
    )
end
JuMP.is_fixed(vref::MyVariableRef) = variable_info(vref).has_fix
function JuMP.fix_value(vref::MyVariableRef)::Float64
    return variable_info(vref).fixed_value
end
function JuMP.fix(vref::MyVariableRef, value; force::Bool = false)
    info = variable_info(vref)
    if !force && (info.has_lb || info.has_ub)
        error(
            "Unable to fix $(vref) to $(value) because it has existing bounds.",
        )
    end
    update_variable_info(
        vref,
        JuMP.VariableInfo(
            false,
            info.lower_bound,
            false,
            info.upper_bound,
            true,
            value,
            info.has_start,
            info.start,
            info.binary,
            info.integer,
        ),
    )
    return
end
function JuMP.unfix(vref::MyVariableRef)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            false,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            info.integer,
        ),
    )
end
function JuMP.start_value(vref::MyVariableRef)::Union{Nothing,Float64}
    return variable_info(vref).start
end
function JuMP.set_start_value(vref::MyVariableRef, start)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            true,
            start,
            info.binary,
            info.integer,
        ),
    )
end
JuMP.is_binary(vref::MyVariableRef) = variable_info(vref).binary
function JuMP.set_binary(vref::MyVariableRef)
    @assert !JuMP.is_integer(vref)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            true,
            info.integer,
        ),
    )
end
function JuMP.unset_binary(vref::MyVariableRef)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            false,
            info.integer,
        ),
    )
end
JuMP.is_integer(vref::MyVariableRef) = variable_info(vref).integer
function JuMP.set_integer(vref::MyVariableRef)
    @assert !JuMP.is_binary(vref)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            true,
        ),
    )
end
function JuMP.unset_integer(vref::MyVariableRef)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            false,
        ),
    )
end

# Constraints
const MyConstraintRef = JuMP.ConstraintRef{MyModel,ConstraintIndex}
JuMP.constraint_type(::MyModel) = MyConstraintRef
function JuMP.add_constraint(
    model::MyModel,
    c::JuMP.AbstractConstraint,
    name::String = "",
)
    model.nextconidx += 1
    index = ConstraintIndex(model.nextconidx)
    cref = JuMP.ConstraintRef(model, index, JuMP.shape(c))
    model.constraints[index] = c
    JuMP.set_name(cref, name)
    return cref
end
function JuMP.delete(model::MyModel, constraint_ref::MyConstraintRef)
    @assert JuMP.is_valid(model, constraint_ref)
    delete!(model.constraints, constraint_ref.index)
    return delete!(model.con_to_name, constraint_ref.index)
end
function JuMP.delete(model::MyModel, con_refs::Vector{<:MyConstraintRef})
    return JuMP.delete.(model, con_refs)
end
function JuMP.is_valid(model::MyModel, constraint_ref::MyConstraintRef)
    return (
        model === constraint_ref.model &&
        constraint_ref.index in keys(model.constraints)
    )
end
function JuMP.constraint_object(cref::MyConstraintRef)
    return cref.model.constraints[cref.index]
end
function JuMP.num_constraints(
    model::MyModel,
    F::Type{<:JuMP.AbstractJuMPScalar},
    S::Type{<:MOI.AbstractSet},
)
    return count(
        con -> con isa JuMP.ScalarConstraint{F,S},
        values(model.constraints),
    )
end
function JuMP.num_constraints(
    model::MyModel,
    ::Type{<:Vector{F}},
    S::Type{<:MOI.AbstractSet},
) where {F<:JuMP.AbstractJuMPScalar}
    return count(
        con -> con isa JuMP.VectorConstraint{F,S},
        values(model.constraints),
    )
end

# Objective
function JuMP.set_objective_function(m::MyModel, f::JuMP.AbstractJuMPScalar)
    return m.objective_function = f
end
function JuMP.set_objective_function(m::MyModel, f::Real)
    return m.objective_function = JuMP.GenericAffExpr{Float64,MyVariableRef}(f)
end
JuMP.objective_sense(model::MyModel) = model.objectivesense
function JuMP.set_objective_sense(model::MyModel, sense)
    return model.objectivesense = sense
end
JuMP.objective_function_type(model::MyModel) = typeof(model.objective_function)
JuMP.objective_function(model::MyModel) = model.objective_function
function JuMP.objective_function(model::MyModel, FT::Type)
    # InexactError should be thrown, this is needed in `objective.jl`
    if !(model.objective_function isa FT)
        throw(
            InexactError(
                :objective_function,
                FT,
                typeof(model.objective_function),
            ),
        )
    end
    return model.objective_function::FT
end

# Names
JuMP.name(vref::MyVariableRef) = vref.model.var_to_name[vref.idx]
function JuMP.set_name(vref::MyVariableRef, name::String)
    vref.model.var_to_name[vref.idx] = name
    return vref.model.name_to_var = nothing
end
function JuMP.variable_by_name(model::MyModel, name::String)
    if model.name_to_var === nothing
        # Inspired from MOI/src/Utilities/model.jl
        model.name_to_var = Dict{String,Int}()
        for (var, var_name) in model.var_to_name
            if haskey(model.name_to_var, var_name)
                # -1 is a special value that means this string does not map to
                # a unique variable name.
                model.name_to_var[var_name] = -1
            else
                model.name_to_var[var_name] = var
            end
        end
    end
    index = get(model.name_to_var, name, nothing)
    if index isa Nothing
        return nothing
    elseif index == -1
        error("Multiple variables have the name $name.")
    else
        return MyVariableRef(model, index)
    end
end
JuMP.name(cref::MyConstraintRef) = cref.model.con_to_name[cref.index]
function JuMP.set_name(cref::MyConstraintRef, name::String)
    cref.model.con_to_name[cref.index] = name
    return cref.model.name_to_con = nothing
end
function JuMP.constraint_by_name(model::MyModel, name::String)
    if model.name_to_con === nothing
        # Inspired from MOI/src/Utilities/model.jl
        model.name_to_con = Dict{String,ConstraintIndex}()
        for (con, con_name) in model.con_to_name
            if haskey(model.name_to_con, con_name)
                # -1 is a special value that means this string does not map to
                # a unique constraint name.
                model.name_to_con[con_name] = ConstraintIndex(-1)
            else
                model.name_to_con[con_name] = con
            end
        end
    end
    index = get(model.name_to_con, name, nothing)
    if index isa Nothing
        return nothing
    elseif index.value == -1
        error("Multiple constraints have the name $name.")
    else
        # We have no information on whether this is a vector constraint
        # or a scalar constraint
        return JuMP.ConstraintRef(model, index, JuMP.ScalarShape())
    end
end

# Show
function JuMP.show_backend_summary(io::IO, model::MyModel) end
function JuMP.show_objective_function_summary(io::IO, model::MyModel)
    return println(
        io,
        "Objective function type: ",
        JuMP.objective_function_type(model),
    )
end
function JuMP.objective_function_string(print_mode, model::MyModel)
    return JuMP.function_string(print_mode, JuMP.objective_function(model))
end
_plural(n) = (isone(n) ? "" : "s")
function JuMP.show_constraints_summary(io::IO, model::MyModel)
    n = length(model.constraints)
    return print(io, "Constraint", _plural(n), ": ", n)
end
function JuMP.constraints_string(print_mode, model::MyModel)
    strings = String[]
    # Sort by creation order, i.e. ConstraintIndex value
    constraints = sort(collect(model.constraints), by = c -> c.first.value)
    for (index, constraint) in constraints
        push!(strings, JuMP.constraint_string(print_mode, constraint))
    end
    return strings
end

end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
x = asindex(1:2:10, :x)
y = asindex(1:5, :y)

@testset "vcat" begin
    @testset "vcat_keys" begin
        for (x, y, xykeys, xyvals, xyname, xyindex) in ()
            @test vcat_keys(x, y) == xykeys
            @test vcat_keys(keys(x), y) == xykeys
            @test vcat_keys(x, keys(x)) == xykeys
            @test vcat_keys(keys(x), keys(y)) == xykeys

            @test vcat_values(x, y) == xyvals
            @test vcat_values(values(x), y) == xyvals
            @test vcat_values(x, values(x)) == xyvals
            @test vcat_values(values(x), values(y)) == xyvals

            @test vcat_names(x, y) == xyname
            @test vcat_names(dimnames(x), y) == xyname
            @test vcat_names(x, dimnames(x)) == xyname
            @test vcat_names(dimnames(x), dimnames(y)) == xyname

            @test vcat(x, y) == xyname
        end
    end
end


text:
using FITSIO, FITSIO.Libcfitsio

function sxt_correct_level1lbt(lbtfile::String)
    f=fits_open_file(lbtfile,1)
    fits_movabs_hdu(f,2)
    colnum_ccdtemp1 = fits_get_colnum(f,"CCD_Temp_1")
    colnum_ccdtemp2 = fits_get_colnum(f,"CCD_Temp_2")
    nrows = fits_get_num_rows(f)
    ccdtemp1 = Array{Float32}(undef,nrows)
    ccdtemp2 = Array{Float32}(undef,nrows)
    fits_read_col(f,colnum_ccdtemp1,1,1,ccdtemp1)
    fits_read_col(f,colnum_ccdtemp2,1,1,ccdtemp2)
    if (maximum(ccdtemp1) - minimum(ccdtemp1) > maximum(ccdtemp2) - minimum(ccdtemp2))
        cp(lbtfile, lbtfile * ".orig"; force=true)
        fits_update_key(f,"TTYPE24","CCD_Temp_2")
        fits_update_key(f,"TTYPE25","CCD_Temp_1")
        print("CCD_Temp_1 range  greater than CCD_Temp_2 range. CCD_Temp_1 updated.\n")
    else
        print("CCD_Temp_1 range  less or equal to CCD_Temp_2 range. CCD_Temp_1 not changed.\n")
    end
    fits_close_file(f)
    return
end

max_stars_repo_name: gulabd/SXTMerger

max_issues_repo_name: gulabd/SXTMerger


text:
"String for unit multiplier in option line"
const unitstrings = Dict{ Float64, String }(
  1e0 => "Hz",
  1e3 => "kHz",
  1e6 => "MHz",
  1e9 => "GHz"
)
"Default unit value to use in option line."
const default_unit_value = 1e9

"Strings for parameter format symbols in option line."
const parameterstrings = Dict{ Symbol, String }(
  :ScatteringParameters   => "S",
  :AdmittanceParameters   => "Y",
  :ImpedanceParameters    => "Z",
  :HybridHParameters      => "H",
  :HybridGParameters      => "G",
)
"Defaoult parameter format symbo to use in option line."
const default_parameter_symbol = :ScatteringParameters

"String for data format symbols in option line."
const formatstrings = Dict{ Symbol, String }(
 :DecibelAngle      => "dB",
 :MagnitudeAngle    => "MA",
 :RealImaginary     => "RI"
)
"Default data format symbol to use in option line."
const default_format_symbol = :MagnitudeAngle

"Default characteristic iimpedance to use in option line."
const default_resistance_string = "50"

"""
    write_comment_line( comment )

Returns a comment line string for a comment.
"""
write_comment_line( comment ) = "!" * comment

"""
    write_frequency_unit_option( [ value ] )

Returns an option line string for a valid frequency unit multiplier.
"""
write_frequency_unit_option( value = default_unit_value ) = unitstrings[ value ]

"""
    write_parameter_option( [ symbol ] )

Returns an option line string for a valid parameter option symbol.
"""
write_parameter_option( symbol = default_parameter_symbol ) = parameterstrings[ symbol ]

"""
    write_format_option( [ symbol ] )

Returns an option line string for a valid data format symbol.
"""
write_format_option( symbol = default_format_symbol ) = formatstrings[ symbol ]

"""
    write_resistance_option( [ resistance ] )

Returns an option line string for the characteristic impedance.
"""
write_resistance_option( resistance = 50 ) = string( resistance )

"""
     write_option_line( [ options ] )

Returns the option line string for an option structure.
"""
function write_option_line( options::Options = Options() )
  "# " *
  write_frequency_unit_option( options.unit ) * " " *
  write_parameter_option( options.parameter ) * " " *
  write_format_option( options.format ) * " " *
  "R " * write_resistance_option( options.resistance )
end

"""
    compstring( z )

Returns a pair with the real and imaginary components of a complex number.
"""
compstring( z ) = string( real( z ) ) * " " * string( imag( z ) )

"""
    comp2ma( z )

Returns a pair with the linear magnitude and angle in degrees of a complex number.
"""
comp2ma( z ) = ( abs( z ), rad2deg( angle( z ) ) )

"""
    comp2da( z )

Returns a pair with the magnitude in dB and angle in degrees of a complex number.
"""
function comp2da( z )
  m, a = comp2ma( z )
  20log10( m ), a
end

"""
    pairstring( pair )

Returns a string combining the components of a pair.
"""
pairstring( pair ) = string( pair[ 1 ] ) * " " * string( pair[ 2 ] )

"""
    mastring( z )

Returns a string in magnitude / angle format of a complex number.
"""
mastring( z ) = pairstring( comp2ma( z ) )

"""
    dastring( z )

Returns a string in dB / angle format of a complex number.
"""
dastring( z ) = pairstring( comp2da( z ) )

"Holds conversion functions for parameter format symbols when writing."
const WriteConversions = Dict{ Symbol, Function }(
  :RealImaginary  => compstring,
  :MagnitudeAngle => mastring,
  :DecibelAngle   => dastring
)

"""
    write_data( data, N [ options ] )

Returns a string for N-port Touchstone formated lines for a data point.
"""
function write_data( data::DataPoint, N::Integer, options::Options = Options() )
  freqstring = string( data.frequency / options.unit )
  if N == 2
    paras = reshape( data.parameter, 1, 4 )
    lines = 1
  else
    paras = data.parameter
    lines = N
  end
  if options.parameter == :ImpedanceParameters
    paras = map( x -> x / options.resistance, paras )
  end
  if options.parameter == :AdmittanceParameters
    paras = map( x -> x * options.resistance, paras )
  end
  if N == 2
    if options.parameter == :HybridGParameters
      paras[ 1 ] *= options.resistance
      paras[ 4 ] /= options.resistance
    end
    if options.parameter == :HybridHParameters
      paras[ 1 ] /= options.resistance
      paras[ 4 ] *= options.resistance
    end
  end
  strings = map( snn -> WriteConversions[ options.format ]( snn ), paras )
  datastrings = ones( String, lines )
  for row in 1:lines
    # V 1.0 not more than 4 pairs ( + freq ) in a line
    if N > 4
      for col in 1:N
        if col > 0 && col % 4 == 0
          strings[ row, col ] *= "\n"
        end
      end
    end
    datastrings[ row ] = replace( join( strings[ row, : ], " " ), "\n " => "\n" )
  end
  datastrings[ 1 ] = freqstring * " " * datastrings[ 1 ]
  datastring = join( datastrings, "\n" )
  if datastring[ end ] == "\n"
    return datastring
  else
    return datastring * "\n"
  end
end

function writeNoiseData( nd::NoiseDataPoint, options::Options = Options(), version = 1 )
  freqstring = string( nd.frequency / options.unit )
  noisefigstring = 20log10( nd.minNoiseFigure )
  reflcoeffstring = mastring( nd.reflCoeff )
  effnoiseresstring = version == 1 ? string( nd.effNoiseRes / options.resistance ) : string( nd.effNoiseRes )
  res = "$freqstring $noisefigstring $reflcoeffstring $effnoiseresstring\n"
  return res
end

"""
    write_touchstone_stream( stream, ts )

Writes formated Touchstone data for a TouchstoneData structure to a stream.
"""
function write_touchstone_stream( stream::IO, ts::TouchstoneData )
  for comment in ts.comments
    write( stream, write_comment_line( comment ) )
    write( stream, "\n" )
  end

  write( stream, write_option_line( ts.options ) )
  write( stream, "\n" )
  #
  # ports = -1
  # if size( ts.data, 1 ) > 0
  #   ports = size( ts.data[ 1 ].parameter, 1 )
  # end
  #
  for datapoint in ts.data
    write( stream, write_data( datapoint, ports( ts ), ts.options ) )
  end
  for ndp in ts.noiseData
    write( stream, writeNoiseData( ndp, ts.options ) )
  end
end

"""
    write_touchstone_string( ts )

Returns a string with formated Touchstone data for a TouchstoneData structure.
"""
function write_touchstone_string( ts::TouchstoneData )
  buffer = IOBuffer()
  write_touchstone_stream( buffer, ts )
  String( take!( copy( buffer ) ) )
end

"""
    write_touchstone_file( filename, ts )

Writes formated Touchstone data for a TouchstoneData structure to a file.
"""
function write_touchstone_file( filename::String, ts::TouchstoneData )
  open( filename, "w" ) do io
    write_touchstone_stream( io, ts )
  end
end


text:
"""
    TypeC()

Lorem ipsum dolor sit amet.
"""
struct TypeC

end


text:
in = "128392-643281"
sta, sto = parse.(Int, split(in, '-'))

ispass(dif) = any(==(0), dif) && all(>=(0), dif)
ispassword(num, pred) = pred(diff(reverse(digits(num))))
count(x->ispassword(x, ispass), sta:sto)

function ispass2(dif)
    all(>=(0), dif) || return false
    reps = findall(==(0), dif)
    isempty(reps) && return false
    for r ∈ reps
        r-1 ∈ reps || r + 1 ∈ reps || return true
    end
    false
end

count(x->ispassword(x, ispass2), sta:sto)


text:
module GenericMRRR

using LinearAlgebra
using LinearAlgebra: eigsortby, sorteig!, checksquare, reflector!
export geigen!, geigvals!

using Requires: @require

include("symtridiag.jl")
include("dqds.jl")

struct AlgorithmFailure <: Exception
    msg::String
end
struct ConvergenceFailure <: Exception
    msg::String
end
struct RepresentationFailure <: Exception
    msg::String
end

"""
    geigen!(A::SymTridiagonal{T},...) where T <: AbstractFloat -> Eigen

Computes full or partial eigen-decomposition of a symmetric tridiagonal matrix.

Signatures and functionality are (intended to be!) identical to
`LinearAlgebra.eigen!(::SymTridiagonal,...)`, except that these handle more general element
types.
"""
function geigen!(A::SymTridiagonal{T}) where T <: AbstractFloat
    λ, V = _st_schur!(A.dv, A.ev, wantV = true)
    LinearAlgebra.Eigen(sorteig!(λ, V, eigsortby)...)
end
function geigen!(A::SymTridiagonal{T}, irange::UnitRange) where T <: AbstractFloat
    λ, V = _st_schur!(A.dv, A.ev, wantV = true,
                     select=IndexedEigvalSelector(first(irange),last(irange)))
    LinearAlgebra.Eigen(sorteig!(λ, V, eigsortby)...)
end
function geigen!(A::SymTridiagonal{T}, vl::Real, vu::Real) where T <: AbstractFloat
    λ, V = _st_schur!(A.dv, A.ev, wantV = true,
                     select=IntervalEigvalSelector(vl,vu))
    LinearAlgebra.Eigen(sorteig!(λ, V, eigsortby)...)
end
"""
    geigvals!(A::SymTridiagonal{T},...) where T <: AbstractFloat -> Vector{T}

Computes some or all eigenvalues of a symmetric tridiagonal matrix.

Signatures and functionality are (intended to be!) identical to
`LinearAlgebra.eigvals!(::SymTridiagonal,...)`, except that these handle more general element
types.
"""
function geigvals!(A::SymTridiagonal{T}) where T <: AbstractFloat
    λ, _ = _st_schur!(A.dv, A.ev, wantV = false)
    return λ
end
function geigvals!(A::SymTridiagonal{T}, irange::UnitRange) where T <: AbstractFloat
    λ, _ = _st_schur!(A.dv, A.ev, wantV = false,
                     select=IndexedEigvalSelector(first(irange),last(irange)))
    return λ
end
function geigvals!(A::SymTridiagonal{T}, vl::Real, vu::Real) where T <: AbstractFloat
    λ, _ = _st_schur!(A.dv, A.ev, wantV = false,
                     select=IntervalEigvalSelector(vl,vu))
    return λ
end

function __init__()
    @require GenericLinearAlgebra="14197337-ba66-59df-a3e3-ca00e7dcff7a" begin
        include("gla.jl")
    end
end

end # module


text:
using QuantumOptics, CollectiveSpins
using PyCall
@pyimport matplotlib.pyplot as plt

const T = [0:0.01:1]
const γ = 1.
const edipole = [0,0,1.]
const D = range(1.65, stop=1.67, length=20)

# Initial state (Bloch state)
const phi = 0.
const theta = pi/2.

for d=D
    println(d)
    systemgeometry = CollectiveSpins.geometry.cube(d)
    system = SpinCollection(systemgeometry, edipole, γ)
    N = length(system.spins)

    td_ind = Float64[]
    td_mf = Float64[]
    td_mpc = Float64[]

    td2_ind = Float64[]
    td2_mf = Float64[]

    # Independent
    state0 = CollectiveSpins.independent.blochstate(phi,theta,N)
    tout, state_ind_t = CollectiveSpins.independent.timeevolution(T, system, state0)

    # Meanfield
    state0 = CollectiveSpins.meanfield.blochstate(phi,theta,N)
    tout, state_mf_t = CollectiveSpins.meanfield.timeevolution(T, system, state0)

    # Meanfield + Correlations
    state0 = CollectiveSpins.mpc.blochstate(phi,theta,N)
    tout, state_mpc_t = CollectiveSpins.mpc.timeevolution(T, system, state0)

    # Quantum: master equation

    function fout(t, rho::Operator)
        i = findfirst(T, t)
        rho_ind = CollectiveSpins.independent.densityoperator(state_ind_t[i])
        rho_mf = CollectiveSpins.meanfield.densityoperator(state_mf_t[i])
        rho_mpc = CollectiveSpins.mpc.densityoperator(state_mpc_t[i])
        push!(td_ind, tracedistance(rho, rho_ind))
        push!(td_mf, tracedistance(rho, rho_mf))
        push!(td_mpc, tracedistance(rho, rho_mpc))
        push!(td2_ind, tracedistance(rho_mpc, rho_ind))
        push!(td2_mf, tracedistance(rho_mpc, rho_mf))
    end

    Ψ₀ = CollectiveSpins.quantum.blochstate(phi,theta,N)
    ρ₀ = Ψ₀⊗dagger(Ψ₀)
    CollectiveSpins.quantum.timeevolution(T, system, ρ₀, fout=fout)
    plt.figure(1)
    plt.plot(d, maximum(td_ind), ".", lw=3.0, color="darkorange", label="Independent")
    plt.plot(d, maximum(td_mf), ".", lw=3.0, color="gray", label="meanfield")
    plt.plot(d, maximum(td_mpc), ".", lw=1.5, color="navy", label="MPC")
    plt.figure(2)
    plt.plot(d, maximum(td2_ind), ".", lw=1.5, color="green", label="MPC")
    plt.plot(d, maximum(td2_mf), ".", lw=1.5, color="red", label="MPC")
end
plt.show()