text:
using ..Atmos
using ..Atmos: thermo_state
using ..Mesh.Topologies
using ..Mesh.Grids
using ..MoistThermodynamics
using LinearAlgebra

"""
    atmos_core_init(bl, currtime)

Initialize the 'AtmosCore' diagnostics group.
"""
function atmos_core_init(dgngrp::DiagnosticsGroup, currtime)
    dg = Settings.dg
    bl = dg.balancelaw
    if isa(bl.moisture, DryModel)
        @warn """
            Diagnostics ($dgngrp.name): cannot be used with `DryModel`
            """
    end

    atmos_collect_onetime(Settings.mpicomm, Settings.dg, Settings.Q)

    return nothing
end

# Simple horizontal averages
function vars_atmos_core_simple(m::AtmosModel, FT)
    @vars begin
        u_core::FT
        v_core::FT
        w_core::FT
        avg_rho_core::FT        # ρ
        rho_core::FT            # ρρ
        qt_core::FT             # q_tot
        thl_core::FT            # θ_liq
        ei_core::FT             # e_int
    end
end
num_atmos_core_simple_vars(m, FT) = varsize(vars_atmos_core_simple(m, FT))
atmos_core_simple_vars(m, array) =
    Vars{vars_atmos_core_simple(m, eltype(array))}(array)

function atmos_core_simple_sums!(
    atmos::AtmosModel,
    state_conservative,
    thermo,
    MH,
    sums,
)
    sums.u_core += MH * state_conservative.ρu[1]
    sums.v_core += MH * state_conservative.ρu[2]
    sums.w_core += MH * state_conservative.ρu[3]
    sums.avg_rho_core += MH * state_conservative.ρ
    sums.rho_core += MH * state_conservative.ρ * state_conservative.ρ
    sums.qt_core += MH * state_conservative.moisture.ρq_tot
    sums.thl_core += MH * thermo.moisture.θ_liq_ice * state_conservative.ρ
    sums.ei_core += MH * thermo.e_int * state_conservative.ρ

    return nothing
end

# Variances and covariances
function vars_atmos_core_ho(m::AtmosModel, FT)
    @vars begin
        var_u_core::FT          # u′u′
        var_v_core::FT          # v′v′
        var_w_core::FT          # w′w′
        var_qt_core::FT         # q_tot′q_tot′
        var_thl_core::FT        # θ_liq_ice′θ_liq_ice′
        var_ei_core::FT         # e_int′e_int′

        cov_w_rho_core::FT      # w′ρ′
        cov_w_qt_core::FT       # w′q_tot′
        cov_w_thl_core::FT      # w′θ_liq_ice′
        cov_w_ei_core::FT       # w′e_int′
        cov_qt_thl_core::FT     # q_tot′θ_liq_ice′
        cov_qt_ei_core::FT      # q_tot′e_int′
    end
end
num_atmos_core_ho_vars(m, FT) = varsize(vars_atmos_core_ho(m, FT))
atmos_core_ho_vars(m, array) = Vars{vars_atmos_core_ho(m, eltype(array))}(array)

function atmos_core_ho_sums!(
    atmos::AtmosModel,
    state_conservative,
    thermo,
    MH,
    ha,
    sums,
)
    u = state_conservative.ρu[1] / state_conservative.ρ
    u′ = u - ha.u_core
    v = state_conservative.ρu[2] / state_conservative.ρ
    v′ = v - ha.v_core
    w = state_conservative.ρu[3] / state_conservative.ρ
    w′ = w - ha.w_core
    q_tot = state_conservative.moisture.ρq_tot / state_conservative.ρ
    q_tot′ = q_tot - ha.qt_core
    θ_liq_ice′ = thermo.moisture.θ_liq_ice - ha.thl_core
    e_int′ = thermo.e_int - ha.ei_core

    sums.var_u_core += MH * u′^2 * state_conservative.ρ
    sums.var_v_core += MH * v′^2 * state_conservative.ρ
    sums.var_w_core += MH * w′^2 * state_conservative.ρ
    sums.var_qt_core += MH * q_tot′^2 * state_conservative.ρ
    sums.var_thl_core += MH * θ_liq_ice′^2 * state_conservative.ρ
    sums.var_ei_core += MH * e_int′^2 * state_conservative.ρ

    sums.cov_w_rho_core +=
        MH *
        w′ *
        (state_conservative.ρ - ha.avg_rho_core) *
        state_conservative.ρ
    sums.cov_w_qt_core += MH * w′ * q_tot′ * state_conservative.ρ
    sums.cov_w_thl_core += MH * w′ * θ_liq_ice′ * state_conservative.ρ
    sums.cov_qt_thl_core += MH * q_tot′ * θ_liq_ice′ * state_conservative.ρ
    sums.cov_qt_ei_core += MH * q_tot′ * e_int′ * state_conservative.ρ
    sums.cov_w_ei_core += MH * w′ * e_int′ * state_conservative.ρ

    return nothing
end

"""
    atmos_core_collect(bl, currtime)

Perform a global grid traversal to compute various diagnostics.
"""
function atmos_core_collect(dgngrp::DiagnosticsGroup, currtime)
    dg = Settings.dg
    bl = dg.balancelaw
    if isa(bl.moisture, DryModel)
        @warn """
            Diagnostics $(dgngrp.name): cannot be used with `DryModel`
            """
    end

    mpicomm = Settings.mpicomm
    Q = Settings.Q
    mpirank = MPI.Comm_rank(mpicomm)
    current_time = string(currtime)

    # extract grid information
    grid = dg.grid
    topology = grid.topology
    N = polynomialorder(grid)
    Nq = N + 1
    Nqk = dimensionality(grid) == 2 ? 1 : Nq
    npoints = Nq * Nq * Nqk
    nrealelem = length(topology.realelems)
    nvertelem = topology.stacksize
    nhorzelem = div(nrealelem, nvertelem)

    # get needed arrays onto the CPU
    if Array ∈ typeof(Q).parameters
        host_state_conservative = Q.realdata
        host_state_auxiliary = dg.state_auxiliary.realdata
        host_vgeo = grid.vgeo
    else
        host_state_conservative = Array(Q.realdata)
        host_state_auxiliary = Array(dg.state_auxiliary.realdata)
        host_vgeo = Array(grid.vgeo)
    end
    FT = eltype(host_state_conservative)

    zvals = AtmosCollected.zvals

    # Visit each node of the state variables array and:
    # - generate and store the thermo variables,
    # - if core condition holds (q_liq > 0 && w > 0)
    #   - count that point in the core fraction for that z
    #   - count the point's weighting towards averaging for that z, and
    #   - accumulate the simple horizontal sums
    #
    core_repdvsr = zeros(FT, Nqk * nvertelem)
    thermo_array =
        [zeros(FT, num_thermo(bl, FT)) for _ in 1:npoints, _ in 1:nrealelem]
    simple_sums = [
        zeros(FT, num_atmos_core_simple_vars(bl, FT))
        for _ in 1:(Nqk * nvertelem)
    ]
    ql_w_gt_0 = [zeros(FT, (Nq * Nq * nhorzelem)) for _ in 1:(Nqk * nvertelem)]
    @visitQ nhorzelem nvertelem Nqk Nq begin
        evk = Nqk * (ev - 1) + k

        state_conservative =
            extract_state_conservative(dg, host_state_conservative, ijk, e)
        state_auxiliary =
            extract_state_auxiliary(dg, host_state_auxiliary, ijk, e)
        MH = host_vgeo[ijk, grid.MHid, e]

        thermo = thermo_vars(bl, thermo_array[ijk, e])
        compute_thermo!(bl, state_conservative, state_auxiliary, thermo)

        if thermo.moisture.q_liq > 0 && state_conservative.ρu[3] > 0
            idx = (Nq * Nq * (eh - 1)) + (Nq * (j - 1)) + i
            ql_w_gt_0[evk][idx] = one(FT)
            core_repdvsr[evk] += MH

            simple = atmos_core_simple_vars(bl, simple_sums[evk])
            atmos_core_simple_sums!(bl, state_conservative, thermo, MH, simple)
        end
    end

    # reduce horizontal sums and core fraction across ranks and compute averages
    simple_avgs = [
        zeros(FT, num_atmos_core_simple_vars(bl, FT))
        for _ in 1:(Nqk * nvertelem)
    ]
    core_frac = zeros(FT, Nqk * nvertelem)
    MPI.Allreduce!(core_repdvsr, +, mpicomm)
    for evk in 1:(Nqk * nvertelem)
        tot_ql_w_gt_0 = MPI.Reduce(sum(ql_w_gt_0[evk]), +, 0, mpicomm)
        tot_horz = MPI.Reduce(length(ql_w_gt_0[evk]), +, 0, mpicomm)

        MPI.Allreduce!(simple_sums[evk], simple_avgs[evk], +, mpicomm)
        simple_avgs[evk] .= simple_avgs[evk] ./ core_repdvsr[evk]

        if mpirank == 0
            core_frac[evk] = tot_ql_w_gt_0 / tot_horz
        end
    end

    # complete density averaging
    simple_varnames = flattenednames(vars_atmos_core_simple(bl, FT))
    for vari in 1:length(simple_varnames)
        for evk in 1:(Nqk * nvertelem)
            simple_ha = atmos_core_simple_vars(bl, simple_avgs[evk])
            avg_rho = simple_ha.avg_rho_core
            if simple_varnames[vari] != "avg_rho_core"
                simple_avgs[evk][vari] /= avg_rho
            end
        end
    end

    # compute the variances and covariances
    ho_sums =
        [zeros(FT, num_atmos_core_ho_vars(bl, FT)) for _ in 1:(Nqk * nvertelem)]
    @visitQ nhorzelem nvertelem Nqk Nq begin
        evk = Nqk * (ev - 1) + k

        state_conservative =
            extract_state_conservative(dg, host_state_conservative, ijk, e)
        thermo = thermo_vars(bl, thermo_array[ijk, e])
        MH = host_vgeo[ijk, grid.MHid, e]

        if thermo.moisture.q_liq > 0 && state_conservative.ρu[3] > 0
            simple_ha = atmos_core_simple_vars(bl, simple_avgs[evk])
            ho = atmos_core_ho_vars(bl, ho_sums[evk])
            atmos_core_ho_sums!(
                bl,
                state_conservative,
                thermo,
                MH,
                simple_ha,
                ho,
            )
        end
    end

    # reduce across ranks and compute averages
    ho_avgs =
        [zeros(FT, num_atmos_core_ho_vars(bl, FT)) for _ in 1:(Nqk * nvertelem)]
    for evk in 1:(Nqk * nvertelem)
        MPI.Reduce!(ho_sums[evk], ho_avgs[evk], +, 0, mpicomm)
        if mpirank == 0
            ho_avgs[evk] .= ho_avgs[evk] ./ core_repdvsr[evk]
        end
    end

    # complete density averaging and prepare output
    if mpirank == 0
        varvals = OrderedDict()
        varvals["core_frac"] = (("z",), core_frac)

        for vari in 1:length(simple_varnames)
            davg = zeros(FT, Nqk * nvertelem)
            for evk in 1:(Nqk * nvertelem)
                davg[evk] = simple_avgs[evk][vari]
            end
            varvals[simple_varnames[vari]] = (("z",), davg)
        end

        ho_varnames = flattenednames(vars_atmos_core_ho(bl, FT))
        for vari in 1:length(ho_varnames)
            davg = zeros(FT, Nqk * nvertelem)
            for evk in 1:(Nqk * nvertelem)
                simple_ha = atmos_core_simple_vars(bl, simple_avgs[evk])
                avg_rho = simple_ha.avg_rho_core
                davg[evk] = ho_avgs[evk][vari] / avg_rho
            end
            varvals[ho_varnames[vari]] = (("z",), davg)
        end

        # write output
        dprefix = @sprintf(
            "%s_%s_%s_num%04d",
            dgngrp.out_prefix,
            dgngrp.name,
            Settings.starttime,
            dgngrp.num
        )
        dfilename = joinpath(Settings.output_dir, dprefix)
        write_data(
            dgngrp.writer,
            dfilename,
            OrderedDict("z" => zvals),
            varvals,
            currtime,
        )
    end

    MPI.Barrier(mpicomm)
    return nothing
end # function collect

function atmos_core_fini(dgngrp::DiagnosticsGroup, currtime) end


text:
#=
  Bayesian Null hypothesis test in Turing.jl
  From https://mhtess.github.io/bdappl/chapters/04-hypothesisTesting.html
  (WebPPL model)
  """
  Let’s see an example: Consider our children helping study from the previous study. Rather than estimate the
  propensity_to_help, we might want to perform a “Null Hypothesis” test: What is the probability that the behavior
  observed in our experiment (that 15 out of 20 kids helped) could be explained by true randomness (i.e.,
  propensity_to_help = 0.5)?

  ...

  Bayesian Null Hypothesis Testing requires you to specify both hypotheses:
  It is not enough to define what the null hypothesis is. The alternative
  hypothesis is intuitively M1:θ≠0.5\mathcal{M}_1: \theta \neq 0.5M1:θ≠0.5,
  but we must be more specific. One standard way of defining the null
  hypothesis is to put say that the parameter is drawn from a uniform
  distribution.
  """

  From the WebPPL model:
  Marginal:
    {"the_better_model":"alternative"} : 0.758
    {"the_better_model":"null"} : 0.242

  This Turing.jl model:
   Distributions of variable the_better_model (num:0)
   alternative model =>    7612  (0.761200)
   null       =>    2388  (0.238800)

=#

using Turing
include("jl_utils.jl")

@model function supermodel(k=15,n=20)
    null_model ~ Dirac(0.5)
    alternative_model ~ Uniform(0,1)

    # Pick one of the models
    null = 1
    alternative = 2
    the_better_model ~ DiscreteUniform(1,2) # Categorical([0.5,0.5])
    p ~ Dirac(the_better_model == null ? null_model : alternative_model)
    k ~ Binomial(n,p)

