text:
using x265_jll;

mutable struct x265_picture
    pts::Clong
    dts::Clong
    userData::Ptr{Cvoid}
    planeR::Ptr{Cuchar}
    planeG::Ptr{Cuchar}
    planeB::Ptr{Cuchar}
    strideR::Cint
    strideG::Cint
    strideB::Cint
    bitDepth::Cint
end

x265_picture(picture::Ptr) = unsafe_load(Ptr{x265_picture}(picture))

function x265_apiver()
    @static if Sys.isapple()
        parts = split(x265_jll.get_libx265_path(), ".")
        return parts[length(parts)-1]
    end

    @static if Sys.islinux()
        parts = split(readlink(x265_jll.get_libx265_path()), ".")
        return last(parts)
    end

    @static if Sys.iswindows()
        error("Not implemented: don't know how to access a shared lib on Windows")
    end
end

param = C_NULL
encoder = C_NULL
picture = C_NULL

# x265 parameters
param = ccall((:x265_param_alloc, libx265), Ptr{Cvoid}, ())
println("typeof(param): ", typeof(param), "; value: $param")

if param == C_NULL
    error("NULL x265_param")
end

# default parameters
ccall(
    (:x265_param_default_preset, libx265),
    Cvoid,
    (Ptr{Cvoid}, Cstring, Cstring),
    param,
    "superfast",
    "zerolatency",
)

# set extra parameters manually
fps = Integer(5)
stat = ccall(
    (:x265_param_parse, libx265),
    Cint,
    (Ptr{Cvoid}, Cstring, Cstring),
    param,
    "fps",
    string(fps),
)
println("x265_param_parse::$stat")

# bRepeatHeaders = 1
stat = ccall(
    (:x265_param_parse, libx265),
    Cint,
    (Ptr{Cvoid}, Cstring, Ptr{Cvoid}),
    param,
    "repeat-headers",
    C_NULL,
)
#stat = ccall((:x265_param_parse, libx265), Cint, (Ptr{Cvoid}, Cstring, Cstring), param, "repeat-headers", "1")
println("x265_param_parse::$stat")

# internalCsp = X265_CSP_I444
stat = ccall(
    (:x265_param_parse, libx265),
    Cint,
    (Ptr{Cvoid}, Cstring, Cstring),
    param,
    "input-csp",
    "i444",
)
println("x265_param_parse::$stat")

# internalBitDepth = 8
# "input-depth", "output-depth" are CLI-ONLY
#stat = ccall((:x265_param_parse, libx265), Cint, (Ptr{Cvoid}, Cstring, Cstring), param, "output-depth", "8")
#println("x265_param_parse::$stat")

image_width = Integer(200)
image_height = Integer(175)
res = string(image_width) * "x" * string(image_height)
stat = ccall(
    (:x265_param_parse, libx265),
    Cint,
    (Ptr{Cvoid}, Cstring, Cstring),
    param,
    "input-res",
    res,
)
println("x265_param_parse::$stat")

# rc.bitrate = bitrate
# bitrate = Integer(1000)
# stat = ccall((:x265_param_parse, libx265), Cint, (Ptr{Cvoid}, Cstring, Cstring), param, "bitrate", string(bitrate))
# println("x265_param_parse::$stat")

# set constant quality rate
crf = Integer(28)
stat = ccall(
    (:x265_param_parse, libx265),
    Cint,
    (Ptr{Cvoid}, Cstring, Cstring),
    param,
    "crf",
    string(crf),
)
println("x265_param_parse::$stat")

# x265 encoder
const encoder_open = "x265_encoder_open_" * x265_apiver()
encoder = ccall((encoder_open, libx265), Ptr{Cvoid}, (Ptr{Cvoid},), param)
println("typeof(encoder): ", typeof(encoder), "; value: $encoder")

if encoder == C_NULL
    error("NULL x265_encoder")
end

# x265 picture
picture = ccall((:x265_picture_alloc, libx265), Ptr{Cvoid}, ())
println("typeof(picture): ", typeof(picture), "; value: $picture")

if picture == C_NULL
    error("NULL x265_picture")
end

ccall((:x265_picture_init, libx265), Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}), param, picture)

j_picture = x265_picture(picture)
display(j_picture)

# a sample 2D plane
width = 200
height = 175

planeB = fill(UInt8(128), width, height)
println("strides:", strides(planeB))

vectorB = Vector{UInt8}(undef, width*height)
vectorB .= 128

# pointer(planeB) gives the same answer as Base.unsafe_convert(Ptr{Cuchar}, Base.cconvert(Ptr{UInt8},Ref(planeB, 1)))

j_picture.planeR = pointer(vectorB)
j_picture.planeG = Base.unsafe_convert(Ptr{Cuchar}, Base.cconvert(Ptr{UInt8},Ref(planeB, 1)))
# j_picture.planeB = Base.unsafe_convert(Ptr{Cuchar}, Base.cconvert(Ptr{UInt8},Ref(vectorB, 1)))
# j_picture.planeB = Base.unsafe_convert(Base.cconvert((vectorB)))
j_picture.planeB = pointer(planeB)
j_picture.strideR = 0
j_picture.strideG = 0
j_picture.strideB = strides(planeB)[2]
j_picture.bitDepth = 8

# sync the Julia structure back to C
unsafe_store!(Ptr{x265_picture}(picture), j_picture)

# check if the structure is getting written to
j_picture2 = x265_picture(picture)
display(j_picture2)

# release memory
ccall((:x265_param_free, libx265), Cvoid, (Ptr{Cvoid},), param)
param = C_NULL

ccall((:x265_encoder_close, libx265), Cvoid, (Ptr{Cvoid},), encoder)
encoder = C_NULL

ccall((:x265_picture_free, libx265), Cvoid, (Ptr{Cvoid},), picture)
picture = C_NULL
;


text:
function VorticityBoundaryConditions!(ω,ψ,Δx,Δy,un,us,ve,vw)
  ω[:,end] .= 2*((ψ[:,end]  - ψ[:,end-1] )/(Δx^2) .- ve/Δx)
  ω[:,1] .= 2*((ψ[:,1]  - ψ[:,2]   )/(Δx^2) .- vw/Δx)
  ω[end,:] .= 2*((ψ[end,:]  - ψ[end-1,:] )/(Δy^2) .+ us/Δy)
  ω[1,:] .= 2*((ψ[1,:]  - ψ[2,:]   )/(Δy^2) .+ un/Δy)
end
function GaussSiedel!(ϕ,Ap,An,As,Ae,Aw,Rp,res; λ=1, maxiter = 1000)
  normRes = 1
  k = 0
  Ny,Nx = size(ϕ)
  while normRes >= 1e-8 && k < maxiter
    k += 1
    for i in 2:Ny-1
      for j in 2:Nx-1
        ϕP = ϕ[i,j]
        ϕE = ϕ[i+0,j+1]
        ϕW = ϕ[i+0,j-1]
        ϕN = ϕ[i-1,j+0]
        ϕS = ϕ[i+1,j+0]
        res[i,j] = Rp[i,j] - (Ap*ϕP
          + An*ϕN
          + As*ϕS
          + Ae*ϕE
          + Aw*ϕW)
        Δϕ = res[i,j]/Ap
        ϕ[i,j] = λ*(ϕ[i,j] + Δϕ) + (1-λ)*ϕ[i,j]
      end
    end
    normRes = norm(res)
  end
  return k
end
function LinearSolve!(A,x,b)
  # Solves the equation Ax = b assuming zero Dirichlet BCs everywhere
  Ny,Nx = size(b)
  Ny,Nx = Ny-2, Nx-2
  x_int = x[2:end-1,2:end-1]
  b_int = b[2:end-1,2:end-1]
  b_vec = reshape(b_int,Ny*Nx)
  # x_int = A\b_vec
  x_vec = reshape(x_int,Ny*Nx)
  cg!(x_vec,A,b_vec, log = true)
  x[2:end-1,2:end-1] .= reshape(x_int,(Ny,Nx))
end
function BuildPoissonMatrix(Ny,Nx,Δx,Δy)
  # This function returns a (Ny*Nx) × (Ny*Nx) matrix in the form of
  # a sparse array, corresponding to the discrete 2D Laplacian operator.
  Ny = Ny-2
  Nx = Nx-2

  Isx = [1:Ny; 1:Ny-1; 2:Ny]
  Jsx = [1:Ny; 2:Ny; 1:Ny-1]

  Isy = [1:Nx; 1:Nx-1; 2:Nx]
  Jsy = [1:Nx; 2:Nx; 1:Nx-1]

  Vsx = [fill(-2,Ny); fill(1, 2Ny-2)]
  Vsy = [fill(-2,Nx); fill(1, 2Nx-2)]
  D²x = sparse(Isx, Jsx, Vsx)
  D²y = sparse(Isy, Jsy, Vsy)
  # D_xx = 1/(Δx^2) .* kron(sparse(I,Nx,Nx), D²x)
  # D_yy = 1/(Δy^2) .* kron(D²y, sparse(I,Ny,Ny))
  D_yy = 1/(Δy^2) .* kron(sparse(I,Nx,Nx), D²x)
  D_xx = 1/(Δx^2) .* kron(D²y, sparse(I,Ny,Ny))
  Lap = D_xx + D_yy
end
function BuildAdvectionDiffusionCoefficients(Re,Δt,Δx,Δy)
  # Time-derivative
  ap = 1/Δt
  # Diffusion
  ap += 2/(Re*Δx^2) + 2/(Re*Δy^2)
  an = -1/(Re*Δy^2)
  aw = -1/(Re*Δx^2)
  as = -1/(Re*Δy^2)
  ae = -1/(Re*Δx^2)
  return ap,an,as,ae,aw
end
function BuildAdvectionDiffusionRHS!(Rp,ϕ,ψ,Δt,Δx,Δy,Ny,Nx,Re)
  # Time derivative
  Rp .= ϕ/Δt

  # Diffusion term (fully implicit)

  # Convection term
  for i in 2:Ny-1
    for j in 2:Nx-1
      ϕE = ϕ[i+0,j+1]; ϕW = ϕ[i+0,j-1]; ϕN = ϕ[i-1,j+0]; ϕS = ϕ[i+1,j+0]
      ψE = ψ[i+0,j+1]; ψW = ψ[i+0,j-1]; ψN = ψ[i-1,j+0]; ψS = ψ[i+1,j+0]

      u    = (ψN - ψS)/(2Δy); v    = -(ψE - ψW)/(2Δx)
      ∂ϕ∂y = (ϕN - ϕS)/(2Δy); ∂ϕ∂x = (ϕE - ϕW)/(2Δx)

      Rp[i,j] += - (u*∂ϕ∂x + v*∂ϕ∂y)
      # Rp[i,j] += (ψE - ψW)/(2Δx) * (ϕN - ϕS)/(2Δy) -
      #            (ψN - ψS)/(2Δy) * (ϕE - ϕW)/(2Δx)
    end
  end
end

function RecordHistory!(ϕ,ϕ_old,ϕ_hist)
  Δϕ = norm(ϕ - ϕ_old)
  ϕ_old .= ϕ
  push!(ϕ_hist,Δϕ)
  return(Δϕ)
