text:
module TestUtils

export gradcheck

using YAAD
import YAAD: Value

make_jacobian(input::AbstractArray, num_out::Int) = (zeros(eltype(input), length(input), num_out), )
make_jacobian(input::Variable, num_out::Int) = make_jacobian(value(input), num_out)
make_jacobian(inputs::Tuple, num_out::Int) = Tuple(first(make_jacobian(each, num_out)) for each in inputs)

zero_like(x::T) where {T <: Number} = zero(T)
zero_like(x) = fill!(similar(x), 0)
zero_like(x::Broadcast.Broadcasted) = zero_like(Broadcast.materialize(x))
zero_like(x::Value) = zero_like(value(x))

function get_numerical_jacobian(f, inputs...; target=inputs, eps=1e-3)
    output_size = length(value(f(inputs...)))
    jacobian = make_jacobian(target, output_size)

    x_tensors = [value(t) for t in target]
    for (x_tensor, d_tensor) in zip(x_tensors, jacobian)
        for (d_idx, x_idx) in enumerate(eachindex(x_tensor))
            orig = x_tensor[x_idx]
            x_tensor[x_idx] = orig - eps
            outa = copy(value(f(inputs...)))
            x_tensor[x_idx] = orig + eps
            outb = copy(value(f(inputs...)))
            x_tensor[x_idx] = orig

            r = (outb - outa) / (2 * eps)
            _set_numerical_jacobian_elem!(d_tensor, d_idx, r)
        end
    end
    return jacobian
end

_set_numerical_jacobian_elem!(d_tensor, d_idx, r::AbstractArray) = d_tensor[d_idx, :] = reshape(r, :)
_set_numerical_jacobian_elem!(d_tensor, d_idx, r::Number) = d_tensor[d_idx, 1] = r

function get_analytical_jacobian(inputs, output::Value)
    jacobian = make_jacobian(inputs, length(value(output)))
    jacobian_reentrant = make_jacobian(inputs, length(value(output)))
    grad_output = zero_like(output)
    reentrant = true
    correct_grad_sizes = true

    _analytical_jacobian!(jacobian, jacobian_reentrant, grad_output, output)

    for (jacobian_x, jacobian_reentrant_x) in zip(jacobian, jacobian_reentrant)
        length(jacobian_x) != 0 && maximum(abs.(jacobian_x - jacobian_reentrant_x)) != 0 && (reentrant = false)
    end

    jacobian, reentrant, correct_grad_sizes
end

function _analytical_jacobian!(jacobian::Tuple, jacobian_reentrant::Tuple, grad_output::AbstractArray, output::Value)
    for i in 1:length(grad_output)
        grad_output = fill!(grad_output, 0)
        grad_output[i] = 1
        for jacobian_c in (jacobian, jacobian_reentrant)
            grads_input = gradient(output, grad_output)
            for (jacobian_x, d_x) in zip(jacobian_c, grads_input)
                jacobian_x[:, i] = reshape(d_x, :)
                # TODO: check gradient size
            end
        end
    end
    jacobian, jacobian_reentrant
end

function _analytical_jacobian!(jacobian::Tuple, jacobian_reentrant::Tuple, grad_output::Number, output::Value)
    grad_output = 1
    for jacobian_c in (jacobian, jacobian_reentrant)
        grads_input = gradient(output, grad_output)
        for (jacobian_x, d_x) in zip(jacobian_c, grads_input)
            jacobian_x[:, 1] = reshape(d_x, :)
            # TODO: check gradient size
        end
    end

    jacobian, jacobian_reentrant
end

_differentiable_outputs(x) = Tuple(o for o in _as_tuple(x) if o isa Value)
_as_tuple(x) = (x, )
_as_tuple(x::Tuple) = x

function gradcheck(f, inputs...; eps=1e-6, atol=1e-5, rtol=1e-3, raise_exception=true)
    any_input_requiring_grad = false
    for inp in inputs
        if inp isa Variable
            any_input_requiring_grad = true
        end
    end

    any_input_requiring_grad || error("gradcheck expects at least one input to require gradient")

    output = _differentiable_outputs(f(inputs...))

    # replace this with assert in the future
    function fail_test(msgs...)
        raise_exception && error(msgs...)
        false
    end

    for (i, o) in enumerate(output)
        fn(inputs...) = _as_tuple(f(inputs...))[i]

        analytical, reentrant, correct_grad_sizes = get_analytical_jacobian(inputs, o)
        numerical = get_numerical_jacobian(fn, inputs...; eps=eps)

        correct_grad_sizes || fail_test("Analytical gradient has incorrect size")

        for (j, (a, n)) in enumerate(zip(analytical, numerical))
            if length(a) != 0 || length(n) != 0
                isapprox(a, n; rtol=rtol, atol=atol) ||
                    fail_test("Jacobian mismatch for output $(i) with respect to input $j, \n",
                              "numerical: $n \nanalytical: $a\n")
            end
        end

        # TODO: reentrant is not actually checked
        reentrant || fail_test("Backward is not reentrant, i.e., running backward with same ",
                               "input and grad_output multiple times gives different values, ",
                               "although analytical gradient matches numerical gradient")
    end

    true
end

end # module TestUtils


text:
module AtmosphericProfilesLibrary

import Dierckx

# Large data-based profiles
include("profiles/Soares.jl")
include("profiles/Nieuwstadt.jl")
include("profiles/Bomex.jl")
include("profiles/LifeCycleTan2018.jl")
include("profiles/Rico.jl")
include("profiles/TRMM_LBA.jl")
include("profiles/ARM_SGP.jl")
include("profiles/GATE_III.jl")
include("profiles/Dycoms_RF01.jl")
include("profiles/Dycoms_RF02.jl")
include("profiles/GABLS.jl")
include("profiles/SP.jl")
include("profiles/DryBubble.jl")

end # module

ext: jl


text:
using Distributions, StatsBase

function InteractionNetworkModels.imcmc_gibbs_update!(
    S_curr::InteractionSequence,
    S_prop::InteractionSequence,
    i::Int,
    model::SIM,
    mcmc::SimMcmcInsertDeleteGibbs
    )

    @inbounds n = length(S_curr[i])
    δ = rand(1:mcmc.ν)
    a,b = (lb(n, δ, model), ub(n, δ))
    # @show n, δ, a, b
    d = rand(a:b)
    m = n + δ - 2*d
    # @show m
    # Set-up views
    ind_del = view(mcmc.ind_del, 1:d)
    ind_add = view(mcmc.ind_add, 1:(δ-d))
    vals = view(mcmc.vals, 1:(δ-d))

    # Sample indexing info and new entries (all in-place)
    StatsBase.seqsample_a!(1:n, ind_del)
    StatsBase.seqsample_a!(1:m, ind_add)
    sample!(model.V, vals)

    delete_insert!(S_prop[i], δ, d, ind_del, ind_add, vals)

    log_ratio = log(b - a + 1) - log(ub(m, δ) - lb(m, δ, model) +1) + (δ - 2*d) * log(length(model.V))

    # Multinomial term
    log_multinom_ratio = log_multinomial_ratio(S_curr, S_prop)

    # @show curr_dist, prop_dist
    @inbounds log_α = (
        -model.γ * (
            model.dist(model.mode, S_prop)-model.dist(model.mode, S_curr)
            ) + log_ratio + log_multinom_ratio
        )
    if log(rand()) < log_α
        # println("accepted")
        copy!(S_curr[i], S_prop[i])
        # @inbounds S_curr[i] = copy(S_prop[i])
        # println("$(i) value proposal $(tmp_prop) was accepted")
        return 1
    else
        copy!(S_prop[i], S_curr[i])
        # @inbounds S_prop[i] = copy(S_curr[i])
        # println("$(i) value proposal $(tmp_prop) was rejected")
        return 0
    end
    # @show S_curr
end

function InteractionNetworkModels.imcmc_gibbs_scan!(
    S_curr::InteractionSequence,
    S_prop::InteractionSequence,
    model::SIM,
    mcmc::SimMcmcInsertDeleteGibbs
    )
    count = 0
    N = length(S_curr)
    for i = 1:N
        count += imcmc_gibbs_update!(S_curr, S_prop, i, model, mcmc)
    end
    return count
end

function InteractionNetworkModels.draw_sample!(
    sample_out::InteractionSequenceSample{Int},
    mcmc::SimMcmcInsertDeleteGibbs,
    model::SIM;
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    init::InteractionSequence{Int}=model.mode
    )

    # Define aliases for pointers to the storage of current vals and proposals
    curr_pointers = mcmc.curr_pointers
    prop_pointers = mcmc.prop_pointers

    S_curr = InteractionSequence{Int}()
    S_prop = InteractionSequence{Int}()
    for i in 1:length(init)
        migrate!(S_curr, curr_pointers, i, 1)
        migrate!(S_prop, prop_pointers, i, 1)
        copy!(S_curr[i], init[i])
        copy!(S_prop[i], init[i])
    end
    # @show S_curr, S_prop, init
    # S_curr = copy(init)
    # S_prop = copy(init)

    ind = 0
    count = 0
    sample_count = 1
    i = 0
    acc_count = 0
    gibbs_scan_count = 0
    gibbs_tot_count = 0
    gibbs_acc_count = 0
    # is_insert = false

    # Bounds for uniform sampling
    lb(x::Vector{Path{Int}}) = max(1, length(x) - 1 )
    ub(x::Vector{Path{Int}}) = min(model.K_outer, length(x) + 1)

    # prob gibbs (effectively with ϵ=1)
    prob_gibbs = 1/(2+1) + mcmc.β

    while sample_count ≤ length(sample_out)
        i += 1
        # Store value
        if (i > burn_in) & (((i-1) % lag)==0)
            sample_out[sample_count] = deepcopy(S_curr)
            sample_count += 1
        end
        # Gibbs scan
        if rand() < prob_gibbs
            # println("Gibbs")
            gibbs_scan_count += 1
            gibbs_tot_count += length(S_curr)
            gibbs_acc_count += imcmc_gibbs_scan!(S_curr, S_prop, model, mcmc) # This enacts the scan, changing curr, and outputs number of accepted moves.
        # Else do insert or delete
        else
            # println("Transdim")
            count += 1
            # If only one interaction we only insert
            if length(S_curr) == 1
                is_insert = true
                ind = rand(1:(length(S_curr)+1))
                log_ratio = imcmc_insert_prop_sample!(
                    S_curr, S_prop,
                    curr_pointers, prop_pointers,
                    ind,
                    mcmc.path_dist
                )
                log_ratio += 0.5 # Adjust log ratio

            # If max mnumber of interactions we only delete
            elseif length(S_curr) == model.K_outer
                is_insert = false
                ind = rand(1:length(S_curr))
                log_ratio = imcmc_delete_prop_sample!(
                    S_curr, S_prop,
                    curr_pointers, prop_pointers,
                    ind,
                    mcmc.path_dist
                )
                log_ratio += 0.5 # Adjust log ratio

            elseif rand() < 0.5  # Coin flip for insert
                is_insert = true
                ind = rand(1:(length(S_curr)+1))
                log_ratio = imcmc_insert_prop_sample!(
                    S_curr, S_prop,
                    curr_pointers, prop_pointers,
                    ind,
                    mcmc.path_dist
                )
            else # Else delete
                is_insert = false
                ind = rand(1:length(S_curr))
                log_ratio = imcmc_delete_prop_sample!(
                    S_curr, S_prop,
                    curr_pointers, prop_pointers,
                    ind,
                    mcmc.path_dist
                )
            end
            # Multinomal term
            log_multinom_ratio = log_multinomial_ratio(S_curr, S_prop)

            # println(S_curr)
            log_α = - model.γ * (
                model.dist(model.mode, S_prop) - model.dist(model.mode, S_curr)
            ) + log_ratio + log_multinom_ratio

            if rand() < exp(log_α)
                if is_insert
                    migrate!(S_curr, curr_pointers, ind, 1)
                    copy!(S_curr[ind], S_prop[ind])
                    # insert!(S_curr, ind, copy(S_prop[ind]))
                else
                    migrate!(curr_pointers, S_curr, 1, ind)
                    # deleteat!(S_curr, ind)
                end
                acc_count += 1
            else
                if is_insert
                    migrate!(prop_pointers, S_prop, 1, ind)
                    # deleteat!(S_prop, ind)
                else
                    migrate!(S_prop, prop_pointers, ind, 1)
                    # insert!(S_prop, ind, copy(S_curr[ind]))
                end
            end
            # sample_out[i] = copy(S_curr)
        end
    end
    # Send storage back
    # @show S_curr, S_prop
    for i in 1:length(S_curr)
        migrate!(curr_pointers, S_curr, 1, 1)
        migrate!(prop_pointers, S_prop, 1, 1)
    end
    return count, acc_count, gibbs_tot_count, gibbs_scan_count, gibbs_acc_count

end

function InteractionNetworkModels.draw_sample(
    mcmc::SimMcmcInsertDeleteGibbs,
    model::SIM;
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    init::InteractionSequence{Int}=model.mode
    )

    sample_out = InteractionSequenceSample{Int}(undef, desired_samples)
    (
        count,
        acc_count,
        gibbs_tot_count,
        gibbs_scan_count,
        gibbs_acc_count
        ) = draw_sample!(sample_out, mcmc, model, burn_in=burn_in, lag=lag, init=init)

    p_measures = Dict(
            "Proportion Gibbs moves" => gibbs_scan_count/(count + gibbs_scan_count),
            "Trans-dimensional move acceptance probability" => acc_count/count,
            "Gibbs move acceptance probability" => gibbs_acc_count/gibbs_tot_count
    )
    output = SimMcmcOutput(
        model,
        sample_out,
        p_measures
        )

    return output

end


text:
function write_array_set_to_disk(data_set::Set{Array{Float64,2}},save_directory::String,file_prefix::String)

  counter = 1
  while (length(data_set)>0)

    # Grab a data array from the set -
    data_array = pop!(data_set)

    # write the array -
    full_path_string = save_directory*"/"*file_prefix*"_T$(counter).dat"
    writedlm(full_path_string,data_array)

    # update -
    counter = counter + 1
  end
end