end

# data
k=15
n=20
println("We observe $(k) successes of $(n) attempts")
model = supermodel(k,n)
# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(15), 10_000)
# chns = sample(model, IS(), 10_000)
chns = sample(model, SMC(), 10_000)

display(chns)

show_var_dist_pct(chns,:the_better_model,["null","alternative model"])

ext: jl

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Test, Countries

@testset "Countries" begin

    @testset "Country" begin

        @testset "from alpha2" begin
            @test Country("AA").idx == 1
            @test Country("AB").idx == 2
            @test Country("BA").idx == 27
            @test Country("ZZ").idx == 26*26
        end

        @testset "from alpha3" begin
            @test Country("GBR") === Country("GB")
            @test Country("FRA") === Country("FR")
            @test_throws Exception Country("ZZZ")
        end

        @testset "from numeric" begin
            @test Country(826) === Country("GB")
            @test Country(250) === Country("FR")
            @test_throws Exception Country(999)
        end

        @testset "from country" begin
            @test Country(Country("GB")) === Country("GB")
        end

        @testset "from name" begin
            @test Country("United Kingdom of Great Britain and Northern Ireland") === Country("GB")
            @test Country("France") === Country("FR")
            @test_throws Exception Country("This Is Not A Country")
        end

        @testset "partial" begin
            @test Country("United Kingdom") === Country("GB")
            @test Country("Britain") === Country("GB")
            @test_throws Exception Country("Korea")
        end

        @testset "case insensitive" begin
            @test Country("de") === Country("DE")
            @test Country("De") === Country("DE")
            @test Country("dE") === Country("DE")
            @test Country("france") === Country("FR")
            @test Country("FRANCE") === Country("FR")
            @test Country("FrAnCe") === Country("FR")
        end

    end

    @testset "properties" begin

        @testset "country_numeric" begin
            @test country_numeric("Norway") == 578
            @test country_numeric("ZZ") == 0
        end

        @testset "country_alpha2" begin
            @test country_alpha2("Portugal") == "PT"
            @test country_alpha2("ZZ") == "ZZ"
        end

        @testset "country_alpha3" begin
            @test country_alpha3("Sweden") == "SWE"
            @test country_alpha3("ZZ") == ""
        end

        @testset "country_name" begin
            @test country_name("BE") == "Belgium"
            @test country_name("ZZ") == ""
        end

        @testset "country_assigned" begin
            @test country_assigned("Italy") == true
            @test country_assigned("ZZ") == false
        end

    end

    @testset "each_country" begin
        cs = collect(each_country())
        @test cs isa Vector{Country}
        @test length(cs) == 249
    end

    @testset "IO" begin

        c = Country("France")
        @test string(c) == "FR"
        @test repr(c) == "Country(\"FR\")"
        @test repr(c, context=:typeinfo=>Country) == "FR"
        @test repr("text/plain", c) == "FR: France"
        @test repr("text/plain", c, context=:compact=>true) == "Country(\"FR\")"

        io = IOBuffer()
        write(io, c)
        seekstart(io)
        @test read(io, Country) === c

    end

    @testset "comparison" begin

        @testset "==" begin
            @test Country("FR") == Country("France")
            @test Country("FR") != Country("ZZ")
        end

        @testset "isequal" begin
            @test isequal(Country("FR"), Country("France"))
            @test !isequal(Country("FR"), Country("ZZ"))
        end

        @testset "hash" begin
            @test hash(Country("FR")) == hash(Country("France"))
            @test hash(Country("FR")) != hash(Country("ZZ"))
            @test Set(Country.(["FR", "France", "ES", "Spain"])) == Set([Country("FR"), Country("ES")])
        end

        @testset "isless" begin
            @test !isless(Country("FR"), Country("FR"))
            @test isless(Country("AA"), Country("FR"))
            @test isless(Country("FR"), Country("ZZ"))
        end

    end

    @testset "new_country" begin

        @test country_assigned("ZL") == false
        new_country(alpha2="ZL", alpha3="ZLD", name="Zedland", numeric=900)
        @test country_assigned("ZL") == true
        @test country_alpha2("Zedland") == "ZL"
        @test country_alpha3("Zedland") == "ZLD"
        @test country_name("ZL") == "Zedland"
        @test country_numeric("ZLD") == 900

        @test_throws Exception new_country(alpha2="ZL", alpha3="ZLD", name="Zedland2", numeric=900)
        @test country_name("ZL") == "Zedland"

        @test_throws Exception new_country(alpha2="ZZ", alpha3="ZLD", name="Zedland2", numeric=999)
        @test country_assigned("ZZ") == false

        @test_throws Exception new_country(alpha2="ZZ", alpha3="ZZZ", name="Zedland", numeric=999)
        @test country_assigned("ZZ") == false

        @test_throws Exception new_country(alpha2="ZZ", alpha3="ZZZ", name="ZedZedZedland", numeric=900)
        @test country_assigned("ZZ") == false

    end

    @testset "alias_country" begin

        @test_throws Exception country_alpha2("my alias")
        alias_country("My Alias", "GBR")
        @test country_alpha2("my alias") === "GB"

        @test_throws Exception alias_country("My Alias", "FRA")
        @test country_alpha2("my alias") == "GB"

        alias_country("My Alias", "GBR")
        @test country_alpha2("my alias") === "GB"

        @test_throws Exception alias_country("", "GBR")
        @test_throws Exception alias_country("Z", "GBR")
        @test_throws Exception alias_country("ZZ", "GBR")
        @test_throws Exception alias_country("ZZZ", "GBR")
        @test_throws Exception alias_country("France", "GBR")

    end

end


text:
function parse_encoded_string(input_string, lower_bound=0, upper_bound=127, upper_string="F", lower_string="B")

    for direction in input_string
        letter = string(direction)
        if letter == upper_string
            upper_bound = upper_bound - Int32(div(upper_bound - lower_bound, 2, RoundUp))
        elseif letter == lower_string
            lower_bound = lower_bound + Int32(div(upper_bound - lower_bound, 2, RoundUp))
        else
            throw("Invalid string type $letter")
        end
        # println("letter: $letter")
        # println((lower_bound, upper_bound))
    end
    @assert lower_bound == upper_bound
    lower_bound
end

function parse_rows(input_string)
    return parse_encoded_string(input_string, 0, 127, "F", "B")
end

function parse_columns(input_string)
    return parse_encoded_string(input_string, 0, 7, "L", "R")
end

function get_seat_id(input_string)
    row_string = input_string[1:7]
    col_string = input_string[8:end]

    row_num = parse_rows(row_string)
    col_num = parse_columns(col_string)

    seat_id = row_num * 8 + col_num
    seat_id
end

# turns out we didn't need this
function id_to_row_column(seat_id)
    column = seat_id % 8
    row = (seat_id - column) / 8
    row, column
end

function main()
    input_file = "day5_input.txt"

    highest_boarding_pass = open(input_file) do file
        file_string = read(file, String)
        all_seats = split(file_string, "\n")
        all_seats_ids = Vector{Int32}()

        for seat_entry in all_seats
            seat_id = get_seat_id(strip(seat_entry))
            push!(all_seats_ids, seat_id)
        end

        set_of_all_seats = Set(all_seats_ids) # easier to test if things are in a set
        for i in minimum(all_seats_ids):maximum(all_seats_ids)
            if !in(i, set_of_all_seats) && in(i-1, set_of_all_seats) && in(i+1, set_of_all_seats)
                println("My seat ID is $i")
            end
        end

        # Debugging
        # sort!(all_seats_ids)
        # println("all_seats_ids: $all_seats_ids")
        maximum(all_seats_ids)
    end
    println("Highest Seat ID: $highest_boarding_pass")
end

#if abspath(PROGRAM_FILE) == @__FILE__
    main()
#end


text:
module MeshIO

using GeometryTypes
using GeometryTypes: raw
using ColorTypes
using Printf
import FileIO

import FileIO: DataFormat, @format_str, Stream, File, filename, stream, skipmagic
import Base.show

include("io/off.jl")
include("io/ply.jl")
include("io/stl.jl")
include("io/obj.jl")
include("io/2dm.jl")

load(fn::File{format}, MeshType=GLNormalMesh) where {format} = open(fn) do s
    skipmagic(s)
    load(s, MeshType)
end
save(fn::File{format}, msh::AbstractMesh) where {format} = open(fn, "w") do s
    save(s, msh)
end

end # module


text:
#=
    A Robust Complex Division in Scilab
    by Michael Baudin, Robert L. Smith
=#