end

ext: jl


text:
using StanSample

ProjDir = @__DIR__

gq = "
  data {
    int<lower=0> N;
    int<lower=0> y[N];
  }
  parameters {
    real<lower=0.00001> Theta;
  }
  model {
    Theta ~ gamma(4, 1000);
    for (n in 1:N)
      y[n] ~ exponential(Theta);
  }
  generated quantities{
    real y_pred;
    y_pred = exponential_rng(Theta);
  }
";

data = Dict(
  "N" => 3,
  "y" => [100, 950, 450]
);

stanmodel = SampleModel("Generate_quantities", gq);

rc = stan_sample(stanmodel; data,
  use_cpp_chains=true, check_num_chains=false,
  num_cpp_chains=2, num_julia_chains=2)

if success(rc)
  # Convert to an MCMCChains.Chains object
  chns = read_samples(stanmodel)

  # Show the same output in DataFrame format
  df = StanSample.read_summary(stanmodel)

  available_chains(stanmodel) |> display
  println()

  gq_df = stan_generate_quantities(stanmodel, 1, "2_3")
end


text:
function ImGui_ImplOpenGL3_SetupRenderState(ctx::Context, draw_data, fb_width::Cint, fb_height::Cint, vertex_array_object::GLuint)
    # setup render state:
    # - alpha-blending enabled
    # - no face culling
    # - no depth testing
    # - scissor enabled
    # - polygon fill
    glEnable(GL_BLEND)
    glBlendEquation(GL_FUNC_ADD)
    glBlendFuncSeparate(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA, GL_ONE, GL_ONE_MINUS_SRC_ALPHA)
    glDisable(GL_CULL_FACE)
    glDisable(GL_DEPTH_TEST)
    glDisable(GL_STENCIL_TEST)
    glEnable(GL_SCISSOR_TEST)

    ctx.GlslVersion > 310 && glDisable(GL_PRIMITIVE_RESTART)

    glPolygonMode(GL_FRONT_AND_BACK, GL_FILL)

    # FIXME: support for glClipControl

    # setup viewport, orthographic projection matrix
    glViewport(0, 0, GLsizei(fb_width), GLsizei(fb_height))
    disp_pos = unsafe_load(draw_data.DisplayPos)
    disp_size = unsafe_load(draw_data.DisplaySize)
    L = disp_pos.x
    R = disp_pos.x + disp_size.x
    T = disp_pos.y
    B = disp_pos.y + disp_size.y
    ortho_projection = GLfloat[2.0/(R-L), 0.0, 0.0, 0.0,
                               0.0, 2.0/(T-B), 0.0, 0.0,
                               0.0, 0.0, -1.0, 0.0,
                               (R+L)/(L-R), (T+B)/(B-T), 0.0, 1.0]

    glUseProgram(ctx.ShaderHandle)
    glUniform1i(ctx.AttribLocationTex, 0)
    glUniformMatrix4fv(ctx.AttribLocationProjMtx, 1, GL_FALSE, ortho_projection)
    ctx.GlslVersion > 330 && glBindSampler(0, 0)

    glBindVertexArray(vertex_array_object)
    glBindBuffer(GL_ARRAY_BUFFER, ctx.VboHandle)
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, ctx.ElementsHandle)
    glEnableVertexAttribArray(ctx.AttribLocationVtxPos)
    glEnableVertexAttribArray(ctx.AttribLocationVtxUV)
    glEnableVertexAttribArray(ctx.AttribLocationVtxColor)
    pos_offset = fieldoffset(ImDrawVert, 1)
    uv_offset = fieldoffset(ImDrawVert, 2)
    col_offset = fieldoffset(ImDrawVert, 3)
    glVertexAttribPointer(ctx.AttribLocationVtxPos[],   2, GL_FLOAT,         GL_FALSE, sizeof(ImDrawVert), Ptr{GLCvoid}(pos_offset))
    glVertexAttribPointer(ctx.AttribLocationVtxUV[],    2, GL_FLOAT,         GL_FALSE, sizeof(ImDrawVert), Ptr{GLCvoid}(uv_offset))
    glVertexAttribPointer(ctx.AttribLocationVtxColor[], 4, GL_UNSIGNED_BYTE, GL_TRUE,  sizeof(ImDrawVert), Ptr{GLCvoid}(col_offset))

    return nothing
end

function ImGui_ImplOpenGL3_RenderDrawData(ctx::Context, draw_data)
    # avoid rendering when minimized, scale coordinates for retina displays
    fb_width = trunc(Cint, unsafe_load(draw_data.DisplaySize.x) * unsafe_load(draw_data.FramebufferScale.x))
    fb_height = trunc(Cint, unsafe_load(draw_data.DisplaySize.y) * unsafe_load(draw_data.FramebufferScale.y))
    (fb_width ≤ 0 || fb_height ≤ 0) && return nothing

    # backup GL state
    last_active_texture = GLint(0); @c glGetIntegerv(GL_ACTIVE_TEXTURE, &last_active_texture)
    glActiveTexture(GL_TEXTURE0)
    last_program = GLint(0); @c glGetIntegerv(GL_CURRENT_PROGRAM, &last_program)
    last_texture = GLint(0); @c glGetIntegerv(GL_TEXTURE_BINDING_2D, &last_texture)
    last_sampler = GLint(0); @c glGetIntegerv(GL_SAMPLER_BINDING, &last_sampler)
    last_array_buffer = GLint(0)
    ctx.GlslVersion ≥ 330 && @c glGetIntegerv(GL_ARRAY_BUFFER_BINDING, &last_array_buffer)
    last_vertex_array_object = GLint(0); @c glGetIntegerv(GL_VERTEX_ARRAY_BINDING, &last_vertex_array_object)
    last_polygon_mode = GLint[0,0]; glGetIntegerv(GL_POLYGON_MODE, last_polygon_mode)
    last_viewport = GLint[0,0,0,0]; glGetIntegerv(GL_VIEWPORT, last_viewport)
    last_scissor_box = GLint[0,0,0,0]; glGetIntegerv(GL_SCISSOR_BOX, last_scissor_box)
    last_blend_src_rgb = GLint(0); @c glGetIntegerv(GL_BLEND_SRC_RGB, &last_blend_src_rgb)
    last_blend_dst_rgb = GLint(0); @c glGetIntegerv(GL_BLEND_DST_RGB, &last_blend_dst_rgb)
    last_blend_src_alpha = GLint(0); @c glGetIntegerv(GL_BLEND_SRC_ALPHA, &last_blend_src_alpha)
    last_blend_dst_alpha = GLint(0); @c glGetIntegerv(GL_BLEND_DST_ALPHA, &last_blend_dst_alpha)
    last_blend_equation_rgb = GLint(0); @c glGetIntegerv(GL_BLEND_EQUATION_RGB, &last_blend_equation_rgb)
    last_blend_equation_alpha = GLint(0); @c glGetIntegerv(GL_BLEND_EQUATION_ALPHA, &last_blend_equation_alpha)
    last_enable_blend = glIsEnabled(GL_BLEND)
    last_enable_cull_face = glIsEnabled(GL_CULL_FACE)
    last_enable_depth_test = glIsEnabled(GL_DEPTH_TEST)
    last_enable_stencil_test = glIsEnabled(GL_STENCIL_TEST)
    last_enable_scissor_test = glIsEnabled(GL_SCISSOR_TEST)
    last_enable_primitive_restart = ctx.GlslVersion ≥ 310 ? glIsEnabled(GL_PRIMITIVE_RESTART) : GL_FALSE

    vertex_array_object = GLuint(0)
    @c glGenVertexArrays(1, &vertex_array_object)
    ImGui_ImplOpenGL3_SetupRenderState(ctx, draw_data, fb_width, fb_height, vertex_array_object)

    # will project scissor/clipping rectangles into framebuffer space
    clip_off = unsafe_load(draw_data.DisplayPos)         # (0,0) unless using multi-viewports
    clip_scale = unsafe_load(draw_data.FramebufferScale) # (1,1) unless using retina display which are often (2,2)

    # render command lists
    data = unsafe_load(draw_data)
    cmd_lists = unsafe_wrap(Vector{Ptr{ImDrawList}}, data.CmdLists, data.CmdListsCount)
    for cmd_list in cmd_lists
        vtx_buffer = cmd_list.VtxBuffer |> unsafe_load
        idx_buffer = cmd_list.IdxBuffer |> unsafe_load
        # glBindBuffer(GL_ARRAY_BUFFER, g_VboHandle[])
        glBufferData(GL_ARRAY_BUFFER, vtx_buffer.Size * sizeof(ImDrawVert), Ptr{GLCvoid}(vtx_buffer.Data), GL_STREAM_DRAW)

        # glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, g_ElementsHandle[])
        glBufferData(GL_ELEMENT_ARRAY_BUFFER, idx_buffer.Size * sizeof(ImDrawIdx), Ptr{GLCvoid}(idx_buffer.Data), GL_STREAM_DRAW)

        cmd_buffer = cmd_list.CmdBuffer |> unsafe_load
        for cmd_i = 0:cmd_buffer.Size-1
            pcmd = cmd_buffer.Data + cmd_i * sizeof(ImDrawCmd)
            elem_count = unsafe_load(pcmd.ElemCount)
            cb_funcptr = unsafe_load(pcmd.UserCallback)
            if cb_funcptr != C_NULL
                # user callback (registered via ImDrawList_AddCallback)
                # ImDrawCallback_ResetRenderState is a special callback value used by the user to request the renderer to reset render state.
                if cb_funcptr == ctx.ImDrawCallback_ResetRenderState
                    ImGui_ImplOpenGL3_SetupRenderState(ctx, draw_data, fb_width, fb_height, vertex_array_object);
                else
                    ccall(cb_funcptr, Cvoid, (Ptr{ImDrawList}, Ptr{ImDrawCmd}), cmd_list, pcmd)
                end
            else
                # project scissor/clipping rectangles into framebuffer space
                rect = unsafe_load(pcmd.ClipRect)
                clip_rect_x = (rect.x - clip_off.x) * clip_scale.x
                clip_rect_y = (rect.y - clip_off.y) * clip_scale.y
                clip_rect_z = (rect.z - clip_off.x) * clip_scale.x
                clip_rect_w = (rect.w - clip_off.y) * clip_scale.y
                if clip_rect_x < fb_width && clip_rect_y < fb_height && clip_rect_z ≥ 0 && clip_rect_w ≥ 0
                    # apply scissor/clipping rectangle
                    ix = trunc(Cint, clip_rect_x)
                    iy = trunc(Cint, fb_height - clip_rect_w)
                    iz = trunc(Cint, clip_rect_z - clip_rect_x)
                    iw = trunc(Cint, clip_rect_w - clip_rect_y)
                    glScissor(ix, iy, iz, iw)
                    # bind texture, draw
                    glBindTexture(GL_TEXTURE_2D, UInt(unsafe_load(pcmd.TextureId)))
                    format = sizeof(ImDrawIdx) == 2 ? GL_UNSIGNED_SHORT : GL_UNSIGNED_INT
                    glDrawElementsBaseVertex(GL_TRIANGLES, GLsizei(elem_count), format, Ptr{Cvoid}(unsafe_load(pcmd.IdxOffset) * sizeof(ImDrawIdx)), GLint(unsafe_load(pcmd.VtxOffset)))
                end
            end
        end
    end
    @c glDeleteVertexArrays(1, &vertex_array_object)

    # restore modified GL state
    glUseProgram(last_program)
    glBindTexture(GL_TEXTURE_2D, last_texture)
    ctx.GlslVersion > 330 && glBindSampler(0, last_sampler)
    glActiveTexture(last_active_texture)
    glBindVertexArray(last_vertex_array_object)
    glBindBuffer(GL_ARRAY_BUFFER, last_array_buffer)
    glBlendEquationSeparate(last_blend_equation_rgb, last_blend_equation_alpha)
    glBlendFuncSeparate(last_blend_src_rgb, last_blend_dst_rgb, last_blend_src_alpha, last_blend_dst_alpha)
    last_enable_blend ? glEnable(GL_BLEND) : glDisable(GL_BLEND)
    last_enable_cull_face ? glEnable(GL_CULL_FACE) : glDisable(GL_CULL_FACE)
    last_enable_depth_test ? glEnable(GL_DEPTH_TEST) : glDisable(GL_DEPTH_TEST)
    last_enable_stencil_test ? glEnable(GL_STENCIL_TEST) : glDisable(GL_STENCIL_TEST)
    last_enable_scissor_test ? glEnable(GL_SCISSOR_TEST) : glDisable(GL_SCISSOR_TEST)
    if ctx.GlslVersion ≥ 310
        last_enable_primitive_restart ? glEnable(GL_PRIMITIVE_RESTART) : glDisable(GL_PRIMITIVE_RESTART)
    end

    glPolygonMode(GL_FRONT_AND_BACK, GLenum(last_polygon_mode[1]))
    glViewport(last_viewport[1], last_viewport[2], GLsizei(last_viewport[3]), GLsizei(last_viewport[4]))
    glScissor(last_scissor_box[1], last_scissor_box[2], GLsizei(last_scissor_box[3]), GLsizei(last_scissor_box[4]))

    return nothing
end


text:
"""
    MDModel <: AbstractFDDObject

Type for component synthesis models employed to solve model detection problems.

If `sysc::MDModel` is a component synthesis model object, the underlying descriptor system model
can be obtained via `sysc.sys` and the dimensions of the control, disturbance, noise and auxiliary
vectors are contained in the integers `sysm.mu`, `sysm.md`, `sysm.mw` and `sysm.ma`, respectively.
"""
struct MDModel <: AbstractFDDObject
    sys::DescriptorStateSpace
    mu::Int
    md::Int
    mw::Int
    ma::Int
    function MDModel(sys::DescriptorStateSpace, mu::Int, md::Int, mw::Int, ma::Int)
        mdmodel_validation(sys, mu, md, mw, ma)
        new(sys, mu, md, mw, ma)
    end
end
"""
    MDModel(sys; mu, md, mw, ma) -> sysc::MDModel

Build for a linear time-invariant descriptor system `sys = (A-λE,B,C,D)`
a component synthesis model object `sysm::MDModel`
to be used in conjunction with the analysis and synthesis functions of model detection filters.

The resulting synthesis model object `sysc` contains the component model with partitioned inputs,
`sysc.sys = (A-λE,[Bu Bd Bw Bv],C,[Du Dd Dw Dv])`, where
`Bu`, `Bd`, `Bw` and `Bv` are formed from the successive columns of `B` and are
 the input matrices from the control inputs `u`, disturbance inputs `d`,
noise inputs `w` and auxiliary inputs `v`, respectively,
and `Du`, `Dd`, `Dw` and `Dv` are formed from the successive columns of `D` and
are the feedthrough matrices from those inputs.
The dimensions of control, disturbance, noise and auxiliary input vectors are contained in
`sysm.mu`, `sysm.md`, `sysm.mw` and `sysm.ma`, respectively.

The information on the partition of the input components
in control, disturbance, noise and auxiliary inputs can be specified using the following keyword arguments:

`mu = nu` specifies the dimension `nu` of the control input vector `u` (default: `nu = 0`)

`md = nd` specifies the dimension `nd` of the disturbance input vector `d` (default: `nd = 0`)

`mw = nw` specifies the dimension `nw` of the noise input vector `w` (default: `nw = 0`)

`ma = na` specifies the dimension `na` of the auxiliary input vector `v` (default: `na = 0`)
"""
function MDModel(sys::DescriptorStateSpace; mu::Int = 0, md::Int = 0, mw::Int = 0, ma::Int = 0)
    MDModel(sys, mu, md, mw, ma)
end
function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, sysc::MDModel)
    summary(io, sysc); println(io)
    display(sysc.sys)
    mu = sysc.mu
    md = sysc.md
    mw = sysc.mw
    maux = sysc.ma
    if mu+md+mw+maux > 0
       println(io, "Input groups:")
       println(io, "Name         Channels")
       mu > 0 &&   println(io,"controls     $(1:mu)")
       md > 0 &&   println(io,"disturbances $(mu+1:mu+md)")
       mw > 0 &&   println(io,"noise        $(mu+md+1:mu+md+mw)")
       maux > 0 && println(io,"aux          $(mu+md+mw+1:mu+md+mw+maux)")
    end
end
"""
    MDMModel <: AbstractFDDObject

Type for multiple synthesis models employed to solve model detection problems.

If `sysm::MDMModel` is the multiple synthesis model object, the underlying vector of descriptor system models
can be obtained via `sysm.sys`, the common dimension of the control vectors is contained in the integer
`sysm.mu`, and the dimensions of the disturbance, noise
and auxiliary input vectors for the `i`-th model `sysm.sys[i]` are contained in the `i`-th
components of the integer vectors `sysm.md`, `sysm.mw` and `sysm.ma`, respectively.
"""
struct MDMModel <: AbstractFDDObject
    sys::Vector{<:DescriptorStateSpace}
    mu::Int
    md::Vector{Int}
    mw::Vector{Int}
    ma::Vector{Int}
    function MDMModel(sys::Vector{<:DescriptorStateSpace}, mu::Int, md::Vector{Int},
                     mw::Vector{Int}, ma::Vector{Int})
        mdmodel_validation(sys, mu, md, mw, ma)
        new(sys, mu, md, mw, ma)
    end
end
"""
    MDMModel(sys; mu, md, mw, ma) -> sysm::MDMModel

Build for a vector of linear time-invariant system models `sys`, with the `i`-th component model
`sys[i] = (Ai-λEi,Bi,Ci,Di)`, a multiple synthesis model object `sysm::MDMModel`
to be used in conjunction with the analysis and synthesis functions of model detection filters.

The resulting multiple synthesis model object `sysm` contains the vector `sysm.sys`
of component models with partitioned inputs,
with the `i`-th model `sysm.sys[i] = (Ai-λEi,[Bui Bdi Bwi Bvi],Ci,[Dui Ddi Dwi Dvi])`, where
`Bui`, `Bdi`, `Bwi` and `Bvi` are formed from the successive columns of `Bi` and are
 the input matrices from the control inputs `u`, disturbance inputs `di`,
noise inputs `wi` and auxiliary inputs `vi`, respectively,
and `Dui`, `Ddi`, `Dwi` and `Dvi` are formed from the successive columns of `Di` and
are the feedthrough matrices from those inputs.
The dimensions of control, disturbance, noise and auxiliary input vectors are contained in
`sysm.mu`, `sysm.md[i]`, `sysm.mw[i]` and `sysm.ma[i]`, respectively.

If `N` is the number of component models, the information on the partition of the input components
in control, disturbance, noise and auxiliary inputs can be specified using the following keyword arguments:

`mu = nu` specifies the dimension `nu` of control input vector `u` (default: `nu = 0`)

`md = nd` specifies the vector `nd` containing the `N` dimensions of disturbance input vectors, such that
          the `i`-th disturbance vector `di` has dimension `nd[i]` (default: `nd = zeros(Int,N)`)

`mw = nw` specifies the vector `nw` containing the `N` dimensions of noise input vectors, such that
          the `i`-th noise vector `wi` has dimension `nw[i]` (default: `nw = zeros(Int,N)`)

`ma = na` specifies the vector `na` containing the `N` dimensions of auxiliary input vectors, such that
          the `i`-th auxiliary vector `vi` has dimension `na[i]` (default: `na = zeros(Int,N)`)
"""
function MDMModel(sys::Vector{<:DescriptorStateSpace}; mu::Int = 0,
                  md::Vector{Int} = fill(0,length(sys)),
                  mw::Vector{Int} = fill(0,length(sys)),
                  ma::Vector{Int} = fill(0,length(sys)))
    MDMModel(sys, mu, md, mw, ma)
end
function mdmodel_validation(sys::DST, mu::Int, md::Int, mw::Int, ma::Int) where {DST <: DescriptorStateSpace}
    m = size(sys,2)
    mu > m && error("number of control inputs $mu exceeds the number of system inputs $m")
    mu < 0 && error("number of control inputs must be nonnegative")
    md < 0 && error("number of disturbance inputs must be nonnegative")
    mw < 0 && error("number of noise inputs must be nonnegative")
    ma < 0 && error("number of auxiliary inputs must be nonnegative")
    mu+md+mw+ma > m && error("the specified total number of inputs exceeds the number of system inputs")
end

function mdmodel_validation(sys::Vector{<:DescriptorStateSpace}, mu::Int, md::Vector{Int},
                            mw::Vector{Int}, ma::Vector{Int})
    N = length(sys)
    N == length(md) || error("improper length of disturbance input vector dimensions")
    N == length(mw) || error("improper length of noise input vector dimensions")
    N == length(ma) || error("improper length of auxliary input vector dimensions")
    m = minimum(size.(sys,2))
    mu > m && error("number of control inputs $mu exceeds the minimum number of system inputs $m")
    mu < 0 && error("number of control inputs must be nonnegative")
    any(md .< 0) && error("number of disturbance inputs must be nonnegative")
    any(mw .< 0) && error("number of noise inputs must be nonnegative")
    any(ma .< 0) && error("number of auxiliary inputs must be nonnegative")
    mi = size.(sys,2)
    any(md+mw+ma-mi .> -mu) && error("specified number of inputs exceeds the number of system inputs")
end
function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, sysm::MDMModel)
    summary(io, sysm); println(io)
    mu = sysm.mu
    for i = 1:length(sysm.sys)
        println(io, "\nModel #$i:")
        display(sysm.sys[i])
        mdi = sysm.md[i]
        mwi = sysm.mw[i]
        mauxi = sysm.ma[i]
        if mu+mdi+mwi+mauxi > 0
            println(io, "Input groups:")
            println(io, "Name         Channels")
            mu > 0 &&    println(io,"controls     $(1:mu)")
            mdi > 0 &&   println(io,"disturbances $(mu+1:mu+mdi)")
            mwi > 0 &&   println(io,"noise        $(mu+mdi+1:mu+mdi+mwi)")
            mauxi > 0 && println(io,"aux          $(mu+mdi+mwi+1:mu+mdi+mwi+mauxi)")
         end

    end
end
"""
    mdmodset(sys; controls, c, disturbances, d, noise, n, aux) -> sysm::MDMModel

Build for a vector of linear time-invariant system models `sys`, with the `i`-th component model
`sys[i] = (Ai-λEi,Bi,Ci,Di)`, a multiple synthesis model object `sysm::MDMModel`
to be used in conjunction with the analysis and synthesis functions of model detection filters.

The information on the partition of the input components in control, disturbance, noise and auxiliary
inputs can be specified using the following keyword arguments:

`controls = inpu` or `c = inpu` specifies the indices `inpu` of the control inputs (default: void)

`disturbances = inpd` or `d = inpd` specifies the indices or a vector of indices `inpd` of the disturbance inputs (default: void)

`noise = inpn` or  `noise = inpn` specifies the indices or a vector of indices `inpn` of the noise inputs (default: void)

`aux = inpa` specifies the indices or a vector of indices `inpa` of the auxiliary inputs (default: void)

The indices of inputs can be specified as integer vectors, integer scalars or integer `UnitRange`s.
For disturbance, noise and auxiliary inputs, vectors of integer vectors or
vectors of integer `UnitRange`s can be used to specify possibly different set of indices for each component model.

The resulting `sysm` contains the vector `sysm.sys` of partitioned systems, with the `i`-th model
`sysm.sys[i] = (Ai-λEi,[Bui Bdi Bwi Bvi],Ci,[Dui Ddi Dwi Dvi])`, where
`Bui`, `Bdi`, `Bwi` and `Bvi` are the input matrices from the control inputs `u`, disturbance inputs `di`,
noise inputs `wi` and auxiliary inputs `vi`, respectively, and `Dui`, `Ddi`, `Dwi` and `Dvi` are the feedthrough matrices from those inputs.
The dimensions of control, disturbance, fault, noise and auxiliary input vectors are contained in
`sysm.mu`, `sysm.md[i]`, `sysm.mw[i]` and `sysm.ma[i]`, respectively.

_Method:_ If `Gi(λ)` is the `p x mi` transfer function matrix of `sys[i]`, then the resulting component system `sysm.sys[i]` has an
equivalent input output form `[Gui(λ) Gdi(λ) Gwi(λ) Gvi(λ)]`, where the following relations define the component matrices:
`Gui(λ) = Gi(λ)*Su`, `Gdi(λ) = Gi(λ)*Sdi`,  `Gwi(λ) = G(λ)*Swi`,  `Gvi(λ) = Gi(λ)*Svi`, with
the selection matrices `Su`, `Sdi`, `Swi` and `Svi` formed from the columns of the `mi`-th order identity matrix
(note that `Su` is the same for all component models).
"""
function mdmodset(sys::Vector{<:DescriptorStateSpace};
                   controls::VRS = Int[], c::VRS = Int[],
                   disturbances::Union{VRS,Vector{Vector{Int}},Vector{UnitRange{Int}}} = fill(Int[],length(sys)),
                   d::Union{VRS,Vector{Vector{Int}},Vector{UnitRange{Int}}} = fill(Int[],length(sys)),
                   noise::Union{VRS,Vector{Vector{Int}},Vector{UnitRange{Int}}} = fill(Int[],length(sys)),
                   n::Union{VRS,Vector{Vector{Int}},Vector{UnitRange{Int}}} = fill(Int[],length(sys)),
                   aux::Union{VRS,Vector{Vector{Int}},Vector{UnitRange{Int}}} = fill(Int[],length(sys)))
    N = length(sys)
    if N > 0
       p = size(sys[1],1)
       Ts = sys[1].Ts
    end
    m = maximum(size.(sys,2))
    sysn = similar(sys,N)
    inpu = sort(unique([controls;c]))
    mu = length(inpu)
    mu == 0 || inpu[mu] <= m || error("selected index/indices of control inputs larger than the minimum number of system inputs $m")
    vecdist = isa(disturbances,Vector)
    vecd = isa(d,Vector)
    vecnoise = isa(noise,Vector)
    vecn = isa(n,Vector)
    vecaux = isa(aux,Vector)
    vecdist && length(disturbances) != N && error("the vector of disturbance input indices must have length $N")
    vecd && length(disturbances) != N && error("the vector of disturbance input indices must have length $N")
    vecnoise && length(noise) != N && error("the vector of noise input indices must have length $N")
    vecn && length(n) != N && error("the vector of noise input indices must have length $N")
    vecaux && length(aux) != N && error("the vector of auxiliary input indices must have length $N")
    md = similar(Vector{Int},N)
    mw = similar(Vector{Int},N)
    ma = similar(Vector{Int},N)
    for i = 1:N
        i == 1 || (Ts = DescriptorSystems.promote_Ts(Ts,sys[i].Ts))
    end
    for i = 1:N
        i == 1 || p == size(sys[i],1) || error("all systems must have the same number of outputs")
        if vecdist && vecd
            inpdi = unique([disturbances[i];d[i]])
        elseif vecn
            inpdi = unique([disturbances;d[i]])
        elseif vecnooise
            inpdi = unique([disturbances[i];d])
        else
            inpdi = unique([disturbances;d])
        end
        isempty(intersect(inpu,inpdi)) || error("control and disturbance inputs must be distinct")
        if vecnoise && vecn
            inpwi = unique([noise[i];n[i]])
        elseif vecn
            inpwi = unique([noise;n[i]])
        elseif vecnooise
            inpwi = unique([noise[i];n])
        else
            inpwi = unique([noise;n])
        end
        inpai = vecaux ? unique([aux[i];Int[]]) : unique([aux;Int[]])
        mdi = length(inpdi)
        mwi = length(inpwi)
        mai = length(inpai)
        mdi == 0 || maximum(inpdi) <= m || error("selected index/indices of disturbance inputs larger than the minimum number of system inputs $m")
        mwi == 0 || maximum(inpwi) <= m || error("selected index/indices of noise inputs larger than the minimum number of system inputs $m")
        mai == 0 || maximum(inpauxi) <= m || error("selected index/indices of auxiliary inputs larger than the minimum number of system inputs $m")
        sysn[i] = dss(dssdata(sys[i][:,[inpu;inpdi;inpwi;inpai]])...;Ts)
        md[i] = mdi
        mw[i] = mwi
        ma[i] = mai
    end
    return MDMModel(sysn; mu, md, mw, ma)
end
"""
    mdmodset(sysc::Vector{<:MDModel}) -> sysm::MDMModel

Build for a vector of component synthesis models `sysc`, a multiple synthesis model object `sysm::MDMModel`
to be used in conjunction with the analysis and synthesis functions of model detection filters.
"""
function mdmodset(sysc::Vector{<:MDModel})
    N = length(sysc)
    sysn = similar(Vector{DescriptorStateSpace},N)
    md = similar(Vector{Int},N)
    mw = similar(Vector{Int},N)
    ma = similar(Vector{Int},N)
    if N > 0
       p = size(sysc[1].sys,1)
       mu = sysc[1].mu
       Ts = sysc[1].sys.Ts
       sysn[1] = sysc[1].sys
       md[1] = sysc[1].md
       mw[1] = sysc[1].mw
       ma[1] = sysc[1].ma
    end
    for i = 2:N
        p == size(sysc[i].sys,1) || error("all component models must have the same number of outputs")
        mu == sysc[i].mu || error("all component models must have the same number of control inputs")
        Ts = DescriptorSystems.promote_Ts(Ts,sysc[i].sys.Ts)
        sysn[i] = sysc[i].sys
        md[i] = sysc[i].md
        mw[i] = sysc[i].mw
        ma[i] = sysc[i].ma
    end
    return MDMModel(sysn; mu, md, mw, ma)
end

"""
    MDFilter <: AbstractFDDObject

Type for model detection filters resulted as solutions of model detection problems.

If `filter::MDFilter` is the model detection filter object,
the underlying `i`-th descriptor system model
can be obtained via `filter.sys[i]` and the dimensions of the partitioned filter input vectors as
`outputs` and `controls`,
can be accessed as the integers contained in `filter.ny` and `filter.mu`, respectively.
"""
struct MDFilter{T} <: AbstractFDDObject where T
    sys::Vector{DescriptorStateSpace{T}}
    ny::Int
    mu::Int
    function MDFilter{T}(sys::Vector{DescriptorStateSpace{T}}, ny::Int, mu::Int) where T
        N = length(sys)
        sysn = similar(sys,N)

        mf = ny+mu
        mf > size(sys[1],2) && error("total number of inputs exceeds the number of filter inputs $m")

        sysn[1] = sys[1][:,1:mf]
        for i = 1:N
            mf > size(sys[i],2) && error("total number of inputs exceeds the number of filter inputs $m")
            sysn[i] = sys[i][:,1:mf]
        end
        new{T}(sysn, ny, mu)
    end
end
# """
#     MDFilter(sys, p, mu) -> Q::MDFilter

# Build for a vector of linear time-invariant descriptor system models `sys[i] = (Ai-λEi,Bi,Ci,Di)`
# with the same number of inputs, a model detection filter object `Q`,
# as determined with the synthesis functions of model detection filters.
# `p` and `mu` are the number of measured outputs and the number of control inputs, respectively.
# It is assumed that each `Bi = [Byi Bui Bvi]` and `Di = [Dyi Dui Dvi]` are partitioned matrices such that
# `Byi` and `Dyi` have `p` columns, and `Bui` and `Dui` have `mu` columns,
# where `Byi` and `Bui` are the input matrices from the measured outputs `y` and
# control inputs `u`, `Dyi` and `Dui` are the feedthrough matrices from the measured outputs `y` and
# control inputs `u`.

# The resulting `Q` contains the vector of partitioned systems
# `Q.sys[i] = (Ai-λEi,[Byi Bdi],Ci,[Dyi Dui])` and the indices of inputs corresponding
# to the measured outputs and control inputs are contained in the associated
# integer vectors `Q.outputs` and `Q.controls`.
# """
# function MDFilter(sys::Vector{DescriptorStateSpace{T}}, p::Int, mu::Int) where T
#     m = size(sys[1],2)
#     p < 0 && error("number of measured outputs must be non-negative")
#     mu < 0 && error("number of control inputs must be non-negative")
#     p+mu > m && error("number of measured outputs and control inputs exceeds the number of filter inputs $m")
#     return MDFilter{T}(sys, Vector(1:p), Vector(p+1:p+mu))
# end
# function MDFilter(sys::Vector{DescriptorStateSpace}, p::Int, mu::Int)
#     m = size(sys[1],2)
#     p < 0 && error("number of measured outputs must be non-negative")
#     mu < 0 && error("number of control inputs must be non-negative")
#     p+mu > m && error("number of measured outputs and control inputs exceeds the number of filter inputs $m")
#     return MDFilter{eltype(sys[1])}(sys, Vector(1:p), Vector(p+1:p+mu))
# end

# function MDFilter(sys::Vector{DescriptorStateSpace{T}}; outputs::VRS = Int[], controls::VRS = Int[]) where T
#     return MDFilter{T}(sys, vec([outputs; Int[]]), vec([controls; Int[]]))
# end
function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, filter::MDFilter)
    summary(io, filter); println(io)
    for i = 1:length(filter.sys)
        println(io, "Filter #$i:")
        display(filter.sys[i])
    end
    p = filter.ny
    mu = filter.mu
    if p+mu > 0
       println(io, "Input groups:")
       println(io,          "Name      Channels")
       p > 0 &&  println(io,"outputs   $(1:p)")
       mu > 0 && println(io,"controls  $(p+1:p+mu)")
   end
end

"""
    MDFilterIF <: AbstractFDDObject

Type for the internal form of model detection filters resulted as solutions of model detection problems.

If `filter::MDFilterIF` is the model detection filter internal form object,
the underlying `(i,j)`-th descriptor system model
can be obtained via `filter.sys[i,j]` and the corresponding
dimensions of control, disturbance, fault, noise and auxiliary input vectors are contained in
`sysm.mu`, `sysm.md[j]`, `sysm.mw[j]` and `sysm.ma[j]`, respectively.
"""
struct MDFilterIF{T} <: AbstractFDDObject where T
    sys::Matrix{DescriptorStateSpace{T}}
    mu::Int
    md::Vector{Int}
    mw::Vector{Int}
    ma::Vector{Int}
    function MDFilterIF{T}(sys::Matrix{DescriptorStateSpace{T}}, mu::Int, md::Vector{Int},
                           mw::Vector{Int}, ma::Vector{Int}) where T
        M, N = size(sys)
        sysn = similar(sys,M,N)
        for j = 1:N
            mj = mu+md[j]+mw[j]+ma[j]
            mj > size(sys[1,j],2) && error("number of inputs exceeds the number of system inputs $m")
            sysn[1,j] = sys[1,j][:,1:mj]
            for i = 2:M
                mj > size(sys[i,j],2) && error("number of inputs exceeds the number of system inputs $m")
                sysn[i,j] = sys[i,j][:,1:mj]
            end
        end
        new{T}(sysn, mu, md, mw, ma)
    end
end
# function mdfilterIF_validation(sys::DescriptorStateSpace{T}, mu::Int, md::Int, mw::Int, ma::Int) where T
#     m = size(sys,2)

#     inpu = unique(controls);
#     inpd = unique(disturbances);
#     inpw = unique(noise);
#     inpaux = unique(aux);

#     isempty(intersect(inpu,inpd)) || error("control and disturbance inputs must be distinct")
#     isempty(intersect(inpd,inpw)) || error("disturbance and noise inputs must be distinct")
#     isempty(intersect(inpw,inpaux)) || error("noise and aux inputs must be distinct")

#     mu+md+mw+ma > m && error("number of inputs exceeds the number of system inputs $m")

#     return inpu, inpd, inpw, inpaux
# end
function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, filter::MDFilterIF)
    summary(io, filter); println(io)
    M, N = size(filter.sys)
    mu = filter.mu
    for j = 1:N
        md = filter.md[j]
        mw = filter.mw[j]
        maux = filter.ma[j]
        for i = 1:M
            println(io, "\n(Filter,Model) = (#$i,#$j):")
            display(filter.sys[i,j])
            if mu+md+mw+maux > 0
               println(io, "Input groups:")
               println(io, "Name         Channels")
               mu > 0 && println(io,"controls     $(1:mu)")
               md > 0 && println(io,"disturbances $(mu+1:mu+md)")
               mw > 0 && println(io,"noise        $(mu+md+1:mu+md+mw)")
               maux > 0 && println(io,"aux         $(mu+md+mw+1:mu+md+mw+maux)")
            end
        end
    end
end
"""
    mdIFeval(sysQ::MDFilter, sysm::MDMModel; minimal = false, atol, atol1 = atol, atol2 = atol, rtol, fast = true) -> sysR::MDFilterIF

Compute the internal form `sysR` of the model detection filter `sysQ` applied to the multiple synthesis model `sysm`.
If the `j`-th component model `sysm.sys[j]` has the partitioned transfer function matrix `Gj(λ) = [ Guj(λ)  Gdj(λ) Gwj(λ) Gvj(λ) ]` in accordance with
the partitioned system inputs as `controls`, `disturbances`, `noise` and `auxiliary` inputs, respectively,
and `Qi(λ) = [ Qyi(λ) Qui(λ) ]` is the partitioned transfer function matrix of the `i`-th filter `sysQ.sys[i]`
in accordance with the partitioned filter inputs as `outputs` and `controls`, then
the transfer function matrix `Rij(λ)` corresponding to the (`i`-th filter,`j`-th model) pair in the resulting internal form `sysR.sys[i,j]` is given by

     Rij(λ) = | Qyi(λ)  Qui(λ) | * | Guj(λ)  Gdj(λ) Gwj(λ) Gvj(λ) |
                                   |  I      0      0      0      |

Minimal descriptor realizations are computed if `minimal = true` and a (possibly) non-minimal
realization is determined if `minimal = false` (default).

The minimal realization computation relies on pencil manipulation algorithms which
employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

If `(Arij-λErij,Brij,Crij,Drij)` is the full order descriptor realization of `sysR.sys[i,j]`, then
the keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of matrices `Arij`, `Brij`, `Crij`, `Drij`, the absolute tolerance for the nonzero elements of `Erij`,
and the relative tolerance for the nonzero elements of `Arij`, `Brij`, `Crij`, `Drij` and `Eirj`.
The default relative tolerance is `nij*ϵ`, where `ϵ` is the working _machine epsilon_
and `nij` is the order of the system `sysR.sys[i,j]`.
The keyword argument `atol` can be used to simultaneously set `atol1 = atol` and `atol2 = atol`.

"""
function mdIFeval(Q::MDFilter, sysm::MDMModel; minimal::Bool = false,
                  atol::Real = 0., atol1::Real = atol, atol2::Real = atol, rtol::Real = 0., fast::Bool = true)
   # compute the internal form as
   #   R[i,j] = Q[i] * [ Gu[j] Gd[j] Gw[j] Gaux[j]]
   #                   [ I     0     0     0      ]
   M = length(Q.sys)
   N = length(sysm.sys)
   T1 = eltype(Q.sys[1])
   sysR = similar(Matrix{DescriptorStateSpace{T1}},M,N)
   p = size(sysm.sys[1],1)
   mu = sysm.mu
   (Q.ny == p && Q.mu == mu) || error("filter Q is incompatible with the given multiple model")
   for j = 1:N
       m = size(sysm.sys[j],2)
       sysej = [sysm.sys[j];  I zeros(mu,m-mu)]
       for i = 1:M
           sysR[i,j] = minimal ? gminreal(Q.sys[i]*sysej; atol1, atol2, rtol, fast) : Q.sys[i]*sysej
       end
   end
   return MDFilterIF{T1}(sysR, sysm.mu, sysm.md, sysm.mw, sysm.ma)
end

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
#in matlab:
#x=rand(512,1);
#[s,f,t,p]=spectrogram(x,ones(1,256),128,256,10);
#save
#
#in julia:
#using MAT
#
#matdata=matread("matlab.mat")
#
#for i in ("x", "f", "t", "p")
#  fid=open("spectrogram_$i.txt","w")
#  print(fid,matdata["$i"])
#  close(fid)
#end

using DSP, Test
using Statistics: mean
using FFTW: fftfreq

@testset "matlab ref" begin
    x0 = vec(readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t'))
    f0 = vec(readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_f.txt"),'\t'))
    t0 = vec(readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_t.txt"),'\t'))
    p0 = readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_p.txt"),'\t')
    spec = spectrogram(x0, 256, 128; fs=10)
    p, f, t = power(spec), freq(spec), time(spec)

    # with real input matlab outputs a 1-sided PSD
    @test p0 ≈ p
    @test f0 ≈ f
    @test t0 ≈ t

    mt_spec = mt_spectrogram(x0, 256, 128; fs=10)
    @test freq(mt_spec) == f
    @test time(mt_spec) == t
    @test power(mt_spec)[:,1] ≈ power(mt_pgram(x0[1:256]; fs=10))

    # in-place full precision version:
    spec_config = MTSpectrogramConfig{Float64}(length(x0), 256, 128; fs=10)
    out = allocate_output(spec_config)
    mt_spec2 = mt_spectrogram!(out, x0, spec_config)
    @test power(mt_spec2) ≈ power(mt_spec)
    @test freq(mt_spec2) == freq(mt_spec)
    @test time(mt_spec2) == time(mt_spec)

    # out-of-place with config:
    r = mt_spectrogram(x0, spec_config)
    @test power(r) ≈ power(mt_spec)
    @test freq(r) == freq(mt_spec)
    @test time(r) == time(mt_spec)

    # in-place without config:
    r = mt_spectrogram!(out, x0, 256, 128; fs=10)
    @test power(r) ≈ power(mt_spec)
    @test freq(r) == freq(mt_spec)
    @test time(r) == time(mt_spec)

    # in-place half precision version:
    spec_config = MTSpectrogramConfig{Float32}(length(x0), 256, 128; fs=10)
    out = allocate_output(spec_config)
    mt_spec3 = mt_spectrogram!(out, x0, spec_config)
    @test power(mt_spec3) ≈ power(mt_spec)
    @test freq(mt_spec3) == freq(mt_spec)
    @test time(mt_spec3) == time(mt_spec)

    mt_spec4 = mt_spectrogram!(out, Float32.(x0), spec_config)
    @test power(mt_spec4) ≈ power(mt_spec)
    @test freq(mt_spec4) == freq(mt_spec)
    @test time(mt_spec4) == time(mt_spec)

    # We can also only pass the window config. Full precision:
    config = MTConfig{Float64}(256; fs=10)
    mt_spec5 = mt_spectrogram(x0, config, 128)
    @test power(mt_spec5) ≈ power(mt_spec)
    @test freq(mt_spec5) == freq(mt_spec)
    @test time(mt_spec5) == time(mt_spec)

    # We can also only pass the window config. Half precision:
    config = MTConfig{Float32}(256; fs=10)
    mt_spec6 = mt_spectrogram(x0, config, 128)
    @test power(mt_spec6) ≈ power(mt_spec)
    @test freq(mt_spec6) == freq(mt_spec)
    @test time(mt_spec6) == time(mt_spec)

    # with Float32 input:
    mt_spec7 = mt_spectrogram(Float32.(x0), config, 128)
    @test power(mt_spec7) ≈ power(mt_spec)
    @test freq(mt_spec7) == freq(mt_spec)
    @test time(mt_spec7) == time(mt_spec)

    @test_throws DimensionMismatch mt_spectrogram!(similar(out, size(out, 1), size(out, 2)+1), x0, spec_config)
    @test_throws DimensionMismatch mt_spectrogram!(out, vcat(x0, x0), spec_config)
end

@testset "0:7" begin
    #Matlab: p = pwelch(0:7, [1, 1, 1, 1, 1, 1, 1, 1], 0, 8, 1, 'twosided')
    data = 0:7
    data0 = Float64[98.0,
                    13.656854249492380,
                    4.0,
                    2.343145750507620,
                    2.0,
                    2.343145750507620,
                    4.0,
                    13.656854249492380]
    @test power(periodogram(data, onesided=false)) ≈ data0
    @test power(welch_pgram(data, length(data), 0, onesided=false)) ≈ data0
    @test power(spectrogram(data, length(data), 0, onesided=false)) ≈ data0
    @test power(periodogram(complex.([data;], [data;]), onesided=false)) ≈ data0*2
    @test power(welch_pgram(complex.([data;], [data;]), length(data), 0, onesided=false)) ≈ data0*2
    @test power(spectrogram(complex.([data;], [data;]), length(data), 0, onesided=false)) ≈ data0*2

    # # ~~~~~~~~ Tests with no window ~~~~~~~~~~~~~~~~~~~
    # Matlab: p = pwelch(0:7, [1, 1], 0, 2, 1, 'twosided')
    expected = Float64[34.5, 0.5]
    @test power(welch_pgram(data, 2, 0; onesided=false)) ≈ expected
    @test mean(power(spectrogram(data, 2, 0; onesided=false)), dims=2) ≈ expected

    # Matlab: p = pwelch(0:7, [1, 1, 1], 0, 3, 1, 'twosided')
    expected = Float64[25.5, 1.0, 1.0]
    @test power(welch_pgram(data, 3, 0; onesided=false)) ≈ expected
    @test mean(power(spectrogram(data, 3, 0; onesided=false)), dims=2) ≈ expected

    # Matlab: p = pwelch(0:7, [1, 1, 1], 1, 3, 1, 'twosided')
    expected = Float64[35.0, 1.0, 1.0]
    @test power(welch_pgram(data, 3, 1; onesided=false)) ≈ expected
    @test mean(power(spectrogram(data, 3, 1; onesided=false)), dims=2) ≈ expected

    # Matlab: p = pwelch(0:7, [1, 1, 1, 1], 1, 4, 1, 'twosided')
    expected = Float64[45, 2, 1, 2]
    @test power(welch_pgram(data, 4, 1; onesided=false)) ≈ expected
    @test mean(power(spectrogram(data, 4, 1; onesided=false)), dims=2) ≈ expected

    # ~~~~~~~~~~~ This one tests periodogram ~~~~~~~~~~~~
    # ~ If functionality of the other arguments has been
    # ~ tested above, we only test here that the correct
    # ~ value of the spectral density is obtained when
    # ~ using a window. More tests to be added if needed
    #Matlab: p = pwelch(0:7, window_func(8), 0, 8, 1, 'twosided')
    cases = (
        (hamming,  Float64[65.461623986801527,
                           20.556791795515764,
                            0.369313143650544,
                            0.022167446610882,
                            0.025502985564107,
                            0.022167446610882,
                            0.369313143650544,
                           20.556791795515764]),
        (bartlett, Float64[62.999999999999993,
                           21.981076052592442,
                            0.285714285714286,
                            0.161781090264695,
                            0.142857142857143,
                            0.161781090264695,
                            0.285714285714286,
                           21.981076052592442])
    )

    @testset for (window1, expected) in cases
        @test power(periodogram(data; window=window1, onesided=false)) ≈ expected
        @test power(welch_pgram(data, length(data), 0; window=window1, onesided=false)) ≈ expected
        @test power(spectrogram(data, length(data), 0; window=window1, onesided=false)) ≈ expected
        @test power(periodogram(data; window=window1(length(data)), onesided=false)) ≈ expected
        @test power(welch_pgram(data, length(data), 0; window=window1(length(data)), onesided=false)) ≈ expected
        @test power(spectrogram(data, length(data), 0; window=window1(length(data)), onesided=false)) ≈ expected
    end

    # Padded periodogram
    # MATLAB: a = periodogram(0:7, [], 32);
    expected = [
                      98
        174.463067389405
        121.968086934209
        65.4971744936088
        27.3137084989848
        12.1737815028909
        10.3755170959439
        10.4034038628775
                       8
        5.25810953219633
        4.47015397150535
        4.89522578856669
        4.68629150101524
        3.69370284475603
         3.1862419983415
        3.61553458569862
                       2
    ]
    @test power(periodogram(data; nfft=32)) ≈ expected
    @test power(welch_pgram(data, length(data), 0; nfft=32)) ≈ expected
    @test power(spectrogram(data, length(data), 0; nfft=32)) ≈ expected

    # Padded periodogram with window
    # MATLAB: a = periodogram(0:7, hamming(8), 32, 1)
    expected = [
          65.4616239868015
          122.101693164395
          98.8444689598445
           69.020252632913
          41.1135835910315
          20.5496474310966
          8.43291449161938
          2.78001620362588
         0.738626287301088
         0.174995741770789
        0.0501563022944516
        0.0327357460012861
        0.0443348932217643
        0.0553999745503552
        0.0561319901616643
        0.0526025934871384
        0.0255029855641069
    ]
    @test power(periodogram(data; window=hamming, nfft=32)) ≈ expected
    @test power(welch_pgram(data, length(data), 0; window=hamming, nfft=32)) ≈ expected
    @test power(spectrogram(data, length(data), 0; window=hamming, nfft=32)) ≈ expected

    # Test fftshift
    p = periodogram(data)
    @test power(p) == power(fftshift(p))
    @test freq(p) ≈ freq(fftshift(p))

    p = periodogram(data; onesided=false)
    @test fftshift(power(p)) == power(fftshift(p))
    @test fftshift(freq(p)) == freq(fftshift(p))
end

@testset "fftshift" begin
    data = 1:100

    p = spectrogram(data)
    @test power(p) == power(fftshift(p))
    @test freq(p) ≈ freq(fftshift(p))

    p = spectrogram(data; onesided=false)
    @test fftshift(power(p), 1) == power(fftshift(p))
    @test fftshift(freq(p)) == freq(fftshift(p))
end

@testset "2D" begin
    data2d = readdlm(joinpath(dirname(@__FILE__), "data", "per2dx.txt"),'\t')
    expectedsum = vec(readdlm(joinpath(dirname(@__FILE__), "data", "per2dsum.txt"),'\t'))
    expectedmean = vec(readdlm(joinpath(dirname(@__FILE__), "data", "per2dmean.txt"),'\t'))
    # 2-d periodgram (radialsum)
    # computed in octave with raPsd2d ((C) E. Ruzanski) replacing nanmean with nansum
    # P = raPsd2d(x,1)'*n^2
    @test power(periodogram(data2d,fs=1, radialsum=true)) ≈ expectedsum

    # 2-d periodgram (radialavg)
    # computed in octave with raPsd2d ((C) E. Ruzanski)
    # P = raPsd2d(x,1)'*n^2
    @test power(periodogram(data2d, fs=1, radialavg=true)) ≈ expectedmean

    # 2-d periodgram 2-d PSD
    @test power(periodogram(data2d, fs=1)) ≈ abs2.(fft(data2d))*1/prod(size(data2d))
    # 2-d periodgram 2-d PSD with padding
    pads = (size(data2d,1)+4,size(data2d,1)+7)
    data2dpad = zeros(Float64,pads...)
    data2dpad[1:size(data2d,1),1:size(data2d,2)] = data2d
    @test power(periodogram(data2d, fs=1, nfft=pads)) ≈ abs2.(fft(data2dpad))*1/prod(size(data2d))
    # 2-d periodgram radial freq
    @test freq(periodogram(data2d, fs=3.3, radialsum=true)) ≈ freq(periodogram(vec(data2d[1,:]), fs=3.3))
    # 2-d periodgram 2-d freq
    f1,f2 = freq(periodogram(data2d, fs=3.3))
    f1d = freq(periodogram(vec(data2d[1,:]), fs=3.3, onesided=false))
    @assert size(data2d,1)==size(data2d,2)
    for j=1:size(data2d,2)
        for i=1:size(data2d,1)
            @test  [f1[i],f2[j]] ≈ [f1d[i],f1d[j]]
        end
    end
    # Test fftshift
    p = periodogram(data2d)
    @test fftshift(power(p)) == power(fftshift(p))
    f = freq(p)
    @test (fftshift(f[1]),fftshift(f[2])) == freq(fftshift(p))
end

@testset "radial" begin
    # 2-d periodgram radial test for a non-square signal sparse in fft space
    n1 = 52
    n2 = 46  # assuming n1>n2
    nf = (22,7) # the non-zero location
    F = (fftfreq(n1,1),fftfreq(n2,1))
    a = [F[1][nf[1]],F[2][nf[2]]]
    FB = Bool[[F[1][i], F[2][j]]==a || [F[1][i], F[2][j]]==-a for i = 1:n1, j = 1:n2]
    ind = findall(FB)
    x = zeros(n1,n2)*0im;
    x[ind] = [1+2im,1-2im]
    y = real(ifft(x))

    fwn = round(Int, sqrt((a[1])^2+(a[2])^2)*n2)
    pe = zeros(n2>>1 + 1)
    pe[fwn+1] = 2*abs2(x[nf...])/n1/n2
    P = periodogram(y,nfft=(n1,n2),radialsum=true)
    @test power(P) ≈ pe
    @test freq(P)[fwn+1] ≈ fwn/n2
end

# Testing STFT function and comparing results with MATLAB
@testset "stft" begin
    fs = 16000
    nfft = 512
    nwin = 400
    nhop = 160
    s = vec(readdlm(joinpath(dirname(@__FILE__), "data", "stft_x.txt"),'\t'))

    Sjl = stft(s, nwin, nwin-nhop; nfft=nfft, fs=fs, window=hanning)
    Sml_re = readdlm(joinpath(dirname(@__FILE__), "data", "stft_S_real.txt"),'\t')
    Sml_im = readdlm(joinpath(dirname(@__FILE__), "data", "stft_S_imag.txt"),'\t')
    Sml = complex.(Sml_re, Sml_im)
    @test Sjl ≈ Sml
end

@testset "fft2oneortwosided!" begin
    n = 10
    floatt = Float64
    for onesided in (true, false),
            nfft in (n, n+2, n+3),
           atype in (floatt, Complex{floatt})
        nout = nout = onesided ? (nfft >> 1)+1 : nfft
        x = zeros(atype, nfft)
        if atype <: Real
            x[1:n] = rand(atype, n)
            xrcfft = rfft(x)
        else
            x[1:n] = rand(floatt, n)+im*rand(floatt, n)
            xrcfft = fft(x)
        end
        xfft = fft(x)
        out = zeros(fftouttype(atype),nout,3)
        if !(onesided == true && atype <: Complex)
            outft = DSP.Periodograms.fft2oneortwosided!(out, xrcfft, nfft, onesided, nout)
        end
        if onesided == true && atype <: Real
            @test out[:,2] ≈ xrcfft
            @test out[:,[1,3]] ≈ xrcfft.*[0 0]
        elseif onesided == false && atype <: Real
            @test out[:,2] ≈ xfft
            @test out[:,[1,3]] ≈ xfft.*[0 0]
        elseif onesided == false && atype <: Complex
            @test out[:,2] ≈ xfft
            @test out[:,[1,3]] ≈ xfft.*[0 0]
        else
            #onesided complex
        end
    end
end

@testset "mt_pgram" begin
    # MATLAB: x = pmtm(stft_x, 4, 5000, 16000, 'unity')
    s = vec(readdlm(joinpath(dirname(@__FILE__), "data", "stft_x.txt"),'\t'))
    mtdata = vec(readdlm(joinpath(dirname(@__FILE__), "data", "mt_pgram.txt")))
    @test power(mt_pgram(s; fs=16000)) ≈ mtdata
    @test power(mt_pgram(s; fs=16000, window=dpss(length(s), 4))) ≈ mtdata

    # error tests
    @test_throws ArgumentError periodogram([1 2 3])
    @test_throws ArgumentError periodogram(rand(2,3), nfft=(3,2))
    @test_throws ArgumentError periodogram([1 2;3 4],radialsum=true, radialavg=true)

    # #124
    q = arraysplit(ones(Float64, 1000),100,10);
    @test map(mean, q) == ones(Float64, 11)

    # test that iterating ArraySplit always yields the same Vector
    q = arraysplit(-10:10,4,2)
    for x in q
        @assert isa(x, Vector)
        @test pointer(x) == pointer(q.buf)
    end

    x = vec(readdlm(joinpath(dirname(@__FILE__), "data", "pmtm_x.txt")))
    @test eltype(x) <: Float64

    # MATLAB code:
    # fft = 2^nextpow2(length(x));
    # [pxx,fx] = pmtm(x,4,nfft,1000,'unity');
    fx = vec(readdlm(joinpath(dirname(@__FILE__), "data", "pmtm_fx.txt")))
    pxx = vec(readdlm(joinpath(dirname(@__FILE__), "data", "pmtm_pxx.txt")))

    fs = 1000
    nw=4
    nfft = nextpow(2, length(x))
    result = mt_pgram(x; fs=fs, nw=nw, nfft=nfft)
    @test freq(result) ≈ fx
    @test power(result) ≈ pxx

    # Test against in-place. Full precision:
    config = MTConfig{Float64}(length(x); fs=fs, nw=nw, nfft=nfft)
    out = allocate_output(config)
    @test eltype(out) == Float64
    result2 = mt_pgram!(out, x, config)
    @test freq(result2) ≈ fx
    @test power(result2) ≈ pxx

    # in-place without config
    r = mt_pgram!(out, x; fs=fs, nw=nw, nfft=nfft)
    @test freq(r) ≈ fx
    @test power(r) ≈ pxx

    # out-of-place with config
    r = mt_pgram(x, config)
    @test freq(r) ≈ fx
    @test power(r) ≈ pxx

    # Lower precision output:
    config = MTConfig{Float32}(length(x); fs=fs, nw=nw, nfft=nfft)
    out = allocate_output(config)
    @test eltype(out) == Float32
    result3 = mt_pgram!(out, x, config)
    @test freq(result3) ≈ fx
    @test power(result3) ≈ pxx

    # with Float32 input:
    result4 = mt_pgram!(out, Float32.(x), config)
    @test freq(result4) ≈ fx
    @test power(result4) ≈ pxx

    @test_throws DimensionMismatch mt_pgram!(similar(out, length(out)+1), x, config)
    @test_throws DimensionMismatch mt_pgram!(out, vcat(x, one(eltype(x))), config)

    y = vec(readdlm(joinpath(dirname(@__FILE__), "data", "pmtm_y.txt")))
    z = x + im*y
    @test eltype(z) <: Complex{Float64}

    # MATLAB code: `[pzz,fz] = pmtm(z,4,nfft,1000, 'unity')`
    fz = vec(readdlm(joinpath(dirname(@__FILE__), "data", "pmtm_fz.txt")))
    pzz = vec(readdlm(joinpath(dirname(@__FILE__), "data", "pmtm_pzz.txt")))
    result = mt_pgram(z; fs=fs, nw=nw, nfft=nfft)
    result_mask = 0 .< freq(result) .< 500
    freqs = freq(result)[result_mask]
    @test freqs ≈ fz[2:length(freqs)+1]
    @test power(result)[result_mask] ≈ pzz[2:length(freqs)+1]

    # Test against in-place. Full precision:
    config = MTConfig{Complex{Float64}}(length(z); fs=fs, nw=nw, nfft=nfft)
    out = allocate_output(config)
    @test eltype(out) == Float64

    result2 = mt_pgram!(out, z, config)
    @test freq(result2) == freq(result)
    @test power(result2) ≈ power(result)

    # Lower precision output:
    config = MTConfig{Complex{Float32}}(length(z); fs=fs, nw=nw, nfft=nfft)
    out = allocate_output(config)
    @test eltype(out) == Float32

    result2 = mt_pgram!(out, z, config)
    @test freq(result2) == freq(result)
    @test power(result2) ≈ power(result)

    # Lower precision computation:
    result3 = mt_pgram!(out, Complex{Float32}.(z), config)
    @test freq(result3) == freq(result)
    @test power(result3) ≈ power(result)

end


text:
using Robotlib, Test, Random
@testset "Calib Force" begin
N         = 100
Rf        = Robotlib.toOrthoNormal(randn(3,3))
mf        = 1
POSES     = cat([Robotlib.toOrthoNormal(randn(3,3)) for i = 1:N]..., dims=3)
forceoffs = zeros(3)
forces    = hcat([Rf'*(POSES[1:3,1:3,i]'*[0, 0, mf*-9.82] - forceoffs) for i = 1:N]...)' |> copy

R,m = Robotlib.Calibration.calibForce(POSES,forces; offset=false, verbose=false)

@test R ≈ Rf
@test m ≈ mf

R = Robotlib.Calibration.calibForce2(POSES,forces,mf)
@test R ≈ Rf

forceoffs = randn(3)
forces    = hcat([Rf'*(POSES[1:3,1:3,i]'*[0, 0, mf*-9.82] - forceoffs) for i = 1:N]...)'

R,m,offs = Robotlib.Calibration.calibForce(POSES,forces; offset=true, verbose=false)

@test R ≈ Rf
@test m ≈ mf
@test offs ≈ forceoffs

end

function savefig2(filename)

    savefig(filename)
    temppath, temp = mktemp()
    pattern = r",\s*axis background/.style={fill={rgb,1:red,\d.\d+;green,\d.\d+;blue,\d.\d+}}"
    hp = r"height = {[\d.]+mm}"
    line = read(filename,String)
    # println(line)
    line = replace(line, pattern => "")
    line = replace(line, r"width = {[\d.]+mm}" => "")
    line = replace(line, r"height = {[\d.]+mm}" => "height = \\figureheight, width = \\figurewidth")
    println(temp, line)

    close(temp)
    mv(temppath, filename, remove_destination=true)
end

# using Plots, DSP, Robotlib, Random, LaTeXStrings
# pgfplots()
# σr = exp10.(range(-4,stop=1,length=1000))
# res = map(σr) do σ
#
#     N         = 100
#     Rf        = expω(1randn(3))
#     mf        = 1
#     POSES     = cat([expω(randn(3)) for i = 1:N]..., dims=3)
#     forceoffs = zeros(3)
#     forces    = hcat([Rf'*(POSES[1:3,1:3,i]'*[0, 0, mf*-9.82] - forceoffs) for i = 1:N]...)'
#     forces  .+= σ*randn(size(forces))
#
#     R1,m = Robotlib.Calibration.calibForce(POSES,forces,1.1mf; offset=false)
#
#     R2 = Robotlib.Calibration.calibForce2(POSES,forces,1.1mf)
#     sum(abs2,R1-Rf), sum(abs2,R2-Rf)
# end
#
# e1,e2 = getindex.(res,1),getindex.(res,2)
# plot(σr,filtfilt(ones(20),[20],[e1 e2]), xscale=:log10, yscale=:log10, lab=["Relaxation" "Cayley"], xlabel=L"$\sigma$ [N]", ylabel=L"$||R_e||_F^2$",legend=:bottomright)
# hline!([3, 3], lab="", l=:dash)
# savefig2("/home/fredrikb/phdthesis/calibpaper/figs/calibforce2.tex")


text:
let
    @test caldate_to_mjd(2000, 1, 1, 12, 0, 0) == 51544.5
end

let
    year, month, day, hour, minute, second, microsecond = mjd_to_caldate(51544.5)

    @test year        == 2000
    @test month       == 1
    @test day         == 1
    @test hour        == 12
    @test minute      == 0
    @test second      == 0
    @test microsecond == 0
end

let
    @test caldate_to_jd(2000, 1, 1, 12, 0, 0) == 2451545.0
end

let
    year, month, day, hour, minute, second, microsecond = jd_to_caldate(2451545.0)
    @test year        == 2000
    @test month       == 1
    @test day         == 1
    @test hour        == 12
    @test minute      == 0
    @test second      == 0
    @test microsecond == 0
end

let
    # MJD_ZERO = SatelliteDynamics.MJD_ZERO
    t = elapsed_from_epoch(MJD_ZERO + 175.500011574074074074073e-5, MJD_ZERO)
    @test isapprox(t, 15163200.000011574074074074073e-5, rtol=1e-6)

    t = elapsed_from_epoch(MJD_ZERO + 365.25*20, MJD_ZERO)
    @test isapprox(t, 631152000, rtol=1e-6)
end

let
    mjd = days_from_elapsed(15163200.000011574074074074073e-5, 0)
    @test isapprox(mjd, 175.500011574074074074073e-5, rtol=1e-3)

    mjd = days_from_elapsed(631152000, 0)
    @test isapprox(mjd, 7305, rtol=1e-3)
end

let
    @test SatelliteDynamics.valid_time_system("GPS") == true
    @test SatelliteDynamics.valid_time_system("TAI") == true
    @test SatelliteDynamics.valid_time_system("TT") == true
    @test SatelliteDynamics.valid_time_system("UTC") == true
    @test SatelliteDynamics.valid_time_system("UT1") == true

    @test SatelliteDynamics.valid_time_system("BLAH") == false
end

let
    jd = caldate_to_jd(2018, 6, 1)

    dutc = -37.0

    # GPS
    @test time_system_offset(jd, 0, "GPS", "GPS") == 0
    @test time_system_offset(jd, 0, "GPS", "TT")  == TT_GPS
    @test time_system_offset(jd, 0, "GPS", "UTC") == -18
    @test isapprox(time_system_offset(jd, 0, "GPS", "UT1"), -17.92267, atol=1e-4)
    @test time_system_offset(jd, 0, "GPS", "TAI") == TAI_GPS

    # TT
    @test time_system_offset(jd, 0, "TT", "GPS") == GPS_TT
    @test time_system_offset(jd, 0, "TT", "TT")  == 0
    @test time_system_offset(jd, 0, "TT", "UTC") == dutc + TAI_TT
    @test isapprox(time_system_offset(jd, 0, "TT", "UT1"), -69.10667, atol=1e-4)
    @test time_system_offset(jd, 0, "TT", "TAI") == TAI_TT

    # UTC
    @test time_system_offset(jd, 0, "UTC", "GPS") == 18
    @test time_system_offset(jd, 0, "UTC", "TT")  == -dutc + TT_TAI
    @test time_system_offset(jd, 0, "UTC", "UTC") == 0.0
    @test isapprox(time_system_offset(jd, 0, "UTC", "UT1"), 0.0769968, atol=1e-4)
    @test time_system_offset(jd, 0, "UTC", "TAI") == -dutc

    # UT1
    @test isapprox(time_system_offset(jd, 0, "UT1", "GPS"), 17.9230032, atol=1e-4)
    @test isapprox(time_system_offset(jd, 0, "UT1", "TT"), 69.1070032, atol=1e-4)
    @test isapprox(time_system_offset(jd, 0, "UT1", "UTC"), -0.0769968, atol=1e-4)
    @test time_system_offset(jd, 0, "UT1", "UT1") == 0
    @test isapprox(time_system_offset(jd, 0, "UT1", "TAI"), 36.9230032, atol=1e-4)

    # TAI
    @test time_system_offset(jd, 0, "TAI", "GPS") == GPS_TAI
    @test time_system_offset(jd, 0, "TAI", "TT")  == TT_TAI
    @test time_system_offset(jd, 0, "TAI", "UTC") == dutc
    @test isapprox(time_system_offset(jd, 0, "TAI", "UT1"), -36.92267, atol=1e-4)
    @test time_system_offset(jd, 0, "TAI", "TAI") == 0
end

let
    epc = Epoch(2018, 12, 20, 0, 0, 0, 0.0, tsys="TAI")
    @test epc.days        == 2458472
    @test epc.seconds     == 43200
    @test epc.nanoseconds == 0.0
    @test epc.tsys        == "TAI"

    epc = Epoch(2018, 12, 20, 0, 0, .5, 1.0001, tsys="TAI")
    @test epc.days        == 2458472
    @test epc.seconds     == 43200
    @test epc.nanoseconds == 500000001.0001
    @test epc.tsys        == "TAI"
end

# String constructor
let
    epc = Epoch("2018-12-20")
    year, month, day, hour, minute, seconds, nanoseconds = caldate(epc)
    @test year        == 2018
    @test month       == 12
    @test day         == 20
    @test hour        == 0
    @test minute      == 0
    @test seconds     == 0
    @test nanoseconds == 0.0
    @test epc.tsys    == "UTC"

    epc = Epoch("2018-12-20T16:22:19.0Z")
    year, month, day, hour, minute, seconds, nanoseconds = caldate(epc)
    @test year        == 2018
    @test month       == 12
    @test day         == 20
    @test hour        == 16
    @test minute      == 22
    @test seconds     == 19
    @test nanoseconds == 0.0
    @test epc.tsys    == "UTC"

    epc = Epoch("2018-12-20T16:22:19.123Z")
    year, month, day, hour, minute, seconds, nanoseconds = caldate(epc)
    @test year        == 2018
    @test month       == 12
    @test day         == 20
    @test hour        == 16
    @test minute      == 22
    @test seconds     == 19
    @test nanoseconds == 123000000
    @test epc.tsys    == "UTC"

    epc = Epoch("2018-12-20T16:22:19.123456789Z")
    year, month, day, hour, minute, seconds, nanoseconds = caldate(epc)
    @test year        == 2018
    @test month       == 12
    @test day         == 20
    @test hour        == 16
    @test minute      == 22
    @test seconds     == 19
    @test nanoseconds == 123456789
    @test epc.tsys    == "UTC"

    epc = Epoch("2018-12-20T16:22:19Z")
    year, month, day, hour, minute, seconds, nanoseconds = caldate(epc)
    @test year        == 2018
    @test month       == 12
    @test day         == 20
    @test hour        == 16
    @test minute      == 22
    @test seconds     == 19
    @test nanoseconds == 0.0
    @test epc.tsys    == "UTC"

    epc = Epoch("20181220T162219Z")
    year, month, day, hour, minute, seconds, nanoseconds = caldate(epc)
    @test year        == 2018
    @test month       == 12
    @test day         == 20
    @test hour        == 16
    @test minute      == 22
    @test seconds     == 19
    @test nanoseconds == 0.0
    @test epc.tsys    == "UTC"

    epc = Epoch("2018-12-01 16:22:19 GPS")
    year, month, day, hour, minute, seconds, nanoseconds = caldate(epc)
    @test year        == 2018
    @test month       == 12
    @test day         == 1
    @test hour        == 16
    @test minute      == 22
    @test seconds     == 19
    @test nanoseconds == 0.0
    @test epc.tsys    == "GPS"

    epc = Epoch("2018-12-01 16:22:19.0 GPS")
    year, month, day, hour, minute, seconds, nanoseconds = caldate(epc)
    @test year        == 2018
    @test month       == 12
    @test day         == 1
    @test hour        == 16
    @test minute      == 22
    @test seconds     == 19
    @test nanoseconds == 0.0
    @test epc.tsys    == "GPS"

    epc = Epoch("2018-12-01 16:22:19.123 GPS")
    year, month, day, hour, minute, seconds, nanoseconds = caldate(epc)
    @test year        == 2018
    @test month       == 12
    @test day         == 1
    @test hour        == 16
    @test minute      == 22
    @test seconds     == 19
    @test nanoseconds == 123000000
    @test epc.tsys    == "GPS"

    epc = Epoch("2018-12-01 16:22:19.123456789 GPS")
    year, month, day, hour, minute, seconds, nanoseconds = caldate(epc)
    @test year        == 2018
    @test month       == 12
    @test day         == 1
    @test hour        == 16
    @test minute      == 22
    @test seconds     == 19
    @test nanoseconds == 123456789
    @test epc.tsys    == "GPS"

end

# Addition
let
    epc = Epoch("2019-01-01 12:00:00 TAI")
    @test epc.days        == 2458485
    @test epc.seconds     == 0
    @test epc.nanoseconds == 0

    epc += 1.0e-9
    @test epc.days        == 2458485
    @test epc.seconds     == 0
    @test epc.nanoseconds == 1

    epc -= 2.0e-9
    @test epc.days        == 2458484
    @test epc.seconds     == 86399
    @test epc.nanoseconds == 999999999

    for i in 1:365*86400
        epc += 1
    end
    @test epc.days        == 2458849
    @test epc.seconds     == 86399
    @test epc.nanoseconds == 999999999

    year, month, day, hour, minute, second, nanoseconds = caldate(epc)
    @test year        == 2020
    @test month       == 1
    @test day         == 1
    @test hour        == 11
    @test minute      == 59
    @test second      == 59
    @test nanoseconds == 999999999
end

# Test long term addition
let
    epc = Epoch(1980, 1, 1, 0, 0, 0, 1)
    year, month, day, hour, minute, second, nanoseconds = caldate(epc)
    @test year        == 1980
    @test month       == 1
    @test day         == 1
    @test hour        == 0
    @test minute      == 0
    @test second      == 0
    @test nanoseconds == 1
    @test epc.tsys    == "TAI"

    for i in 1:40*365.25*86400
        epc += 1
    end

    year, month, day, hour, minute, second, nanoseconds = caldate(epc)
    @test year        == 2020
    @test month       == 1
    @test day         == 1
    @test hour        == 0
    @test minute      == 0
    @test second      == 0
    @test nanoseconds == 1
    @test epc.tsys    == "TAI"
end

# Direct conversion
let
    epc = Epoch("2018-12-20T16:22:19.123456789Z")
    year, month, day, hour, minute, second, nanoseconds = caldate(epc)
    @test year        == 2018
    @test month       == 12
    @test day         == 20
    @test hour        == 16
    @test minute      == 22
    @test second      == 19
    @test nanoseconds == 123456789.0
end

# Conversion into different time system
let
    epc = Epoch("2018-12-20T16:22:19.123456789Z")
    year, month, day, hour, minute, second, nanoseconds = caldate(epc, tsys="TAI")
    @test year        == 2018
    @test month       == 12
    @test day         == 20
    @test hour        == 16
    @test minute      == 22
    @test second      == 19 + 37
    @test nanoseconds == 123456789.0

    epc = Epoch("2018-12-20T16:22:19.123456789Z")
    year, month, day, hour, minute, second, nanoseconds = caldate(epc, tsys="GPS")
    @test year        == 2018
    @test month       == 12
    @test day         == 20
    @test hour        == 16
    @test minute      == 22
    @test second      == 19 + 37 - 19
    @test nanoseconds == 123456789.0
end

let
    epc = Epoch(2000, 1, 1)
    @test jd(epc) == 2451544.5
    @test jd(epc, tsys="UTC") < 2451544.5
end

let
    epc = Epoch(2000, 1, 1, 0)
    @test mjd(epc) == 51544
    @test mjd(epc, tsys="UTC") < 51544
end

let
    epc = Epoch(2000, 1, 1)
    @test day_of_year(epc) == 1

    epc = Epoch(2000, 1, 1, 12)
    @test day_of_year(epc) == 1.5

    epc = Epoch(2000, 12, 31)
    @test day_of_year(epc) == 366

    epc = Epoch(2001, 1, 1)
    @test day_of_year(epc) == 1

    epc = Epoch(2001, 12, 31)
    @test day_of_year(epc) == 365
end

let
    epc = Epoch(2000, 1, 1)
    @test isapprox(gmst(epc, use_degrees=true), 99.835, atol=1e-3)

    epc = Epoch(2000, 1, 1)
    @test isapprox(gmst(epc, use_degrees=false), 1.742, atol=1e-3)
end

let
    epc = Epoch(2000, 1, 1)
    @test isapprox(gast(epc, use_degrees=true), 99.832, atol=1e-3)

    epc = Epoch(2000, 1, 1)
    @test isapprox(gast(epc, use_degrees=false), 1.742, atol=1e-3)
end

# Arithmetic Comparisons
let
    @test Epoch(2000, 1, 1, 12, 23, 59, 123456789) == Epoch(2000, 1, 1, 12, 23, 59, 123456789)
    @test Epoch(2000, 1, 1, 12, 23, 59, 123456789) != Epoch(2000, 1, 1, 12, 23, 59, 123456788)
end

let
    @test !(Epoch(2000, 1, 1, 12, 23, 59, 123456789) < Epoch(2000, 1, 1, 12, 23, 59, 123456789))
    @test !(Epoch(2000, 1, 1, 12, 23, 59, 123456789) <= Epoch(2000, 1, 1, 12, 23, 59, 123456788))
    @test Epoch(2000, 1, 1, 12, 23, 59, 123456788) < Epoch(2000, 1, 1, 12, 23, 59, 123456789)
    @test Epoch(2000, 1, 1, 12, 23, 59, 123456788) <= Epoch(2000, 1, 1, 12, 23, 59, 123456789)
end

let
    @test !(Epoch(2000, 1, 1, 12, 23, 59, 123456789) > Epoch(2000, 1, 1, 12, 23, 59, 123456789))
    @test !(Epoch(2000, 1, 1, 12, 23, 59, 123456788) >= Epoch(2000, 1, 1, 12, 23, 59, 123456789))
    @test Epoch(2000, 1, 1, 12, 23, 59, 123456789) > Epoch(2000, 1, 1, 12, 23, 59, 123456788)
    @test Epoch(2000, 1, 1, 12, 23, 59, 123456789) >= Epoch(2000, 1, 1, 12, 23, 59, 123456788)
end