function calculate_mean_array_from_set(data_set::Set{Array{Float64,2}},number_of_rows::Int,number_of_cols::Int)

    # Initialize the mean -
    mean_array = zeros(number_of_rows,number_of_cols)

    counter = 0
    while (length(data_set)>0)

      # Grab a data array from the set -
      data_array = pop!(data_set)

      # iterate -
      for row_index = 1:number_of_rows
        for col_index = 1:number_of_cols

          # calculate the mean value -
          mean_array[row_index,col_index] = (data_array[row_index,col_index]+counter*mean_array[row_index,col_index])/(counter+1)

        end
      end

      # update the counter -
      counter = counter + 1
    end

    # return -
    return mean_array
end

function rand_normal(mean, stdev)

    if stdev <= 0.0
        error("standard deviation must be positive")
    end

    u1 = rand()
    u2 = rand()
    r = sqrt(-2.0*log(u1))

    theta = 2.0*pi*u2
    return mean + stdev*r*sin(theta)
end


text:
# Base.promote_rule(PeriodicFunctionMatrix, PeriodicSymbolicMatrix) = PeriodicFunctionMatrix{:c,T}
# Base.promote_rule(PeriodicFunctionMatrix, HarmonicArray) = PeriodicFunctionMatrix
function promote_period(PM1,args...; ndigits = 4)
   nlim = 2^ndigits
   if typeof(PM1) <: AbstractVecOrMat
      period = nothing
      isconst = true
   else
      period = PM1.period
      isconst = isconstant(PM1)
   end
   for a in args
       (typeof(a) <: AbstractVecOrMat || isconstant(a)) && continue
       isconst && (isconst = false; period = a.period; continue)
       peri = a.period
       r = rationalize(period/peri)
       num = numerator(r)
       den = denominator(r)
       (num <= nlim && den <= nlim) || error("incommensurate periods")
       period = period*den
   end
   return period
end
function promote_Ts(PM1,args...)
    Ts = PM1.Ts
    isconst = isconstant(PM1)
    for a in args
        isconstant(a) && continue
        isconst && (isconst = false; Ts = a.Ts; continue)
        Ts ≈ a.Ts || error("incompatible sampling times")
    end
    return Ts
end
function promote_nperiod(PM1,args...)
   if typeof(PM1) <: AbstractVecOrMat
      nperiod = nothing
      isconst = true
   else
      nperiod = PM1.nperiod
      isconst = isconstant(PM1)
   end
   for a in args
       (typeof(a) <: AbstractVecOrMat || isconstant(a)) && continue
       isconst && (isconst = false; period = a.nperiod; continue)
       nperiod = gcd(a.nperiod,nperiod)
       nperiod == 1 && break
   end
   return nperiod
end
function promote_cpmtype(PM1,args...)
   # promote continuous-time periodic matrix type
   if typeof(PM1) <: AbstractVecOrMat
      pmtype = nothing
      isconst = true
   else
      pmtype = typeof(PM1).name.wrapper
      isconst = isconstant(PM1)
      iscontinuous(PM1) || error("only continuous-time or constant matrices supported")
   end
   for a in args
       (typeof(a) <: AbstractVecOrMat || isconstant(a)) && continue
       iscontinuous(a) || error("only continuous-time or constant matrices supported")
       isconst && (isconst = false; pmtype = typeof(a).name.wrapper; continue)
       typeof(a).name.wrapper == pmtype && continue
       pmtype = PeriodicFunctionMatrix
       break
   end
   return isnothing(pmtype) ? PeriodicFunctionMatrix : pmtype
end
function promote_dpmtype(PM1,args...)
   # promote discrete-time periodic matrix type
   if typeof(PM1) <: AbstractVecOrMat
      pmtype = nothing
      isconst = true
   else
      pmtype = typeof(PM1).name.wrapper
      isconst = isconstant(PM1)
      iscontinuous(PM1) && error("only discrete-time or constant matrices supported")
   end
   for a in args
       (typeof(a) <: AbstractVecOrMat || isconstant(a)) && continue
       iscontinuous(a) && error("only discrete-time or constant matrices supported")
       isconst && (isconst = false; pmtype = typeof(a).name.wrapper; continue)
       typeof(a).name.wrapper == pmtype && continue
       pmtype = PeriodicMatrix
       break
   end
   return isnothing(pmtype) ? PeriodicMatrix : pmtype
end

struct PeriodicStateSpace{PM} <: AbstractPeriodicStateSpace
    A::PM
    B::PM
    C::PM
    D::PM
    period::Float64
end
"""
    PeriodicStateSpace(A::PM, B::PM, C::PM, D::PM) -> psys::PeriodicStateSpace{PM}

Construct a `PeriodicStateSpace` object from a quadruple of periodic matrix objects.

The periodic matrix objects `A`, `B`, `C`, `D` specifies the periodic matrices
`A(t)`, `B(t)`, `C(t)` and `D(t)` of a linear periodic time-varying state space model
in the continuous-time form

     dx(t)/dt = A(t)x(t) + B(t)u(t) ,
     y(t)     = C(t)x(t) + D(t)u(t) ,

or in the discrete-time form

     x(t+1)  = A(t)x(t) + B(t)u(t) ,
     y(t)    = C(t)x(t) + D(t)u(t) ,

where `x(t)`, `u(t)` and `y(t)` are the system state vector,
system input vector and system output vector, respectively,
and `t` is the continuous or discrete time variable.
The system matrices satisfy `A(t) = A(t+T₁)`, `B(t) = B(t+T₂)`, `C(t) = C(t+T₃)`, `D(t) = D(t+T₄)`,
i.e., are periodic with periods `T₁`, `T₂`, `T₃` and `T₄`, respectively.
The different periods must be commensurate (i.e., their ratios must be rational numbers with
numerators and denominators up to at most 4 decimal digits).
All periodic matrix objects must have the same type `PM`, where
`PM` stays for one of the supported periodic matrix types, i.e.,
[`PeriodicMatrix`](@ref), [`PeriodicArray`](@ref), [`PeriodicFunctionMatrix`](@ref), [`PeriodicSymbolicMatrix`](@ref),
[`HarmonicArray`](@ref), [`FourierFunctionMatrix`](@ref) or [`PeriodicTimeSeriesMatrix`](@ref).
"""
function PeriodicStateSpace(A::PFM1, B::PFM2, C::PFM3, D::PFM4) where {PFM1 <: PeriodicFunctionMatrix, PFM2 <: PeriodicFunctionMatrix, PFM3 <: PeriodicFunctionMatrix, PFM4 <: PeriodicFunctionMatrix}
    period = ps_validation(A, B, C, D)
    T = promote_type(eltype(A),eltype(B),eltype(C),eltype(D))
    PeriodicStateSpace{PeriodicFunctionMatrix{:c,T}}((period == A.period && T == eltype(A)) ? A : PeriodicFunctionMatrix{:c,T}(A,period),
                                                     (period == B.period && T == eltype(B)) ? B : PeriodicFunctionMatrix{:c,T}(B,period),
                                                     (period == C.period && T == eltype(C)) ? C : PeriodicFunctionMatrix{:c,T}(C,period),
                                                     (period == D.period && T == eltype(D)) ? D : PeriodicFunctionMatrix{:c,T}(D,period),
                                                     Float64(period))
end
function PeriodicStateSpace(A::FFM1, B::FFM2, C::FFM3, D::FFM4) where {FFM1 <: FourierFunctionMatrix, FFM2 <: FourierFunctionMatrix, FFM3 <: FourierFunctionMatrix, FFM4 <: FourierFunctionMatrix}
    period = ps_validation(A, B, C, D)
    T = promote_type(eltype(A),eltype(B),eltype(C),eltype(D))
    PeriodicStateSpace{FourierFunctionMatrix{:c,T}}((period == A.period && T == eltype(A)) ? A : FourierFunctionMatrix{:c,T}(A,period),
                                                     (period == B.period && T == eltype(B)) ? B : FourierFunctionMatrix{:c,T}(B,period),
                                                     (period == C.period && T == eltype(C)) ? C : FourierFunctionMatrix{:c,T}(C,period),
                                                     (period == D.period && T == eltype(D)) ? D : FourierFunctionMatrix{:c,T}(D,period),
                                                     Float64(period))
end

function PeriodicStateSpace(A::PSM, B::PSM, C::PSM, D::PSM) where {T,PSM <: PeriodicSymbolicMatrix{:c,T}}
    period = ps_validation(A, B, C, D)
    PeriodicStateSpace{PSM}(period == A.period ? A : PeriodicSymbolicMatrix{:c,T}(A,period),
                            period == B.period ? B : PeriodicSymbolicMatrix{:c,T}(B,period),
                            period == C.period ? C : PeriodicSymbolicMatrix{:c,T}(C,period),
                            period == D.period ? D : PeriodicSymbolicMatrix{:c,T}(D,period),
                            Float64(period))
end
function PeriodicStateSpace(A::PHR1, B::PHR2, C::PHR3, D::PHR4) where {PHR1 <: HarmonicArray, PHR2 <: HarmonicArray, PHR3 <: HarmonicArray, PHR4 <: HarmonicArray}
    period = ps_validation(A, B, C, D)
    T = promote_type(eltype(A),eltype(B),eltype(C),eltype(D))
    PeriodicStateSpace{HarmonicArray{:c,T}}((period == A.period && T == eltype(A)) ? A : HarmonicArray{:c,T}(A,period),
                                            (period == B.period && T == eltype(B)) ? B : HarmonicArray{:c,T}(B,period),
                                            (period == C.period && T == eltype(C)) ? C : HarmonicArray{:c,T}(C,period),
                                            (period == D.period && T == eltype(D)) ? D : HarmonicArray{:c,T}(D,period),
                                            Float64(period))
end
function PeriodicStateSpace(A::PTSM1, B::PTSM2, C::PTSM3, D::PTSM4) where {PTSM1 <: PeriodicTimeSeriesMatrix, PTSM2 <: PeriodicTimeSeriesMatrix, PTSM3 <: PeriodicTimeSeriesMatrix, PTSM4 <: PeriodicTimeSeriesMatrix}
    period = ps_validation(A, B, C, D)
    T = promote_type(eltype(A),eltype(B),eltype(C),eltype(D))
    PeriodicStateSpace{PeriodicTimeSeriesMatrix{:c,T}}((period == A.period && T == eltype(A)) ? A : PeriodicTimeSeriesMatrix{:c,T}(A,period),
                                            (period == B.period && T == eltype(B)) ? B : PeriodicTimeSeriesMatrix{:c,T}(B,period),
                                            (period == C.period && T == eltype(C)) ? C : PeriodicTimeSeriesMatrix{:c,T}(C,period),
                                            (period == D.period && T == eltype(D)) ? D : PeriodicTimeSeriesMatrix{:c,T}(D,period),
                                            Float64(period))
end

function ps_validation(A::PM1, B::PM2, C::PM3, D::PM4) where {T1, T2, T3, T4, PM1 <: AbstractPeriodicArray{:c,T1}, PM2 <: AbstractPeriodicArray{:c,T2}, PM3 <: AbstractPeriodicArray{:c,T3}, PM4 <: AbstractPeriodicArray{:c,T4}}
    nx = size(A,1)
    nx == size(A,2) || error("matrix A(t) is not square")
    (ny, nu) = size(D)
    size(B,1) == nx ||  error("B(t) must have the same row size as A(t)")
    size(C,2) == nx ||  error("C(t) must have the same column size as A(t)")
    nu == size(B,2) ||  error("D(t) must have the same column size as B(t)")
    ny == size(C,1) ||  error("D(t) must have the same row size as C(t)")
    # validate sampling time
    period = promote_period(A, B, C, D)
    return period
end
function PeriodicStateSpace(A::PM1, B::PM2, C::PM3, D::PM4) where {PM1 <: PeriodicMatrix, PM2 <: PeriodicMatrix, PM3 <: PeriodicMatrix, PM4 <: PeriodicMatrix}
    period = ps_validation(A, B, C, D)
    promote_Ts(A,B,C,D)
    T = promote_type(eltype(A),eltype(B),eltype(C),eltype(D))
    PeriodicStateSpace{PeriodicMatrix{:d,T}}((period == A.period && T == eltype(A)) ? A : PeriodicMatrix{:d,T}(A,period),
                                             (period == B.period && T == eltype(B)) ? B : PeriodicMatrix{:d,T}(B,period),
                                             (period == C.period && T == eltype(C)) ? C : PeriodicMatrix{:d,T}(C,period),
                                             (period == D.period && T == eltype(D)) ? D : PeriodicMatrix{:d,T}(D,period),
                                             Float64(period))
end

function ps_validation(A::PM1, B::PM2, C::PM3, D::PM4) where {T1, T2, T3, T4, PM1 <: PeriodicMatrix{:d,T1}, PM2 <: PeriodicMatrix{:d,T2}, PM3 <: PeriodicMatrix{:d,T3}, PM4 <: PeriodicMatrix{:d,T4}}
    period = promote_period(A, B, C, D)
    max(length(A),length(B),length(C),length(D)) == 0 && (return period)

    # Validate dimensions
    ny = size(D,1)
    nu = size(D,2)
    any(ny .!= ny[1]) && error("all matrices D[i] must have the same row size")
    any(nu .!= nu[1]) && error("all matrices D[i] must have the same column size")

    ndx, nx = size(A)
    if all(ndx .== ndx[1]) && all(nx .== nx[1])
       # constant dimensions of A, B, C, D
       ndx[1] == nx[1] || error("A[i] must be square matrices")
       any(size(B, 1) .!= ndx[1]) && error("all matrices A[i] and B[i] must have the same row size")
       any(size(C, 2) .!= nx[1]) &&  error("all matrices A[i] and C[i] must have the same column size")
       any(size(C, 1) .!= ny[1]) && error("all matrices C[i] and D[i] must have the same row size")
       any(size(B, 2) .!= nu[1]) && error("all matrices B[i] and D[i] must have the same column size")
    else
       N = A.dperiod
       (N != B.dperiod || N != C.dperiod) && error("the number of component matrices of A, B and C must be the same ")
       for i = 1:N
           i == N && ndx[i] != nx[1] && error("the number of columns of A[i+1] must be equal to the number of rows of A[i]")
           i != N && ndx[i] != nx[i+1] && error("the number of columns of A[i+1] must be equal to the number of rows of A[i]")
       end
       any(size(B, 1) .!= ndx) && error("all matrices A[i] and B[i] must have the same row size")
       any(size(C, 2) .!= nx) &&  error("all matrices A[i] and C[i] must have the same column size")
       any(size(C, 1) .!= ny[1]) && error("all matrices C[i] and D[i] must have the same row size")
       any(size(B, 2) .!= nu[1]) && error("all matrices B[i] and D[i] must have the same column size")
    end
    return period
end
function PeriodicStateSpace(A::PA1, B::PA2, C::PA3, D::PA4) where {PA1 <: PeriodicArray, PA2 <: PeriodicArray, PA3 <: PeriodicArray, PA4 <: PeriodicArray}
    period = ps_validation(A, B, C, D)
    promote_Ts(A,B,C,D)
    T = promote_type(eltype(A),eltype(B),eltype(C),eltype(D))
    PeriodicStateSpace{PeriodicArray{:d,T}}((period == A.period && T == eltype(A)) ? A : PeriodicArray{:d,T}(A,period),
                                            (period == B.period && T == eltype(B)) ? B : PeriodicArray{:d,T}(B,period),
                                            (period == C.period && T == eltype(C)) ? C : PeriodicArray{:d,T}(C,period),
                                            (period == D.period && T == eltype(D)) ? D : PeriodicArray{:d,T}(D,period),
                                            Float64(period))
end
function ps_validation(A::PA1, B::PA2, C::PA3, D::PA4) where {T1, T2, T3, T4, PA1 <: PeriodicArray{:d,T1}, PA2 <: PeriodicArray{:d,T2}, PA3 <: PeriodicArray{:d,T3}, PA4 <: PeriodicArray{:d,T4}}
    period = promote_period(A, B, C, D)
    max(length(A),length(B),length(C),length(D)) == 0 && (return period)

    # Validate dimensions
    ny = size(D,1)
    nu = size(D,2)
    nx = size(A,1)
    nx == size(A,2) || error("matrix A(t) is not square")

    any(size(B, 1) .!= nx) && error("A and B must have the same row size")
    any(size(C, 2) .!= nx) &&  error("A and C must have the same column size")
    any(size(C, 1) .!= ny) && error("C and D must have the same row size")
    any(size(B, 2) .!= nu) && error("B and D must have the same column size")

    return period
end

#  conversions
function Base.convert(::Type{PeriodicStateSpace{PM}}, psys::PeriodicStateSpace) where {PM <: AbstractPeriodicArray}
   convert(PM,psys.A)
   convert(PM,psys.B)
   convert(PM,psys.C)
   convert(PM,psys.D)
   PeriodicStateSpace(convert(PM,psys.A), convert(PM,psys.B), convert(PM,psys.C), convert(PM,psys.D))
end
# properties
Base.size(sys::PeriodicStateSpace) = maximum.(size(sys.D))
function Base.size(sys::PeriodicStateSpace, d::Integer)
    return d <= 2 ? size(sys)[d] : 1
end
Base.length(sys::PeriodicStateSpace) = 1
Base.eltype(sys::PeriodicStateSpace) = eltype(sys.A)

function Base.getproperty(sys::PeriodicStateSpace, d::Symbol)
    if d === :nx
        return size(getfield(sys, :A), 2)
    elseif d === :ny
        return size(getfield(sys, :C), 1)
    elseif d === :nu
        return size(getfield(sys, :B), 2)
    else
        getfield(sys, d)
    end
end
Base.propertynames(sys::PeriodicStateSpace) =
    (:nx, :ny, :nu, fieldnames(typeof(sys))...)

# display sys
Base.print(io::IO, sys::PeriodicStateSpace) = show(io, sys)
Base.show(io::IO, sys::PeriodicStateSpace{PM}) where
    {PM <: Union{PeriodicMatrix,PeriodicArray,PeriodicTimeSeriesMatrix,PeriodicSymbolicMatrix,PeriodicFunctionMatrix,HarmonicArray}} =
    show(io, MIME("text/plain"), sys)

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, sys::PeriodicStateSpace{<:FourierFunctionMatrix})
    summary(io, sys); println(io)
    n = size(sys.A,1)
    p, m = size(sys.D)
    T = eltype(sys)
    if n > 0
       nperiod = sys.A.nperiod
       println(io, "\nState matrix A::$T($n×$n): subperiod: $(sys.A.period/nperiod)    #subperiods: $nperiod ")
       isconstant(sys.A) ? show(io, mime, sys.A.M(0)) : show(io, mime, sys.A.M)
       if m > 0
          nperiod = sys.B.nperiod
          println(io, "\n\nInput matrix B::$T($n×$m): $(sys.B.period/nperiod)    #subperiods: $nperiod ")
          isconstant(sys.B) ? show(io, mime, sys.B.M(0)) : show(io, mime, sys.B.M)
       else
          println(io, "\n\nEmpty input matrix B.")
       end

       if p > 0
          nperiod = sys.C.nperiod
          println(io, "\n\nOutput matrix C::$T($p×$n): $(sys.C.period/nperiod)    #subperiods: $nperiod ")
          isconstant(sys.C) ? show(io, mime, sys.C.M(0)) : show(io, mime, sys.C.M)
       else
          println(io, "\n\nEmpty output matrix C.")
       end
       if m > 0 && p > 0
          nperiod = sys.D.nperiod
          println(io, "\n\nFeedthrough matrix D::$T($p×$m): $(sys.D.period/nperiod)    #subperiods: $nperiod ")
          isconstant(sys.D) ? show(io, mime, sys.D.M(0)) : show(io, mime, sys.D.M)
       else
          println(io, "\n\nEmpty feedthrough matrix D.")
       end
       println(io, "\n\nContinuous-time periodic state-space model.")
    elseif m > 0 && p > 0
       nperiod = sys.D.nperiod
       println(io, "\nFeedthrough matrix D::$T($p×$m): $(sys.D.period/nperiod)    #subperiods: $nperiod ")
       isconstant(sys.D) ? show(io, mime, sys.D.M(0)) : show(io, mime, sys.D.M)
       println(io, "\n\nTime-varying gain.")
    else
       println(io, "\nEmpty state-space model.")
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, sys::PeriodicStateSpace{<:HarmonicArray})
    summary(io, sys); println(io)
    n = sys.nx
    p, m = size(sys)
    period = sys.period
    if n > 0
       nharmonics, nperiod = size(sys.A.values,3)-1, sys.A.nperiod
       subperiod = period/nperiod
       println(io, "\nState matrix A: subperiod: $subperiod    #subperiods: $nperiod   #harmonics: $nharmonics")
       ts = 2*pi*nperiod/period
       if nharmonics >= 0
          println("Constant term")
          display(mime, real(sys.A.values[:,:,1]))
       end
       for i in [1:min(4,nharmonics); nharmonics-2:nharmonics]
         i <= 0 && break
         println("cos($(i*ts)t) factor")
          display(mime, real(sys.A.values[:,:,i+1]))
          println("sin($(i*ts)t) factor")
          display(mime, imag(sys.A.values[:,:,i+1]))
          i == 4 && nharmonics > 4 && println("   ⋮")
       end
       if m > 0
          nharmonics, nperiod = size(sys.B.values,3)-1, sys.B.nperiod
          subperiod = period/nperiod
          println(io, "\nInput matrix B: subperiod: $subperiod    #subperiods: $nperiod   #harmonics: $nharmonics")
          ts = 2*pi*nperiod/period
          if nharmonics >= 0
             println("Constant term")
             display(mime, real(sys.B.values[:,:,1]))
          end
          for i in [1:min(4,nharmonics); nharmonics-2:nharmonics]
             i <= 0 && break
             println("cos($(i*ts)t) factor")
             display(mime, real(sys.B.values[:,:,i+1]))
             println("sin($(i*ts)t) factor")
             display(mime, imag(sys.B.values[:,:,i+1]))
             i == 4 && nharmonics > 4 && println("   ⋮")
          end
       else
          println(io, "\n\nEmpty input matrix B.")
       end
       if p > 0
          nharmonics, nperiod = size(sys.C.values,3)-1, sys.C.nperiod
          subperiod = period/nperiod
          println(io, "\nOutput matrix C: subperiod: $subperiod    #subperiods: $nperiod   #harmonics: $nharmonics")
          ts = 2*pi*nperiod/period
          if nharmonics >= 0
            println("Constant term")
            display(mime, real(sys.C.values[:,:,1]))
         end
         for i in [1:min(4,nharmonics); nharmonics-2:nharmonics]
            i <= 0 && break
            println("cos($(i*ts)t) factor")
            display(mime, real(sys.C.values[:,:,i+1]))
            println("sin($(i*ts)t) factor")
            display(mime, imag(sys.C.values[:,:,i+1]))
            i == 4 && nharmonics > 4 && println("   ⋮")
         end
      else
          println(io, "\n\nEmpty output matrix C.")
       end
       if m > 0 && p > 0
          nharmonics, nperiod = size(sys.D.values,3)-1, sys.D.nperiod
          subperiod = period/nperiod
          println(io, "\nFeedthrough matrix D: subperiod: $subperiod    #subperiods: $nperiod   #harmonics: $nharmonics")
          ts = 2*pi*nperiod/period
          if nharmonics >= 0
            println("Constant term")
            display(mime, real(sys.D.values[:,:,1]))
          end
          for i in [1:min(4,nharmonics); nharmonics-2:nharmonics]
             i <= 0 && break
             println("cos($(i*ts)t) factor")
             display(mime, real(sys.D.values[:,:,i+1]))
             println("sin($(i*ts)t) factor")
             display(mime, imag(sys.C.values[:,:,i+1]))
             i == 4 && nharmonics > 4 && println("   ⋮")
          end
       else
          println(io, "\n\nEmpty feedthrough matrix D.")
       end
       println(io, "\n\nPeriod:      $(sys.period) second(s).")
       println(io, "Periodic continuous-time state-space model.")
    elseif m > 0 && p > 0
       nharmonics, nperiod = size(sys.D.values,3)-1, sys.D.nperiod
       subperiod = period/nperiod
       println(io, "\nFeedthrough matrix D: subperiod: $subperiod    #subperiods: $nperiod   #harmonics: $nharmonics")
       ts = 2*pi*nperiod/period
       if nharmonics >= 0
          println("Constant term")
          display(mime, real(sys.D.values[:,:,1]))
       end
       for i in [1:min(4,nharmonics); nharmonics-2:nharmonics]
          i <= 0 && break
          println("cos($(i*ts)t) factor")
          display(mime, real(sys.D.values[:,:,i+1]))
          println("sin($(i*ts)t) factor")
          display(mime, imag(sys.D.values[:,:,i+1]))
          i == 4 && nharmonics > 4 && println("⋮")
       end
       println(io, "\n\nPeriod:      $(sys.period) second(s).")
       println(io, "Time-varying gains.")
    else
       println(io, "\nEmpty Periodic continuous-time state-space model.")
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, sys::PeriodicStateSpace{<:PeriodicFunctionMatrix})
    summary(io, sys); println(io)
    n = size(sys.A,1)
    p, m = size(sys.D)
    T = eltype(sys)
    if n > 0
       nperiod = sys.A.nperiod
       println(io, "\nState matrix A::$T($n×$n): subperiod: $(sys.A.period/nperiod)    #subperiods: $nperiod ")
       isconstant(sys.A) ? show(io, mime, sys.A.f(0)) : show(io, mime, sys.A.f)
       if m > 0
          nperiod = sys.B.nperiod
          println(io, "\n\nInput matrix B::$T($n×$m): $(sys.B.period/nperiod)    #subperiods: $nperiod ")
          isconstant(sys.B) ? show(io, mime, sys.B.f(0)) : show(io, mime, sys.B.f)
       else
          println(io, "\n\nEmpty input matrix B.")
       end

       if p > 0
          nperiod = sys.C.nperiod
          println(io, "\n\nOutput matrix C::$T($p×$n): $(sys.C.period/nperiod)    #subperiods: $nperiod ")
          isconstant(sys.C) ? show(io, mime, sys.C.f(0)) : show(io, mime, sys.C.f)
       else
          println(io, "\n\nEmpty output matrix C.")
       end
       if m > 0 && p > 0
          nperiod = sys.D.nperiod
          println(io, "\n\nFeedthrough matrix D::$T($p×$m): $(sys.D.period/nperiod)    #subperiods: $nperiod ")
          isconstant(sys.D) ? show(io, mime, sys.D.f(0)) : show(io, mime, sys.D.f)
       else
          println(io, "\n\nEmpty feedthrough matrix D.")
       end
       println(io, "\n\nPeriod:      $(sys.period) second(s).")
       println(io, "Continuous-time periodic state-space model.")
    elseif m > 0 && p > 0
       nperiod = sys.D.nperiod
       println(io, "\nFeedthrough matrix D::$T($p×$m): $(sys.D.period/nperiod)    #subperiods: $nperiod ")
       isconstant(sys.D) ? show(io, mime, sys.D.f(0)) : show(io, mime, sys.D.f)
       println(io, "Period:      $(sys.period) second(s).")
       println(io, "\n\nTime-varying gain.")
    else
       println(io, "\nEmpty state-space model.")
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, sys::PeriodicStateSpace{<:PeriodicSymbolicMatrix})
    summary(io, sys); println(io)
    n = size(sys.A,1)
    p, m = size(sys.D)
    if n > 0
       nperiod = sys.A.nperiod
       println(io, "\nState matrix A: subperiod: $(sys.A.period/nperiod)    #subperiods: $nperiod ")
       show(io, mime, sys.A.F)
       if m > 0
          nperiod = sys.B.nperiod
          println(io, "\n\nInput matrix B: $(sys.B.period/nperiod)    #subperiods: $nperiod ")
          show(io, mime, sys.B.F)
       else
          println(io, "\n\nEmpty input matrix B.")
       end

       if p > 0
          nperiod = sys.C.nperiod
          println(io, "\n\nOutput matrix C: $(sys.C.period/nperiod)    #subperiods: $nperiod ")
          show(io, mime, sys.C.F)
       else
          println(io, "\n\nEmpty output matrix C.")
       end
       if m > 0 && p > 0
          nperiod = sys.D.nperiod
          println(io, "\n\nFeedthrough matrix D: $(sys.D.period/nperiod)    #subperiods: $nperiod ")
          show(io, mime, sys.D.F)
       else
          println(io, "\n\nEmpty feedthrough matrix D.")
       end
       println(io, "\n\nPeriod:      $(sys.period) second(s).")
       println(io, "Continuous-time periodic state-space model.")
    elseif m > 0 && p > 0
       nperiod = sys.D.nperiod
       println(io, "\nFeedthrough matrix D: $(sys.D.period/nperiod)    #subperiods: $nperiod ")
       show(io, mime, sys.D.F)
       println(io, "\n\nPeriod:      $(sys.period) second(s).")
       println(io, "Time-varying gain.")
    else
       println(io, "\nEmpty state-space model.")
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, sys::PeriodicStateSpace{<:PeriodicTimeSeriesMatrix})
    summary(io, sys); println(io)
    n = sys.nx
    N = length(sys.A)
    p, m = size(sys)
    period = sys.period
    if n > 0
       dperiod, nperiod = length(sys.A), sys.A.nperiod
       subperiod = period/nperiod
       Ts = subperiod/dperiod
       println(io, "\nState matrix A: subperiod: $subperiod    #subperiods: $nperiod ")
       println(io, "time values: t[1:$dperiod] = [$(0*Ts)  $(Ts) ... $((dperiod-1)*Ts)]")
       for i in [1:min(4,dperiod); dperiod-2:dperiod]
          i <= 0 && break
          println("t[$i] = $((i-1)*Ts)")
          display(mime, sys.A.values[i])
          i == 4 && dperiod > 4 && println("   ⋮")
       end
       if m > 0
          dperiod, nperiod = length(sys.B), sys.B.nperiod
          subperiod = period/nperiod
          Ts = subperiod/dperiod
          println(io, "\n\nInput matrix B: subperiod: $subperiod    #subperiods: $nperiod ")
          println(io, "time values: t[1:$dperiod] = [$(0*Ts)  $(Ts) ... $((dperiod-1)*Ts)]")
          for i in [1:min(4,dperiod); dperiod-2:dperiod]
              i <= 0 && break
              println("t[$i] = $((i-1)*Ts)")
              display(mime, sys.B.values[i])
              i == 4 && dperiod > 4 && println("   ⋮")
          end
       else
          println(io, "\n\nEmpty input matrix B.")
       end
       if p > 0
          dperiod, nperiod = length(sys.C), sys.C.nperiod
          subperiod = period/nperiod
          Ts = subperiod/dperiod
          println(io, "\n\nOutput matrix C: subperiod: $subperiod    #subperiods: $nperiod ")
          println(io, "time values: t[1:$dperiod] = [$(0*Ts)  $(Ts) ... $((dperiod-1)*Ts)]")
          for i in [1:min(4,dperiod); dperiod-2:dperiod]
              i <= 0 && break
              println("t[$i] = $((i-1)*Ts)")
              display(mime, sys.C.values[i])
              i == 4 && dperiod > 4 && println("   ⋮")
          end
       else
          println(io, "\n\nEmpty output matrix C.")
       end
       if m > 0 && p > 0
          dperiod, nperiod = length(sys.D), sys.D.nperiod
          subperiod = period/nperiod
          Ts = subperiod/dperiod
          println(io, "\n\nFeedthrough matrix D: subperiod: $subperiod    #subperiods: $nperiod ")
          println(io, "time values: t[1:$dperiod] = [$(0*Ts)  $(Ts) ... $((dperiod-1)*Ts)]")
          for i in [1:min(4,dperiod); dperiod-2:dperiod]
              i <= 0 && break
              println("t[$i] = $((i-1)*Ts)")
              display(mime, sys.D.values[i])
              i == 4 && dperiod > 4 && println("   ⋮")
          end
       else
           println(io, "\n\nEmpty feedthrough matrix D.")
       end
       println(io, "\n\nPeriod:      $(sys.period) second(s).")
       println(io, "Periodic continuous-time state-space model.")
    elseif m > 0 && p > 0
       dperiod, nperiod = length(sys.D), sys.D.nperiod
       subperiod = period/nperiod
       Ts = subperiod/dperiod
       println(io, "\n\nFeedthrough matrix D: subperiod: $subperiod    #subperiods: $nperiod ")
       println(io, "time values: t[1:$dperiod] = [$(0*Ts)  $(Ts) ... $((dperiod-1)*Ts)]")
       for i in [1:min(4,dperiod); dperiod-2:dperiod]
           i <= 0 && break
           println("t[$i] = $((i-1)*Ts)")
           display(mime, sys.D.values[i])
           i == 4 && dperiod > 4 && println("   ⋮")
       end
       println(io, "\n\nPeriod:      $(sys.period) second(s).")
       println(io, "Time-varying gains.")
    else
       println(io, "\nEmpty Periodic continuous-time state-space model.")
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, sys::PeriodicStateSpace{<:PeriodicMatrix})
    summary(io, sys); println(io)
    n = sys.nx
    N = sys.A.dperiod
    p, m = size(sys)
    period = sys.period
    Ts = sys.A.Ts
    if any(n .> 0)
       dperiod, nperiod = sys.A.dperiod, sys.A.nperiod
       subperiod = period/nperiod
       println(io, "\nState matrix A: subperiod: $subperiod    #subperiods: $nperiod ")
       println(io, "time values: t[1:$dperiod]")
       for i in [1:min(4,dperiod); dperiod-2:dperiod]
           i <= 0 && break
           println("t[$i] = $(i*Ts)")
           display(mime, sys.A.M[i])
           i == 4 && dperiod > 4 && println("   ⋮")
       end
       if m > 0
          dperiod, nperiod = sys.B.dperiod, sys.B.nperiod
          subperiod = period/nperiod
          println(io, "\n\nInput matrix B: subperiod: $subperiod    #subperiods: $nperiod ")
          println(io, "time values: t[1:$dperiod]")
          for i in [1:min(4,dperiod); dperiod-2:dperiod]
              i <= 0 && break
              println("t[$i] = $(i*Ts)")
              display(mime, sys.B.M[i])
              i == 4 && dperiod > 4 && println("   ⋮")
          end
       else
          println(io, "\n\nEmpty input matrix B.")
       end
       if p > 0
          dperiod, nperiod = sys.C.dperiod, sys.C.nperiod
          subperiod = period/nperiod
          println(io, "\n\nOutput matrix C: subperiod: $subperiod    #subperiods: $nperiod ")
          println(io, "time values: t[1:$dperiod]")
          for i in [1:min(4,dperiod); dperiod-2:dperiod]
              i <= 0 && break
              println("t[$i] = $(i*Ts)")
              display(mime, sys.C.M[i])
              i == 4 && dperiod > 4 && println("   ⋮")
          end
       else
          println(io, "\n\nEmpty output matrix C.")
       end
       if m > 0 && p > 0
          dperiod, nperiod = sys.D.dperiod, sys.D.nperiod
          subperiod = period/nperiod
          println(io, "\n\nFeedthrough matrix D: subperiod: $subperiod    #subperiods: $nperiod ")
          println(io, "time values: t[1:$dperiod]")
          for i in [1:min(4,dperiod); dperiod-2:dperiod]
              i <= 0 && break
              println("t[$i] = $(i*Ts)")
              display(mime, sys.D.M[i])
              i == 4 && dperiod > 4 && println("   ⋮")
          end
       else
           println(io, "\n\nEmpty feedthrough matrix D.")
       end
       println(io, "\n\nPeriod:      $(sys.period) second(s).")
       println(io,     "Sample time: $Ts second(s).")
       println(io, "Periodic discrete-time state-space model.")
    elseif m > 0 && p > 0
       dperiod, nperiod = sys.D.dperiod, sys.D.nperiod
       subperiod = period/nperiod
       println(io, "\n\nFeedthrough matrix D: subperiod: $subperiod    #subperiods: $nperiod ")
       println(io, "time values: t[1:$dperiod]")
       for i in [1:min(4,dperiod); dperiod-2:dperiod]
           i <= 0 && break
           println("t[$i] = $(i*Ts)")
           display(mime, sys.D.M[i])
           i == 4 && dperiod > 4 && println("   ⋮")
       end
       println(io, "\n\nPeriod:      $(sys.period) second(s).")
       println(io,     "Sample time: $Ts second(s).")
       println(io, "Time-varying gains.")
    else
       println(io, "\nEmpty Periodic discrete-time state-space model.")
    end
end
function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, sys::PeriodicStateSpace{<:PeriodicArray})
    summary(io, sys); println(io)
    n = sys.nx
    N = sys.A.dperiod
    p, m = size(sys)
    period = sys.period
    Ts = sys.A.Ts
    if any(n .> 0)
        dperiod, nperiod = sys.A.dperiod, sys.A.nperiod
        subperiod = period/nperiod
        println(io, "\nState matrix A: subperiod: $subperiod    #subperiods: $nperiod ")
        println(io, "time values: t[1:$dperiod]")
        for i in [1:min(4,dperiod); dperiod-2:dperiod]
            i <= 0 && break
            println("t[$i] = $(i*Ts)")
            display(mime, sys.A.M[:,:,i])
            i == 4 && dperiod > 4 && println("   ⋮")
         end
        if m > 0
           dperiod, nperiod = sys.B.dperiod, sys.B.nperiod
           subperiod = period/nperiod
           println(io, "\n\nInput matrix B: subperiod: $subperiod    #subperiods: $nperiod ")
           println(io, "time values: t[1:$dperiod]")
           for i in [1:min(4,dperiod); dperiod-2:dperiod]
               i <= 0 && break
               println("t[$i] = $(i*Ts)")
               display(mime, sys.B.M[:,:,i])
               i == 4 && dperiod > 4 && println("   ⋮")
            end
        else
           println(io, "\n\nEmpty input matrix B.")
        end
        if p > 0
           dperiod, nperiod = sys.C.dperiod, sys.C.nperiod
           subperiod = period/nperiod
           println(io, "\n\nOutput matrix C: subperiod: $subperiod    #subperiods: $nperiod ")
           println(io, "time values: t[1:$dperiod]")
           for i in [1:min(4,dperiod); dperiod-2:dperiod]
               i <= 0 && break
               println("t[$i] = $(i*Ts)")
               display(mime, sys.C.M[:,:,i])
               i == 4 && dperiod > 4 && println("   ⋮")
            end
        else
           println(io, "\n\nEmpty output matrix C.")
        end
        if m > 0 && p > 0
           dperiod, nperiod = sys.D.dperiod, sys.D.nperiod
           subperiod = period/nperiod
           println(io, "\n\nFeedthrough matrix D: subperiod: $subperiod    #subperiods: $nperiod ")
           println(io, "time values: t[1:$dperiod]")
           for i in [1:min(4,dperiod); dperiod-2:dperiod]
               i <= 0 && break
               println("t[$i] = $(i*Ts)")
               display(mime, sys.D.M[:,:,i])
               i == 4 && dperiod > 4 && println("   ⋮")
            end
        else
           println(io, "\n\nEmpty feedthrough matrix D.")
        end
        println(io, "\n\nPeriod:      $(sys.period) second(s).")
        println(io,     "Sample time: $Ts second(s).")
        println(io, "Periodic discrete-time state-space model.")
    elseif m > 0 && p > 0
        dperiod, nperiod = sys.D.dperiod, sys.D.nperiod
        subperiod = period/nperiod
        println(io, "\n\nFeedthrough matrix D: subperiod: $subperiod    #subperiods: $nperiod ")
        println(io, "time values: t[1:$dperiod]")
        for i in [1:min(4,dperiod); dperiod-2:dperiod]
            i <= 0 && break
            println("t[$i] = $(i*Ts)")
            display(mime, sys.D.M[:,:,i])
            i == 4 && dperiod > 4 && println("   ⋮")
         end
        println(io, "\n\nPeriod:      $(sys.period) second(s).")
        println(io,     "Sample time: $Ts second(s).")
        println(io, "Time-varying gains.")
    else
        println(io, "\nEmpty Periodic discrete-time state-space model.")
    end
end

# """
#     PeriodicDiscreteStateSpace{T}(A::Vector{Matrix{T}}, E::Vector{Union{Matrix{T},UniformScaling}},
#                             B::Vector{Matrix{T}}, C::Vector{Matrix{T}}, D::Vector{Matrix{T}},
#                             Ts::Real) where T <: Number

# Construct a periodic discrete-time descriptor state-space model from a quintuple of
# `N`-dimensional vectors of matrices `(A,E,B,C,D)` and a sampling time `Ts`.
# `N` is the period of the system.

# If `SYS::PeriodicDiscreteStateSpace{T}` is a periodic descriptor system model object
# defined by the 5-tuple `(A,E,B,C,D)`, then:

# `SYS.A[i]` is the `ndx[i] × nx[i]` state matrix `A[i]` with elements of type `T`.
# The dimensions must satisfy `ndx[i] = nx[i+1]`, for `i = 1, ..., N-1` and `ndx[N] = nx[1]`.

# `SYS.E[i]` is the `ndx[i] × ndx[i]` descriptor matrix `E[i]` with elements of type `T`.
#  For a standard periodic discrete-time system, for all `i` `SYS.E[i] = I`, the `UniformScaling` of type `Bool`.

# `SYS.B[i]` is the `ndx[i] × nu[i]` system input matrix `B[i]` with elements of type `T`.
# All input dimensions must be the same.

# `SYS.C[i]` is the `nx[i] × ny[i]` system output matrix `C[i]` with elements of type `T`.

# `SYS.D` is the `ny × nu` system feedthrough matrix `D` with elements of type `T`.
# All output dimensions must be the same.

# `SYS.Ts` is the real sampling time `Ts`, where `Ts > 0` or `Ts = -1`, which
# indicates a discrete-time system with an unspecified sampling time.

# The vector of dimensions `nx`, ndx, `ny`, `nu` can be obtained as `SYS.nx`, `SYS.ndx`, `SYS.ny` and `SYS.nu`, respectively,
# the period N can be obtained as  `SYS.period`.
# """
# struct PeriodicDiscreteDescriptorStateSpace{T} <: AbstractPeriodicStateSpace
#     A::Vector{Matrix{T}}
#     E::Union{Vector{Matrix{T}},UniformScaling{Bool}}
#     B::Vector{Matrix{T}}
#     C::Vector{Matrix{T}}
#     D::Vector{Matrix{T}}
#     Ts::Float64
#     function PeriodicDiscreteDescriptorStateSpace{T}(A::Vector{Matrix{T}}, E::Union{Vector{Matrix{T}},UniformScaling{Bool}},
#                 B::Vector{Matrix{T}}, C::Vector{Matrix{T}}, D::Vector{Matrix{T}},  Ts::Real) where T
#         ps_validation(A, E, B, C, D, Ts)
#         new{T}(A, E, B, C, D, Float64(Ts))
#     end
# end
# function ps_validation(A::Vector{Matrix{T}}, E::Union{Vector{Matrix{T}},UniformScaling{Bool}},
#          B::Vector{Matrix{T}}, C::Vector{Matrix{T}}, D::Vector{Matrix{T}},  Ts::Real) where T
#     N = length(A)
#     desc = (E != I)
#     ((desc && (N != length(E))) || N != length(B) || N != length(C) || N != length(D)) &&
#         error("all vectors of state-space matrices must have the same length")
#     # Validate sampling time
#     Ts > 0 || Ts == -1 || error("Ts must be either a positive number or -1 (unspecified)")
#     N == 0 && (return N)

#     # Validate dimensions
#     ny = size.(D,1)
#     nu = size.(D,2)
#     # any(ny .!= ny[1]) && error("all matrices D[i] must have the same row size")
#     # any(nu .!= nu[1]) && error("all matrices D[i] must have the same column size")
#     ndx = size.(A,1)
#     nx = size.(A,2)

#     for i = 1:N
#         if desc
#            ndx[i] == LinearAlgebra.checksquare(E[i]) || error("all matrices A[i] and E[i] must have the same row size")
#         end
#         i == N && ndx[i] != nx[1] && error("the number of columns of A[i+1] must be equal to the number of rows of A[i]")
#         i != N && ndx[i] != nx[i+1] && error("the number of columns of A[i+1] must be equal to the number of rows of A[i]")
#     end
#     any(size.(B, 1) .!= ndx) && error("all matrices A[i] and B[i] must have the same row size")
#     any(size.(C, 2) .!= nx) &&  error("all matrices A[i] and C[i] must have the same column size")
#     any(size.(C, 1) .!= ny) && error("all matrices C[i] and D[i] must have the same row size")
#     any(size.(B, 2) .!= nu) && error("all matrices B[i] and D[i] must have the same column size")

#     return N
# end

# function Base.getproperty(sys::PeriodicDiscreteDescriptorStateSpace, d::Symbol)
#     if d === :nx
#         return size.(sys.A,2)
#     elseif d === :ndx
#         return size.(sys.A,1)
#     elseif d === :ny
#         return size.(sys.C,1)
#     elseif d === :nu
#         return size.(sys.B,2)
#     elseif d === :period
#         return length(sys.A)
#     else
#         getfield(sys, d)
#     end
# end
# Base.propertynames(sys::PeriodicDiscreteDescriptorStateSpace) =
#     (:nx, :ndx, :ny, :nu, :period, fieldnames(typeof(sys))...)
# """
#     size(sys) -> (p,m)
#     size(sys,1) -> p
#     size(sys,2) -> m

# Return the number of outputs `p` and the number of inputs `m` of a descriptor system `sys`.
# """
# function Base.size(sys::PeriodicDiscreteDescriptorStateSpace)
#     return sys.period == 0 ? (0,0) : size(sys.D[1])
# end
# function Base.size(sys::PeriodicDiscreteDescriptorStateSpace, d::Integer)
#     return d <= 2 ? size(sys)[d] : 1
# end
# function Base.length(sys::PeriodicDiscreteDescriptorStateSpace)
#     return 1
# end
# function Base.eltype(sys::PeriodicDiscreteDescriptorStateSpace)
#     return sys.period == 0 ? Float64 : eltype(sys.A[1])
# end

# # display sys
# Base.print(io::IO, sys::PeriodicDiscreteDescriptorStateSpace) = show(io, sys)
# Base.show(io::IO, sys::PeriodicDiscreteDescriptorStateSpace) = show(io, MIME("text/plain"), sys)

# function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, sys::PeriodicDiscreteDescriptorStateSpace)
#     summary(io, sys); println(io)
#     n = size.(sys.A,1)
#     N = sys.period
#     p, m = size(sys)
#     desc = (sys.E != I)
#     if any(n .> 0)
#        println(io, "\nState matrix A:")
#        for i = 1:N
#            println("@time: $i")
#            display(mime, sys.A[i])
#        end
#        if desc
#           println(io, "\n\nDescriptor matrix E:")
#           for i = 1:N
#               println("@time: $i")
#               display(mime, sys.E[i])
#           end
#        end
#        if m > 0
#           println(io, "\n\nInput matrix B:")
#           for i = 1:N
#               println("@time: $i")
#               display(mime, sys.B[i])
#           end
#        else
#           println(io, "\n\nEmpty input matrix B.")
#        end
#        if p > 0
#           println(io, "\n\nOutput matrix C:")
#           for i = 1:N
#             println("@time: $i")
#             display(mime, sys.C[i])
#           end
#        else
#           println(io, "\n\nEmpty output matrix C.")
#        end
#        if m > 0 && p > 0
#           println(io, "\n\nFeedthrough matrix D:")
#           for i = 1:N
#               println("@time: $i")
#               display(mime, sys.D[i])
#            end
#        else
#            println(io, "\n\nEmpty feedthrough matrix D.")
#        end
#        sys.Ts < 0 && println(io, "\n\nSample time: unspecified.")
#        sys.Ts > 0 && println(io, "\n\nSample time: $(sys.Ts) second(s).")
#        desc ? println(io, "Periodic discrete-time descriptor state-space model.") :
#               println(io, "Periodic discrete-time state-space model.")
#     elseif m > 0 && p > 0
#        println(io, "\nFeedthrough matrix D:")
#        show(io, mime, sys.D)
#        println(io, "\n\nStatic gains.")
#     else
#        println(io, "\nEmpty Periodic discrete-time tate-space model.")
#     end
# end

# function dss(psys::PeriodicDiscreteDescriptorStateSpace{T}; kstart::Int = 1, compacted::Bool = false,
#              fast::Bool = true, atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol,
#              rtol::Real = eps(real(float(one(T))))*iszero(max(atol1,atol2))) where T
#     indp(j,n) = mod(j-1,n)+1
#     desc = (psys.E != I)
#     K = psys.period
#     Ts = psys.Ts
#     LTs = Ts > 0 ? Ts*K : -1
#     #[ap,bp,cp,dp,ep] = pdsdata(psys);
#     # generate the stacked lifted system of (Grasselli and Longhi, 1991)
#     kind = indp.(kstart:kstart+K-1,K)
#     ni = psys.nx[kind]
#     mui = psys.ndx[kind]
#     m = psys.nu[kind]
#     n = sum(ni)
#     mukm1 = sum(view(mui,1:K-1))
#     nukm1 = sum(view(ni,1:K-1))
#     KK = kind[K]
#     kindm1 = view(kind,1:K-1)
#     A = [zeros(T, mui[K], nukm1) psys.A[KK]; blockdiag(view(psys.A,kindm1)...) zeros(T, mukm1, ni[K])] +
#         [zeros(T, mui[K], n); zeros(T, mukm1, ni[1]) desc ? -blockdiag(view(psys.E,kindm1)...) : -I]
#     E = [desc ? psys.E[KK] : I zeros(T, mui[K], n-ni[1]); zeros(T, mukm1, n)]
#     B = [ zeros(T, mui[K], sum(view(m,1:K-1))) psys.B[KK]; blockdiag(view(psys.B,kindm1)...) zeros(T, mukm1, m[K])]
#     C = blockdiag(view(psys.C,kind)...)
#     D = blockdiag(view(psys.D,kind)...)
#     if compacted
#        A, E, B, C, D = lsminreal(A, E, B, C, D; fast, atol1, atol2, rtol, contr = false, obs = false, noseig = true)
#     end
#     return dss(A, E, B, C, D; Ts = LTs)
# end

# function psreduc_fast(S::Vector{Matrix{T1}}, T::Vector{Matrix{T1}}; atol::Real = 0) where T1
#     # PSREDUC_FAST  Finds a matrix pair having the same finite and infinite
#     #               eigenvalues as a given periodic pair.
#     #
#     # [A,E] = PSREDUC_FAST(S,T,tol) computes a matrix pair (A,E) having
#     # the same finite and infinite eigenvalues as the periodic pair (S,T).
#     #

#     #   A. Varga 30-3-2004.
#     #   Revised: .
#     #
#     #   Reference:
#     #   A. Varga & P. Van Dooren
#     #      Computing the zeros of periodic descriptor systems.
#     #      Systems and Control Letters, vol. 50, 2003.

#     K = length(S)
#     if K == 1
#        return S[1], T[1]
#     else
#        m = size.(S,1)
#        n = size.(S,2)
#        if sum(m) > sum(n)
#            # build the dual pair
#         #    [S,T]=celltranspose(S,T);
#         #    [S{1:K}] = deal(S{K:-1:1});
#         #    [T{1:K-1}] = deal(T{K-1:-1:1});
#         #    m = size.(S,1)
#         #    n = size.(S,2)
#        end

#        si = S[1];  ti = -T[1];
#        tolr = atol
#        for i = 1:K-2
#            F = qr([ ti; S[i+1] ], ColumnNorm())
#            println("F = $(size(F.R))")
#            nr = minimum(size(F.R))
#            # compute rank of r
#            ss = abs.(diag(F.R[1:nr,1:nr]))
#            atol == 0 && ( tolr = (m[i]+m[i+1]) * maximum(ss) * eps())
#            rankr = count(ss .> tolr)

#            si = F.Q'*[si; zeros(T1,m[i+1],n[1])]; si=si[rankr+1:end,:]
#            ti = F.Q'*[ zeros(T1,size(ti,1),n[i+2]); -T[i+1]]; ti = ti[rankr+1:end,:]
#        end
#        a = [ zeros(T1,m[K],n[1]) S[K]; si ti]
#        e = [ T[K] zeros(T1,m[K],n[K]); zeros(size(si)...) zeros(size(ti)...)]
#        return a, e
#     end
# end

# function Base.getindex(sys::DST, inds...) where DST <: DescriptorStateSpace
#     size(inds, 1) != 2 &&
#         error("Must specify 2 indices to index descriptor state-space models")
#     rows, cols = index2range(inds...)
#     return DescriptorStateSpace{eltype(sys)}(copy(sys.A), copy(sys.E), sys.B[:, cols], sys.C[rows, :], sys.D[rows, cols], sys.Ts)
# end
# index2range(ind1, ind2) = (index2range(ind1), index2range(ind2))
# index2range(ind::T) where {T<:Number} = ind:ind
# index2range(ind::T) where {T<:AbstractArray} = ind
# index2range(ind::Colon) = ind
# function Base.lastindex(sys::DST, dim::Int) where DST <: DescriptorStateSpace
#     lastindex(sys.D,dim)
# end
# # Basic Operations
# function ==(sys1::DST1, sys2::DST2) where {DST1<:DescriptorStateSpace, DST2<:DescriptorStateSpace}
#     # fieldnames(DST1) == fieldnames(DST2) || (return false)
#     return all(getfield(sys1, f) == getfield(sys2, f) for f in fieldnames(DST1))
# end

# function isapprox(sys1::DST1, sys2::DST2; atol = zero(real(eltype(sys1))),
#                   rtol = rtol::Real =  ((max(size(sys1.A)...))+1)*eps(real(float(one(real(eltype(sys1))))))*iszero(atol)) where
#                   {DST1<:DescriptorStateSpace,DST2<:DescriptorStateSpace}
#     #fieldnames(DST1) == fieldnames(DST2) || (return false)
#     return all(isapprox(getfield(sys1, f), getfield(sys2, f); atol = atol, rtol = rtol) for f in fieldnames(DST1))
# end

# Base.copy(sys::DescriptorStateSpace{T}) where T = DescriptorStateSpace{T}(copy(sys.A), copy(sys.E), copy(sys.B), copy(sys.C), copy(sys.D), sys.Ts)

# # sum sys1+sys2
# function +(sys1::DescriptorStateSpace{T1}, sys2::DescriptorStateSpace{T2}) where {T1,T2}
#     sys1.nu == 1 && sys1.ny == 1 && (sys2.ny > 1 || sys2.nu > 1) && (return ones(T1,sys2.ny,1)*sys1*ones(T1,1,sys2.nu) + sys2)
#     sys2.nu == 1 && sys2.ny == 1 && (sys1.nu > 1 || sys1.ny > 1) && (return sys1 + ones(T1,sys1.ny,1)*sys2*ones(T1,1,sys1.nu))
#     #Ensure systems have same dimensions and sampling times
#     size(sys1) == size(sys2) || error("The systems have different shapes.")
#     Ts = promote_Ts(sys1.Ts,sys2.Ts)

#     T = promote_type(T1, T2)
#     n1 = sys1.nx
#     n2 = sys2.nx
#     A = [sys1.A  zeros(T,n1,n2);
#          zeros(T,n2,n1) sys2.A]
#     B = [sys1.B ; sys2.B]
#     C = [sys1.C sys2.C;]
#     D = [sys1.D + sys2.D;]
#     if sys1.E == I && sys2.E == I
#         E = I
#     else
#         E = [sys1.E  zeros(T,n1,n2);
#              zeros(T,n2,n1) sys2.E]
#     end

#     return DescriptorStateSpace{T}(A, E, B, C, D, Ts)
# end
# # difference sys1-sys2
# function -(sys1::DescriptorStateSpace{T1}, sys2::DescriptorStateSpace{T2}) where {T1,T2}
#     sys1.nu == 1 && sys1.ny == 1 && (sys2.ny > 1 || sys2.nu > 1) && (return ones(T1,sys2.ny,1)*sys1*ones(T1,1,sys2.nu) - sys2)
#     sys2.nu == 1 && sys2.ny == 1 && (sys1.nu > 1 || sys1.ny > 1) && (return sys1 - ones(T1,sys1.ny,1)*sys2*ones(T1,1,sys1.nu))
#     #Ensure systems have same dimensions and sampling times
#     size(sys1) == size(sys2) || error("The systems have different shapes.")
#     Ts = promote_Ts(sys1.Ts,sys2.Ts)

#     T = promote_type(T1, T2)
#     n1 = sys1.nx
#     n2 = sys2.nx

#     A = [sys1.A  zeros(T,n1,n2);
#          zeros(T,n2,n1) sys2.A]
#     B = [sys1.B ; sys2.B]
#     C = [sys1.C -sys2.C;]
#     D = [sys1.D - sys2.D;]
#     if sys1.E == I && sys2.E == I
#         E = I
#     else
#         E = [sys1.E  zeros(T,n1,n2);
#              zeros(T,n2,n1) sys2.E]
#     end

#     return DescriptorStateSpace{T}(A, E, B, C, D, Ts)
# end
# # negation -sys
# function -(sys::DescriptorStateSpace{T}) where T
#     return DescriptorStateSpace{T}(sys.A, sys.E, sys.B, -sys.C, -sys.D, sys.Ts)
# end
# # sys+mat and mat+sys
# function +(sys::DescriptorStateSpace{T1}, mat::VecOrMat{T2}) where {T1,T2}
#     p, m = typeof(mat) <: Vector ? (length(mat),1) : size(mat)
#     size(sys) == (p, m) || error("The input-output dimensions of system does not match the shape of matrix.")
#     T = promote_type(T1, T2)
#     return DescriptorStateSpace{T}(copy_oftype(sys.A,T), sys.E == I ? I : copy_oftype(sys.E,T),
#                                    copy_oftype(sys.B,T), copy_oftype(sys.C,T),
#                                    copy_oftype(sys.D,T) + mat, sys.Ts)
# end
# +(mat::VecOrMat{T1}, sys::DescriptorStateSpace{T2}) where {T1,T2} = +(sys,mat)
# # sys+I and I+sys
# function +(sys::DescriptorStateSpace{T1}, mat::UniformScaling{T2}) where {T1,T2}
#     size(sys,1) == size(sys,2) || error("The system must have the same number of inputs and outputs")
#     T = promote_type(T1, T2)
#     return DescriptorStateSpace{T}(copy_oftype(sys.A,T), sys.E == I ? I : copy_oftype(sys.E,T),
#                                    copy_oftype(sys.B,T), copy_oftype(sys.C,T),
#                                    copy_oftype(sys.D,T) + mat, sys.Ts)
# end

# # sys-mat and mat-sys
# function -(sys::DescriptorStateSpace{T1}, mat::VecOrMat{T2}) where {T1,T2}
#     p, m = typeof(mat) <: Vector ? (length(mat),1) : size(mat)
#     size(sys) == (p, m) || error("The input-output dimensions of system does not match the shape of matrix.")
#     T = promote_type(T1, T2)
#     return DescriptorStateSpace{T}(copy_oftype(sys.A,T), sys.E == I ? I : copy_oftype(sys.E,T),
#                                    copy_oftype(sys.B,T), copy_oftype(sys.C,T),
#                                    copy_oftype(sys.D,T) - mat, sys.Ts)
# end
# -(mat::VecOrMat{T1}, sys::DescriptorStateSpace{T2}) where {T1,T2} = +(-sys,mat)

# # sys-I and I-sys
# function -(sys::DescriptorStateSpace{T1}, mat::UniformScaling{T2}) where {T1,T2}
#     size(sys,1) == size(sys,2) || error("The system must have the same number of inputs and outputs")
#     T = promote_type(T1, T2)
#     return DescriptorStateSpace{T}(copy_oftype(sys.A,T), sys.E == I ? I : copy_oftype(sys.E,T),
#                                    copy_oftype(sys.B,T), copy_oftype(sys.C,T),
#                                    copy_oftype(sys.D,T) - mat, sys.Ts)
# end
# +(mat::UniformScaling{T1}, sys::DescriptorStateSpace{T2}) where {T1,T2} = +(sys,mat)
# -(mat::UniformScaling{T1}, sys::DescriptorStateSpace{T2}) where {T1,T2} = +(-sys,mat)
# # sys ± n and n ± sys
# function +(sys::DescriptorStateSpace{T1}, n::Number) where T1
#     T = promote_type(T1, eltype(n))
#     return DescriptorStateSpace{T}(copy_oftype(sys.A,T), sys.E == I ? I : copy_oftype(sys.E,T),
#                                    copy_oftype(sys.B,T), copy_oftype(sys.C,T),
#                                    copy_oftype(sys.D,T) .+ n, sys.Ts)
# end

# +(n::Number, sys::DescriptorStateSpace{T1}) where T1 = +(sys, n)
# -(n::Number, sys::DescriptorStateSpace{T1}) where T1 = +(-sys, n)
# -(sys::DescriptorStateSpace{T1},n::Number) where T1 = +(sys, -n)

# # multiplication sys1*sys2
# function *(sys1::DescriptorStateSpace{T1}, sys2::DescriptorStateSpace{T2}) where {T1,T2}
#     sys1.nu == 1 && sys1.ny == 1 && sys2.ny > 1 && (return dsdiag(sys1,sys2.ny)*sys2)
#     sys2.nu == 1 && sys2.ny == 1 && sys1.nu > 1 && (return sys1*dsdiag(sys2,sys1.nu))
#     sys1.nu == sys2.ny || error("sys1 must have same number of inputs as sys2 has outputs")
#     Ts = promote_Ts(sys1.Ts, sys2.Ts)
#     T = promote_type(T1, T2)
#     n1 = sys1.nx
#     n2 = sys2.nx

#     A = [sys1.A    sys1.B*sys2.C;
#          zeros(T,n2,n1)   sys2.A]
#     B = [sys1.B*sys2.D ; sys2.B]
#     C = [sys1.C   sys1.D*sys2.C;]
#     D = [sys1.D*sys2.D;]
#     if sys1.E == I && sys2.E == I
#         E = I
#     else
#         E = [sys1.E  zeros(T,n1,n2);
#              zeros(T,n2,n1) sys2.E]
#     end

#     return DescriptorStateSpace{T}(A, E, B, C, D, Ts)
# end
# # sys*mat
# function *(sys::DescriptorStateSpace{T1}, mat::VecOrMat{T2}) where {T1,T2}
#     p, m = typeof(mat) <: Vector ? (length(mat),1) : size(mat)
#     sys.nu == p || error("The input dimension of system does not match the number of rows of the matrix.")
#     T = promote_type(T1, T2)
#     return DescriptorStateSpace{T}(copy_oftype(sys.A,T),
#                                    sys.E == I ? I : copy_oftype(sys.E,T),
#                                    copy_oftype(sys.B,T)*to_matrix(T,mat), copy_oftype(sys.C,T),
#                                    copy_oftype(sys.D,T)*to_matrix(T,mat), sys.Ts)
# end
# # mat*sys
# function *(mat::VecOrMat{T1}, sys::DescriptorStateSpace{T2}) where {T1 <: Number,T2}
#     p, m = typeof(mat) <: Vector ? (length(mat),1) : size(mat)
#     sys.ny == m || error("The output dimension of system does not match the number of columns of the matrix.")
#     T = promote_type(T1, T2)
#     return DescriptorStateSpace{T}(copy_oftype(sys.A,T),
#                                    sys.E == I ? I : copy_oftype(sys.E,T),
#                                    copy_oftype(sys.B,T), to_matrix(T,mat)*copy_oftype(sys.C,T),
#                                    to_matrix(T,mat)*copy_oftype(sys.D,T), sys.Ts)
# end
# # sI*sys
# function *(s::Union{UniformScaling,Number}, sys::DescriptorStateSpace{T}) where T
#     T1 = promote_type(eltype(s),T)
#     return DescriptorStateSpace{T1}(copy_oftype(sys.A,T1), sys.E == I ? I : copy_oftype(sys.E,T1), copy_oftype(sys.B,T1),
#                                    lmul!(s,copy_oftype(sys.C,T1)), lmul!(s,copy_oftype(sys.D,T1)), sys.Ts)

# end
# # sys*sI
# function *(sys::DescriptorStateSpace{T},s::Union{UniformScaling,Number}) where T
#     T1 = promote_type(eltype(s),T)
#     return DescriptorStateSpace{T1}(copy_oftype(sys.A,T1), sys.E == I ? I : copy_oftype(sys.E,T1), rmul!(copy_oftype(sys.B,T1),s),
#                                    copy_oftype(sys.C,T1), rmul!(copy_oftype(sys.D,T1),s), sys.Ts)

# end

# # right division sys1/sys2 = sys1*inv(sys2)
# /(n::Union{UniformScaling,Number}, sys::DescriptorStateSpace) = n*inv(sys)
# /(sys::DescriptorStateSpace, n::Number) = sys*(1/n)
# /(sys::DescriptorStateSpace, n::UniformScaling) = sys*(1/n.λ)

# # left division sys1\sys2 = inv(sys1)*sys2
# \(n::Number, sys::DescriptorStateSpace) = (1/n)*sys
# \(n::UniformScaling, sys::DescriptorStateSpace) = (1/n.λ)*sys
# \(sys::DescriptorStateSpace, n::Union{UniformScaling,Number}) = inv(sys)*n

# # promotions
# Base.promote_rule(::Type{DST}, ::Type{P}) where {DST <: DescriptorStateSpace, P <: Number } = DST
# Base.promote_rule(::Type{DST}, ::Type{P}) where {DST <: DescriptorStateSpace, P <: UniformScaling } = DST
# Base.promote_rule(::Type{DST}, ::Type{P}) where {DST <: DescriptorStateSpace, P <: Matrix{<:Number} } = DST
# #Base.promote_rule(::Type{AbstractVecOrMat{<:Number}}, ::Type{DST}) where {DST <: DescriptorStateSpace} = DST
# #Base.promote_rule(::Type{DST}, ::Type{Union{AbstractVecOrMat{<:Number}, Number, UniformScaling}}) where {DST <: DescriptorStateSpace} = DST
# #Base.promote_rule(::Type{Union{AbstractVecOrMat{<:Number}, Number, UniformScaling}}, ::Type{DST}) where {DST <: DescriptorStateSpace} = DST
# # Base.promote_rule(::Type{DST}, ::Type{UniformScaling}) where {DST <: DescriptorStateSpace, P <: Number} = DST
# # Base.promote_rule(::Type{P}, ::Type{DST}) where {DST <: DescriptorStateSpace, P <: Number} = DST

# # conversions
# function Base.convert(::Type{DST}, p::Number) where {T, DST <: DescriptorStateSpace{T}}
#     T1 = promote_type(eltype(p),T)
#     dss(T1(p), Ts = sampling_time(DST))
# end
# sampling_time(::Type{DST}) where {DST <: DescriptorStateSpace} = 0.


text:
"""A node which lags its input by a fixed number of knots."""
struct Lag{T} <: NodeOp{T}
    n::Int64
end

struct LagState{T} <: NodeEvaluationState
    # A buffer of the last number of values required.
    value_buffer::CircularBuffer{T}
end

function create_evaluation_state(::Tuple{Node}, op::Lag{T}) where {T}
    return LagState(CircularBuffer{T}(op.n))
end

function run_node!(
    op::Lag{T},
    state::LagState{T},
    ::DateTime,  # time_start
    ::DateTime,  # time_end
    input::Block{T},
) where {T}
    if isempty(input)
        # The input is empty, so regardless of the state we cannot lag onto anything.
        # The state is unmodified, and return the empty block.
        return input
    end

    w = op.n  # The window length.
    ni = length(input)  # The input length.
    ns = length(state.value_buffer)  # The state buffer length.

    # The output length.
    #   - at most the input size
    #   - at most total number of knots available in this batch: max((ns + ni) - w, 0)
    no = min(ni, max((ns + ni) - w, 0))

    # The number of elements to take from the state buffer.
    #   - at most the window size
    #   - at most the number of available elements in the buffer
    #   - at most the output size
    m = min(w, ns, no)

    times = @inbounds @view(input.times[(1 + ni - no):end])
    values = if m == 0
        @inbounds @view(input.values[1:no])
    else
        # Merge values from the state and in the Block.
        result = _allocate_values(T, no)
        for i in 1:m
            @inbounds result[i] = popfirst!(state.value_buffer)
        end
        @inbounds result[(m + 1):end] = @view(input.values[1:(no - m)])
        result
    end

    # Update the state with the remaining values in the input.
    @inbounds append!(state.value_buffer, @view(input.values[(1 + no - m):end]))

    return Block(unchecked, times, values)
end

"""
    lag(x::Node, n::Integer)

Construct a node which takes values from `x`, but lags them by `n` knots.

This means that we do not introduce any new timestamps that do not appear in `x`, however
we will not emit knots for the first `n` values that appear when evaluating `x`.
"""
function lag(x::Node, n::Integer)
    return if n == 0
        # Optimisation.
        x
    elseif n < 0
        throw(ArgumentError("Cannot lag by $n."))
    elseif _is_constant(x)
        # Constant nodes shouldn't be affected by lagging.
        x
    else
        obtain_node((x,), Lag{value_type(x)}(n))
    end
end

"""
    diff(x::Node[, n=1])

Compute the `n`-knot difference of `x`, i.e. `x - lag(x, n)`.
"""
function Base.diff(x::Node, n::Integer=1)
    n > 0 || throw(ArgumentError("n must be positive, got $n"))
    return x - lag(x, n)
end


text:
using PartonDensity
using Test

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


text:
using GLMakie
function flower(n; npetals = 8)
    n = div(n, npetals)
    x = mapreduce(hcat, (1:npetals) .* (2π/npetals)) do θ
        ct = cos(θ)
        st = sin(θ)

        x0 = tanh.(randn(1, n) .- 1) .+ 4.0 .+ 0.05.* randn(1, n)
        y0 = randn(1, n) .* 0.3

        x₁ = x0 * cos(θ) .- y0 * sin(θ)
        x₂ = x0 * sin(θ) .+ y0 * cos(θ)
        vcat(x₁, x₂)
    end
    _y = mapreduce(i -> fill(i, n), vcat, 1:npetals)
    y = zeros(npetals, length(_y))
    foreach(i -> y[_y[i], i] = 1, 1:length(_y))
    Float32.(x), Float32.(y)
end

x, y = flower(900)
scatter(x[1,:], x[2,:], color = mapslices(argmax, y, dims = 1)[:])

#######
#   Define a Tracked Array for operator overloading AD
#######
struct TrackedArray{T,N,V<:AbstractArray{T,N}} <: AbstractArray{T,N}
    value::V
    grad::Union{Nothing,V}
    tape::Vector{Any}
end

TrackedArray(a::AbstractArray) = TrackedArray(a, similar(a) .= 0, [])
TrackedMatrix{T,V} = TrackedArray{T,2,V} where {T,V<:AbstractMatrix{T}}
TrackedVector{T,V} = TrackedArray{T,1,V} where {T,V<:AbstractVector{T}}
Base.size(a::TrackedArray) = size(a.value)
Base.show(io::IO, ::MIME"text/plain", a::TrackedArray) = show(io, a)
Base.show(io::IO, a::TrackedArray) = print(io, "TrackedArray($(size(a.value)))")
value(A::TrackedArray) = A.value
resetgrad!(A::TrackedArray) = (A.grad .= 0; empty!(A.tape))
value(A) = A
track(A) = TrackedArray(A)
track(a::Number) = TrackedArray(reshape([a], 1, 1))

function accum!(A::TrackedArray)
    isempty(A.tape) && return(A.grad)
    A.grad .= sum(g(accum!(r)) for (r, g) in A.tape)
    empty!(A.tape)
    A.grad
end

#######
#   Define AD rules for few operations appearing in FFNN
#######
import Base: +, *
import Base.Broadcast: broadcasted
function *(A::TrackedMatrix, B::TrackedMatrix)
    a, b = value.((A, B))
    C = track(a * b)
    push!(A.tape, (C, Δ -> Δ * b'))
    push!(B.tape, (C, Δ -> a' * Δ))
    C
end

function *(A::TrackedMatrix, B::AbstractMatrix)
    a, b = value.((A, B))
    C = track(a * b)
    push!(A.tape, (C, Δ -> Δ * b'))
    C
end

function broadcasted(::typeof(+), A::TrackedMatrix, B::TrackedVector)
    C = track(value(A) .+ value(B))
    push!(A.tape, (C, Δ -> Δ))
    push!(B.tape, (C, Δ -> sum(Δ, dims = 2)[:]))
    C
end

function σ(x::Real)
    t = @fastmath exp(-abs(x))
    y = ifelse(x ≥ 0, inv(1 + t), t / (1 + t))
    ifelse(x > 40, one(y), ifelse(x < -80, zero(y), y))
end

broadcasted(::typeof(identity), A::TrackedArray) = A

function broadcasted(::typeof(σ), A::TrackedArray)
    Ω = σ.(value(A))
    C = track(Ω)
    push!(A.tape, (C, Δ -> Δ .* Ω .* (1 .- Ω)))
    C
end

function mse(A::TrackedMatrix, B::AbstractMatrix)
    n = size(A, 2)
    a = value(A)
    c = similar(a, 1, 1)
    c .= sum((a .- B).^2)/2n
    C = track(c)
    push!(A.tape, (C, Δ -> Δ .* (a .- B) ./ n))
    C
end

mse(x::AbstractMatrix, y::AbstractMatrix) = sum((x - y).^2) / (2*size(x,2))

#######
#   Define a Dense layer
#######
struct Dense{F,W,B}
    σ::F
    w::W
    b::B
end

Base.show(io::IO, m::Dense) = print(io, "Dense($(size(m.w,2)) → $(size(m.w,1)))")
Dense(i::Int, o::Int, σ = identity) = Dense(σ, randn(Float32, o, i), randn(Float32, o))
track(m::Dense) = Dense(m.σ, track(m.w), track(m.b))
track(m::ComposedFunction) = track(m.outer) ∘ track(m.inner)
(m::Dense)(x) = m.σ.(m.w * x .+ m.b)
params(m::ComposedFunction) = vcat(params(m.outer), params(m.inner))
params(m::Dense) = [m.w, m.b]

#######
#   Let's try to actually train a model
#######
x, y = flower(900)
function initmodel()
    m₁ = track(Dense(2, 20, σ))
    m₂ = track(Dense(20, 20, σ))
    m₃ = track(Dense(20, size(y,1)))
    m = m₃ ∘ m₂ ∘ m₁
end
m = initmodel()
m(x) |> value

######
#   Let's try to learn the parameters
######
α = 0.01
ps = params(m)
@elapsed for i in 1:10000
    foreach(resetgrad!, ps)
    loss = mse(m(x), y)
    fill!(loss.grad, 1)
    foreach(accum!, ps)
    foreach(x -> x.value .-= α .* x.grad, ps)
    mod(i,250) == 0 && println("loss after $(i) iterations = ", sum(value(loss)))
end

all(mapslices(argmax, value(m(x)), dims = 1)[:] .== mapslices(argmax, y, dims = 1)[:])
scatter(x[1,:], x[2,:], color = mapslices(argmax, value(m(x)), dims = 1)[:])

######
#   Let's try to move the computation to GPU
######
using CUDA
gpu(x::AbstractArray) = CuArray(x)
gpu(x::TrackedArray) = TrackedArray(CuArray(value(x)))
gpu(m::Dense) = Dense(m.σ, gpu(m.w), gpu(m.b))
gpu(m::ComposedFunction) = gpu(m.outer) ∘ gpu(m.inner)

gx, gy = gpu(x), gpu(y)
m = gpu(m)
ps = params(m)
@elapsed for i in 1:10000
    foreach(resetgrad!, ps)
    loss = mse(m(gx), gy)
    fill!(loss.grad, 1)
    foreach(accum!, ps)
    foreach(x -> x.value .-= α .* x.grad, ps)
    mod(i,250) == 0 && println("loss after $(i) iterations = ", sum(value(loss)))
end

#######
#   Why we see a small speed-up? The problem is small
#######
using BenchmarkTools
p = randn(Float32, 20, 2)
@benchmark $(p) * $(x)
gp = gpu(p)
@benchmark $(gp) * $(gx)

######
#   Let's verify the gradients
######
using FiniteDifferences
ps = [m₃.w, m₃.b, m₂.w, m₂.b, m₁.w, m₁.b]
map(ps) do p
    foreach(resetgrad!, ps)
    loss = mse(m(x), y)
    fill!(loss.grad, 1)
    foreach(accum!, ps)
    accum!(p)
    θ = deepcopy(value(p))
    Δθ = deepcopy(p.grad)
    f = θ -> begin
        p.value .= θ
        value(mse(m(x), y))
    end
    sum(abs2.(grad(central_fdm(5, 1), f, θ)[1] - Δθ))
end


text:
"""
url to remote JuliaPOMDP organization repo
"""
const REMOTE_URL = "https://github.com/JuliaPOMDP/"

"""
url to JuliaPOMDP registry
"""
const POMDP_REGISTRY = "https://github.com/JuliaPOMDP/Registry"

# TODO (max): would it be better to have a dict of form: string => url for solvers?
"""
Set containing string names of officially supported solvers and utility packages
(e.g. `MCTS`, `SARSOP`, `POMDPToolbox`, etc).
If you have a validated solver that supports the POMDPs.jl API,
contact the developers to add your solver to this list.
"""
const NATIVE_PACKAGES = Set{String}(
                        ["POMDPModels",
                         "ParticleFilters",

                         "POMDPModelTools",
                         "POMDPPolicies",
                         "BeliefUpdaters",
                         "POMDPSimulators",
                         "RLInterface",

                         "DiscreteValueIteration",
                         "MCTS",
                         "QMDP",
                         "FIB",
                         "POMCP",
                         "BasicPOMCP",
                         "ARDESPOT",
                         "DESPOT",
                         "MCVI",
                         "GenerativeModels",
                         "POMDPModels",
                         "POMDPXFiles",
                         "POMDPFiles",
                         "ParticleFilters",
                         "TabularTDLearning",
                         "POMDPReinforce",
                         "POMCPOW",
                         "AEMS",
                         "IncrementalPruning",
                         "LocalApproximationValueIteration",
                         "DeepQLearning"
                        ])

# Packages registered on METADATA
const REGISTERED_PACKAGES = Set{String}(["ParticleFilters",
                                         "POMDPModels",
                                        ])

const NON_NATIVE_PACKAGES = Set{String}(
                             ["SARSOP",
                             "POMDPSolve"])

const SUPPORTED_PACKAGES = union(NATIVE_PACKAGES, NON_NATIVE_PACKAGES)

const EXPORTED_TYPES = [MDP,
                        POMDP,
                        Policy,
                        Simulator,
                        Solver,
                        Updater,
                        Union{MDP, POMDP}]


text:
using DataDeps

register(DataDep(
  "optical-recognition-handwritten-digits",
  "http://archive.ics.uci.edu/ml/datasets/Optical+Recognition+of+Handwritten+Digits",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/optdigits/optdigits.tra",
  "e1b683cc211604fe8fd8c4417e6a69f31380e0c61d4af22e93cc21e9257ffedd",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=65,
        feature_indices=1:64,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
@symbol_func function b(cur_reactor::AbstractReactor)
  cur_P = P_W(cur_reactor)

  cur_n = cur_reactor.n_bar

  if isa(cur_n, SymEngine.Basic)
    cur_P = subs(
      cur_P, symbols(:n_bar) => calc_n_bar(cur_reactor)
    )
  end

  cur_b = 1.23

  cur_b += 0.074 * log( cur_P )

  cur_b
end


text:
module MacroHelpers

using MacroTools

"""
    ensure_symbol(input)

This function ensures that the given argument has a Symbol type
"""
ensure_symbol(input::Symbol) = input
ensure_symbol(input) = error("$(input) is not a symbol.")

"""
    bottom_type(type)

This function returns `T` expression for the following input expressions:
1. typeof(T)
2. Type{ <: T }
3. Type{ T }
4. T

# Arguments
- `type`: Type expression to be lowered

# See also: [`extract_fformtype`](@ref), [`upper_type`](@ref)
"""
function bottom_type(type)
    @capture(type, (typeof(T_)) | (Type{ <: T_ }) | (Type{ T_ }) | (T_)) ||
        error("Expression $(type) doesnt seem to be a valid type expression.")
    return T
end

"""
    upper_type(type)

This function returns `Type{ <: T }` expression for the following input expressions:
1. typeof(T)
2. Type{ <: T }
3. Type{ T }
4. T

# Arguments
- `type`: Type expression to be extended

# See also: [`extract_fformtype`](@ref), [`bottom_type`](@ref)
"""
function upper_type(type)
    if @capture(type, typeof(T_))
        return :(typeof($(ensure_symbol(T))))
    else
        return :(Type{ <: $(bottom_type(type)) })
    end
end

"""
    proxy_type(proxy, type)

Returns a type wrapped with a proxy type in a form of `ProxyType{ <: Type }`.

# Arguments
- `proxy`: Proxy type used to wrap `type`
- `type`: Type to be wrapped
"""
function proxy_type(proxy, type::Symbol)
    return :($(proxy){ <: $(type) })
end

function proxy_type(proxy, type::Expr)
    if @capture(type, NTuple{N_, T_})
        return :(NTuple{ $N, <: $(proxy_type(proxy, T)) })
    elseif @capture(type, Tuple{ T__ })
        return :(Tuple{ $(map(t -> :( <: $(proxy_type(proxy, t))), T)...) })
    else
        return :($(proxy){ <: $(type) })
    end
end

"""
    rearranged_tuple(name::Symbol, length::Int, swap::Tuple{Int, Int})
"""
function rearranged_tuple(name::Symbol, length::Int, swap::Tuple{Int, Int})
    args = map(i -> :($(name)[$i]), 1:length)
    i, j = first(swap), last(swap)
    args[i], args[j] = args[j], args[i]
    return Expr(:tuple, args...)
end

"""
    @proxy_methods(proxy_type, proxy_getter, proxy_methods)

Generates proxy methods for a specified `proxy_type` using `proxy_getter`. For example:
```julia
@proxy_methods Message getdata [
    Distributions.mean,
    Distributions.var
]
```

generates:
```julia
Distributions.mean(proxy::Message) = Distributions.mean(getdata(proxy))
Distributions.var(proxy::Message)  = Distributions.mean(getdata(proxy))
```
"""
macro proxy_methods(proxy_type, proxy_getter, proxy_methods)

    @capture(proxy_methods, [ methods__ ]) || error("Invalid specification of proxy methods, should be an array of methods")

    output      = Expr(:block)
    output.args = map(method -> :(($method)(proxy::$(proxy_type)) = ($method)($(proxy_getter)(proxy))), methods)

    return esc(output)
end

function expression_convert_eltype(eltype::Type{T}, expr::Expr) where T
    if @capture(expr, f_(args__))
        return :(ReactiveMP.convert_eltype($f, $T, $expr))
    elseif @capture(expr, (elems__, ))
        entries = map(elems) do elem
            @capture(elem, (name_ = value_)) || error("Invalid expression specification in expression_convert_eltype() function: $expr. Expression should be in the form of a constructor call or tuple of (name = value) elements.")
            return (name, value)
        end
        return Expr(:tuple, map((entry) -> :($(first(entry)) = $(expression_convert_eltype(eltype, last(entry)))), entries)...)
    end
    error("Invalid expression specification in expression_convert_eltype() function: $expr. Expression should be in the form of a constructor call or tuple of (name = value) elements.")
end

end


text:
using Random
using LinearAlgebra: norm
using Plots

@testset "VortexFilament constructor" begin
    v1 = rand(3)
    v2 = rand(3)
    v3 = rand(3)

    @test VortexFilament(1.0,[v1,v2]) == VortexFilament(1.0,[v1,v2],[Segment(v1,v2),Segment(v2,v1)],[1,2],Int64[])
    @test VortexFilament(1.0,[v1,v2,v3]) == VortexFilament(1.0,[v1,v2,v3],[Segment(v1,v2),Segment(v2,v3),Segment(v3,v1)],[1,2,3],Int64[])
    @test VortexFilament(1.0,[v1,v2,v3],[2]) == VortexFilament(1.0,[v1,v2,v3],[Segment(v1,v2),Segment(v2,v3),Segment(v3,v1)],[1,3],[2])
end

@testset "isinf" begin
    v1 = [Inf,-1,0]
    v2 = [0,-1,0]
    v3 = [0,1,0]
    v4 = [Inf,1,0]
    vf = VortexFilament(1,[v1,v2,v3,v4])
    @test isinf(vf)

    v1 = [-Inf,1,0]
    v2 = [Inf,1,0]
    vf = VortexFilament(1,[v1,v2])
    @test isinf(vf)
end

@testset "issemiinf" begin
    v1 = rand(3)
    v2 = [Inf,rand(Float64),rand(Float64)]
    vf = VortexFilament(1,[v1,v2])
    @test issemiinf(vf)

    v1 = [rand(Float64),-Inf,rand(Float64)]
    v2 = rand(3)
    v3 = rand(3)
    vf = VortexFilament(1,[v1,v2,v3])
    @test issemiinf(vf)

    v1 = [Inf,-1,0]
    v2 = rand(3)
    v3 = rand(3)
    v4 = [Inf,1,0]
    vf = VortexFilament(1,[v1,v2,v3,v4])
    @test issemiinf(vf) == false
end

@testset "_closestpointonfinitesegmentaxis" begin
    v1 = [0,0,0]
    v2 = [1,0,0]
    @test VortexFilaments._closestpointonfinitesegmentaxis(Segment(v1,v2),[2,1,0]) == ([2,0,0],2)
    @test VortexFilaments._closestpointonfinitesegmentaxis(Segment(v1,v2),[-2,1,0]) == ([-2,0,0],-2)
end

@testset "infdir" begin
    v1 = rand(3)
    v2 = rand(3)
    v3 = rand(3)
    vtwoinf = [rand(Float64), Inf, Inf]
    vplusinfx = [Inf,rand(Float64),rand(Float64)]
    vplusinfy = [rand(Float64),Inf,rand(Float64)]
    vplusinfz = [rand(Float64),rand(Float64),Inf]
    vmininfx = [-Inf,rand(Float64),rand(Float64)]
    vmininfy = [rand(Float64),-Inf,rand(Float64)]
    vmininfz = [rand(Float64),rand(Float64),-Inf]
    @test_throws ErrorException VortexFilaments.infdir(VortexFilament(1.0,[v1,v2,vtwoinf]))
    @test_throws ErrorException VortexFilaments.infdir(VortexFilament(1.0,[vplusinfx,v2,vplusinfy]))
    @test VortexFilaments.infdir(VortexFilament(1.0,[v1,v2,vplusinfx])) == 1
    @test VortexFilaments.infdir(VortexFilament(1.0,[v1,v2,vplusinfy])) == 2
    @test VortexFilaments.infdir(VortexFilament(1.0,[v1,v2,vplusinfz])) == 3
    @test VortexFilaments.infdir(VortexFilament(1.0,[v1,v2,vmininfx])) == 1
    @test VortexFilaments.infdir(VortexFilament(1.0,[v1,v2,vmininfy])) == 2
    @test VortexFilaments.infdir(VortexFilament(1.0,[v1,v2,vmininfz])) == 3
    @test VortexFilaments.infdir(VortexFilament(1.0,[vplusinfx,v1,v2,vmininfx])) == 1
    @test VortexFilaments.infdir(VortexFilament(1.0,[vplusinfx,v1,v2,vplusinfx])) == 1
end

@testset "inducevelocity" begin
    v0 = [0,0,0]
    vplusinfx = [Inf,rand(Float64),rand(Float64)]
    vmininfx = [-Inf,rand(Float64),rand(Float64)]
    vplusinfx2 = [Inf,rand(Float64),rand(Float64)]
    vplusinfy = [rand(Float64),Inf,rand(Float64)]
    v1 = [-1,-1,0]
    v2 = [-1,1,0]
    v3 = [1,1,1]
    v4 = [1,-1,1]
    r = rand(Float64) + 1
    rvec = rand(3) .+ 1

    # Special cases for evaluation point
    s = Segment(v0,v3)
    @test inducevelocity(s,r*v0+(r-1)*v3) == [0.0,0.0,0.0] # evaluation point on segment
    @test inducevelocity(s,v0) == [0.0,0.0,0.0] # evaluation point on vertex 1
    @test inducevelocity(s,v3) == [0.0,0.0,0.0] # evaluation point on vertex 2
    @test inducevelocity(s,v0+(1+r)*v3) == [0.0,0.0,0.0] # evaluation point on segment axis but not on segment itself
    @test inducevelocity(Segment([0,1,2],[Inf,1,2]),[1,1,2]) == [0.0,0.0,0.0] # evaluation point on semi-infinite segment
    @test inducevelocity(Segment([0,1,2],[Inf,1,2]),[-1,1,2]) == [0.0,0.0,0.0] # evaluation point on semi-infinite segment axis but not on segment itself
    @test inducevelocity(Segment([-Inf,1,2],[Inf,1,2]),[1,1,2]) == [0.0,0.0,0.0] # evaluation point on infinite segment

    # Special segments
    s = Segment(v0,vplusinfx) # semi-infinite segment
    @test isapprox(norm(inducevelocity(s,[0.0,rvec[2],rvec[3]])), 1/(4π*norm([rvec[2],rvec[3]])), atol=1e-3)
    s = Segment(vmininfx,vplusinfx) # infinite segment
    centerpoint = 0.5*(vplusinfx + vmininfx)
    @test isapprox(norm(inducevelocity(s,rvec)), 1/(2π*norm([rvec[2]-centerpoint[2],rvec[3]-centerpoint[3]])), atol=1e-3)
    s = Segment(vmininfx,vplusinfy) # infinite segment with different infinite coordinates
    @test_throws ErrorException inducevelocity(s,rvec)
    s = Segment(vplusinfx,vplusinfx2) # infinite segment with non-opposite infinite vertices
    @test_throws ErrorException inducevelocity(s,rvec)

    # Horseshoe vortex
    b = rand(Float64)+1
    Γ = rand(Float64)
    v1 = [Inf,-b/2,0]
    v2 = [0,-b/2,0]
    v3 = [0,b/2,0]
    v4 = [Inf,b/2,0]
    vf = VortexFilament(Γ,[v1,v2,v3,v4]) # with this order of vertices, the circulation should be positive to have a negative downwash
    downwash(Γ,b,y) = -Γ/(4π)*b/((b/2)^2-y^2)
    yrange = range(-b/2,b/2,length=10)
    xevals = [[0.0,y,0.0] for y in yrange]
    @test isapprox((v->v[3]).(inducevelocity.(Ref(vf),xevals[2:end-1])), downwash.(Γ,b,yrange[2:end-1]), atol=1e-2)

    # Infinite vortex filament
    Γ = rand(Float64)
    v1 = [-Inf,0,0]
    v2 = [Inf,0,0]
    vf = VortexFilament(Γ,[v1,v2])
    @test isapprox(inducevelocity(vf,[0,1,0]), [0.0,0.0,Γ/(2π)], atol=1e-2)

    # Normal example
    v1 = [-1,-1,0]
    v2 = [-1,1,0]
    v3 = [1,1,1]
    v4 = [1,-1,1]

    vf = VortexFilament(1.0,[v1,v2,v3,v4])
    @test isapprox(inducevelocity(vf,[-1,0,1]),[0.15667; 0.0; -0.12386]; atol=1e-5)

    v1 = rand(3)
    v2 = rand(3)
    vtwodim = rand(2)
    @test_throws DimensionMismatch inducevelocity(VortexFilament(1.0,[v1,vtwodim,v2]),rand(3))

end

@testset "plot" begin
    vertices = [[0.0,0.0,0.0], [0.0,1.0,0.0], [1.0,1.0,0.0], [1.0,0.0,0.0]]
    vf = VortexFilament(1.0,vertices)
    plot(vf)

    vertices = [[0,0.0,0.0], [Inf,0.0,0.0]]
    vf = VortexFilament(1.0,vertices) # semi-infinite vortex filament
    plot(vf,xlims=[0,2])

    vertices = [[-Inf,0.0,0.0], [Inf,0.0,0.0]]
    vf = VortexFilament(1.0,vertices) # infinite vortex filament
    plot(vf,xlims=[-2,2])
end


text:
@testset "720.longest-word-in-dictionary.jl" begin
    @test longest_common_word(["w", "wo", "wor", "worl", "world"]) == "world"
    @test longest_common_word(["a", "banana", "app", "appl", "ap", "apply", "apple"]) ∈
        ("apply", "apple")
end


text:
# Model tests, to ensure that the loaded models are running
# with a decent accuracy.

# Tests files have been taken from ONNX.
# You need PyCall and numpy installed to run these tests.

using Keras
using PyCall
@pyimport numpy

function read_input(filename)
    f = numpy.load(open(filename))
    return convert(Array{Float64, N} where N, reshape(permutedims(f, (5,4,3,2,1)), 224,224,3,1))
end

function read_output(filename)
    f = numpy.load(open(filename))
    return reshape(permutedims(f, (3,2,1)), 1000,1)
end

# ResNet test
loss(x, y) = sum((x .- y).^2)
netloss = 0
for x=0:2
    output = Keras.load("resnet.json", "resnet.h5", read_input("inputs$x.npy"))
    op_expected = read_output("outputs$x.npy")
    netloss += loss(output, op_expected)
end

println("Average net loss is ", netloss/3)


text: # defs2.jl : any definition that Clang did not get correctly


text:
using Maracas
# abstract MaracasTestSet <: Base.Test.AbstractTestSet
is_magenta(ts::MaracasTestSet)=occursin(Base.text_colors[:magenta], ts.description)
is_blue(ts::MaracasTestSet)=occursin(Base.text_colors[:blue], ts.description)
is_cyan(ts::MaracasTestSet)=occursin(Base.text_colors[:cyan], ts.description)

@describe "it is a test suite" begin
    @it "has specs" begin
        a_spec = @it("is made with macro '@it'", begin end)
        @test occursin("[Spec]", a_spec.description)
    end
    @it "has tests" begin
        a_test = @unit("made with macro '@unit'", begin end)
        @test occursin("[Test]", a_test.description)
    end

    @unit "test suite title is magenta by default" begin
        nested_describe = @describe("you can document your code with your tests", begin end)
        @test is_magenta(nested_describe)
    end

    @unit "spec title is cyan by default" begin
        @test is_cyan(@it("is cyan", begin end))
    end
    @unit "test title is blue by default" begin
        @test is_blue(@unit("in blue", begin end))
    end

    @unit "'it' is prepended to specs" begin
        description = "had a spec description not starting with it"
        a_spec = @it("had a spec description not starting with it", begin end)
        @test occursin(string("it ", description), a_spec.description)
    end

    @skip @describe "a whole describe can be skipped" begin
        @it "should not be executed" begin
            @test false
        end
    end
    @skip @it "can skip '@it' with @skip" begin
        @test false
    end
    @skip @unit "'@unit' can be skipped with @skip" begin
        @test false
    end
    @skip @testset "@skip works with default testset but it's not shown in results" begin
        @test false
    end
    @unit "@skip can also skip @test assertions" begin
        @skip @test false
        @skip @test_throws false
        @skip @test_skip false
        @skip @test_broken false
    end
end


text: Base.fill!(a::Image, c::Color) = Base.fill!(a.pic, c)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using ThermProps
using Base.Test

# write your own tests here
LiquidDensity("1,1,1-Trichloroethane",298.15)
LiquidDensity("1,1,1-Trichloroethane1",298.15)
LiquidDensity("1,1,1-Trichloroethane",600.00)
LiquidDensity()

GasViscosity("1,3-Butadiene",298.15)
GasViscosity("1,3-Butadiene1",298.15)
GasViscosity("1,3-Butadiene",800.0)
GasViscosity()

LiquidViscosity("Carbon monoxide",71.00)
LiquidViscosity("Carbon monoxide1",298.15)
LiquidViscosity("Carbon monoxide",600.00)
LiquidViscosity()


text:
"""
Removes asset node (sne) from the network
"""
function drop_node!(sne::SimNE,model::ABM)
    set_down!(sne)
    g = model.ntw_graph
    log_info(model.ticks," All neighbours of $(sne.id) are: $(all_neighbors(model.ntw_graph,get_address(sne.id,g))) ")
    for nb in all_neighbors(model.ntw_graph,get_address(sne.id,g))
       sne_nb = getindex(model,get_eid(nb,model))
       link_down!(sne_nb,sne.id,model)
    end

    #delete pending queries to control
    sne.requested_ctl = Dict()

    #it simulates control detects sne down:
    aid = get_control_agent(sne.id,model)
    a = getindex(model,abs(aid))
    controlled_sne_down!(a,sne.id,model)

    #soft remove
    model.ntw_graph = soft_remove_vertex(g,get_address(sne.id,g))
    push!(model.dropped_nodes,(model.ticks,sne.id))
end

"""
   Function to deteriorate a network element
"""
function deteriorate!(sne::SimNE,model::ABM)
    state = get_state(sne)
    if state.up
        # if sne.id == 5
        #     state.rul = state.rul - 3
        # else
            state.rul = lineal_d(sne.maintenance.deterioration_parameter,state.rul,1)
        # end

        if state.rul == 0
            a = getindex(model,abs(sne.controller_id))
            drop_node!(sne,model)
            schedule_event!(a,CTL_Event(2),model.ticks+1,[sne.id])
            schedule_event!(a,CTL_Event(1),model.ticks+1,[-1*sne.id])
        end
    end
end

function init_condition!(sne::SimNE,model::ABM)
    #initialise condition: we generate entire time series to failure
    # TODO: Review
    #simulated sensor functions
    # funs = [
    #     (exp_f,(1.0,0.05),exp_ts,(),exp_c),
    #     (weibull_f,(1.0,1.0),wb_ts,(6.0),wb_c),
    #     (log_f,(50,0.1),log_ts,(),log_c)
    #     ]

    # ttf = 100
    # Δᵩ = 0.05
    # downtime = 2 #time steps
    # sne.condition_ts = generate_sensor_series(ttf,model.N,Δᵩ,0.05,downtime,funs)
    # vrul = generate_rul_series(ttf,Δᵩ,model.N,downtime)
    # sne.rul = vrul#reshape(vrul,length(vrul),1)
    Random.seed!(model.seed)
    #get_random(model.seed,1)
    # TODO check heterogenous assets with different expected rul
    state = get_state(sne)
    #randomly initialize condition of sne
    # state.rul = rand((sne.maintenance.eul-30):sne.maintenance.eul,nv(model.ntw_graph))[sne.id]

    # state.rul = rand(50:100,nv(model.ntw_graph))[sne.id]
    #TODO allow for different degradation profiles to be passed by user, e.g. time-dependent where the deterioration function changes at a given time-step
    state.rul = sne.id in model.init_sne_params.ids ? model.init_sne_params.ruls[first(indexin(sne.id,model.init_sne_params.ids))] : Int.(round.(rand(Uniform(50,100),nv(model.ntw_graph))))[sne.id]

    #set maitenance due time
    state.maintenance_due = model.ticks + state.rul
    set_state!(sne,state)
end