function cdiv(x::Complex{DoubleFloat{T}},
             y::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    a, b = reim(x)
    c, d = reim(y)
    if abs(d) <= abs(c)
        e, f = robustdiv(a, b, c, d)
    else
        e, f = robustdiv(b, a, d, c)
        f = -f
    end
    z = Complex{DoubleFloat{T}}(e, f)
    return z
end

function robustdiv(a::DF, b::DF, c::DF, d::DF) where {T, DF<:DoubleFloat{T}}
    r = d / c
    t = inv( fma(r, d, c) )
    e = robustcomp(a,b,c,d,r,t)
    a = -a
    f = robustcomp(b,a,c,d,r,t)
    return e, f
end

function robustcomp(a,b,c,d,r,t)
    if !iszero(r)
        if !iszero(b) & !iszero(r)
            e = fma(b, r, a) * t
        else
            e = a * t + (b * t) * r
        end
    else
        e = fma(d, (b/c), a) * t
    end
    return e
end


text: # This file is a part of UltraNest.jl, licensed under the MIT License (MIT).


text:
"""
    getSTRdiff(datafile, haplotype_size; showhaps=false, show="full_pedigree")

STRdiff is a program to extract and phase genotypes in the
region surrounding a de novo STR event found in the CEPH pedigrees. All
possible trios from a pedigree are phased and a consensus chromosome from the
children bearing the de novo STR is used to identify any difference in the
STR allele length between grandparent and parent.  A table of differences,
based on local chromosomal phasing, and the sex of the transmitting parent
is returned.

Input file format example: (single spaced, tab-delimited, header optional)\n
#PEDIGREE       VCFFILE                IID     STRLOC         STRTYPE(bp)\n
1463.ped        1463\\_phase_in.vcf.gz    2188    2:174994186    2\n
1463.ped        1463\\_phase_in.vcf.gz    2188    1:2581321      4\n

This version requires that the pedigree has two parents and grandparents
and that the de novo event occurred in one of the CEPH parents.

An example command: getSTRdiff("myGenotypes.vcf.gz", 50000)

**Important:**\n
The pedigree file **must** have the grandparents and parent samples
in the order described below.

Paternal trio

line1: paternal grandfather

line2: paternal grandmother

line3: father

Maternal trio

line4: maternal grandfather

line5: maternal grandmother

line6: mother

then all the children.

For example (single spaced, tab-delimited),

1358    8542    0    0    1    0\n
1358    8544    0    0    2    0\n
1358    8538    8542    8544    1    0\n
1358    8541    0    0    1    0\n
1358    8543    0    0    2    0\n
1358    8540    8541    8543    2    0\n
1358    8534    8538    8540    1    0\n

additional children ...

The haplotype carrying the novel STR allele in the offspring
is compared with the haplotypes of the grandparents to
determine which grandparental chromosome most likely gave rise
to the new STR allele.

It is strongly recommended that a given STR be analyzed at different
haplotype lengths and that only runs with a unique solution score of >10
be considered for additional analysis.

Use the repSTRdiff function to automate multiple runs.

During testing, we documented a small number cases where beagle produced
offspring with different haplotypes depending on the input order of the
children in the pedigree file.  This beagle-related issue appeared only
at a small number of STR loci.

Dependencies: bcftools, bgzip
"""
function getSTRdiff(infile, hapsize::Int64; kwargs...)

kw = Dict(kwargs)
sh = get(kw, :show, nothing)
s_haps = get(kw, :showhaps, false)

fname = splitext(infile)
fn = fname[1]
loc_list = String[]

if (hapsize < 10000) || (hapsize > 500000)
    error("hapsizes, bp on each side of the STR, should be between 10000 and 500000.")
end

OUT = open("$fn.$hapsize.out", "w")
header =  "#ped\tiid\tloc\ttranmitting_grandparent\ttransmitting_grandparent_sex\ttransmitting_grandparent_str\ttransmitting_grandparent_str_size\tparent_str\tparent_str_size\tstr_allele_change\tsize_change_bp\toffspring_with_de_novo_allele\tgp_match_probability\tuniq_soln_score\thapsize\tstr_type"
println(OUT, header)
close(OUT)

open(infile) do file

    for line in eachline(file)

        if occursin(r"^#", line)
            continue
        end

        l = '-'^80

        println("\n$l")
        printstyled("START RUN\nINFO: Processing $line\nINFO: Haplotype size = $hapsize\n", color=:cyan, bold=true)
        println("$l")

        z = preprocessDenovoSTRs(line, window=hapsize)

        datin = split(line, '\t')
        if length(datin) != 5
            error("expecting 5 fields for each STR location,\nfound line: $datin")
        end
        pedname = splitext(datin[1])
        ped = pedname[1]
        shf = datin[1]
        did = datin[3]
        expbp = datin[5]

        sex = hashfile("$shf", 2, [5])
        fathers = hashfile("$shf", 2, [3])

        f1 = findmax(countmap(collect(values(fathers))))[2][1]
        mothers = hashfile("$shf", 2, [4])
        m1 = findmax(countmap(collect(values(mothers))))[2][1]
        othpar = ""

        if did == f1
            othpar = m1
        elseif did == m1
            othpar = f1
        end

        println("INFO: Parent with de novo STR is $did")
        println("INFO: Other parent id is $othpar")

        fa = z[1][1]; fb = z[1][2]; site = z[1][3]; iid = z[1][4]
        ta = z[1][5]; tb = z[1][6]; mutSTR = z[1][7]; par_genos = z[1][8]

        (p, q, r, u, mutCon, fa_strs, num_iht_mut) = STR_iht(fa, fb, site, iid, mutSTR, othpar; showhaps=s_haps)

        (tp, tq, tr, tu, mutCon_empty, t_strs, num_iht_mut_null) = STR_iht(ta, tb, site, iid, mutSTR, othpar; showhaps=s_haps)

#println("%%%%%$p  $q  $r  $u  $mutCon  $fa_strs $num_iht_mut \n\n")

        if (mutCon != nothing) && ( size(tr, 1) != length(mutCon) )
            error("Variant site $site has failed for hapsize $hapsize.
Try a different haplotype size or consider removing this site
if it fails at multiple hapsizes.")
            exit()
        end

        if num_iht_mut == -9

            na = "na"; rr = -9
            OUT = open("$fn.$hapsize.out", "a+")
            println(OUT, "$ped\t$iid\t$site\t$na\t$na\t$na\t$na\t$na\t$na\t$na\t$na\t$num_iht_mut\t$na\t$rr\t$hapsize\t$expbp")
            close(OUT)
            println("\n$l\n$l")
            continue

        else

            match_p = []

            for i in 3:6
                p = sum(tr[!,i] .== mutCon)/length(mutCon)
                push!(match_p, p)
            end

            parent_idx = rand(findall( match_p .== maximum(match_p) ) ) + 2

            match_p_sorted = sort(match_p, rev=true)

            mdiff = round(match_p_sorted[1] - match_p_sorted[2], digits=2) * 100

            uss = Int64(floor(mdiff))

            rflag = ""

            if length(findall( match_p .== maximum(match_p) ) ) > 1
                rflag = "rflag"
                uss = 0
            end

            trans_p = names(tr,parent_idx)  #transmitting parent id

            trans_ps = trans_p[1]  #convert to plain string for patmatch

            #Test transmitting parent against parents genotypes and get
            #the alleles.  Test homs then hets

            psize = []
            pg = []

            if occursin(par_genos[1], trans_ps) && (par_genos[2] == par_genos[3])

                push!(pg, parse(Int, par_genos[2]))
                push!(psize, get(t_strs, pg[1], 0))

            elseif occursin(par_genos[4], trans_ps) && (par_genos[5] == par_genos[6])

                push!(pg, parse(Int, par_genos[5]))
                push!(psize, get(t_strs, pg[1], 0))

            elseif occursin(par_genos[1], trans_ps) && (par_genos[2] != par_genos[3])

                pg_1 = parse(Int, par_genos[2])
                pg_2 = parse(Int, par_genos[3])
                push!(pg, pg_1, pg_2)
                push!(psize, get(t_strs, pg_1, 0))
                push!(psize, get(t_strs, pg_2, 0))

            elseif occursin(par_genos[4], trans_ps) && (par_genos[5] != par_genos[6])

                pg_1 = parse(Int, par_genos[5])
                pg_2 = parse(Int, par_genos[6])
                push!(pg, pg_1, pg_2)

                push!(psize, get(t_strs, pg_1, 0))
                push!(psize, get(t_strs, pg_2, 0))

            else

                println("ERROR: unable to uniquely determine all phased STR
alleles transmitting the de novo event, please vary the hapsize value.")
                rr = -9
                na = "na"
                println("$par_genos $trans_ps")
                OUT = open("$fn.$hapsize.out", "a+")
                println(OUT, "$ped\t$iid\t$site\t$na\t$na\t$na\t$na\t$na\t$na\t$na\t$na\t$na\t$na\t$rr\t$hapsize\t$expbp")
                close(OUT)
                continue

            end

            dsize = get(fa_strs, mutSTR, 0)  #de novo size

            if length(psize) > 1

                dsize_diff1 = dsize - psize[1]
                dsize_diff2 = dsize - psize[2]
                str_s_diff1 = mutSTR - pg[1]
                str_s_diff2 = mutSTR - pg[2]

                if abs(dsize_diff1) == abs(dsize_diff2)
                    dsize_diff = join(["\u00B1", abs(dsize_diff1)], "")
                elseif abs(dsize_diff1) != abs(dsize_diff2)
                    dsize_diff = NaN
                end

                if abs(str_s_diff1) == abs(str_s_diff2)
                    str_s_diff = join(["\u00B1", abs(str_s_diff1)], "")
                elseif abs(str_s_diff1) != abs(str_s_diff2)
                    str_s_diff = NaN
                end

            else

                dsize_diff = dsize - psize[1]
                str_s_diff = mutSTR - pg[1]

            end

            trans_pid = trans_ps[1:end-1]
            trans_ps_sex = get(sex, trans_pid, 0)
            tps = trans_ps_sex[1]

            function namesex(x)
                tps = ""
                if x == "2"
                    tps = "female"
                elseif x == "1"
                    tps = "male"
                else
                    tps = "unknown"
                end
                return tps
            end

            if sh == "full_table"
                println("TABLE: predicted inheritance patterns given the inferred haplotypes ...")
                println(q, "\n")
            elseif sh == "full_pedigree"
                println("TABLE: haplotypes for the pedigree ...")
                println(r, "\n", tr)
            end

            mprobs = round.(match_p, digits=6)
            tps_sex = namesex(tps)
            println()
            println("Grandparent-offspring match probabilites: ", round.(match_p, digits=4) )
            OUT=open("matchprobs.txt", "a+")
                println(OUT, hapsize, "\t", round.(match_p, digits=4), "\t", uss)
            close(OUT)
            println("Grandparent transmitting the de novo allele: $trans_ps  $rflag")
            println("Transmitting grandparent's sex: $tps_sex")
            println("Unique solution score: $uss")
            println("Number of de novo alleles found in phased offspring chromosomes: $num_iht_mut")

            ppg = join(pg, ",")
            ppsize = join(psize, ",")

            println("Grandparental STR allele and size: $ppg and $ppsize bp")
            println("Parental de novo STR allele and size: $mutSTR and $dsize bp")
            println("De novo STR allele size difference (bp): $dsize_diff bp")

            OUT = open("$fn.$hapsize.out", "a+")
            println(OUT, "$ped\t$iid\t$site\t$trans_ps\t$tps_sex\t$ppg\t$ppsize\t$mutSTR\t$dsize\t$str_s_diff\t$dsize_diff\t$num_iht_mut\t$mprobs\t$uss\t$hapsize\t$expbp")
            close(OUT)

            println("$l\n$l")
        end

        run(`bash -c "rm -f *fam*"`)
        run(`bash -c "rm -f *trio*"`)

    end

end

end

max_stars_repo_name: ScottWatkins/STRdiff

max_issues_repo_name: ScottWatkins/STRdiff

max_forks_repo_name: ScottWatkins/STRdiff


text:
using Plots
using Plots.PlotMeasures
using LambertW
import PlotlyBase

plotly()
Plots.PlotlyBackend()

Plots.scalefontsizes()
Plots.scalefontsizes(1.3)

R = 0:0.01:1;

## No-repair reliability
reliability_single = R;
reliability_dr = @. 1 - (1 - R)^2;
#reliability_3r = @. 1 - (1 - R)^3;
reliability_tmr = @. 3 * R^2 - 2 * R^3

plot(
    R,
    [reliability_single, reliability_dr, reliability_tmr],
    label = ["Απλό εξάρτημα" "Διπλός πλεονασμός" "Τριπλός πλεονασμός με ψηφοφορία"],
    # title = "Ανθεκτικότητα σύνθετου συστήματος σε μόνιμες βλάβες",
    legend = :outerright,
    xlabel = "Αξιοπιστία εξαρτήματος",
    ylabel = "Αξιοπιστία συστήματος",
    lw = 2,
)
savefig("reliability_norepair.pdf")

## No-repair reliability with vulnerable voter
Rv = 0:0.01:1;
reliability_tmrv = reliability_tmr * permutedims(Rv)

heatmap(
    R,
    Rv,
    reliability_tmrv,
    xlabel = "Αξιοπιστία R εξαρτήματος",
    ylabel = "Αξιοπιστία Rv ψηφοφόρου",
    colorbar_title = "Αξιοπιστία συστήματος",
    aspect_ratio = 1,
    c = :inferno
)
savefig("reliability_norepair_voter.pdf")

## Ρrepair reliability
λ = 10 .^ (LinRange(-11, -2, 1000))

λ_single = λ;
λ_dr = @. -lambertw(- (2λ * exp(-λ) - λ^2 * exp(-2λ)));
λ_tmr = @. -lambertw(- (3(λ * exp(-λ))^2 - 2(λ * exp(-λ))^3));

plot(
    λ,
    [λ_single, λ_single, λ_dr, λ_tmr],
    label = permutedims([
        "Απλό εξάρτημα"
        "Διπλός πλεονασμός, ψυχρός"
        "Διπλός πλεονασμός, ενεργός"
        "Τριπλός πλεονασμός με ψηφοφορία"]),
    #title = "Ανθεκτικότητα σύνθετου συστήματος σε παροδικά σφάλματα",
    legend = (0.6,0.35),
    xlabel = "Ρυθμός αποτυχίας εξαρτήματος",
    ylabel = "Ρυθμός αποτυχίας συστήματος",
    yaxis = :log,
    xaxis = :log,
    left_margin = -150px,
    lw = 2,
    linecolor=[palette(:tab10)[1] palette(:tab10)[1] palette(:tab10)[2] palette(:tab10)[3]]
)
savefig("reliability_repair.pdf")


text:
#A simple DSV (delimiter-seperated values) file reader
#-------------------------------------------------------------------------------
module DSVFiles
#==Implements a basic character-delimited text reader.
Uses somewhat impractical high-level tools (readline/split) to provide a more
traditional handling of character-delimited files.

#TODO: Doesn't enforce that CSV data MUST HAVE a "," between values.
==#

import Base: open, read, write, close, readline, eof

#==Constants
===============================================================================#
const splitter_default = [' ', '\t']
const splitter_csv = [' ', '\t', ',']

#==Main data structures
===============================================================================#
mutable struct DSVReader
	s::IO
	splitter
	linebuf::Vector
end

#Default splitter value:
(RT::Type{T})(s::IO, splitter = splitter_default) where T<:DSVReader =
	RT(s, splitter, String[])

#==Helper functions
===============================================================================#
#Advance DSVReader linebuf, if necessary:
#TODO: Re-implement with more efficient, lower level functions
function refreshbuf(r::DSVReader)
	while length(r.linebuf) < 1
		if eof(r.s); throw(EOFError()); end
		linebuf = split(chomp(readline(r.s)), r.splitter)

		for token in linebuf
			if token != ""
				push!(r.linebuf, token)
			end
		end
	end
end

#==Open/read/close functions
===============================================================================#

open(RT::Type{T}, path::String, splitter = splitter_default) where T<:DSVReader =
	RT(open(path, "r"), splitter)

#Read in entire text file as string
function read(RT::Type{T}, path::String) where T<:DSVReader
	open(RT, path) do reader
		return read(reader.s, String)
	end
end

#Read in next token as String:
function read(r::DSVReader, ::Type{DT}) where DT<:AbstractString
	refreshbuf(r)
	v = popfirst!(r.linebuf)
	try
		#Read ahead if necessary (Update state for eof()):
		refreshbuf(r)
	catch
	end
	return v
end

#Read in next token & interpret as of type DT:
function read(r::DSVReader, ::Type{DT}) where DT
	return parse(DT, read(r, String))
end

#Read in next token & interpret as most appropriate type:
function read(r::DSVReader, ::Type{Any})
	return parse(read(r, String))
end

close(r::DSVReader) = close(r.s)

#==Support functions
===============================================================================#
eof(r::DSVReader) = (eof(r.s) && length(r.linebuf) < 1)

function readline(r::DSVReader)
	linebuf = []
	return readline(r.s)
end

function read(r::DSVReader, String)
	linebuf = []
	return read(r.s, String)
end

#==High-level interface
===============================================================================#

#Default delimiter-seperated value:
open_dsv(path::String, splitter = splitter_default) =
	open(DSVReader, path, splitter = splitter)
open_csv(path::String) =
	open(DSVReader, path, splitter = splitter_csv)

end
#Last line


text:
"""
    ESh(p1::Float64,p2::Float64)

Produces the effect size `h` for two proportions. `p1` and `p2`
are two proportions. `h` is 2*asin(sqrt(p1))-2*asin(sqrt(p2)).
"""
function ESh(p1::Float64,p2::Float64)
    return 2*asin(sqrt(p1)) - 2*asin(sqrt(p2))
end

max_stars_repo_name: mwsohn/pwr

max_issues_repo_name: mwsohn/pwr

max_forks_repo_name: mwsohn/pwr


text:
"Holds the information for the various methods' tableaus."
abstract Tableau{T}

"Holds the tableau of a Runge-Kutta method."
abstract TableauRK{T} <: Tableau{T}

"Holds the tableau of an implicit Runge-Kutta method."
abstract TableauIRK{T} <: TableauRK{T}

"Holds the tableau of a partitioned Runge-Kutta method."
abstract TableauPRK{T} <: TableauRK{T}

@define HeaderTableauRK begin
    name::Symbol
    o::Int
    s::Int
    a::Matrix{T}
    b::Vector{T}
    c::Vector{T}
end

@define HeaderTableauPRK begin
    name::Symbol
    o::Int
    s::Int
    a_q::Matrix{T}
    a_p::Matrix{T}
    b_q::Vector{T}
    b_p::Vector{T}
    c_q::Vector{T}
    c_p::Vector{T}
end

Base.hash(tab::TableauRK, h::UInt) = hash(tab.o, hash(tab.a, hash(tab.b, hash(tab.c, hash(:TableauRK, h)))))
Base.:(==){T1, T2}(tab1::TableauRK{T1}, tab2::TableauRK{T2}) = (tab1.o == tab2.o
                                                             && tab1.s == tab2.s
                                                             && tab1.a == tab2.a
                                                             && tab1.b == tab2.b
                                                             && tab1.c == tab2.c)

Base.isequal{T1, T2}(tab1::TableauRK{T1}, tab2::TableauRK{T2}) = (tab1 == tab2 && T1 == T2 && typeof(tab1) == typeof(tab2))

"Print Runge-Kutta tableau to standard output."
function showTableau{T}(tab::TableauRK{T})
    println("Runge-Kutta Method ", tab.name, "with ", tab.s, " stages and order ", tab.o)
    println("  a = ", tab.a)
    println("  b = ", tab.b)
    println("  c = ", tab.c)
end

"Reads and parses Tableau metadata from file."
function readTableauRKHeaderFromFile(file)
    f = open(file, "r")
    header = readline(f)
    close(f)

    if header[1] == '#'
        header = split(header[2:end])
    else
        header = ()
    end

    if length(header) ≥ 1
        O = parse(Int, header[1])
    else
        O = 0
    end

    if length(header) ≥ 2
        S = parse(Int, header[2])
    else
        S = 0
    end

    if length(header) ≥ 3
        T = eval(parse(header[3]))
    else
        T = Float64
    end

    return O, S, T
end

"Write Runge-Kutta tableau to file."
function writeTableauToFile{T}(dir::AbstractString, tab::TableauRK{T})
    # tab_array = zeros(T, S+1, S+1)
    # tab_array[1:S, 2:S+1] = tab.a
    # tab_array[S+1, 2:S+1] = tab.b
    # tab_array[1:S, 1] = tab.c
    # tab_array[S+1, 1] = tab.order

    tab_array = zeros(T, tab.s, tab.s+2)
    tab_array[1:tab.s, 1] = tab.c
    tab_array[1:tab.s, 2] = tab.b
    tab_array[1:tab.s, 3:tab.s+2] = tab.a

    header = string("# ", tab.o, " ", tab.s, " ", T, "\n")
    file   = string(dir, "/", tab.name, ".tsv")

    info("Writing Runge-Kutta tableau ", tab.name, " with ", tab.s, " stages and order ", tab.o, " to file\n      ", file)

    f = open(file, "w")
    write(f, header)
    writedlm(f, float(tab_array))
    close(f)
    # TODO Write data in original format (e.g., Rational).
end

# TODO function writeTableauToFile{Name, T}(dir::AbstractString, tab::TableauPRK{Name, T})
# TODO function writeTableauToFile{Name, T}(dir::AbstractString, tab::TableauSPARK{Name, T})
# TODO function writeTableauToFile{Name, T}(dir::AbstractString, tab::TableauGLM{Name, T})

"Holds the tableau of an explicit Runge-Kutta method."
immutable TableauERK{T} <: TableauRK{T}
    @HeaderTableauRK

    function TableauERK(name,o,s,a,b,c)
        @assert T <: Real
        @assert isa(name, Symbol)
        @assert isa(s, Integer)
        @assert isa(o, Integer)
        @assert s > 0 "Number of stages must be > 0"
        @assert s==size(a,1)==size(a,2)==length(b)==length(c)
        @assert c[1] == 0
        @assert istrilstrict(a)
        @assert !(s==1 && a[1,1] ≠ 0)
        new(name,o,s,a,b,c)
    end
end

function TableauERK{T}(name::Symbol, order::Int, a::Matrix{T}, b::Vector{T}, c::Vector{T})
    TableauERK{T}(name, order, length(c), a, b, c)
end

"Read explicit Runge-Kutta tableau from file."
function readTableauERKFromFile(dir::AbstractString, name::AbstractString)
    file = string(dir, "/", name, ".tsv")

#    run(`cat $file`)

    o, s, T = readTableauRKHeaderFromFile(file)

    # TODO Read data in original format (e.g., Rational).
    #      For this we need to save tableaus as jld or hdf5.
#    tab_array = readdlm(file, T)
    tab_array = readdlm(file)

    if s == 0
        s = size(tab_array, 1)
    end

    @assert s == size(tab_array, 1) == size(tab_array, 2)-2

    c = tab_array[1:s, 1]
    b = tab_array[1:s, 2]
    a = tab_array[1:s, 3:s+2]

    info("Reading explicit Runge-Kutta tableau ", name, " with ", s, " stages and order ", o, " from file\n      ", file)

    TableauERK(Symbol(name), o, a, b, c)
end

"Holds the tableau of a diagonally implicit Runge-Kutta method."
immutable TableauDIRK{T} <: TableauIRK{T}
    @HeaderTableauRK

    function TableauDIRK(name,o,s,a,b,c)
        @assert T <: Real
        @assert isa(name, Symbol)
        @assert isa(s, Integer)
        @assert isa(o, Integer)
        @assert s > 0
        @assert s==size(a,1)==size(a,2)==length(b)==length(c)
        @assert istril(a)
        @assert !(s==1 && a[1,1] ≠ 0)

        if s > 1 && istrilstrict(a)
            warn("Initializing TableauDIRK with explicit tableau ", name, ".\n",
                 "You might want to use TableauERK instead.")
        end

        new(name,o,s,a,b,c)
    end
end

function TableauDIRK{T}(name::Symbol, order::Int, a::Matrix{T}, b::Vector{T}, c::Vector{T})
    TableauDIRK{T}(name, order, length(c), a, b, c)
end

# TODO function readTableauDIRKFromFile(dir::AbstractString, name::AbstractString)

"Holds the tableau of a fully implicit Runge-Kutta method."
immutable TableauFIRK{T} <: TableauIRK{T}
    @HeaderTableauRK

    function TableauFIRK(name,o,s,a,b,c)
        @assert T <: Real
        @assert isa(name, Symbol)
        @assert isa(s, Integer)
        @assert isa(o, Integer)
        @assert s > 0 "Number of stages must be > 0"
        @assert s==size(a,1)==size(a,2)==length(b)==length(c)

        if (s > 1 && istrilstrict(a)) || (s==1 && a[1,1] == 0)
            warn("Initializing TableauFIRK with explicit tableau ", name, ".\n",
                 "You might want to use TableauERK instead.")
        elseif s > 1 && istril(a)
            warn("Initializing TableauFIRK with diagonally implicit tableau ", name, ".\n",
                 "You might want to use TableauDIRK instead.")
        end

        new(name,o,s,a,b,c)
    end
end

function TableauFIRK{T}(name::Symbol, order::Int, a::Matrix{T}, b::Vector{T}, c::Vector{T})
    TableauFIRK{T}(name, order, length(c), a, b, c)
end

# TODO function readTableauFIRKFromFile(dir::AbstractString, name::AbstractString)

"Holds the tableau of a singly implicit Runge-Kutta method."
immutable TableauSIRK{T} <: TableauIRK{T}
    # TODO
end

function TableauSIRK{T}(name::Symbol, order::Int, a::Matrix{T}, b::Vector{T}, c::Vector{T})
    TableauSIRK{T}(name, order, length(c), a, b, c)
end

"Holds the tableau of an explicit partitioned Runge-Kutta method."
# TODO Need explicit and implicit version?
immutable TableauEPRK{T} <: Tableau{T}
    @HeaderTableauPRK

    function TableauEPRK(name, o, s, a_q, a_p, b_q, b_p, c_q, c_p)
        @assert T <: Real
        @assert isa(name, Symbol)
        @assert isa(s, Integer)
        @assert isa(o, Integer)
        @assert s > 0 "Number of stages must be > 0"
        @assert s==size(a_q,1)==size(a_q,2)==length(b_q)==length(c_q)
        @assert s==size(a_p,1)==size(a_p,2)==length(b_p)==length(c_p)
        # TODO check that both tableaus are lower triangular and that only one element
        #      a_q[i,i] or a_p[i,i] is non-zero for all i.
        new(name, o, s, a_q, a_p, b_q, b_p, c_q, c_p)
    end
end

function TableauEPRK{T}(name::Symbol, order::Int,
                        a_q::Matrix{T}, a_p::Matrix{T},
                        b_q::Vector{T}, b_p::Vector{T},
                        c_q::Vector{T}, c_p::Vector{T})
    @assert length(c_q) == length(c_p)
    TableauEPRK{T}(name, order, length(c_q), a_q, a_p, b_q, b_p, c_q, c_p)
end

function TableauEPRK{T}(name::Symbol, order::Int, tab::TableauERK{T})
    TableauEPRK{T}(name, order, tab.s, tab.a, tab.a, tab.b, tab.b, tab.c, tab.c)
end

function TableauEPRK{T}(name::Symbol, order::Int, tab_q::TableauRK{T}, tab_p::TableauRK{T})
    @assert tab_q.s == tab_p.s
    TableauEPRK{T}(name, order, tab_q.s, tab_q.a, tab_p.a, tab_q.b, tab_p.b, tab_q.c, tab_p.c)
end

# TODO function readTableauPRKFromFile(dir::AbstractString, name::AbstractString)

"Holds the tableau of a implicit partitioned Runge-Kutta method."
immutable TableauIPRK{T} <: Tableau{T}
    @HeaderTableauPRK

    function TableauIPRK(name, o, s, a_q, a_p, b_q, b_p, c_q, c_p)
        @assert T <: Real
        @assert isa(name, Symbol)
        @assert isa(s, Integer)
        @assert isa(o, Integer)
        @assert s > 0 "Number of stages must be > 0"
        @assert s==size(a_q,1)==size(a_q,2)==length(b_q)==length(c_q)
        @assert s==size(a_p,1)==size(a_p,2)==length(b_p)==length(c_p)

        new(name, o, s, a_q, a_p, b_q, b_p, c_q, c_p)
    end
end

function TableauIPRK{T}(name::Symbol, order::Int,
                        a_q::Matrix{T}, a_p::Matrix{T},
                        b_q::Vector{T}, b_p::Vector{T},
                        c_q::Vector{T}, c_p::Vector{T})
    @assert length(c_q) == length(c_p)
    TableauIPRK{T}(name, order, length(c_q), a_q, a_p, b_q, b_p, c_q, c_p)
end

function TableauIPRK{T}(name::Symbol, order::Int, tab_q::TableauRK{T}, tab_p::TableauRK{T})
    @assert tab_q.s == tab_p.s
    TableauIPRK{T}(name, order, tab_q.s, tab_q.a, tab_p.a, tab_q.b, tab_p.b, tab_q.c, tab_p.c)
end

# TODO function readTableauIPRKFromFile(dir::AbstractString, name::AbstractString)

"Holds the tableau of a spezialized additive Runge-Kutta method."
immutable TableauSARK{T} <: Tableau{T}
    name::Symbol
    o::Int
    s::Int
    r::Int

    a_q::Matrix{T}
    α_q::Matrix{T}

    a_q̃::Matrix{T}
    α_q̃::Matrix{T}

    b_q::Vector{T}
    β_q::Vector{T}

    c_q::Vector{T}
    c_λ::Vector{T}

    ω_q::Matrix{T}
    ω_λ::Matrix{T}

    function TableauSARK(name, o, s, r,
                         a_q, α_q, a_q̃, α_q̃,
                         b_q, β_q, c_q, c_λ,
                         ω_q, ω_λ)
        # TODO Make ω_q, ω_λ optional arguments.
        @assert T <: Real
        @assert isa(name, Symbol)
        @assert isa(s, Integer)
        @assert isa(r, Integer)
        @assert isa(o, Integer)
        @assert s > 0 "Number of stages s must be > 0"
        @assert r > 0 "Number of stages r must be > 0"
        @assert s==size(a_q,1)==size(a_q,2)==length(b_q)==length(c_q)
        @assert s==size(α_q,1)==length(β_q)
        @assert r==size(α_q,2)==size(α_p,2)
        @assert r==length(c_λ)
        @assert r==size(a_qᵠ,1)==size(α_qᵠ,1)==size(α_qᵠ,2)
        @assert s==size(a_qᵠ,2)
        # TODO Add assertions on ω_q, ω_λ to be (S-1)x(S) or (R-1)x(R) if set.
        new(name, o, s, r, a_q, α_q, a_q̃, α_q̃, b_q, β_q, c_q, c_λ, ω_q, ω_λ)
    end
end

# TODO Add external constructor for TableauSARK.

# TODO function readTableauSARKFromFile(dir::AbstractString, name::AbstractString)

"Holds the tableau of a spezialized partitioned additive Runge-Kutta method."
immutable TableauSPARK{T} <: Tableau{T}
    name::Symbol
    o::Int
    s::Int
    r::Int

    a_q::Matrix{T}
    a_p::Matrix{T}
    α_q::Matrix{T}
    α_p::Matrix{T}

    a_q̃::Matrix{T}
    a_p̃::Matrix{T}
    α_q̃::Matrix{T}
    α_p̃::Matrix{T}

    b_q::Vector{T}
    b_p::Vector{T}
    β_q::Vector{T}
    β_p::Vector{T}

    c_q::Vector{T}
    c_p::Vector{T}
    c_λ::Vector{T}

    ω_q::Matrix{T}
    ω_p::Matrix{T}
    ω_λ::Matrix{T}

    function TableauSPARK(name, o, s, r,
                          a_q, a_p, α_q, α_p, a_q̃, a_p̃, α_q̃, α_p̃,
                          b_q, b_p, β_q, β_p,
                          c_q, c_p, c_λ,
                          ω_q, ω_p, ω_λ)
        # TODO Make ω_q, ω_p, ω_λ optional arguments.
        @assert T <: Real
        @assert isa(name, Symbol)
        @assert isa(s, Integer)
        @assert isa(r, Integer)
        @assert isa(o, Integer)
        @assert s > 0 "Number of stages s must be > 0"
        @assert r > 0 "Number of stages r must be > 0"
        @assert s==size(a_q,1)==size(a_q,2)==length(b_q)==length(c_q)
        @assert s==size(a_p,1)==size(a_p,2)==length(b_p)==length(c_p)
        @assert s==size(α_q,1)==size(α_p,1)==length(β_q)==length(β_p)
        @assert r==size(α_q,2)==size(α_p,2)
        @assert r==length(c_λ)
        @assert r==size(a_qᵠ,1)==size(a_pᵠ,1)
        @assert r==size(α_qᵠ,1)==size(α_qᵠ,2)
        @assert r==size(α_pᵠ,1)==size(α_pᵠ,2)
        @assert s==size(a_qᵠ,2)==size(a_pᵠ,2)
        # TODO Add assertions on ω_q, ω_p, ω_λ to be (S-1)x(S) or (R-1)x(R) if set.
        new(name, o, s, r,
            a_q, a_p, α_q, α_p, a_q̃, a_p̃, α_q̃, α_p̃,
            b_q, b_p, β_q, β_p,
            c_q, c_p, c_λ,
            ω_q, ω_p, ω_λ)
    end
end

# TODO Add external constructor for TableauSPARK.

# TODO function readTableauSPARKFromFile(dir::AbstractString, name::AbstractString)

"Holds the tableau of a general linear method."
immutable TableauGLM{T} <: Tableau{T}
    name::Symbol
    o::Int
    s::Int
    r::Int

    a::Matrix{T}
    b::Matrix{T}
    u::Matrix{T}
    v::Matrix{T}
    c::Vector{T}

    function TableauGLM(name, o, s, r, a, b, u, v, c)
        @assert T <: Real
        @assert isa(name, Symbol)
        @assert isa(s, Integer)
        @assert isa(r, Integer)
        @assert isa(o, Integer)
        @assert s > 0 "Number of stages s must be > 0"
        @assert r > 0 "Number of stages r must be > 0"
        @assert s==length(c)
        @assert s==size(a,1)==size(a,2)==length(u,1)==length(b,2)
        @assert r==size(v,1)==size(v,2)==length(u,2)==length(b,1)
        new(name, o, s, r, a, b, u, v, c)
    end
end

# TODO Add external constructor for TableauGLM.

# TODO function readTableauGLMFromFile(dir::AbstractString, name::AbstractString)

# TODO Add TableauSGALM. # Super General Additive Linear Method
# TODO Add TableauTSRK.  # Two Step Runge Kutta Method
# TODO Add TableauSplitting.


text:
function update_constraints_by_comparison!(secret_word,candidate_word, let_in_pos, let_not_in_pos, let_not_in_word)
    for i in eachindex(candidate_word)
        if secret_word[i] == candidate_word[i]
            let_in_pos[i]=candidate_word[i]
        end
    end
    for i in eachindex(candidate_word)
        if candidate_word[i] ∉ values(let_in_pos) && occursin(candidate_word[i], secret_word)
            let_not_in_pos[candidate_word[i]] = push!(get(let_not_in_pos,candidate_word[i],Int[]),i)
        end
    end
    for i in eachindex(candidate_word)
        if !occursin(candidate_word[i],secret_word)
          push!(let_not_in_word,candidate_word[i])
        end
    end
end

function firstwordscore_simple(words,trials)
    score=Dict(zip(words,[(0,0,0) for i in eachindex(words)]))
    p = Progress(trials)
    for  i=1:trials
        selected_word = rand(words)
        Threads.@threads for w in words
            w == selected_word && continue #No point attempting with same word
            i = steps_to_solve(words,selected_word,w)
            if i < 6
                score[w] = score[w] .+ (i,1,0)
            else
                score[w] = score[w] .+ (0,0,1)
            end
        end
        next!(p)
    end
    score
end

function firstwordscore_prob(words,freqs,trials)
    score=Dict(zip(words,[(0,0,0) for i in eachindex(words)]))
    p = Progress(trials)
    weights = FrequencyWeights(freqs)
    w_f_dict = Dict(zip(words,freqs))
    for  i=1:trials
        selected_word = sample(words,weights)
        Threads.@threads for w in words
            w == selected_word && continue #No point attempting with same word
            i = steps_to_solve_prob(words, w_f_dict ,selected_word,w)
            if i < 6
                score[w] = score[w] .+ (i,1,0)
            else
                score[w] = score[w] .+ (0,0,1)
            end
        end
        next!(p)
    end
    score
end

"""
    gen_word_reducer_by_target(target_word, guess_word)

Returns a closure that captures the constraints arounds and updates
the `word_set` given and word attempt. This propagation is done automatically
for a given `target_word` which the word that the puzzle solver is trying to guess.

The colsure function has the following signature

    reduce_wordset_by_word_target!(word_set, word_attempt)

where the `word_attempt` is a guess at the `target_word`.
"""
function gen_word_reducer_by_target(target_word, guess_word)
    let_in_pos = Dict{Int,Char}()
    let_not_in_pos = Dict{Char,Vector{Int}}()
    let_not_in_word = Set{Char}()
    w = guess_word
    function reduce_wordset_by_word_target!(word_set, word_attempt)
        update_constraints_by_comparison!(target_word,word_attempt,let_in_pos,let_not_in_pos,let_not_in_word)
        filter!(w->all(w[r[1]]==r[2] for r in let_in_pos), word_set)
        filter!(w->all(occursin(s[1],w) && all(w[p]!=s[1] for p in s[2]) for s in let_not_in_pos), word_set)
        filter!(w->all(!occursin(c,w) for c in let_not_in_word), word_set)
    end
end

function steps_to_solve_prob(word_set, w_f_dict, target_word, guess_word)
    word_set_reducer! = gen_word_reducer_by_target(target_word, guess_word)
    word_attempt = guess_word
    word_set = copy(word_set)
    for i=1:6
        word_attempt == target_word && return i
        i==6  && return i
        word_set = word_set_reducer!(word_set,word_attempt)
        freqs =  [w_f_dict[w] for w in word_set]
        word_attempt = sample(word_set,FrequencyWeights(freqs))
    end
end

function steps_to_solve(word_set, target_word, guess_word)
    word_set_reducer! = gen_word_reducer_by_target(target_word, guess_word)
    word_attempt = guess_word
    word_set = copy(word_set)
    for i=1:6
        word_attempt == target_word && return i
        i==6  && return i
        word_set = word_set_reducer!(word_set,word_attempt)
        word_attempt = rand(word_set)
    end
end

"""
    calc_score_probs(words,scores)
Calcualtes the depth and safety scores for any array or and their assocaited dictionary
of that contains the expermental accumelation of success and failures.
"""
function calc_score_probs(words,scores)
    depth_score = zeros(length(words))
    safety_score = zeros(length(words))
    for i in eachindex(words)
        depth_score[i] = scores[words[i]][1]/scores[words[i]][2]
        safety_score[i] = scores[words[i]][3]/(scores[words[i]][2]+scores[words[i]][3])
    end
    depth_score, safety_score
end


text:
function test_addExp()
	# sequence of arguments should not matter:
	@test addExp(3.0, 5.0) ≈ addExp(5.0, 3.0)
	@test addExp(3.0, 5.0) ≈ log(exp(3)+exp(5))
end

test_addExp()


text:
using ARFFFiles, DataFrames
using Random

"""
    Normalize(vector)
Devuelve un vector de datos normalizados
"""
function Normalize(vector)
    max = maximum(vector)
    min = minimum(vector)
    diff = (max-min)
    reescale = diff ≈ 0 ? 1 : diff

    return map( x -> (x - min)/reescale, vector)
end
"""
    DataLabel(file)
Dado un fichero .arff devuleve:
1. matriz de datos normalizados, cada columna una categoría
2. Las etiquetas

Se supone que:
1. Class es la última categoría y contiene las categorías.
2. Que  el resto son numéricas
"""
function DataLabelArff(file, class_name)
    Random.seed!(0)
    df = ARFFFiles.load(DataFrame, file)
    labels = df[!, "$class_name"]
    attributes  = names(df)
    pop!(attributes) # Elimina la clase
    data = [
            Normalize(df[!,"$name"])
            for name in attributes
    ]
    data = reduce(hcat, data)

    # desordenamos los datos
    index = shuffle(1:length(labels))
    data = data[index, :]
    labels = labels[index]

    return data, labels
end


text:
module ManifoldLearning

    import Base: show, summary
    import SparseArrays: AbstractSparseArray, SparseMatrixCSC, spzeros, spdiagm, findnz
    import Statistics: mean, std
    import StatsBase: StatsBase, fit, standardize
    import MultivariateStats: outdim, projection,  KernelPCA, transform, transform!,
                              principalvars, dmat2gram, gram2dmat, pairwise
    import LinearAlgebra: eigvals, mul!, svd, qr, Symmetric, eigen, eigen!, tr, rmul!, I, norm, Diagonal
    import LightGraphs: neighbors, nv, add_edge!, connected_components, vertices,
                        dijkstra_shortest_paths, induced_subgraph, weights
    import SimpleWeightedGraphs: SimpleWeightedGraph

    export

    # Transformation types
    AbstractDimensionalityReduction,
    Isomap,             # Type: Isomap model
    HLLE,               # Type: Hessian Eigenmaps model
    LLE,                # Type: Locally Linear Embedding model
    LTSA,               # Type: Local Tangent Space Alignment model
    LEM,                # Type: Laplacian Eigenmaps model
    DiffMap,            # Type: Diffusion maps model

    ## common interface
    outdim,             # the output dimension of the transformation
    fit,                # perform the manifold learning
    transform,          # the projection matrix
    eigvals,            # eigenvalues from the spectral analysis
    neighbors,          # the number of nearest neighbors used for aproximate local subspace
    vertices,           # vertices of largest connected component
    projection          # the projection matrix (deprecated)

    include("interface.jl")
    include("utils.jl")
    include("nearestneighbors.jl")
    include("isomap.jl")
    include("hlle.jl")
    include("lle.jl")
    include("ltsa.jl")
    include("lem.jl")
    include("diffmaps.jl")

    # deprecated functions
    @deprecate transform(DimensionalityReduction, X; k=k, d=d) fit(DimensionalityReduction, X; k=k, maxoutdim=d)
    @deprecate projection(DimensionalityReductionModel) transform(DimensionalityReductionModel)
end


text:
"""
    glmcover1(sys1, sys2; fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol) -> (sysx, sysy, info)

Determine for the proper descriptor systems `sys1 = (A1-λE1,B1,C1,D1)` and
`sys2 = (A2-λE2,B2,C2,D2)` with the transfer function matrices `X1(λ)` and `X2(λ)`,
respectively, using a left minimum dynamic cover of Type 1 based
order reduction, the descriptor systems `sysx` and `sysy` with the
transfer function matrices `X(λ)` and `Y(λ)`, respectively, such that

    X(λ) = X1(λ) + Y(λ)*X2(λ) ,

and `sysx` has order less than the order of `sys1`.

The call with

    glmcover1(sys, p1; fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol) -> (sysx, sysy, info)

uses the compound descriptor system `sys = (A-λE,B,[C1; C2],[D1; D2])`,
where `C1` and `D1` have `p1` rows, to define
the proper descriptor systems `sys1 = (A-λE,B,C1,D1)` and `sys2 = (A-λE,B,C2,D2)`
(i.e., `A1-λE1 = A2-λE2 = A-λE` and `B1 = B2 = B`).

The resulting descriptor systems `sysx` and `sysy` have observable realizations
of the form `sysx = (Ao-λEo,Bo1,Co,D1)` and `sysy = (Ao-λEo,Bo2,Co,0)`,
where the pencil `[Ao-λEo; Co]` is in a (observability) staircase form,
with `νl[i] x νl[i+1]` full row rank diagonal blocks, for `i = 1, ..., nl`,
with `νl[nl+1] := p1`.

The resulting named triple `info` contains `(stdim, tcond, fnorm) `,
where `info.stdim = νl` is a vector which contains the column dimensions of the
blocks of the staircase form `[Ao-λEo; Co]`, `info.tcond` is the maximum of
the Frobenius-norm condition numbers of the employed non-orthogonal
transformation matrices, and `info.fnorm` is
the Frobenius-norm of the (internally) employed output-injection gain to reduce the order.
Large values of  `info.tcond` or `info.fnorm` indicate possible loss of
numerical stability of computations.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively,
the absolute tolerance for the nonzero elements of `A1`, `B1`, `C1`, `D1`, `A2`, `B2`, `C2`, `D2`,
the absolute tolerance for the nonzero elements of `E1` and `E2`,
and the relative tolerance for the nonzero elements of
`A1`, `B1`, `C1`, `D1`, `A2`, `B2`, `C2`, `D2`, `E1` and `E2`.
The default relative tolerance is `n*ϵ`, where `ϵ` is the working machine epsilon
and `n` is the maximal order of the systems `sys1` and `sys2`. The keyword argument `atol` can be used
to simultaneously set `atol1 = atol`, `atol2 = atol`.

The rank determinations in the performed reductions
 are based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

_Note:_ `glmcover1` also works for arbitrary descriptor system `sys1`,
if `sys2` is proper. For an improper system `sys1`, the order
reduction is performed only for the proper part of `sys1`, while the
polynomial part of `sys1` is included without modification in the
resulting realization of `sysx`. In this case, `info.stdim = νl` contains
the information corresponding to the proper part of `sysx`.

_Method:_ The dual of method  of [1] is used to compute Type 1 minimum dynamic covers
for standard systems and the dual of method of [2] for proper descriptor systems.
The resulting McMillan degree of `sysx` is the least achievable one
provided the realization of `sys2` is maximally controllable
(i.e., the pair `(A2+F*C2-λE2,B2+F*D2)` is controllable for any `F`).

_References:_

[1] A. Varga, Reliable algorithms for computing minimal dynamic covers,
    Proc. CDC'03, Maui, Hawaii, 2003.

[2] A. Varga. Reliable algorithms for computing minimal dynamic covers for
    descriptor systems. Proc. MTNS Symposium, Leuven, Belgium, 2004.
"""
function glmcover1(sys1::DescriptorStateSpace, sys2::DescriptorStateSpace; kwargs...)
   size(sys1,2) == size(sys2,2) || throw(DimensionMismatch("sys1 and sys2 must have the same number of inputs"))
   sysx, sysy, info = grmcover1(gdual(sys1), gdual(sys2); kwargs...)
   info.stdim[:] = reverse(info.stdim)
   return gdual(sysx,rev = true), gdual(sysy,rev = true), info
end
function glmcover1(sys::DescriptorStateSpace, p1::Int; kwargs...)
   p = size(sys,1);
   (p1 <= p && p1 >= 0) || throw(DimensionMismatch("p1 must be at most $p, got $p1"))
   sysx, sysy, info = grmcover1(gdual(sys,rev = true), p1; kwargs...)
   info.stdim[:] = reverse(info.stdim)
   return gdual(sysx,rev = true), gdual(sysy,rev = true), info
end
"""
    glmcover2(sys1, sys2; fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol) -> (sysx, sysy, info)

Determine for the proper descriptor systems `sys1 = (A1-λE1,B1,C1,D1)` and
`sys2 = (A2-λE2,B2,C2,D2)` with the transfer function matrices `X1(λ)` and `X2(λ)`,
respectively, using a left minimum dynamic cover of Type 2 based
order reduction, the descriptor systems `sysx` and `sysy` with the
transfer function matrices `X(λ)` and `Y(λ)`, respectively, such that

    X(λ) = X1(λ) + Y(λ)*X2(λ) ,

and `sysx` has order less than the order of `sys1`.

The call with

    glmcover2(sys, p1; fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol) -> (sysx, sysy, info)

uses the compound descriptor system `sys = (A-λE,B, [C1; C2],[D1; D2])`,
where `C1` and `D1` have `p1` rows and `E` is invertible, to define
the proper descriptor systems `sys1 = (A-λE,B,C1,D1)` and `sys2 = (A-λE,B,C2,D2)`
(i.e., `A1-λE1 = A2-λE2 = A-λE` and `B1 = B2 = B`).

The resulting descriptor systems `sysx` and `sysy` have observable realizations
of the form `sysx = (Ao-λEo,Bo1,Co,Do1)` and `sysy = (Ao-λEo,Bo2,Co,Do2)`,
where the pencil `[Ao-λEo; Co]` is in a (observability) staircase form,
with `νl[i] x νl[i+1]` full row rank diagonal blocks, for `i = 1, ..., nl`,
with `νl[nl+1] := p1`.

The resulting named triple `info` contains `(stdim, tcond, fnorm, gnorm) `,
where `info.stdim = νl` is a vector which contains the column dimensions of the
blocks of the staircase form `[Ao-λEo; Co]`, `info.tcond` is the maximum of
the Frobenius-norm condition numbers of the employed non-orthogonal
transformation matrices, `info.fnorm` is
the Frobenius-norm of the (internally) employed output-injection gain to reduce the order, and
`info.gnorm` is the Frobenius-norm of the (internally) employed output-feedforward gain.
Large values of  `info.tcond`,`info.fnorm` or `info.gnorm` indicate possible loss of
numerical stability of computations.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively,
the absolute tolerance for the nonzero elements of `A1`, `B1`, `C1`, `D1`, `A2`, `B2`, `C2`, `D2`,
the absolute tolerance for the nonzero elements of `E1` and `E2`,
and the relative tolerance for the nonzero elements of
`A1`, `B1`, `C1`, `D1`, `A2`, `B2`, `C2`, `D2`, `E1` and `E2`.
The default relative tolerance is `n*ϵ`, where `ϵ` is the working machine epsilon
and `n` is the maximal order of the systems `sys1` and `sys2`. The keyword argument `atol` can be used
to simultaneously set `atol1 = atol`, `atol2 = atol`.

The rank determinations in the performed reductions
 are based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

_Note:_ `glmcover2` also works for arbitrary descriptor system `sys1`,
if `sys2` is proper. For an improper system `sys1`, the order
reduction is performed only for the proper part of `sys1`, while the
polynomial part of `sys1` is included without modification in the
resulting realization of `sysx`. In this case, `info.stdim = νl` contains
the information corresponding to the proper part of `sysx`.

_Method:_ The dual of method  of [1] is used to compute Type 2 minimum dynamic covers
for standard systems and the dual of method of [2] for proper descriptor systems.
The resulting McMillan degree of `sysx` is the least achievable one
provided the realization of `sys2` is maximally controllable
(i.e., the pair `(A2+F*C2-λE2,B2+F*D2)` is controllable for any `F`).

References:

[1] A. Varga, Reliable algorithms for computing minimal dynamic covers,
    Proc. CDC'03, Maui, Hawaii, 2003.

[2] A. Varga. Reliable algorithms for computing minimal dynamic covers for
    descriptor systems. Proc. MTNS Symposium, Leuven, Belgium, 2004.
"""
function glmcover2(sys1::DescriptorStateSpace, sys2::DescriptorStateSpace; kwargs...)
   size(sys1,2) == size(sys2,2) || throw(DimensionMismatch("sys1 and sys2 must have the same number of inputs"))
   sysx, sysy, info = grmcover2(gdual(sys1), gdual(sys2); kwargs...)
   info.stdim[:] = reverse(info.stdim)
   return gdual(sysx,rev = true), gdual(sysy,rev = true), info
end
function glmcover2(sys::DescriptorStateSpace, p1::Int; kwargs...)
   p = size(sys,1);
   (p1 <= p && p1 >= 0) || throw(DimensionMismatch("p1 must be at most $p, got $p1"))
   sysx, sysy, info = grmcover2(gdual(sys,rev = true), p1; kwargs...)
   info.stdim[:] = reverse(info.stdim)
   return gdual(sysx,rev = true), gdual(sysy,rev = true), info
end
"""
    grmcover1(sys1, sys2; fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol) -> (sysx, sysy, info)

Determine for the proper descriptor systems `sys1 = (A1-λE1,B1,C1,D1)` and
`sys2 = (A2-λE2,B2,C2,D2)` with the transfer function matrices `X1(λ)` and `X2(λ)`,
respectively, using a right minimum dynamic cover of Type 1 based
order reduction, the descriptor systems `sysx` and `sysy` with the
transfer function matrices `X(λ)` and `Y(λ)`, respectively, such that

    X(λ) = X1(λ) + X2(λ)*Y(λ) ,

and `sysx` has order less than the order of `sys1`.

The call with

    grmcover1(sys, m1; fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol) -> (sysx, sysy, info)

uses the compound descriptor system `sys = (A-λE,[B1 B2],C,[D1 D2])`,
where `B1` and `D1` have `m1` columns, to define
the proper descriptor systems `sys1 = (A-λE,B1,C,D1)` and `sys2 = (A-λE,B2,C,D2)`
(i.e., `A1-λE1 = A2-λE2 = A-λE` and `C1 = C2 = C`).

The resulting descriptor systems `sysx` and `sysy` have controllable realizations
of the form `sysx = (Ar-λEr,Br,Cr1,D1)` and `sysy = (Ar-λEr,Br,Cr2,0)`,
where the pencil `[Br Ar-λEr]` is in a (controllability) staircase form,
with `νr[i] x νr[i-1]` full row rank diagonal blocks, for `i = 1, ..., nr`,
with `νr[0] := m1`.

The resulting named triple `info` contains `(stdim, tcond, fnorm) `,
where `info.stdim = νr` is a vector which contains the row dimensions of the
blocks of the staircase form `[Br Ar-λEr]`, `info.tcond` is the maximum of
the Frobenius-norm condition numbers of the employed non-orthogonal
transformation matrices, and `info.fnorm` is
the Frobenius-norm of the (internally) employed state-feedback to reduce the order.
Large values of  `info.tcond` or `info.fnorm` indicate possible loss of
numerical stability of computations.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively,
the absolute tolerance for the nonzero elements of `A1`, `B1`, `C1`, `D1`, `A2`, `B2`, `C2`, `D2`,
the absolute tolerance for the nonzero elements of `E1` and `E2`,
and the relative tolerance for the nonzero elements of
`A1`, `B1`, `C1`, `D1`, `A2`, `B2`, `C2`, `D2`, `E1` and `E2`.
The default relative tolerance is `n*ϵ`, where `ϵ` is the working machine epsilon
and `n` is the maximal order of the systems `sys1` and `sys2`. The keyword argument `atol` can be used
to simultaneously set `atol1 = atol`, `atol2 = atol`.

The rank determinations in the performed reductions
 are based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

_Note:_ `grmcover1` also works for arbitrary descriptor system `sys1`,
if `sys2` is proper. For an improper system `sys1`, the order
reduction is performed only for the proper part of `sys1`, while the
polynomial part of `sys1` is included without modification in the
resulting realization of `sysx`. In this case, `info.stdim = νr` contains
the information corresponding to the proper part of `sysx`.

_Method:_ The method  of [1] is used to compute Type 1 minimum dynamic covers
for standard systems and the method of [2] for proper descriptor systems.
The resulting order (McMillan degree) of `sysx` is the least achievable one
provided the realization of `sys2` is maximally observable
(i.e., the pair `(A2+B2*F-λE2,C2+D2*F)` is observable for any `F`).

References:

[1] A. Varga, Reliable algorithms for computing minimal dynamic covers,
    Proc. CDC'03, Maui, Hawaii, 2003.

[2] A. Varga. Reliable algorithms for computing minimal dynamic covers for
    descriptor systems. Proc. MTNS Symposium, Leuven, Belgium, 2004.
"""
function grmcover1(sys::DescriptorStateSpace{T}, m1::Int; atol::Real = zero(real(T)),
                   atol1::Real = atol, atol2::Real = atol,
                   rtol::Real = ((size(sys.A,1)+1)*eps(real(float(one(T)))))*iszero(max(atol1,atol2)),
                   fast::Bool = true) where T

   p, m = size(sys);
   (m1 <= m && m1 >= 0) || throw(DimensionMismatch("m1 must be at most $m, got $m1"))
   m2 = m-m1
   Ts = sys.Ts;

   n = size(sys.A,1);
   isys1 = 1:m1; isys2 = m1+1:m;

   # handle simple cases
   if n == 0
      sysx = sys[:,isys1];
      info = (stdim = Int[], tcond = 1, fnorm = 0)
      sysy = dss(zeros(T,m2,m1),Ts = Ts)
      return sysx, sysy, info
   end

   if  m1 == 0
      sysx = dss(zeros(T,p,m1),Ts = Ts);
      info = (stdim = Int[], tcond = 1, fnorm = 0)
      sysy = dss(zeros(T,m2,m1),Ts = Ts)
      return sysx, sysy, info
   end

   T1 = T <: BlasFloat ? T : promote_type(Float64,T)
   At, Et, Bt, Ct, D = dssdata(T1,sys)
   sstype = (Et == I)
   sstype || istriu(Et) || MatrixPencils._qrE!(At, Et, nothing, Bt; withQ = false)
   ONE = one(T1)
   ZERO = zero(T1)

   # handle improper case
   if !sstype && LinearAlgebra.LAPACK.trcon!('1','U','N',triu(Et))  < n*eps(real(T1))
      isproper(sys[:,isys2], atol1 = atol1, atol2 = atol2, rtol = rtol) ||
             error("The system sys2 must be proper")
      At, Et, Bt, Ct, _, _, _, blkdims = fiblkdiag(At, Et, Bt, Ct; fast = fast, finite_infinite = true,
                                                   atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = false, withZ = false)
      n1 = blkdims[1]; i1 = 1:n1; i2 = n1+1:n
      sysi1 = dss(lsminreal(At[i2,i2], Et[i2,i2], Bt[i2,isys1], Ct[:,i2], zeros(T,p,m1),
                  fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)[1:5]...,Ts = Ts)
      γ = iszero(Ts) ? ZERO : ONE
      DCgain = Ct[:,i2]*((lmul!(γ,view(Et,i2,i2))-view(At,i2,i2))\view(Bt,i2,isys2))

      sysfm = dss(At[i1,i1], Et[i1,i1], Bt[i1,:], Ct[:,i1], D + [zeros(T,p,m1) DCgain], Ts = Ts)
      sysfx, sysy, info = grmcover1(sysfm, m1, fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
      return sysfx+sysi1, sysy, info
   end

   # reduce to the special controllability form of [1] or [2]
   if sstype
       _, tau, ncont, nuc = sklf_right2!(At, Bt, m1, Ct; fast = fast, atol1 = atol1, atol2 = atol2,
                                         rtol = rtol, withQ = false)
   else
       _, _, tau, ncont, nuc = sklf_rightfin2!(At, Et, Bt, m1, Ct; fast = fast, atol1 = atol1, atol2 = atol2,
                                               rtol = rtol, withQ = false, withZ = false)
   end
   if ncont == 0
       sysx = dss(D[:,isys1],Ts = Ts);
       sysx.Ts = Ts;
       info = (stdim = [], tcond = 1, fnorm = 0)
       sysy = dss(zeros(T,m2,m1),Ts = Ts)
       return sysx, sysy, info
   end

   if  m2 == 0
       if sstype
          sysx = dss(At,I,Bt,Ct,D,Ts = Ts);
       else
          sysx = dss(At,Et,Bt,Ct,D,Ts = Ts);
       end
       info = (stdim = tau[tau .> 0], tcond = 1, fnorm = 0)
       sysy = dss(zeros(T,m2,m1),Ts = Ts)
       return sysx, sysy, info
   end

   # perform permutation for Type I cover
   p = length(tau)
   cind = Int(p/2)
   ind = zeros(Int,ncont);
   iodd = 1:2:p
   ieven = 2:2:p
   ind1 = tau[iodd]; n1 = sum(ind1)
   ind2 = tau[ieven]; n2 = sum(ind2)
   i2 = 1; i1 = 1; ioff1 = 0; ioff2 = 0;
   for i = 1:cind
       i1i = ind1[i]
       i2i = ind2[i]
       ioff1 +=  i1i
       i2m = i2+i2i-1
       ind[(i2:i2m) .+ n1] =  Vector(i2:i2m) .+ ioff1
       i2 = i2m+1
       i1m = i1+i1i-1
       ind[i1:i1m] =  Vector(i1:i1m) .+ ioff2
       ioff2 +=  i2i
       i1 = i1m+1
   end
   At = At[ind,ind]; Bt = Bt[ind,:]; Ct = Ct[:,ind];
   sstype || (Et = Et[ind,ind])

   # anihilate lower left blocks of At and Et
   tcond = real(T1)(ncont); # the Frobenius condition number of identity matrix
   nlowc1 = n1;             nlowc2 = n1+n2-ind2[cind];
   nlowr1 = n1+n2;          nlowr2 = nlowr1;
   for i = cind:-1:2
      nlowc1 = nlowc1-ind1[i]; nlowc2 = nlowc2-ind2[i-1];
      nlowr1 = nlowr1-ind2[i]; nlowr2 = nlowr2-ind2[i];
      ic1 = nlowc1+1:n1; ic2 = nlowc2+1:nlowc2+ind2[i-1];
      ir2 =  nlowr2+1:nlowr2+ind2[i];
      x = -At[ir2,ic2]\At[ir2,ic1];
      tcond = max(tcond,ncont+norm(x));
      ir = 1:nlowr1;
      mul!(view(At,ir,ic1), view(At,ir,ic2), x, ONE, ONE)
      At[ir2,ic1] .= ZERO
      mul!(view(Ct,:,ic1), view(Ct,:,ic2), x, ONE, ONE)
      if sstype
         x = -x;
         mul!(view(At,ic2,:), x, view(At,ic1,:), ONE, ONE)
         mul!(view(Bt,ic2,:), x, view(Bt,ic1,:), ONE, ONE)
      else
         mul!(view(Et,ir,ic1), view(Et,ir,ic2), x, ONE, ONE)
         y = -Et[ic2,ic1]/Et[ic1,ic1];
         tcond = max(tcond,ncont+norm(y));
         mul!(view(Et,ic2,:), y, view(Et,ic1,:), ONE, ONE)
         mul!(view(At,ic2,:), y, view(At,ic1,:), ONE, ONE)
         mul!(view(Bt,ic2,:), y, view(Bt,ic1,:), ONE, ONE)
      end
   end

   # form the reduced system
   ic = 1:n1; i3 = n1+1:n1+ind2[1];
   f2 = -Bt[i3,isys2]\At[i3,ic];
   mul!(view(At,ic,ic), view(Bt,ic,isys2), f2, ONE, ONE)
   mul!(view(Ct,:,ic), view(D,:,isys2), f2, ONE, ONE)
   if sstype
      sysx = dss(At[ic,ic], I, Bt[ic,isys1], Ct[:,ic], D[:,isys1], Ts = Ts)
      sysy = dss(At[ic,ic], I, Bt[ic,isys1], f2, zeros(T1,m2,m1), Ts = Ts)
   else
      sysx = dss(At[ic,ic], Et[ic,ic], Bt[ic,isys1], Ct[:,ic], D[:,isys1],Ts = Ts)
      sysy = dss(At[ic,ic], Et[ic,ic], Bt[ic,isys1], f2, zeros(T1,m2,m1),Ts = Ts)
   end
   info = (stdim = ind1, tcond = tcond, fnorm = norm(f2))

   return sysx, sysy, info
   # end GRMCOVER1
end
function grmcover1(sys1::DescriptorStateSpace, sys2::DescriptorStateSpace; kwargs...)
    size(sys1,1) == size(sys2,1) || throw(DimensionMismatch("sys1 and sys2 must have the same number of outputs"))
    return grmcover1(gir([sys1 sys2]; kwargs...), size(sys1,2); kwargs...)
end
"""
    grmcover2(sys1, sys2; fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol) -> (sysx, sysy, info)

Determine for the proper descriptor systems `sys1 = (A1-λE1,B1,C1,D1)` and
`sys2 = (A2-λE2,B2,C2,D2)` with the transfer function matrices `X1(λ)` and `X2(λ)`,
respectively, using a right minimum dynamic cover of Type 2 based
order reduction, the descriptor systems `sysx` and `sysy` with the
transfer function matrices `X(λ)` and `Y(λ)`, respectively, such that

    X(λ) = X1(λ) + X2(λ)*Y(λ) ,

and `sysx` has order less than the order of `sys1`.

The call with

    grmcover2(sys, m1; fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol) -> (sysx, sysy, info)

uses the compound descriptor system `sys = (A-λE,[B1 B2],C,[D1 D2])`,
where `B1` and `D1` haves `m1` columns, to define
the proper descriptor systems `sys1 = (A-λE,B1,C,D1)` and `sys2 = (A-λE,B2,C,D2)`
(i.e., `A1-λE1 = A2-λE2 =: A-λE` and `C1 = C2 =: C`).

The resulting descriptor systems `sysx` and `sysy` have controllable realizations
of the form `sysx = (Ar-λEr,Br,Cr1,Dr1)` and `sysy = (Ar-λEr,Br,Cr2,Dr2)`,
where the pencil `[Br Ar-λEr]` is in a (controllability) staircase form,
with `νr[i] x νr[i-1]` full row rank diagonal blocks, for `i = 1, ..., nr`,
with `νr[0] := m1`.

The resulting named triple `info` contains `(stdim, tcond, fnorm, gnorm) `,
where `info.stdim = νr` is a vector which contains the row dimensions of the
blocks of the staircase form `[Br Ar-λEr]`, `info.tcond` is the maximum of
the Frobenius-norm condition numbers of the employed non-orthogonal
transformation matrices, `info.fnorm` is
the Frobenius-norm of the (internally) employed state-feedback gain to reduce the order,
`info.gnorm` is
the Frobenius-norm of the (internally) employed feedforward gain to reduce the order.
Large values of  `info.tcond`, `info.fnorm`  or  `info.gnorm` indicate possible loss of
numerical stability of computations.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively,
the absolute tolerance for the nonzero elements of `A1`, `B1`, `C1`, `D1`, `A2`, `B2`, `C2`, `D2`,
the absolute tolerance for the nonzero elements of `E1` and `E2`,
and the relative tolerance for the nonzero elements of
   `A1`, `B1`, `C1`, `D1`, `A2`, `B2`, `C2`, `D2`, `E1` and `E2`.
The default relative tolerance is `n*ϵ`, where `ϵ` is the working machine epsilon
and `n` is the maximal order of the systems `sys1` and `sys2`. The keyword argument `atol` can be used
to simultaneously set `atol1 = atol`, `atol2 = atol`.

The rank determinations in the performed reductions
 are based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

_Note:_ `grmcover2` also works for arbitrary descriptor system `sys1`,
if `sys2` is proper. For an improper system `sys1`, the order
reduction is performed only for the proper part of `sys1`, while the
polynomial part of `sys1` is included without modification in the
resulting realization of `sysx`. In this case, `info.stdim = νr` contains
the information corresponding to the proper part of `sysx`.

_Method:_ The method  of [1] is used to compute Type 2 minimum dynamic covers
for standard systems and the method of [2] for proper descriptor systems.
The resulting McMillan degree of `sysx` is the least achievable one
provided the realization of `sys2` is maximally observable
(i.e., the pair `(A2+B2*F-λE2,C2+D2*F)` is observable for any `F`).

References:

[1] A. Varga, Reliable algorithms for computing minimal dynamic covers,
    Proc. CDC'03, Maui, Hawaii, 2003.

[2] A. Varga. Reliable algorithms for computing minimal dynamic covers for
    descriptor systems. Proc. MTNS Symposium, Leuven, Belgium, 2004.

"""
function grmcover2(sys::DescriptorStateSpace{T},m1::Int; atol::Real = zero(real(T)),
   atol1::Real = atol, atol2::Real = atol,
   rtol::Real = ((size(sys.A,1)+1)*eps(real(float(one(T)))))*iszero(max(atol1,atol2)),
   fast::Bool = true) where T

   p, m = size(sys);
   (m1 <= m && m1 >= 0) || throw(DimensionMismatch("m1 must be at most $m, got $m1"))
   m2 = m-m1
   Ts = sys.Ts;

   # exchange the roles of B1 and B2
   n = size(sys.A,1);
   isys1 = m1+1:m; isys2 = 1:m1;

   # handle symple cases
   if n == 0
       sysx = sys[:,isys2];
       info = (stdim = Int[], tcond = 1, fnorm = 0, gnorm = 0)
       sysy = dss(zeros(T,m2,m1),Ts = Ts)
       return sysx, sysy, info
   end

   if  m1 == 0
      sysx = dss(zeros(T,p,m1),Ts = Ts);
      info = (stdim = Int[], tcond = 1, fnorm = 0, gnorm = 0)
      sysy = dss(zeros(T,m2,m1),Ts = Ts)
      return sysx, sysy, info
   end

   T1 = T <: BlasFloat ? T : promote_type(Float64,T)
   At, Et, Bt, Ct, D = dssdata(T1,sys)
   sstype = (Et == I)
   ONE = one(T1)
   ZERO = zero(T1)

   sstype || istriu(Et) || MatrixPencils._qrE!(At, Et, nothing, Bt; withQ = false)

   # handle improper case
   if !sstype && LinearAlgebra.LAPACK.trcon!('1','U','N',triu(Et))  < n*eps(real(T1))
      isproper(sys[:,m1+1:m], atol1 = atol1, atol2 = atol2, rtol = rtol) ||
             error("The system sys2 must be proper")
      # sysf, sysi = gsdec(sys, job = "finite", fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
      # sysfm = sysf + [zeros(p,m1) dcgain(sysi[:,m1+1:m]) ];
      # sysfx, sysy, info = grmcover2(sysfm, m1, fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
      # # sysf.D[:,:] = [sysf.D[:,1:m1] sysf.D[:,m1+1:m]+evalfr(sysi[:,m1+1:m],rand())]
      # # sysfx, sysy, info = grmcover2(sysf, m1, fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
      # return sysfx+gminreal(sysi[:,1:m1],fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol), sysy, info
      At, Et, Bt, Ct, _, _, _, blkdims = fiblkdiag(At, Et, Bt, Ct; fast = fast, finite_infinite = true,
                                                   atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = false, withZ = false)
      n1 = blkdims[1]; i1 = 1:n1; i2 = n1+1:n
      sysi1 = dss(lsminreal(At[i2,i2], Et[i2,i2], Bt[i2,1:m1], Ct[:,i2], zeros(T,p,m1),
                  fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)[1:5]...,Ts = Ts)
      γ = iszero(Ts) ? ZERO : ONE
      DCgain = Ct[:,i2]*((lmul!(γ,view(Et,i2,i2))-view(At,i2,i2))\view(Bt,i2,m1+1:m))

      sysfm = dss(At[i1,i1], Et[i1,i1], Bt[i1,:], Ct[:,i1], D + [zeros(T,p,m1) DCgain], Ts = Ts)
      sysfx, sysy, info = grmcover2(sysfm, m1, fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
      return sysfx+sysi1, sysy, info
   end

   ind21 = [Vector(isys1); Vector(isys2)]
   Bt = Bt[:,ind21]

   # reduce to the special controllability form of [1]
   if sstype
      _, tau, ncont, nuc = sklf_right2!(At, Bt, m2, Ct; fast = fast, atol1 = atol1, atol2 = atol2,
                                        rtol = rtol, withQ = false)
   else
      _, _, tau, ncont, nuc = sklf_rightfin2!(At, Et, Bt, m2, Ct; fast = fast, atol1 = atol1, atol2 = atol2,
                                              rtol = rtol, withQ = false, withZ = false)
   end

   if ncont == 0
       sysx = dss(D[:,isys2],Ts = Ts);
       sysx.Ts = Ts;
       info = (stdim = [], tcond = 1, fnorm = 0, gnorm = 0)
       sysy = dss(zeros(T,m2,m1),Ts = Ts)
       return sysx, sysy, info
   end

   if  m2 == 0
      if sstype
         sysx = dss(At,I,Bt,Ct,D,Ts = Ts);
      else
         sysx = dss(At,Et,Bt,Ct,D,Ts = Ts);
      end
      info = (stdim = tau[tau .> 0], tcond = 1, fnorm = 0, gnorm = 0)
      sysy = dss(zeros(T,m2,m1),Ts = Ts)
      return sysx, sysy, info
   end

   # perform permutation for type II cover
   p = length(tau)
   cind = Int(p/2)
   ind = zeros(Int,ncont);
   iodd = 1:2:p
   ieven = 2:2:p
   ind1 = tau[iodd]; n1 = sum(ind1)
   ind2 = tau[ieven]; n2 = sum(ind2)
   i2 = 1; i1 = 1; ioff1 = 0; ioff2 = 0;
  for i = 1:cind
      i1i = ind1[i]
      i2i = ind2[i]
      ioff1 +=  i1i
      i2m = i2+i2i-1
      ind[i2:i2m] =  Vector(i2:i2m) .+ ioff1
      i2 = i2m+1
      i1m = i1+i1i-1
      ind[(i1:i1m) .+ n2] =  Vector(i1:i1m) .+ ioff2
      ioff2 +=  i2i
      i1 = i1m+1
   end
   At = At[ind,ind]; Bt = Bt[ind,:]; Ct = Ct[:,ind]
   sstype || (Et = Et[ind,ind])
   ONE = one(T1)
   ZERO = zero(T1)

   # anihilate lower left blocks of At and Et
   tcond = real(T1)(ncont) # the Frobenius condition number of identity matrix
   nlowc1 = n2-ind2[cind];  nlowc2 = n1+n2-ind1[cind];
   nlowr1 = n1+n2;          nlowr2 = nlowr1-ind1[cind];
   for i = cind:-1:1
      if i == 1
         ice1 = 1:n2; ir1 = n2+1:n2+ind1[i];
      else
         nlowc1 = nlowc1-ind2[i-1]; nlowc2 = nlowc2-ind1[i-1];
         nlowr1 = nlowr1-ind1[i]; nlowr2 = nlowr2-ind1[i-1];
         ic1 = nlowc1+1:n2; ic2 = nlowc2+1:nlowc2+ind1[i-1];
         ir1 =  nlowr1+1:nlowr1+ind1[i];
         ice1 = nlowc1+ind2[i-1]+1:n2;
      end
      if ~sstype
         y = -Et[ir1,ice1]/Et[ice1,ice1];
         tcond = max(tcond,ncont+norm(y));
         mul!(view(Et,ir1,:), y, view(Et,ice1,:), ONE, ONE)
         mul!(view(At,ir1,:), y, view(At,ice1,:), ONE, ONE)
         mul!(view(Bt,ir1,:), y, view(Bt,ice1,:), ONE, ONE)
      end
      i == 1 && break
      x = -At[ir1,ic2]\At[ir1,ic1];
      tcond = max(tcond,ncont+norm(x));
      ir = 1:nlowr1;
      mul!(view(At,ir,ic1), view(At,ir,ic2), x, ONE, ONE)
      At[ir1,ic1] .= ZERO
      mul!(view(Ct,:,ic1), view(Ct,:,ic2), x, ONE, ONE)
      if sstype
         x = -x;
         mul!(view(At,ic2,:), x, view(At,ic1,:), ONE, ONE)
         mul!(view(Bt,ic2,:), x, view(Bt,ic1,:), ONE, ONE)
      else
         mul!(view(Et,ir,ic1), view(Et,ir,ic2), x, ONE, ONE)
      end
   end

   # form the reduced system
   ic = 1:n2; i3 = n2+1:n2+ind1[1];
   f2 = -Bt[i3,1:m2]\At[i3,ic];
   g = -Bt[i3,1:m2]\Bt[i3,m2+1:m];
   mul!(view(Ct,:,ic), view(D,:,isys1), f2, ONE, ONE)
   mul!(view(D,:,isys2), view(D,:,isys1), g, ONE, ONE)
   if sstype
      sysx = dss(At[ic,ic], I, Bt[ic,m2+1:m], Ct[:,ic], D[:,isys2], Ts = Ts)
      sysy = dss(At[ic,ic], I, Bt[ic,m2+1:m], f2, g, Ts = Ts)
   else
      sysx = dss(At[ic,ic], Et[ic,ic], Bt[ic,m2+1:m], Ct[:,ic], D[:,isys2], Ts = Ts)
      sysy = dss(At[ic,ic], Et[ic,ic], Bt[ic,m2+1:m], f2, g, Ts = Ts)
   end
   info = (stdim = ind2[ind2 .> 0], tcond = tcond, fnorm = norm(f2), gnorm = norm(g))

   return sysx, sysy, info

   # end GRMCOVER2
end

function grmcover2(sys1::DescriptorStateSpace, sys2::DescriptorStateSpace; kwargs...)
    size(sys1,1) == size(sys2,1) || throw(DimensionMismatch("sys1 and sys2 must have the same number of outputs"))
    return grmcover2(gir([sys1 sys2]; kwargs...), size(sys1,2);kwargs...)
end


text:
module ModelingStructs

using Requires

include("utils.jl")
include("MStruct.jl")
include("set_get.jl")
include("attributes.jl")
include("similar_convert.jl")
include("show.jl")

function __init__()
    @require Optim = "429524aa-4258-5aef-a3af-852621145aeb" include("../patches/optim.jl")
end

export ModelingStruct, MStruct

end # module

max_stars_repo_name: jonniedie/NamedViewVectors


text:
using RecipesBase

struct AcrobotParams{T} <:Any where {T<:Real}
    m1::T 	# mass of first link
    m2::T 	# mass of second link
    l1::T 	# length of first link
	l2::T 	# length of second link
	lc1::T 	# position of center mass of link1
	lc2::T 	# position of center mass of link2
	i1::T  	# link1 moment of inertia
	i2::T  	# link2 moment of inertia
	g::T   	# gravity force (not directional)
	fmag::T # max force magnitude

    AcrobotParams() = new{Float64}(1., 1., 1., 1., 0.5, 0.5, 1., 1., 9.8, 1.)
    AcrobotParams(T::Type) = new{T}(1., 1., 1., 1., 0.5, 0.5, 1., 1., 9.8, 1.)
	AcrobotParams(T::Type, m1, m2, l1, l2, g) = new{T}(m1, m2, l1, l2, 0.5*l1, 0.5*l2, 1., 1., g, 1.)
end

mutable struct Acrobot{T,TA} <: AbstractEnvironment where {T<:Real,TA<:Real}
    state::Array{T,1}
    reward::T
    done::Bool
    params::AcrobotParams{T}
	dt::T
	integShritte::T  # larger this is more accurate Runge-Kutta approximation
    t::T # time
	deterministic_start::Bool

    function Acrobot(T::Type=Float64, TA::Type=Int, deterministic_start::Bool=false)
        new{T,TA}(zeros(T,4), T(0.0), false,
           AcrobotParams(T), T(0.2), T(10.), T(0.), deterministic_start
       )
    end

	function Acrobot(::Type{T}, ::Type{TA}, deterministic_start::Bool, params::AcrobotParams{T}) where {T, TA}
        new{T,TA}(zeros(T,4), T(0.0), false,
           params, T(0.2), T(10.), T(0.), deterministic_start
       )
    end
end

function get_gamma(env::Acrobot{T})::T where {T}
	return T(1.0)
end

function randomize_env(::Type{Acrobot{T,TA}}, rng::AbstractRNG, deterministic_start::Bool=false) where {T, TA}
	m1 = rand(rng, Uniform(0.8, 1.2)) 	# mass of link 1
	l1 = rand(rng, Uniform(0.8, 1.2))   # length of link 1
	m2 = rand(rng, Uniform(0.8, 1.2))  	# mass of link 2
	l2 = rand(rng, Uniform(0.8, 1.2))   # length of link 2
	g = rand(rng, Uniform(8., 10.))     # gravity
	params = AcrobotParams(T, m1, m2, l1, l2, g)
	return Acrobot(T, TA, deterministic_start, params)
end

function get_state_desc(env::Acrobot{T})::Array{T,2} where {T<:Real}
    ranges = zeros(T, (4,2))
    ranges[1,:] .= [-π, π]           # theta1 range
	ranges[2,:] .= [-π, π]           # theta2 range
	ranges[3,:] .= [-4. * π, 4. * π]       # theta1Dot range
	ranges[4,:] .= [-9. * π, 9. * π]       # theta2Dot range
    return ranges
end

function get_action_desc(env::Acrobot{T,TA})::Int where{T<:Real,TA<:Integer}
    return 3
end

function get_action_desc(env::Acrobot{T,TA})::Array{TA,1} where{T<:Real,TA<:AbstractFloat}
	ranges = zeros(TA, (1,2))
	ranges[1, 1] = -1.
	ranges[1, 2] = 1.
	return ranges
end

function rk_helper!(buff::Array{T, 1}, s::Array{T, 1}, tau::T, params::AcrobotParams{T}) where {T}
	m1, m2, l1, l2, lc1, lc2, i1, i2, g = params.m1, params.m2, params.l1, params.l2, params.lc1, params.lc2, params.i1, params.i2, params.g

	d1 = m1*lc1^2 + m2*(l1^2 + lc2^2 + 2. * l1*lc2*cos(s[2])) + i1 + i2
	d2 = m2*(lc2^2 + l1*lc2*cos(s[2])) + i2

	phi2 = m2*lc2*g*cos(s[1] + s[2] - (π / 2.))
	phi1 = (-m2*l1*lc2*s[4]^2 * sin(s[2]) - 2. * m2*l1*lc2*s[4]*s[3] * sin(s[2]) + (m1*lc1 + m2*l1)*g*cos(s[1] - (π / 2.)) + phi2)

	newa2 = ((1. / (m2*lc2^2 + i2 - (d2^2) / d1)) * (tau + (d2 / d1)*phi1 - m2*l1*lc2*s[3]^2 * sin(s[2]) - phi2))
	newa1 = ((-1. / d1) * (d2*newa2 + phi1))
	buff[1] = s[3]
	buff[2] = s[4]
	buff[3] = newa1
	buff[4] = newa2
end

function Acrobot_sim!(state::Array{T,1}, u::T, params::AcrobotParams{T}, dt, integShritte) where {T}
	theta1, theta2, theta1Dot, theta2Dot = state
	hilf = deepcopy(state)
	s0_dot = zeros(T, 4)
	s1_dot = zeros(T, 4)
	s2_dot = zeros(T, 4)
	s3_dot = zeros(T, 4)
	ss = zeros(T, 4)
	s1 = zeros(T, 4)
	s2 = zeros(T, 4)
	s3 = zeros(T, 4)

	h = dt / integShritte

	for i in 1:integShritte
		rk_helper!(s0_dot, hilf, u, params)
		@. s1 = hilf + (h / 2.) * s0_dot

		rk_helper!(s1_dot, s1, u, params)
		@. s2 = hilf + (h / 2.) * s1_dot

		rk_helper!(s2_dot, s2, u, params)
		@. s3 = hilf + (h / 2.) * s2_dot

		rk_helper!(s3_dot, s3, u, params)
		@. hilf = hilf + (h / 6.) * (s0_dot + 2. * (s1_dot + s2_dot) + s3_dot)
	end

	@. ss = hilf

	theta1 = mod(ss[1] + π, 2. * π) - π
	theta2 = mod(ss[2] + π, 2. * π) - π

	theta1Dot = clamp(ss[3], -4. * π, 4. * π)
	theta2Dot = clamp(ss[4], -9. * π, 9. * π)

	state .= theta1, theta2, theta1Dot, theta2Dot

end

function step!(env::Acrobot{T}, action::Int, rng::AbstractRNG)::T where {T<:Real}
    if action <= 0 || action > 3
        error("Action needs to be an integer in [1, 3]")
    end
    u::T = 0.
	u = (T(action) - 2.) * env.params.fmag

	Acrobot_sim!(env.state, u, env.params, env.dt, env.integShritte)
	env.t += env.dt

    env.reward = -0.1
	env.done = terminalCheck(env)
	if env.done
		env.reward = 0.
	end

	return env.reward
end

function terminalCheck(env::Acrobot)::Bool
	elbowY = -env.params.l1*cos(env.state[1])
	handY = elbowY - env.params.l2*cos(env.state[1] + env.state[2])
	anglecond = handY > env.params.l1
	timecond = env.t ≥ 400.
	done = anglecond | timecond

	return done
end

function step!(env::Acrobot{T}, action::T, rng::AbstractRNG)::T where {T<:AbstractFloat}
	u = clamp(action, -1., 1.) * env.params.fmag

	Acrobot_sim!(env.state, u, env.params, env.dt, env.integShritte)
	env.t += env.dt

    env.reward = -0.1
	env.done = terminalCheck(env)
	if env.done
		env.reward = 0.
	end
end

function get_state(env::Acrobot{T})::Array{T,1} where {T<:Real}
    return env.state
end

function get_reward(env::Acrobot)
    return env.reward
end

function is_terminal(env::Acrobot)::Bool
    return env.done
end

function env_reset!(env::Acrobot, rng)
	fill!(env.state, 0.)
	if !env.deterministic_start
		env.state[1:2] = rand(rng, Uniform(-π * 5. / 180., π * 5. / 180.), 2)
	end
    env.t = 0.
	env.done = false
end

function clone(env::Acrobot{T,TA}) where {T<:Real,TA<:Union{Int,<:AbstractFloat}}
	env2 = Acrobot(T,TA, env.deterministic_start)
	env2.state .= deepcopy(env.state)
	env2.reward = deepcopy(env.reward)
	env2.done = deepcopy(env.done)
	env2.t = deepcopy(env.t)
	env2.params = deepcopy(env.params)
	env2.dt = env.dt
	env2.integShritte = env.integShritte
	return env2
end

@recipe function f(env::Acrobot)
	theta1, theta2 = env.state[1:2]
	l1, l2 = env.params.l1, env.params.l2
	maxlen = (l1+l2)*1.2
	x1,y1 = 0., 0.
	x2 = x1 - sin(theta1)*l1
	y2 = y1 - cos(theta1)*l1
	x3 = x2 - sin(theta1 + theta2)*l2
	y3 = y2 - cos(theta1 + theta2)*l2

	legend := false
	xlims := (-maxlen, maxlen)
	ylims := (-maxlen, maxlen)
	grid := false
	ticks := nothing
	foreground_color := :white
	aspect_ratio := 1.

	# target line
	@series begin
		linecolor := :red
		linewidth := 5

		[0., 0.], [0., maxlen]
	end

	# arms
	@series begin
		linecolor := :black
		linewidth := 10

		[x1, x2, x3], [y1, y2, y3]
	end

end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs